## ## survivalROCsim.R ##---------------------------------------------------------------------- ## ## PURPOSE: simulate data and create C/D ROC estimate using NNE ## ## AUTHOR: (originally created by) P. Heagerty, (edited by) P. Saha ## ## DATE: 12 December 2006 ## ##----------------------------------------------------------------------# ## PART 1 ## TRUTH SET UP AND COMPUTATION ## For the computation of true (FP, TP), R library mvtnorm NEED TO BE ## DOWNLOADED AND INSTALLED. See download.packages() and ## install.packages() for instruction how to download and install new ## packages. library(mvtnorm) rho <- -0.70 FPvalues <- c(1:99)/100 TargetTime <- 1.0 TargetTP <- rep( NA, length(FPvalues) ) TargetThreshold <- rep( NA, length(FPvalues) ) Cmat <- matrix( c(1,rho,rho,1), ncol=2, nrow=2 ) for( j in 1:length(TargetTP) ) { ddd <- 0.10 NewHigh <- -3.5 NewLow <- 3.5 iter <- 0 TOL <- 1e-6 converge <- F TheFnxTarget <- FPvalues[j] MAXITS <- 50 while( iter= TheFnxTarget ) NewHigh <- TheMid if( TheFnx <= TheFnxTarget ) NewLow <- TheMid delta <- abs( NewHigh - NewLow ) if( delta < TOL ) converge <- T iter <- iter+1 } TheMid <- (NewHigh+NewLow)/2 if( itercTime ] <- cTime[ y>cTime ] survival.status[ y>cTime ] <- 0 ## ## # ## # NNE estimator ## # ## fitNNE <- survivalROC.C( Stime=10+survival.time, status=survival.status, marker= x, predict.time = 10 + 1.0, span = 0.25*nobs^(-0.20) ) ## ## ### Interpolate to save at FPvalues selected ## TPout <- rep( NA, length(FPvalues) ) for( j in 1:length(TPout) ){ eval.x <- FPvalues[j] nX <- length( fitNNE$FP ) ## ## Note: reverse order for FP ## LowIndex <- min( c(1:nX)[ fitNNE$FP <= eval.x ] ) HighIndex <- max( c(1:nX)[ fitNNE$FP >= eval.x ] ) delta.x <- fitNNE$FP[HighIndex] - fitNNE$FP[LowIndex] delta.y <- fitNNE$TP[HighIndex] - fitNNE$TP[LowIndex] delta.2.eval <- eval.x - fitNNE$FP[LowIndex] target.y <- fitNNE$TP[LowIndex] if( delta.x > 0 ) target.y <- target.y + (delta.2.eval/delta.x)*delta.y TPout[ j ] <- target.y } SIM.TP = cbind(SIM.TP, TPout) } temp=apply(SIM.TP, 1, mean) All.mean.TP = cbind(All.mean.TP, temp) } ttt <- TargetTP out <- cbind( ttt, All.mean.TP[,1:NROW(sampsi)]-ttt ) print( round( out, 4 ) ) ## ----------------------------------------------------## TP.simu2=c(0,All.mean.TP[,1],1) TP.simu4=c(0,All.mean.TP[,2],1) TP.simu8=c(0,All.mean.TP[,3],1) TP.simu16=c(0,All.mean.TP[,4],1) TP.simu32=c(0,All.mean.TP[,5],1) TP.true=TP.true FP=c(0:100)/100 ## PART 3 ## ROC curves par(mfrow=c(2,3)) ## fig 1, n=200 plot(FP, TP.true, type="l", xlim=c(0,1), ylim=c(0,1), xlab="FP", ylab="TP", main="n = 200") lines(FP, TP.simu2, col="red") ## fig 2, n=400 plot(FP, TP.true, type="l", xlim=c(0,1), ylim=c(0,1), xlab="FP", ylab="TP", main="n = 400") lines(FP, TP.simu4, col="orange") ## fig 3, n=800 plot(FP, TP.true, type="l", xlim=c(0,1), ylim=c(0,1), xlab="FP", ylab="TP", main="n = 800") lines(FP, TP.simu8, col="magenta") ## fig 4, n=1600 plot(FP, TP.true, type="l", xlim=c(0,1), ylim=c(0,1), xlab="FP", ylab="TP", main="n = 1600") lines(FP, TP.simu16, col="navy") ## fig 5, n=3200 plot(FP, TP.true, type="l", xlim=c(0,1), ylim=c(0,1), xlab="FP", ylab="TP", main="n = 3200") lines(FP, TP.simu32, col="green") ## Legend in place of fig 6 plot(FP, TP.true, type="n", xlab="", ylab="", xaxt="n", yaxt="n" ) legend(0.3,0.7, legend=c("True ROC", "Simulated ROC", " n = 200", " n = 400", " n = 800", " n = 1600", " n = 3200"), col=c("black", "white", "red", "orange", "magenta", "navy", "green"), bty="n", lty=1) dev.off() ## end-of-file...