### R code from vignette source 'stroke2-s.rnw'

###################################################
### code chunk number 1: stroke2-s.rnw:4-5
###################################################
library(Epi)


###################################################
### code chunk number 2: stroke2-s.rnw:9-11
###################################################
load( file="../data/from-exc-stroke1.Rdata" )
str( Lst )


###################################################
### code chunk number 3: stroke2-s.rnw:20-25
###################################################
  library( survival )
  mc <- coxph( Surv(dox-dstr,D) ~ sex, data=stroke )
  summary( mc )
  mL <- coxph( Surv(lex.dur,lex.Xst==1) ~ sex, data=Lst )
  summary( mL )


###################################################
### code chunk number 4: stroke2-s.rnw:32-34
###################################################
  mLa <- coxph( Surv(lex.dur,lex.Xst==1) ~ sex + age, data=Lst )
  summary( mLa )


###################################################
### code chunk number 5: stroke2-s.rnw:41-42
###################################################
plot( survfit( Surv(dox-dstr,as.numeric(D)) ~ interaction(sex,age<75), data=stroke ) )


###################################################
### code chunk number 6: 75-KM
###################################################
plot( survfit( Surv(lex.dur,lex.Xst==1) ~ interaction(sex,age<75), 
               data=Lst ),
               col=c("red","blue"), lwd=3 )


###################################################
### code chunk number 7: stroke2-s.rnw:61-62
###################################################
with( Lst, table( interaction(sex,age<75) ) )


###################################################
### code chunk number 8: stroke2-s.rnw:67-72
###################################################
length( unique(Lst$lex.dur[Lst$lex.Xst==1]) )
sLst <- splitLexis( Lst, breaks=seq(0,10,0.05), "Tfs" )
str( sLst )
summary( Lst )
summary( sLst )


###################################################
### code chunk number 9: stroke2-s.rnw:78-80
###################################################
subset( Lst, lex.id %in% 54:55 )
subset( sLst, lex.id %in% 54:55 )


###################################################
### code chunk number 10: stroke2-s.rnw:86-90
###################################################
mCs <- coxph( Surv(lex.dur,lex.Xst==1) ~ sex + age, data=Lst )
ci.lin( mLa )
mC <- coxph( Surv(Tfs,Tfs+lex.dur,lex.Xst==1) ~ sex + age, data=sLst )
ci.lin( mC )


###################################################
### code chunk number 11: stroke2-s.rnw:96-100
###################################################
system.time(
mP <- glm( lex.Xst ~ factor( Tfs ) + sex + age + offset(log(lex.dur)),
                     family=poisson, data=sLst )
           )


###################################################
### code chunk number 12: stroke2-s.rnw:104-105
###################################################
coef( mP )


###################################################
### code chunk number 13: stroke2-s.rnw:109-111
###################################################
ci.lin( mP, subset=c("sex","age"), Exp=TRUE )
ci.lin( mC, Exp=TRUE ) 


###################################################
### code chunk number 14: stroke2-s.rnw:125-132
###################################################
sLst$Tfs.m <- timeBand( sLst, "Tfs", "middle" )
library( splines )
kn <- c(0.05,0.2,0.7,1.5,3) 
Bk <- c(0,4.8) 
mS <- glm( lex.Xst ~ 
           ns( Tfs.m, knots=kn, Bo=Bk ) + sex + age + offset(log(lex.dur)), 
           family=poisson, data=sLst )


###################################################
### code chunk number 15: stroke2-s.rnw:135-138
###################################################
ci.lin( mC )
ci.lin( mP, subset=c("sex","age") )
ci.lin( mS, subset=c("sex","age") )


###################################################
### code chunk number 16: stroke2-s.rnw:144-148
###################################################
t.pt <- seq(0,5,0.01)
CM <- cbind( 1, ns( t.pt, knots=kn, Bo=Bk ), 0, 60 )
hz <- ci.lin( mS, ctr.mat=CM, Exp=TRUE )[,5:7] * 1000
matplot( t.pt, hz, type="l", lwd=c(3,1,1), ylim=c(0,1), lty=1, col="black" )


###################################################
### code chunk number 17: stroke2-s.rnw:154-158
###################################################
nd <- data.frame( Tfs.m=t.pt, sex=0, age=60, lex.dur=1000 )
prhz <- predict( mS, newdata=nd, type="link", se.fit=T )
str( prhz )
prhz <- exp( cbind( prhz$fit, prhz$se.fit ) %*% ci.mat() )


###################################################
### code chunk number 18: stroke2-s.rnw:161-162
###################################################
matplot( hz, prhz )


###################################################
### code chunk number 19: Cox-Pois
###################################################
t.pt <- t.pt[-1]
CM <- cbind( 1, ns( t.pt-0.005, knots=kn, Bo=Bk ), 0, 60 )
Hz <- ci.cum( mS, ctr.mat=CM, intl=0.01 )
matplot( t.pt, exp(-Hz)[,-4], 
         type="l", lwd=c(3,1,1), ylim=c(0,1), lty=1, col="black" )


###################################################
### code chunk number 20: Cox-Pois2
###################################################
matplot( t.pt, exp(-Hz)[,-4], 
         type="l", lwd=c(3,1,1), ylim=c(0,1), lty=1, col="black" )
lines( survfit(mC,newdata=data.frame(sex=0,age=60)), 
       conf.int=TRUE, col="red" )
# overplot the estimate with a thicker line:
lines( survfit(mC,newdata=data.frame(sex=0,age=60)), 
       conf.int=FALSE, col="red", lwd=3 )