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

###################################################
### code chunk number 1: renal-s.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: renal-s.rnw:56-61
###################################################
library( Epi ) ; clear()
load( url("http://BendixCarstensen.com/AdvCoh/courses/Frias-2016/renal.Rda") )
# load( "renal.Rda" )
str( renal )
head( renal )


###################################################
### code chunk number 3: renal-s.rnw:75-83
###################################################
Lr <- Lexis( entry = list( per=doe,
                           age=doe-dob,
                           tfi=0 ),
              exit = list( per=dox ),
       exit.status = factor( event>0, labels=c("NRA","ESRD") ),
              data = renal )
str( Lr )
summary( Lr )


###################################################
### code chunk number 4: Lexis-ups
###################################################
plot( Lr, col="black", lwd=3 )


###################################################
### code chunk number 5: Lexis-fancy
###################################################
# pdf( "lexis-fancy.pdf", height=80/5+1, width=40/5+1 )
# x11( height=80/5+1, width=40/5+1 )
par( mai=c(3,3,1,1)/4, mgp=c(3,1,0)/1.6 )
plot( Lr, 1:2, col=c("blue","red")[Lr$sex], lwd=3, grid=0:20*5,
      xlab="Calendar time", ylab="Age",
      xlim=c(1970,2010), ylim=c(0,80), xaxs="i", yaxs="i", las=1 )
# dev.off()


###################################################
### code chunk number 6: renal-s.rnw:116-120
###################################################
library( survival )
mc <- coxph( Surv( lex.dur, lex.Xst=="ESRD" ) ~
             I(age/10) + sex, data=Lr )
summary( mc )


###################################################
### code chunk number 7: renal-s.rnw:144-149
###################################################
Lc <- cutLexis( Lr, cut = Lr$dor, # where to cut follow up
              timescale = "per",  # the timescale that "dor" refers to
              new.state = "Rem",  # name of the new state
       precursor.states = "NRA" ) # which states are less severe
summary( Lc )


###################################################
### code chunk number 8: renal-s.rnw:157-158 (eval = FALSE)
###################################################
## boxes( Lc )


###################################################
### code chunk number 9: boxes1
###################################################
boxes( Lc, boxpos=TRUE, scale.R=100, show.BE=TRUE )


###################################################
### code chunk number 10: Lexis-rem
###################################################
par( mai=c(3,3,1,1)/4, mgp=c(3,1,0)/1.6 )
plot( Lc, col=c("red","limegreen")[(Lc$lex.Cst=="Rem")+1],
      xlab="Calendar time", ylab="Age",
      lwd=3, grid=0:20*5, xlim=c(1970,2010), ylim=c(0,80), xaxs="i", yaxs="i", las=1 )
points( Lc, pch=c(NA,16)[(Lc$lex.Xst=="ESRD")+1],
            col=c("red","limegreen")[(Lc$lex.Cst=="Rem")+1])
points( Lc, pch=c(NA,1)[(Lc$lex.Xst=="ESRD")+1],
            col="black", lwd=2 )


###################################################
### code chunk number 11: renal-s.rnw:194-197
###################################################
m1 <- coxph( Surv( tfi, tfi+lex.dur, lex.Xst=="ESRD" ) ~
             sex + I((doe-dob-50)/10) + (lex.Cst=="Rem"), data=Lc )
summary( m1 )


###################################################
### code chunk number 12: renal-s.rnw:224-227
###################################################
sLc <- splitLexis( Lc, "tfi", breaks=seq(0,30,1/12) )
summary( Lc )
summary(sLc )


###################################################
### code chunk number 13: renal-s.rnw:241-248
###################################################
library( splines )
mp <- glm( lex.Xst=="ESRD" ~ ns( tfi, df=4 ) +
                  sex + I((doe-dob-40)/10) +  (lex.Cst=="Rem"),
           offset = log(lex.dur),
           family = poisson, 
             data = sLc )
summary( mp )


###################################################
### code chunk number 14: renal-s.rnw:257-265
###################################################
t.kn <- with( subset( sLc, lex.Xst=="ESRD"),
              quantile( tfi+lex.dur, 0:4/5 ) )
mp <- glm( lex.Xst=="ESRD" ~ Ns( tfi, knots=t.kn ) +
                             sex + I((doe-dob-40)/10) +  (lex.Cst=="Rem"),
           offset = log(lex.dur),
           family = poisson, 
             data = sLc )
summary( mp )


###################################################
### code chunk number 15: renal-s.rnw:270-273
###################################################
ci.lin( mp )
ci.exp( mp )
ci.exp( mp, subset=c("sex","dob","Cst"), pval=TRUE )


###################################################
### code chunk number 16: renal-s.rnw:276-279
###################################################
ci.exp( m1 )
ci.exp( mp, subset=c("sex","dob","Cst") )
ci.exp( mp, subset=c("sex","dob","Cst") ) / ci.exp( m1 )


###################################################
### code chunk number 17: renal-s.rnw:284-285
###################################################
termplot( mp, terms=1 )


###################################################
### code chunk number 18: renal-s.rnw:295-301
###################################################
mP <- glm( lex.Xst=="ESRD" ~ -1 + Ns( tfi, knots=t.kn, intercept=TRUE ) +
                  sex + I((doe-dob-40)/10) + (lex.Cst=="Rem"),
           offset = log(lex.dur/100),
           family = poisson, 
             data = sLc )
Termplot( mP, terms=1 )


###################################################
### code chunk number 19: renal-s.rnw:314-315
###################################################
summary( sLc )


###################################################
### code chunk number 20: renal-s.rnw:329-330 (eval = FALSE)
###################################################
## pmax( per-dor, 0, na.rm=TRUE )


###################################################
### code chunk number 21: time-since-rem
###################################################
sLc <- transform( sLc, tfr = pmax( (per-dor)/10, 0, na.rm=TRUE ) )
mPx <- glm( lex.Xst=="ESRD" ~ -1 + Ns( tfi, knots=t.kn, intercept=TRUE ) +
              sex + I((age-tfi-40)/10) + (lex.Cst=="Rem") + tfr,
            offset = log(lex.dur/100),
            family = poisson, 
              data = sLc )
round( ci.exp( mPx ), 3 )
Termplot( mPx, terms=1 )


###################################################
### code chunk number 22: rem-inc
###################################################
mr <- glm( lex.Xst=="Rem" ~ ns( tfi, knots=t.kn ) + sex,
           offset = log(lex.dur),
           family = poisson,
             data = subset( sLc, lex.Cst=="NRA" ) )
ci.exp( mr, pval=TRUE )


###################################################
### code chunk number 23: renal-s.rnw:400-414
###################################################
inL <- subset( sLc, select=1:11 )[NULL,]
str( inL )
timeScales(inL)
inL[1,"lex.id"] <- 1
inL[1,"per"] <- 2000
inL[1,"age"] <- 50
inL[1,"tfi"] <- 0
inL[1,"lex.Cst"] <- "NRA"
inL[1,"lex.Xst"] <- NA
inL[1,"lex.dur"] <- NA
inL[1,"sex"] <- "M"
inL[1,"doe"] <- 2000
inL[1,"dob"] <- 1950
inL


###################################################
### code chunk number 24: renal-s.rnw:422-425
###################################################
Tr <- list( "NRA" = list( "Rem"  = mr,
                          "ESRD" = mp ),
            "Rem" = list( "ESRD" = mp ) )


###################################################
### code chunk number 25: first-sim
###################################################
( iL <- simLexis( Tr, inL, N=10 ) )
summary( iL )


###################################################
### code chunk number 26: 10000-sim
###################################################
system.time(
sM <- simLexis( Tr, inL, N=10000 ) )
summary( sM )


###################################################
### code chunk number 27: nState
###################################################
nSt <- nState( sM, at=seq(0,10,0.1), from=50, time.scale="age" )
head( nSt )


###################################################
### code chunk number 28: pState
###################################################
pp <- pState( nSt, perm=1:3 )
head( pp )
tail( pp )


###################################################
### code chunk number 29: plot-pp
###################################################
plot( pp )


###################################################
### code chunk number 30: add-state
###################################################
xM <- transform( sM, lex.Xst = factor( ifelse( lex.Xst=="ESRD" & lex.Cst=="Rem",
                                               "ESRD(Rem)",
                                               as.character(lex.Xst) ),
                                       levels=c("NRA","Rem","ESRD(Rem)","ESRD") ),
                     lex.Cst = factor( as.character(lex.Cst),
                                       levels=c("NRA","Rem","ESRD(Rem)","ESRD") ) )
summary( sM )
summary( xM )
boxes( xM, boxpos=TRUE, show.BE=TRUE, scale.R=100 )


###################################################
### code chunk number 31: new-nState
###################################################
xSt <- nState( xM, at=seq(0,10,0.1), from=50, time.scale="age" )
xp <- pState( xSt, perm=1:4 )
head( xp )
plot( xp, col=rev(c("pink","limegreen","forestgreen","red")), xlab="Age" )
lines( as.numeric(rownames(xp)), xp[,"Rem"], lwd=4 )


###################################################
### code chunk number 32: redo
###################################################
Lc <- cutLexis( Lr, cut = Lr$dor, # where to cut follow up
              timescale = "per",  # the timescale that "dor" refers to
              new.state = "Rem",  # name of the new state
       precursor.states = "NRA",  # which states are less severe
              new.scale = "tfr",  # define a new timescale as time since Rem
           split.states = TRUE )  # subdivide non-precursor states
str( Lc )
# source("/home/bendix/stat/R/lib.src/Epi/pkg/R/summary.Lexis.r")
# summary( Lc, S=T, scale=100 )
summary( Lc )
boxes( Lc, boxpos=list(x=c(20,80,20,80),y=c(80,80,20,20)), 
           scale.R=100, show.BE=TRUE )
sLc <- splitLexis( Lc, "tfi", breaks=seq(0,30,1/12) )
summary(  Lc )
summary( sLc )
head( subset( sLc, lex.id==2 )[,1:8], 8 )
tail( subset( sLc, lex.id==2 )[,1:8], 3 )
( fl <- levels(Lc)[3:4] )
mp <- glm( lex.Xst %in% fl ~ ns( tfi, df=4 ) +
                             sex + I((age-tfi-40)/10) + (lex.Cst=="Rem"),
           offset = log(lex.dur/100),
           family = poisson, 
             data = sLc )
# the timescale tfr must be given some value for time before Rem
sLc$tfr <- pmax( 0, sLc$tfr, na.rm=TRUE )
head( subset( sLc, lex.id==2 )[,1:8], 8 )
mr <- glm( lex.Xst=="Rem" ~ ns( tfi, df=4 ) + sex,
           offset = log(lex.dur),
           family = poisson,
             data = subset( sLc, lex.Cst=="NRA" ) )
ci.exp( mr, pval=TRUE )
inL <- subset( sLc, select=1:10 )[NULL,]
str( inL )
timeScales(inL)
inL[1,"lex.id"] <- 1
inL[1,"per"] <- 2000
inL[1,"age"] <- 50
inL[1,"tfi"] <- 0
inL[1,"lex.Cst"] <- "NRA"
inL[1,"lex.Xst"] <- NA
inL[1,"lex.dur"] <- NA
inL[1,"sex"] <- "M"
inL
Tr <- list( "NRA" = list( "Rem"  = mr,
                          "ESRD" = mp ),
            "Rem" = list( "ESRD(Rem)" = mp ) )
( iL <- simLexis( Tr, inL, N=10 ) )
summary( iL )

system.time(
sM <- simLexis( Tr, inL, N=10000, t.range=25, n.int=251 ) )
summary( sM )

nSt <- nState( sM, at=seq(0,24,0.1), from=50, time.scale="age" )
head( nSt )
pp <- pState( nSt, perm=c(1,2,4,3) )
head( pp )
tail( pp )
plot( pp )

# Two colors and the corresponding pale ones for the dead states
clr <- c("limegreen","orange")
col2rgb(clr)
cl4 <- cbind(col2rgb(clr),col2rgb(clr)/2+255/2)[,c(1,2,4,3)]
cl4 <- rgb( t(cl4), max=255 )
# Nicer plot
plot( pp, col=cl4, xlab="Age" )
lines( as.numeric(rownames(pp)), pp[,2], lwd=2 )