riskRegression (version 1.4.3)

influenceTest: Influence test [Experimental!!]

Description

Compare two estimates using their influence function

Usage

influenceTest(predict1, predict2, type)

Arguments

predict1

a list containing an estimate and its influence function

predict2

same as predict1 but for another model

type

the type of estimate

Examples

# NOT RUN {
warperNoStrata <- function(n, ...){
  newdata <- data.frame(X1=0,X6=0)
  time <- 1
  
m <- lvm()
regression(m) <- y ~ X6 + X1
distribution(m,~X1) <- binomial.lvm()
distribution(m,~cens) <- coxWeibull.lvm(scale=1)
distribution(m,~y) <- a <- coxWeibull.lvm(scale=0.3)
eventTime(m) <- eventtime ~ min(y=1,cens=0)
d <- as.data.table(sim(m,n))
setkey(d, eventtime)
  
m.cox <- coxph(Surv(eventtime, status) ~ X1+X6, 
               data = d, y = TRUE, x = TRUE)
survNoStrata <- predictCox(m.cox, type = "survival",
                           newdata = newdata, times = time, iid = TRUE, se = TRUE)
  
mStrata.cox <- coxph(Surv(eventtime, status) ~ strata(X1)+X6, 
                     data = d, y = TRUE, x = TRUE)
survStrata <- predictCox(mStrata.cox, type = "survival", 
                         newdata = newdata, times = time, iid = TRUE, se = TRUE)
  
res <- influenceTest(survStrata, survNoStrata, type = "survival")
return(res)
}

warperStrata <- function(n, ...){
  newdata <- data.frame(X1=0,X6=0)
time <- 1
  
m <- lvm()
regression(m) <- y ~ X6
regression(m) <- s ~ exp(-2*X1)
distribution(m,~X1) <- binomial.lvm()
distribution(m,~cens) <- coxWeibull.lvm(scale=1)
distribution(m,~y) <- a <- coxWeibull.lvm(scale=0.3,shape=~s)
eventTime(m) <- eventtime ~ min(y=1,cens=0)
d <- as.data.table(sim(m,n))
  setkey(d, eventtime)
  
m.cox <- coxph(Surv(eventtime, status) ~ X1+X6, data = d, y = TRUE, x = TRUE)
survNoStrata <- predictCox(m.cox, type = "survival", 
                           newdata = newdata, times = time, iid = TRUE, se = TRUE)
  
mStrata.cox <- coxph(Surv(eventtime, status) ~ strata(X1)+X6, data = d, y = TRUE, x = TRUE)
survStrata <- predictCox(mStrata.cox, type = "survival", 
                         newdata = newdata, times = time, iid = TRUE, se = TRUE)
res <- influenceTest(survStrata, survNoStrata, type = "survival")
return(res)
}

n <- 500
resNoStrata <- pbsapply(1:200, warperNoStrata, n = n)

sd(unlist(resNoStrata["delta",]))
quantile(unlist(resNoStrata["se.delta",]))

mean(unlist(resNoStrata["p.value",])<0.05)

resStrata <- pbsapply(1:100, warperStrata, n = n)

sd(unlist(resStrata["delta",]))
quantile(unlist(resStrata["se.delta",]))

mean(unlist(resStrata["p.value",])<0.05)
# }