### R code from vignette source 'cmpr-DMlate.rnw'

###################################################
### code chunk number 1: cmpr-DMlate.rnw:2-7
###################################################
options( width=90,
#        prompt=" ", continue=" ", # Absence of prompts makes it easier for
                                   # Students to copy from the final pdf document
         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: cmpr-DMlate.rnw:27-34
###################################################
library(survival)
library(Epi)
library(popEpi)
# popEpi::splitMulti returns a data.frame rather than a data.table
options("popEpi.datatable" = FALSE)
library(tidyverse)
clear()


###################################################
### code chunk number 3: cmpr-DMlate.rnw:62-67
###################################################
data(DMlate)
# set.seed(1952)
# DMlate <- DMlate[sample(1:nrow(DMlate), 2000),]
str(DMlate)
head(DMlate)


###################################################
### code chunk number 4: DMhists
###################################################
par(mfrow=c(2,2))
hist(DMlate$dobth, breaks = seq(1898, 2010, 1   ), col = "black")
hist(DMlate$dodm , breaks = seq(1995, 2010, 1/12), col = "black")
    abline(v = 1995:2020, col = "red")
hist(DMlate$dodth, breaks = seq(1995, 2010, 1/12), col = "black")
    abline(v = 1995:2020, col = "red")
hist(DMlate$doins, breaks = seq(1995, 2010, 1/12), col = "black")
    abline(v = 1995:2020, col = "red")


###################################################
### code chunk number 5: cmpr-DMlate.rnw:89-97
###################################################
Ldm <- Lexis(entry = list(per = dodm,
                          age = dodm - dobth,
                          tfd = 0),
              exit = list(per = dox),
       exit.status = factor(!is.na(dodth),
                            labels = c("DM", "Dead")),
              data = DMlate)
summary(Ldm)


###################################################
### code chunk number 6: cmpr-DMlate.rnw:100-106
###################################################
Ldm <- sortLexis(Ldm)
Cdm <- cutLexis(Ldm,
                cut = Ldm$doins,
          timescale = "per",
          new.state = "Ins")
summary(Cdm)


###################################################
### code chunk number 7: boxCR
###################################################
Adm <- filter(Cdm, lex.Cst == "DM")
Adm <- subset(Cdm, lex.Cst == "DM")
summary(Adm)
par(mfrow=c(1,2))
boxes(Cdm, boxpos = TRUE, scale.R = 100, show.BE = TRUE)
boxes(Adm, boxpos = TRUE, scale.R = 100, show.BE = TRUE)


###################################################
### code chunk number 8: cmpr-DMlate.rnw:149-158
###################################################
levels(Adm$lex.Xst)
m3 <- survfit(Surv(tfd,
                   tfd + lex.dur,
                   lex.Xst) ~ 1,
                id = lex.id,
              data = Adm)
names(m3)
m3$states
head(cbind(time = m3$time, m3$pstate), 20)


###################################################
### code chunk number 9: cmpr-DMlate.rnw:169-171 (eval = FALSE)
###################################################
## m3$transitions
## summary(Adm)


###################################################
### code chunk number 10: surv2 (eval = FALSE)
###################################################
## par(mfrow=c(1, 2))
## matplot(m3$time, m3$pstate,
##         type = "s", lty = 1, lwd = 4,
##         col = c("ForestGreen", "red", "black"),
##         xlim = c(0, 15), xaxs = "i",
##         ylim = c(0, 1), yaxs = "i" )
## stackedCIF(m3, lwd = 3, xlim = c(0, 15), xaxs = "i", yaxs = "i" )
## text(rep(12, 3), c(0.85, 0.45, 0.15), levels(Cdm))
## box()


###################################################
### code chunk number 11: surv2
###################################################
par(mfrow = c(1, 2))
matshade(m3$time, cbind(m3$pstate,
                        m3$lower,
                        m3$upper)[, c(1, 4, 7, 2, 5, 8, 3, 6, 9)],
         plot = TRUE, lty = 1, lwd = 2,
         col = clr <- c("ForestGreen","red","black"),
         xlim=c(0,15), xaxs="i",
         ylim = c(0,1), yaxs = "i")
mat2pol(m3$pstate, perm = 3:1, x = m3$time, col = clr[3:1])
text(rep(12, 3), c(0.8, 0.5, 0.2), levels(Cdm), col = "white")


###################################################
### code chunk number 12: surv3
###################################################
m2 <- survfit(Surv(tfd,
                   tfd + lex.dur,
                   lex.Xst == "Ins" ) ~ 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"),
        xlim=c(0,15), xaxs="i",
        yaxs = "i", xlab = "time since DM", ylab = "" )
  lines(m2$time, 1 - m2$surv, lwd = 3, col = "red" )
mat2pol(m3$pstate, c(3,2,1), x = m3$time, yaxs = "i",
        col = c("black","red","transparent"),
        xlim=c(0,15), xaxs="i",
        yaxs = "i", xlab = "time since DM", ylab = "" )
  lines(M2$time, 1 - M2$surv, lwd = 3, col = "black" )


###################################################
### code chunk number 13: cmpr-DMlate.rnw:333-336
###################################################
Sdm <- splitMulti(Adm, tfd = seq(0, 20, 0.1))
summary(Adm)
summary(Sdm)


###################################################
### code chunk number 14: cmpr-DMlate.rnw:344-348
###################################################
round(cbind(
with(subset(Sdm, lex.Xst == "Ins" ), quantile(tfd + lex.dur, 0:10/10)),
with(subset(Sdm, lex.Xst == "Dead"), quantile(tfd + lex.dur, 0:10/10))),
3)


###################################################
### code chunk number 15: cmpr-DMlate.rnw:352-356
###################################################
okn <- c(0,0.5,5,6)
dkn <- c(0,2.0,5,9)
 Ins.glm <- glm.Lexis(Sdm, ~ Ns(tfd, knots = okn), from = "DM", to = "Ins" )
Dead.glm <- glm.Lexis(Sdm, ~ Ns(tfd, knots = dkn), from = "DM", to = "Dead")


###################################################
### code chunk number 16: cmpr-DMlate.rnw:362-366
###################################################
int <- 0.01
nd <- data.frame(tfd = seq(0, 15, int))
l.glm <- ci.pred( Ins.glm, nd)
m.glm <- ci.pred(Dead.glm, nd)


###################################################
### code chunk number 17: Ins-mort
###################################################
matshade(nd$tfd,
         cbind(l.glm, m.glm) * 100,
         plot = TRUE,
         log = "y", ylim = c(0.5, 50),
         col = rep(c("red","black"), 2), lwd = 3,
         xlab = "Time since DM (years)",
         ylab = "Rates per 100 PY")


###################################################
### code chunk number 18: int-ill
###################################################
t <- seq(0, 15, 0.01)
nd <- data.frame(tfd = t)
mu <- ci.pred(Dead.glm, nd)[,1]
head(cbind(t, mu))
plot(t, mu, type="l", lwd = 3,
     xlim = c(0, 7), xaxs = "i",
     ylim = c(0, max(mu)), yaxs = "i")
polygon(t[c(1:501,501:1)], c(mu[1:501], rep(0, 501)),
        col = "gray", border = "transparent")


###################################################
### code chunk number 19: cmpr-DMlate.rnw:426-429
###################################################
mid <- function(x) x[-1] - diff(x) / 2
(x <- c(1:5, 7, 10, 17))
mid(x)


###################################################
### code chunk number 20: cmpr-DMlate.rnw:437-438
###################################################
sum(diff(t[1:501]) * mid(mu[1:501]))


###################################################
### code chunk number 21: cmpr-DMlate.rnw:453-455
###################################################
Mu <- c(0, cumsum(diff(t) * mid(mu)))
head(cbind(t, Mu))


###################################################
### code chunk number 22: cmpr-DMlate.rnw:503-507
###################################################
cR <- ci.Crisk(mods = list(Ins =  Ins.glm,
                          Dead = Dead.glm),
                 nd = nd)
str(cR)


###################################################
### code chunk number 23: crisk
###################################################
matshade(as.numeric(dimnames(cR$Crisk)[[1]]),
         cbind(cR$Crisk[,1,],
               cR$Crisk[,2,],
               cR$Crisk[,3,]), plot = TRUE,
         lwd = 2, col = c("limegreen","red","black"))


###################################################
### code chunk number 24: strisk
###################################################
str(cR$Crisk)
mat2pol(cR$Crisk[,3:1,1], col = c("forestgreen","red","black")[3:1])


###################################################
### code chunk number 25: cmpr-DMlate.rnw:537-542
###################################################
matshade(as.numeric(dimnames(cR$Srisk)[[1]]),
         cbind(cR$Srisk[,1,],
               cR$Srisk[,2,]), plot = TRUE,
         lwd = 2, col = c("black","red"),
         ylim = 0:1, yaxs = "i")


###################################################
### code chunk number 26: cmpr-DMlate.rnw:575-577
###################################################
str(cR$Stime)
round(cR$Stime[1:3 * 500 + 1,"Surv",], 1)


###################################################
### code chunk number 27: cmpr-DMlate.rnw:582-584
###################################################
(mY <- matrix(1:3 * 5, 3, 3)) # using the recycling rule
round(cR$Stime[1:3*500+1,"Surv",] / mY * 100, 1)


###################################################
### code chunk number 28: cmpr-DMlate.rnw:589-596
###################################################
time <- as.numeric(dimnames(cR$Stime)[[1]])
matshade(time, cR$Stime[,"Surv",] /
               cbind(time,
                     time,
                     time) * 100,
         plot=TRUE,
         ylim = 0:1*100, yaxs = "i", xaxs = "i")