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

###################################################
### code chunk number 1: lung-sol.rnw:2-6
###################################################
options( width=90,
         prompt=" ", continue=" ",
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )


###################################################
### code chunk number 2: lung-sol.rnw:22-26
###################################################
memory.size( 3000 )
library( Epi )
library( survival )
sessionInfo()


###################################################
### code chunk number 3: lung-sol.rnw:33-36
###################################################
data( lung )
str( lung )
lung[1:10,]


###################################################
### code chunk number 4: lung-sol.rnw:42-43
###################################################
table( lung$status )


###################################################
### code chunk number 5: lung-sol.rnw:48-49
###################################################
addmargins( table( table( lung$time ) ) )


###################################################
### code chunk number 6: lung-sol.rnw:62-67
###################################################
system.time(
m0.cox <- coxph( Surv( time, status==2 ) ~ age + factor( sex ),
                method="breslow", eps=10^-8, iter.max=25, data=lung )
            )
summary( m0.cox )


###################################################
### code chunk number 7: lung-sol.rnw:71-75
###################################################
Lung <- Lexis( exit = list( tfe=time ),
               exit.status = factor(status,labels=c("Alive","Dead")),
               data = lung )
summary( Lung )


###################################################
### code chunk number 8: lung-sol.rnw:81-85
###################################################
mL.cox <- coxph( Surv( tfe, tfe+lex.dur, lex.Xst=="Dead" ) ~
                 age + factor( sex ),
                 method="breslow", eps=10^-8, iter.max=25, data=Lung )
cbind( coef(m0.cox), coef(mL.cox) )


###################################################
### code chunk number 9: lung-sol.rnw:97-102
###################################################
Lung.s <- splitLexis( Lung,
                      breaks=c(0,sort(unique(Lung$time))),
                      time.scale="tfe" )
summary( Lung.s )
subset( Lung.s, lex.id==96 )


###################################################
### code chunk number 10: lung-sol.rnw:106-111
###################################################
system.time(
mLs.cox <- coxph( Surv( tfe, tfe+lex.dur, lex.Xst=="Dead" ) ~
                  age + factor( sex ),
                  method="breslow", eps=10^-8, iter.max=25, data=Lung.s )
            )


###################################################
### code chunk number 11: lung-sol.rnw:114-115
###################################################
cbind( coef(m0.cox), coef(mL.cox), coef(mLs.cox) )


###################################################
### code chunk number 12: lung-sol.rnw:121-122
###################################################
nlevels( factor( Lung.s$tfe ) )


###################################################
### code chunk number 13: lung-sol.rnw:128-135
###################################################
system.time(
mLs.pois.fc <- glm( lex.Xst=="Dead" ~ factor( tfe ) +
                              age + factor( sex ),
                              offset = log(lex.dur),
                    family=poisson, data=Lung.s, eps=10^-8, maxit=25 )
            )
length( coef(mLs.pois.fc) )


###################################################
### code chunk number 14: lung-sol.rnw:139-140
###################################################
rbind( coef( m0.cox), coef( mLs.pois.fc )[188-1:0] )


###################################################
### code chunk number 15: lung-sol.rnw:155-157
###################################################
t.kn <- c(0,25,100,500,1000)
dim( Ns(Lung.s$tfe,knots=t.kn) )


###################################################
### code chunk number 16: lung-sol.rnw:161-169
###################################################
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 )
            )
ci.exp( mLs.pois.sp )
ci.exp( mLs.pois.sp, subset=c("age","sex") )


###################################################
### code chunk number 17: lung-sol.rnw:180-189
###################################################
ests <-
rbind( ci.exp(m0.cox),
       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)  ,],
              ests[c(1,3,5)+1,] )
rownames( cmp ) <- c("Cox","Poisson-factor","Poisson-spline")
colnames( cmp )[c(1,4)] <- c("age","sex")
round( cmp, 5 )


###################################################
### code chunk number 18: lung-sol.rnw:192-196
###################################################
round(
rbind( ci.lin(m0.cox)[,2],
       ci.lin(mLs.pois.fc,subset=c("age","sex"))[,2],
       ci.lin(mLs.pois.sp,subset=c("age","sex"))[,2] ), 6 )


###################################################
### code chunk number 19: lung-sol.rnw:213-214
###################################################
CM <- cbind( 1, Ns( seq(0,1000,10), knots=t.kn ), 60, 1 )


###################################################
### code chunk number 20: lung-sol.rnw:217-218
###################################################
lambda <- ci.exp( mLs.pois.sp, ctr.mat=CM )


###################################################
### code chunk number 21: lung-sol.rnw:224-226
###################################################
Lambda <- ci.cum( mLs.pois.sp, ctr.mat=CM, intl=10 )
Lambda <- rbind( 0, Lambda )


###################################################
### code chunk number 22: lung-sol.rnw:231-234
###################################################
sf <-  survfit( m0.cox,
                newdata=data.frame( sex=factor(2,levels=1:2),
                                    age=c(60) ) )


###################################################
### code chunk number 23: haz-surv
###################################################
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( 0:100*10, 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,5),
         xlab="Days since diagnosis",
         ylab="Mortality rate per year")

# Then the survival curves by the two methods
# Here is the Breslow-estimator; note
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), xaxs="i", lty=1)
matlines( 0:101*10, exp(-Lambda[,1:3]), lwd=c(4,1,1), col="black", lty=1 )