R 3.3.1 --------------------------------------------- Program: KMCox.R Folder: /home/bendix/teach/AdvCoh/talks/Mainz2016-09 Started: Thursday 08. September 2016, 12:14:18 --------------------------------------------- > ### R code from vignette source 'KMCox.rnw' > ### Encoding: UTF-8 > > ################################################### > ### code chunk number 1: KMCox.rnw:201-251 > ################################################### > library( survival ) > library( Epi ) Attaching package: ‘Epi’ The following object is masked from ‘package:base’: merge.data.frame > Lung <- Lexis( exit = list( tfe=time ), + exit.status = factor(status,labels=c("Alive","Dead")), + data = lung ) NOTE: entry.status has been set to "Alive" for all. NOTE: entry is assumed to be 0 on the tfe timescale. > mL.cox <- coxph( Surv( tfe, tfe+lex.dur, lex.Xst=="Dead" ) ~ + age + factor( sex ), + method="breslow", eps=10^-8, iter.max=25, data=Lung ) > Lung.s <- splitLexis( Lung, + breaks=c(0,sort(unique(Lung$time))), + time.scale="tfe" ) > Lung.S <- splitLexis( Lung, + breaks=c(0,sort(unique(Lung$time[Lung$lex.Xst=="Dead"]))), + time.scale="tfe" ) > summary( Lung.s ) Transitions: To From Alive Dead Records: Events: Risk time: Persons: Alive 19857 165 20022 165 69593 228 > summary( Lung.S ) Transitions: To From Alive Dead Records: Events: Risk time: Persons: Alive 15916 165 16081 165 69593 228 > subset( Lung.s, lex.id==96 )[,1:11] lex.id tfe lex.dur lex.Cst lex.Xst inst time status age sex ph.ecog 9235 96 0 5 Alive Alive 12 30 2 72 1 2 9236 96 5 6 Alive Alive 12 30 2 72 1 2 9237 96 11 1 Alive Alive 12 30 2 72 1 2 9238 96 12 1 Alive Alive 12 30 2 72 1 2 9239 96 13 2 Alive Alive 12 30 2 72 1 2 9240 96 15 11 Alive Alive 12 30 2 72 1 2 9241 96 26 4 Alive Dead 12 30 2 72 1 2 > nlevels( factor( Lung.s$tfe ) ) [1] 186 > system.time( + mLs.pois.fc <- glm( lex.Xst=="Dead" ~ - 1 + factor( tfe ) + + age + factor( sex ), + offset = log(lex.dur), + family=poisson, data=Lung.s, eps=10^-8, maxit=25 ) + ) user system elapsed 10.931 0.011 10.942 > length( coef(mLs.pois.fc) ) [1] 188 > system.time( + mLS.pois.fc <- glm( lex.Xst=="Dead" ~ - 1 + factor( tfe ) + + age + factor( sex ), + offset = log(lex.dur), + family=poisson, data=Lung.S, eps=10^-8, maxit=25 ) + ) user system elapsed 3.385 0.012 3.398 > length( coef(mLS.pois.fc) ) [1] 142 > > t.kn <- c(0,25,100,500,1000) > dim( Ns(Lung.s$tfe,knots=t.kn) ) [1] 20022 4 > system.time( + mLs.pois.sp <- glm( lex.Xst=="Dead" ~ Ns( tfe, knots=t.kn ) + + age + factor( sex ), + offset = log(lex.dur), + family=poisson, data=Lung.s, eps=10^-8, maxit=25 ) + ) user system elapsed 0.179 0.000 0.179 > ests <- + rbind( ci.exp(mL.cox), + ci.exp(mLs.pois.fc,subset=c("age","sex")), + ci.exp(mLS.pois.fc,subset=c("age","sex")), + ci.exp(mLs.pois.sp,subset=c("age","sex")) ) > cmp <- cbind( ests[c(1,3,5,7) ,], + ests[c(1,3,5,7)+1,] ) > rownames( cmp ) <- c("Cox","Poisson-factor","Poisson-factor (D)","Poisson-spline") > colnames( cmp )[c(1,4)] <- c("age","sex") > > > ################################################### > ### code chunk number 2: KMCox.rnw:253-254 > ################################################### > round( cmp, 7 ) age 2.5% 97.5% sex 2.5% 97.5% Cox 1.017158 0.9989388 1.035710 0.5989574 0.4313720 0.8316487 Poisson-factor 1.017158 0.9989388 1.035710 0.5989574 0.4313720 0.8316487 Poisson-factor (D) 1.017332 0.9991211 1.035874 0.5984794 0.4310150 0.8310094 Poisson-spline 1.016189 0.9980329 1.034676 0.5998287 0.4319932 0.8328707 > > > ################################################### > ### code chunk number 3: KMCox.rnw:269-273 (eval = FALSE) > ################################################### > ## mLs.pois.sp <- glm( lex.Xst=="Dead" ~ Ns( tfe, knots=t.kn ) + > ## age + factor( sex ), > ## offset = log(lex.dur), > ## family=poisson, data=Lung.s, eps=10^-8, maxit=25 ) > > > ################################################### > ### code chunk number 4: KMCox.rnw:275-279 > ################################################### > CM <- cbind( 1, Ns( seq(10,1000,10)-5, knots=t.kn ), 60, 1 ) > lambda <- ci.exp( mLs.pois.sp, ctr.mat=CM ) > Lambda <- ci.cum( mLs.pois.sp, ctr.mat=CM, intl=10 )[,-4] > survP <- exp(-rbind(0,Lambda)) > > > ################################################### > ### code chunk number 5: KMCox.rnw:281-358 > ################################################### > sf <- survfit( mL.cox, + newdata=data.frame( sex=factor(2,levels=1:2), + age=c(60) ) ) > > length( coef(mLS.pois.fc) ) [1] 142 > br <- ci.exp( mLS.pois.fc, ctr.mat=cbind(diag(139),0,60,1) )*365.25 > bt <- as.numeric( gsub( "factor\\(tfe)", "", names(coef(mLS.pois.fc))[1:139] ) ) > > pdf( "./graph/KMCox1.pdf", width=8, height=6 ) > par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, oma=c(0,0,0,0), + las=1, bty="n" ) > plot( 1, 1, type="n", xaxt="n", yaxt="n",xlab="",ylab="" ) > plot( sf, lwd=c(4,1,1), col="red", conf.int=T, mark.time=F, + xlab="Days since diagnosis", + ylab="Survival", xlim=c(0,900), yaxs="i", lty=1) > dev.off() null device 1 > > pdf( "./graph/KMCox2.pdf", width=8, height=6 ) > par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, oma=c(0,0,0,0), + las=1, bty="n" ) > matplot( 1:100*10-5, lambda * 365.25, + type="l", lwd=c(4,1,1), lty=1, col="black", log="y", + xlim=c(0,900), xaxs="i", ylim=c(1/10,10), + xlab="Days since diagnosis", + ylab="Mortality rate per year") > plot( sf, lwd=c(4,1,1), col="red", conf.int=T, mark.time=F, + xlab="Days since diagnosis", + ylab="Survival", xlim=c(0,900), yaxs="i", lty=1) > dev.off() null device 1 > > pdf( "./graph/KMCox3.pdf", width=8, height=6 ) > par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, oma=c(0,0,0,0), + las=1, bty="n" ) > matplot( 1:100*10-5, lambda * 365.25, + type="l", lwd=c(4,1,1), lty=1, col="black", log="y", + xlim=c(0,900), xaxs="i", ylim=c(1/10,10), + xlab="Days since diagnosis", + ylab="Mortality rate per year") > plot( sf, lwd=c(4,1,1), col="red", conf.int=T, mark.time=F, + xlab="Days since diagnosis", + ylab="Survival", xlim=c(0,900), yaxs="i", lty=1) > matlines( 0:100*10, survP, lwd=c(4,1,1), col="black", lty=1 ) > dev.off() null device 1 > > pdf( "./graph/KMCox4.pdf", width=8, height=6 ) > par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, oma=c(0,0,0,0), + las=1, bty="n" ) > matplot( 1:100*10-5, lambda * 365.25, + type="l", lwd=c(4,1,1), lty=1, col="black", log="y", + xlim=c(0,900), xaxs="i", ylim=c(1/10,10), + xlab="Days since diagnosis", + ylab="Mortality rate per year") > lines( bt, br[,1], col="red", lwd=2, type="s" ) > abline( h=0.1, lwd=2, col="white") > plot( sf, lwd=c(4,1,1), col="red", conf.int=T, mark.time=F, + xlab="Days since diagnosis", + ylab="Survival", xlim=c(0,900), yaxs="i", lty=1) > matlines( 0:100*10, survP, lwd=c(4,1,1), col="black", lty=1 ) > dev.off() null device 1 > > pdf( "./graph/KMCox5.pdf", width=8, height=6 ) > par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, oma=c(0,0,0,0), + las=1, bty="n" ) > matplot( 1:100*10-5, lambda * 365.25, + type="l", lwd=c(4,1,1), lty=1, col="black", log="y", + xlim=c(0,900), xaxs="i", ylim=c(1/10,10), + xlab="Days since diagnosis", + ylab="Mortality rate per year") > lines( bt, br[,1], col="red", lwd=2, type="s" ) > abline( h=0.1, lwd=2, col="white") > plot( sf, lwd=c(4,1,1), col="red", conf.int=T, mark.time=F, + xlab="Days since diagnosis", + ylab="Survival", xlim=c(0,900), yaxs="i", lty=1) > matlines( 0:100*10, survP, lwd=c(4,1,1), col="black", lty=1 ) > lines( sf, lwd=c(4,1,1), col="red", conf.int=T, mark.time=F ) > dev.off() null device 1 > > > > --------------------------------------------- Program: KMCox.R Folder: /home/bendix/teach/AdvCoh/talks/Mainz2016-09 Ended: Thursday 08. September 2016, 12:14:34 Elapsed: 00:00:16 --------------------------------------------- > proc.time() user system elapsed 16.171 0.082 16.245