### R code from vignette source 'surv'
### Encoding: UTF-8

###################################################
### code chunk number 1: surv.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: surv.rnw:116-121
###################################################
library( Epi )
library( popEpi )
library( survival )
# sessionInfo()
data( lung )


###################################################
### code chunk number 3: surv.rnw:125-135
###################################################
lt <- with( lung, addmargins( 
                  table( cut(time/30,seq(0,36,3)), 
                         factor(status,labels=c("Alive","Dead"))) ) )
head( lt )
last <- nrow(lt)
lt[,3] <- nrow(lung) - c(0,cumsum(lt[1:(last-1),3]))
N <- lt[-last,3]
D <- lt[-last,2]
L <- lt[-last,1]
cbind(N,D,L)


###################################################
### code chunk number 4: ltab
###################################################
pi <- D/(N-L/2)
nn <- length(D)
aS <- cumprod( c(1,1-pi[-nn]) )
lm <- D/(N-L/2-D/2)
mS <- cumprod( c(1,exp(-lm[-nn])) )
round( cbind( N, D, L, pi, aS, mS ), 3 )


###################################################
### code chunk number 5: KM1
###################################################
plot( km0 <- survfit( Surv(time,status) ~ 1, data=lung ), yaxs="i" )
km0
medsurv <- summary(km0)$table[7:9]
abline( h=0.5 )
points( medsurv, rep(0.5,3) )
# add the actuarial estimators from above
matlines( seq(0,33,3)*30, cbind(aS,mS), 
          type="l", col=c("gray","black"),
          lty=c("solid","11"), lend="butt", lwd=5, ylim=c(0,1) )
lines( km0, lwd=2, conf.int=FALSE )


###################################################
### code chunk number 6: KMsex
###################################################
lung <- transform( lung, sex = factor(sex,labels=c("M","W")) )
plot( svsx <- survfit( Surv(time,status) ~ sex, data = lung ),
      col=c("blue","red"), conf.int=TRUE, yaxs="i" )
lines( svsx, col=c("blue","red"), lwd=2 )
svsx


###################################################
### code chunk number 7: surv.rnw:332-333
###################################################
survdiff( Surv(time,status) ~ sex, data = lung )


###################################################
### code chunk number 8: surv.rnw:383-385
###################################################
c0 <- coxph( Surv(time,status) ~ sex + age, data = lung )
summary(c0)


###################################################
### code chunk number 9: coxph-pred
###################################################
plot( survfit( c0, newdata=data.frame( sex=c("M","W"), age=60 ) ),   
               col=c("blue","red") )


###################################################
### code chunk number 10: cph-non
###################################################
plot( prm0 <- survfit( c0 ) )


###################################################
### code chunk number 11: surv.rnw:444-445
###################################################
with( lung, c( mean(as.numeric(sex)), mean(age) ) )


###################################################
### code chunk number 12: surv.rnw:470-473
###################################################
nd.lung <- lung[,c("sex","age")]
prlu <- survfit( c0, newdata=nd.lung )
lusurv <- apply( prlu$surv, MARGIN=1, mean )


###################################################
### code chunk number 13: smms
###################################################
plot( prm0, conf.int=FALSE, col=gray(0.5), lwd=2 )
lines( prlu$time, lusurv, type="s", lwd=2 )

xnd.lung <- data.frame( sex=lung$sex, age=rep(c(40,90),c(128,100)) )
( mx <- mean(xnd.lung$age) )
( m0 <- mean(lung$age) )
xnd.lung$age <- xnd.lung$age - mx + m0
mean( xnd.lung$age )
xprlu <- survfit( c0, newdata=xnd.lung )
xlusurv <- apply( xprlu$surv, 1, mean )
lines( xprlu$time, xlusurv, type="s", lty="12", lwd=3 )


###################################################
### code chunk number 14: surv.rnw:599-604
###################################################
library( Epi )
library( popEpi )
data( DMlate )
head( DMlate )
DMlate <- subset( DMlate, dox>dodm )


###################################################
### code chunk number 15: surv.rnw:613-620
###################################################
aDM <- coxph( Surv(dodm-dobth,dox-dobth,!is.na(dodth)) ~ sex, 
              data=DMlate )
dDM <- coxph( Surv(           dox-dodm ,!is.na(dodth)) ~ sex, 
              data=DMlate )
pDM <- coxph( Surv(dodm      ,dox      ,!is.na(dodth)) ~ sex, 
              data=DMlate )
round( rbind( ci.exp(aDM), ci.exp(dDM), ci.exp(pDM) ), 3 )


###################################################
### code chunk number 16: tscmp
###################################################
par( mfrow=c(1,3) )
plot( survfit( aDM, newdata=data.frame(sex=c("M","F")) ),
      col=c("blue","red"), xlab="Current age", lwd=2 )
plot( survfit( dDM, newdata=data.frame(sex=c("M","F")) ),
      col=c("blue","red"), xlab="DM duration", lwd=2 )
plot( survfit( pDM, newdata=data.frame(sex=c("M","F")) ),
      col=c("blue","red"), xlab="Calendar time", lwd=2, 
      xlim=c(1995,2010) )


###################################################
### code chunk number 17: surv.rnw:728-733
###################################################
lung$sex <- factor( lung$sex, labels=c("M","W") )
c0 <- coxph( Surv(time,status) ~ sex + age, data = lung )
p0 <- glm( cbind(status==2,time) ~ sex + age, 
           family=poisreg, data = lung )
round( cbind( ci.exp(c0), ci.exp(p0)[-1,] ), 3 )


###################################################
### code chunk number 18: surv.rnw:753-757
###################################################
c0 <- coxph( Surv(time,status) ~ sex + I((age-60)/10), data = lung )
p0 <- glm( cbind(status==2,time/365.25) ~ sex + I((age-60)/10), 
           family=poisreg, data = lung )
round( cbind( rbind(NA,ci.exp(c0)), ci.exp(p0) ), 3 )


###################################################
### code chunk number 19: surv.rnw:800-806
###################################################
Ll <- Lexis( exit = list(tfd=time),
      exit.status = factor(status,labels=c("Alive","Dead")),
             data = lung )
Sl <- splitMulti( Ll, tfd=seq(0,1100,20) )
summary( Sl )
Sl[1:6+13,1:9]


###################################################
### code chunk number 20: surv.rnw:828-830
###################################################
ps <- gam.Lexis( Sl, formula = ~ s(tfd) + sex + I((age-60)/10) )
round( cbind( ci.exp(c0), ci.exp(ps,subset=c("sex","age"))), 3 )


###################################################
### code chunk number 21: cp-surv
###################################################
nd <- data.frame( tfd=seq(0,1000,10)+5, sex="W", age=60 )
sv.w <- ci.surv( ps, nd, int=10 )
plot( survfit(c0,newdata=data.frame(sex="W",age=60)) )
matshade( nd$tfd-5, sv.w, col="red", lwd=2 )


###################################################
### code chunk number 22: cp-surv2
###################################################
plot( survfit(c0,newdata=data.frame(sex=c("W","M"),age=60)),
      conf.int=TRUE, col=c("red","blue") )
nd <- transform( nd,  sex="M" )
sv.m <- ci.surv( ps, nd, int=10 )
matshade( nd$tfd-5, cbind(sv.w,sv.m), col=c("red","blue"), lwd=2 )


###################################################
### code chunk number 23: surv.rnw:945-951
###################################################
lung$time <- lung$time + runif(nrow(lung),-1,1) 
table(table(lung$time))
Ll <- Lexis( exit = list(tfd=time),
      exit.status = factor(status,labels=c("Alive","Dead")),
             data = lung )
cX <- coxph( Surv(time,lex.Xst=="Dead") ~ sex + age, data=Ll )


###################################################
### code chunk number 24: surv.rnw:968-978
###################################################
Xl <- splitMulti( Ll, tfd=c(0,sort(unique(Ll$time))) )
summary( Xl )
system.time(
pX <- glm( cbind(lex.Xst=="Dead",lex.dur) ~ factor(tfd) + sex + age, 
           family=poisreg, data=Xl )
           )
# save( pX, file='pX.Rda')
# load( file='pX.Rda')
length( coef(pX) )
cbind( ci.exp(cX), ci.exp(pX,subset=c("sex","age")) )


###################################################
### code chunk number 25: surv.rnw:1040-1056
###################################################
Ldm <- Lexis( entry = list( tfd=0,
                            age=dodm-dobth,
                            per=dodm ),
               exit = list( age=dox-dobth ),
        exit.status = factor( !is.na(dodth), labels=c("Alive","Dead") ),
               data = DMlate )
Cdm <- cutLexis( data = Ldm, 
                  cut = Ldm$dooad, 
            timescale = "per", 
            new.state = "OAD", 
     precursor.states = "Alive" )
oo <- options( digits=5 )
subset( Ldm, lex.id %in% c(2,38) )[,c(1:7,10,12)]
subset( Cdm, lex.id %in% c(2,38) )[,c(1:7,10,12)]
options( oo )
summary( Cdm, time=TRUE )


###################################################
### code chunk number 26: surv.rnw:1098-1105
###################################################
ma <- coxph( Surv(tfd,tfd+lex.dur,lex.Xst=="Dead") ~ sex + lex.Cst, 
             data=Cdm )
md <- coxph( Surv(age,age+lex.dur,lex.Xst=="Dead") ~ sex + lex.Cst, 
             data=Cdm )
mp <- coxph( Surv(per,per+lex.dur,lex.Xst=="Dead") ~ sex + lex.Cst, 
             data=Cdm )
round( cbind(ci.exp(ma),ci.exp(md),ci.exp(mp)), 3 )


###################################################
### code chunk number 27: surv.rnw:1125-1127
###################################################
md <- coxph.Lexis( Cdm, formula = tfd ~ sex + lex.Cst )
summary( md )


###################################################
### code chunk number 28: tdep
###################################################
 plot( survfit( md, newdata=data.frame(sex="M",lex.Cst="Alive") ) )
lines( survfit( md, newdata=data.frame(sex="M",lex.Cst="OAD") ), col="red" )


###################################################
### code chunk number 29: box1
###################################################
Adm <- subset( Cdm, lex.Cst=="Alive" )
summary( Adm )
boxes( Adm, boxpos=TRUE, scale.R=100, show.BE=TRUE )


###################################################
### code chunk number 30: surv.rnw:1242-1248
###################################################
levels(Adm$lex.Xst)
m3 <- survfit( Surv(tfd,tfd+lex.dur,lex.Xst) ~ 1, 
               data=Adm, id=lex.id )
matplot( m3$time, m3$pstate, 
         type="s", lty=1, lwd=4,  
         col=c("ForestGreen","red","black") )


###################################################
### code chunk number 31: surv2
###################################################
par( mfrow=c(1,2) )
plot( m3, col=c("red","black","ForestGreen"), lwd=3, ylim=c(0,1) )
stackedCIF( m3, col=c("red","black"), lwd=3 )


###################################################
### code chunk number 32: surv3
###################################################
m2 <- survfit( Surv(tfd,tfd+lex.dur,lex.Xst=="OAD") ~ 1, data=Adm )
M2 <- survfit( Surv(tfd,tfd+lex.dur,lex.Xst=="Dead") ~ 1, data=Adm )
par( mfrow=c(1,2) )
mat2pol( m3$pstate, c(2,3,1), x=m3$time, col=c("red","black","transparent") )
  lines( m2$time, 1-m2$surv, lwd=3, col="red" )
mat2pol( m3$pstate, c(3,2,1), x=m3$time, col=c("black","red","transparent") )
  lines( M2$time, 1-M2$surv, lwd=3, col="black" )


###################################################
### code chunk number 33: surv.rnw:1371-1375
###################################################
Sdm <- splitMulti( Adm, tfd=seq(0,20,0.1) )
summary(Sdm)
gl <- gam.Lexis( Sdm, ~ s(tfd), from="Alive", to="OAD" )
gm <- gam.Lexis( Sdm, ~ s(tfd), from="Alive", to="Dead" )


###################################################
### code chunk number 34: surv.rnw:1378-1382
###################################################
int <- 0.01
nd <- data.frame( tfd=seq(int,15,int)-int/2 )
mrt <- ci.pred( gm, nd )[,1]
lam <- ci.pred( gl, nd )[,1]


###################################################
### code chunk number 35: smpl
###################################################
Lam <- cumsum( lam*int )         # cumulative incidence of OAD
Mrt <- cumsum( mrt*int )         # cumulative mortality
Srv <- exp( -( Mrt+Lam ) )           # survival in ALive
Rdr <- cumsum( lam * Srv * int ) # cumulative risk of OAD                   
Rdt <- cumsum( mrt * Srv * int ) # cumulative risk of Death 


###################################################
### code chunk number 36: comp
###################################################
plot( m3, col=c("red","black","ForestGreen"), lwd=1 )
matlines( nd$tfd, cbind( Srv, Rdr, Rdt ), 
          type="l", lty="21", lwd=3, lend="butt", 
          col=c("ForestGreen","red","black") )


###################################################
### code chunk number 37: surv.rnw:1437-1438 (eval = FALSE)
###################################################
## vignette("simLexis","Epi")