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

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


###################################################
### code chunk number 2: stroke1-s.rnw:38-43
###################################################
stroke <- read.table( url("http://BendixCarstensen.com/AdvCoh/Scot-2014/data/stroke.csv"),
                      sep=";", header=TRUE, na.strings="." )
stroke$id <- 1:nrow(stroke)
str( stroke )
head( stroke )


###################################################
### code chunk number 3: stroke1-s.rnw:50-54
###################################################
stroke <- transform( stroke, dstr=as.Date(dstr,format="%d.%m.%Y"),
                             died=as.Date(died,format="%d.%m.%Y") )
str( stroke )
stroke <- cal.yr(stroke)


###################################################
### code chunk number 4: stroke1-s.rnw:68-72
###################################################
stroke <- transform( stroke, dox = pmin( died, 1996, na.rm=TRUE ) )
subset( stroke, died>1996 )
with( stroke, table( died>1996 ) )
stroke <- transform( stroke, D = ( dox < 1996 ) )


###################################################
### code chunk number 5: KM
###################################################
library( survival )
sst <- with( stroke, Surv( dox-dstr, D ) ~ 1 )
survfit( sst )
plot( survfit( sst ) )


###################################################
### code chunk number 6: KM0
###################################################
plot( survfit( sst ) )
sst0 <- with( subset(stroke,dox>dstr), Surv( dox-dstr, D ) ~ 1 )
lines( survfit( sst0 ), col="red" )


###################################################
### code chunk number 7: stroke1-s.rnw:110-111
###################################################
stroke <- subset( stroke, dox>dstr )


###################################################
### code chunk number 8: stroke1-s.rnw:117-119
###################################################
with( stroke, table( dgn, D ) )
( sdiag <- survfit(  Surv( dox-dstr, D ) ~ dgn, data=stroke ) )


###################################################
### code chunk number 9: KM-dgn
###################################################
plot( sdiag, col=1:4, lwd=3, mark.time=F )
legend( "bottomleft", legend=levels(stroke$dgn),
         col=1:4, lwd=3, bty="n", text.col=1:4 )


###################################################
### code chunk number 10: lch
###################################################
par( mfrow=c(1,3), mar=c(3,3,1,1) )
plot( survfit( Surv(dox-dstr,D) ~ dgn , data=stroke ), col=1:4, fun="cloglog",
      xlim=c(0.002,5.5),ylim=c(-3.5,0.5),lwd=2)
legend("bottomright", legend=levels(stroke$dgn), col=1:4, lty=1, lwd=3, bty="n" )
plot( survfit( Surv(dox-dstr,D) ~ diab, data=stroke ), col=1:2, fun="cloglog",
      xlim=c(0.002,5.5),ylim=c(-3.5,0.5),lwd=2)
legend("bottomright", legend=levels(factor(stroke$diab)),col=1:2,lty=1,lwd=3,bty="n" )
plot( survfit( Surv(dox-dstr,D) ~ sex, data=stroke), col=c("red","blue"), fun="cloglog",
      xlim=c(0.002,5.5), ylim=c(-3.5,0.5), lwd=2 )
legend("bottomright", legend=c("F","M"), col=c("red","blue"), lty=1, lwd=3, bty="n" )


###################################################
### code chunk number 11: KM-sex
###################################################
plot(survfit( Surv(dox-dstr,D) ~ sex, data=stroke),
              col=c("red","blue") )
survdiff( Surv(dox-dstr,D) ~ sex, data=stroke)


###################################################
### code chunk number 12: stroke1-s.rnw:193-197
###################################################
Lst <- Lexis( data=stroke, entry=list(Per=dstr,Age=age,Tfs=dstr-dstr),
                            exit=list(Per=dox),
                     exit.status=as.numeric(stroke$D) )
head( Lst )


###################################################
### code chunk number 13: stroke1-s.rnw:203-204
###################################################
summary( Lst )


###################################################
### code chunk number 14: Lexis
###################################################
plot( Lst )


###################################################
### code chunk number 15: stroke1-s.rnw:226-233
###################################################
pdf( "../graph/stroke1-LexisX.pdf", height=10+1, width=3+1 )
par( mai=c(3,3,1,1)/4, mgp=c(3,1,0)/1.6 )
plot(Lst,xlim=1980+c(0,30),ylim=c(0,100),
         col=c("red","blue")[Lst$sex+1],grid=0:20*5,xaxs="i",yaxs="i")
points( subset(Lst,lex.Xst==TRUE),pch=16,cex=0.6,
        col=c("red","blue")[Lst$sex+1][Lst$lex.Xst==TRUE])
dev.off()


###################################################
### code chunk number 16: stroke1-s.rnw:255-257
###################################################
with( stroke, survfit( Surv( dox-dstr, D ) ~ sex ) )
with( Lst, survfit( Surv( lex.dur, lex.Xst ) ~ sex ) )


###################################################
### code chunk number 17: stroke1-s.rnw:265-266
###################################################
save( stroke, Lst, file="../data/from-exc-stroke1.Rdata" )