### R code from vignette source 'KMCox.rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: KMCox.rnw:201-251
###################################################
library( survival ) 
library( Epi )
Lung <- Lexis( exit = list( tfe=time ),
               exit.status = factor(status,labels=c("Alive","Dead")),
               data = lung )
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 )
summary( Lung.S )
subset( Lung.s, lex.id==96 )[,1:11]
nlevels( factor( Lung.s$tfe ) )
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 )
            )
length( coef(mLs.pois.fc) )
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 )
            )
length( coef(mLS.pois.fc) )

t.kn <- c(0,25,100,500,1000)
dim( Ns(Lung.s$tfe,knots=t.kn) )
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 )
            )
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 )


###################################################
### 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) )
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()

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()

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()

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()

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()