# simulate an artificial data frame
# with survival response and two predictors
set.seed(280180)
N <- 300
X1 <- rnorm(N,mean=10,sd=5)
X2 <- rbinom(N,1,.4)
linPred <- .00001+0.2*X1+2.3*X2
T <- sapply(linPred,function(lp){rexp(n=1,exp(-lp))})
C <- rexp(n=300,exp(-mean(linPred)))
dat <- data.frame(time=pmin(T,C),status=as.numeric(T<=C),X1=X1,X2=X2)
# fit some candidate Cox models and compute the Kaplan-Meier estimate
Models <- list("Kaplan.Meier"=survfit(Surv(time,status)~1,data=dat),
"Cox.X1"=coxph(Surv(time,status)~X1,data=dat),
"Cox.X2"=coxph(Surv(time,status)~X2,data=dat),
"Cox.X1.X2"=coxph(Surv(time,status)~X1+X2,data=dat))
# compute the .632+ estimate of the generalization error
set.seed(17100)
PredError.632plus <- pec(object=Models,
formula=Surv(time,status)~X1+X2,
data=dat,
exact=TRUE,
cens.model="marginal",
splitMethod="boot632plus",
B=10,
keep.matrix=TRUE,
verbose=TRUE)
# plot the .632+ estimates of the generalization error
plot(PredError.632plus,xlim=c(0,30))
# plot the bootstrapped curves, .632+ estimates of the generalization error
# and Apparent error for the Cox model 'Cox.X1' with the 'Cox.X2' model
# as benchmark
plot(PredError.632plus,xlim=c(0,30), models="Cox.X1", special=TRUE, special.bench="Cox.X2", special.benchcol=2, special.addprederr="AppErr")
Run the code above in your browser using DataLab