### R code from vignette source 'comp.rnw'

###################################################
### code chunk number 1: comp.rnw:7-14 (eval = FALSE)
###################################################
## options(width = 90,
##         show.signif.stars = FALSE)
## par(mar = c(3, 3, 1, 1),
##     mgp = c(3, 1, 0) / 1.6,
##     las = 1,
##    lend = "butt",
##     bty = "n")


###################################################
### code chunk number 2: comp.rnw:16-24
###################################################
options(width = 90,
        show.signif.stars = FALSE,
        SweaveHooks=list(fig = function()
                         par(mar = c(3, 3, 1, 1),
                             mgp = c(3, 1, 0) / 1.6,
                             las = 1,
                            lend = "butt",
                             bty = "n")))


###################################################
### code chunk number 3: comp.rnw:26-32
###################################################
library(Epi)
library(popEpi)
library(survival)
library(tidyverse)
library(mgcv)
clear()


###################################################
### code chunk number 4: comp.rnw:35-37
###################################################
anfang <<- now()
cat("Start time:", format(anfang, "%F, %T"), "\n")


###################################################
### code chunk number 5: comp.rnw:39-45
###################################################
vers <-
data.frame(R = substr(R.version.string, 11, 15),
         Epi = as.character(packageVersion(   "Epi")),
      popEpi = as.character(packageVersion("popEpi")))
names(vers) <- paste(" ", names(vers))
print(vers, row.names = FALSE)


###################################################
### code chunk number 6: comp.rnw:52-53 (eval = FALSE)
###################################################
## vignette("crisk", package = "Epi")


###################################################
### code chunk number 7: comp.rnw:56-57 (eval = FALSE)
###################################################
## vignette(package = "Epi")


###################################################
### code chunk number 8: comp.rnw:300-312
###################################################
data(DMlate)
dl <- mutate(DMlate, dofin = pmin(dodth, dooad, doins, dox, na.rm = TRUE),
                     xstat = factor(case_when(dofin == dodth ~ "Dead",
                                              dofin == dooad ~ "OAD",
                                              dofin == doins ~ "Ins",
                                                        TRUE ~ "DM"),
                                    levels = c("DM", "OAD", "Ins", "Dead")))
str(dl)
dL <- Lexis(exit = list(tfd = dofin - dodm),
     exit.status = xstat,
            data = dl)
summary(dL)


###################################################
### code chunk number 9: comp.rnw:330-332
###################################################
Sdm <- splitLexis(dL, time.scale = "tfd", breaks = seq(0, 20, 1/12))
summary(Sdm)


###################################################
### code chunk number 10: boxes4
###################################################
boxes(Relevel(dL, c(1, 4, 2, 3)),
      boxpos  = list(x = c(15, 85, 75, 15),
                     y = c(85, 85, 30, 15)),
      scale.R = 100,
      show.BE = TRUE )


###################################################
### code chunk number 11: comp.rnw:361-365 (eval = FALSE)
###################################################
## system.time(
## mD <- gamLexis(Sdm, ~ s(tfd, k = 5), to = 'Dead'))
## mO <- gamLexis(Sdm, ~ s(tfd, k = 5), to = 'OAD' )
## mI <- gamLexis(Sdm, ~ s(tfd, k = 5), to = 'Ins' )


###################################################
### code chunk number 12: comp.rnw:369-373
###################################################
system.time(
mD <- glmLexis(Sdm, ~ Ns(tfd, knots = c(0, 1, 3, 6)), to = 'Dead'))
mO <- glmLexis(Sdm, ~ Ns(tfd, knots = c(0, 1, 3, 6)), to = 'OAD' )
mI <- glmLexis(Sdm, ~ Ns(tfd, knots = c(0, 1, 3, 6)), to = 'Ins' )


###################################################
### code chunk number 13: comp.rnw:386-390
###################################################
nd <- data.frame(tfd = seq(0, 10, 1/20))
rownames(nd) <- nd$tfd
str(nd)
head(nd)


###################################################
### code chunk number 14: rates
###################################################
matshade(nd$tfd, cbind(ci.pred(mD, nd),
                       ci.pred(mI, nd),
                       ci.pred(mO, nd)) * 1000,
         col = c("black", "red", "blue"), log = "y", lwd = 3, plot = TRUE,
         xlab = "Time since DM diagnosis (years)",
         ylab = "Rates per 1000 PY", ylim = c(0.05, 500), yaxt = "n")
axis(side = 2, at = ll <- outer(c(1,2,5), -2:3, function(x,y) x * 10^y),
               labels = formatC(ll, digits = 4), las = 1)
axis(side = 2, at = outer(c(1.5,2:9), -2:3, function(x,y) x * 10^y),
               labels = NA, tcl = -0.3)
text(0, 0.5*0.6^c(1,2,0),
     c("Dead","Ins","OAD"),
     col = c("black","red","blue"), adj = 0, font = 2)


###################################################
### code chunk number 15: rates-l
###################################################
matshade(nd$tfd, cbind(ci.pred(mD, nd),
                       ci.pred(mI, nd),
                       ci.pred(mO, nd)) * 1000,
         col = c("black", "red", "blue"), lwd = 3, plot = TRUE,
         xlab = "Time since DM diagnosis (years)",
         ylab = "Rates per 1000 PY", ylim = c(0, 500), yaxs = "i")
text(8, 500 - c(2, 3, 1) * 20,
     c("Dead","Ins","OAD"),
     col = c("black","red","blue"), adj = 0, font = 2)


###################################################
### code chunk number 16: comp.rnw:451-474
###################################################
# utility function to compute midpoints between sucessive values in a vector
mp <- function(x) x[-1] - diff(x) / 2
#
int <- 1 / 20
# rates at midpoints of intervals
lD <- mp(ci.pred(mD, nd)[, 1])
lI <- mp(ci.pred(mI, nd)[, 1])
lO <- mp(ci.pred(mO, nd)[, 1])
#
# cumulative rates and survival function at right border of the intervals
LD <- cumsum(lD) * int
LI <- cumsum(lI) * int
LO <- cumsum(lO) * int
# survival function, formula (1.1)
Sv <- exp(- LD - LI - LO)
#
# when integrating to get the cumulative *risks* we use the average
# of the survival function at the two endpoints
# (adding 1 as the first), formula (1.2)
Sv <- c(1, Sv)
rD <- c(0, cumsum(lD * mp(Sv)) * int)
rI <- c(0, cumsum(lI * mp(Sv)) * int)
rO <- c(0, cumsum(lO * mp(Sv)) * int)


###################################################
### code chunk number 17: comp.rnw:480-484
###################################################
summary(rD + rI + rO + Sv)
oo <- options(digits = 20)
cbind(summary(Sv + rD + rI + rO))
options(oo)


###################################################
### code chunk number 18: stack
###################################################
zz <- mat2pol(cbind(rD, rI, rO, Sv), x = nd$tfd,
              xlim = c(0,10), xaxs = "i", yaxs = "i", las = 1,
              xlab = "Time since DM diagnosis (years)",
              ylab = "Probability",
               col =  c("black","red","blue","forestgreen"))
text(9, mp(zz["9", ]), c("Dead", "Ins", "OAD"," DM"), col = "white")
box(col = "white", lwd = 3)


###################################################
### code chunk number 19: comp.rnw:518-523
###################################################
Sj <- c(sjA = sum(Sv * int),
        sjD = sum(rD * int),
        sjI = sum(rI * int),
        sjO = sum(rO * int))
c(Sj, sum(Sj))


###################################################
### code chunk number 20: comp.rnw:575-577
###################################################
head(cbind(ci.pred(mI, nd),
           ci.exp (mI, nd)))


###################################################
### code chunk number 21: comp.rnw:583-585
###################################################
str(ci.lin(mI, nd, sample = 4))
head(cbind(ci.pred(mI, nd), exp(ci.lin(mI, nd, sample = 4))))


###################################################
### code chunk number 22: comp.rnw:631-639
###################################################
system.time(
res <- ci.Crisk(list(OAD = mO,
                     Ins = mI,
                    Dead = mD),
                            nd = data.frame(tfd = seq(0, 10, 1/20)),
                            nB = 10000,
                          perm = 4:1))
str(res)


###################################################
### code chunk number 23: comp.rnw:674-682
###################################################
system.time(
rsm <- ci.Crisk(list(OAD = mO,
                     Ins = mI,
                    Dead = mD),
                            nd = data.frame(tfd = seq(0, 10, 1/20)),
                            nB = 10000,
                       sim.res = 'rates'))
str(rsm)


###################################################
### code chunk number 24: comp.rnw:690-698
###################################################
system.time(
csm <- ci.Crisk(list(OAD = mO,
                     Ins = mI,
                    Dead = mD),
                            nd = data.frame(tfd = seq(0, 10, 1/20)),
                            nB = 10000,
                       sim.res = 'crisk'))
str(csm)


###################################################
### code chunk number 25: comp.rnw:711-718
###################################################
system.time(
Brates <- aperm(apply(rsm,
                      1:2,
                      quantile,
                      probs = c(.5, .025, .975)),
                c(2, 3, 1)))
str(Brates)


###################################################
### code chunk number 26: comp.rnw:721-728
###################################################
system.time(
apply(rsm,
      1:2,
      quantile,
      probs = c(.5, .025, .975)) |>
aperm(c(2, 3, 1)) -> Brates)
str(Brates)


###################################################
### code chunk number 27: rates-ci
###################################################
matshade(nd$tfd, cbind(ci.pred(mD, nd),
                       ci.pred(mI, nd),
                       ci.pred(mO, nd)) * 1000,
         ylim = c(0.1,500), yaxt = "n",
         ylab = "Rates per 1000 PY",
         xlab = "Time since DM diagnosis (years)",
         col = c("black", "red", "blue"), log = "y", lwd = 3, plot = TRUE)
matlines(nd$tfd,
         cbind(Brates[, "Dead", ],
               Brates[, "Ins" , ],
               Brates[, "OAD" , ]) * 1000,
         col = c("white", "black", "black"), lty = 3, lwd = c(3,1,1))
axis(side = 2, at = ll <- outer(c(1,2,5), -2:3, function(x, y) x * 10^y),
               labels = formatC(ll, digits = 4), las = 1)
axis(side = 2, at = outer(c(1.5, 2:9), -2:3, function(x, y) x * 10^y),
               labels = NA, tcl = -0.3)
text(0, 0.5 * 0.6^c(1,2,0),
     c("Dead", "Ins", "OAD"),
     col = c("black", "red", "blue"), adj = 0, font = 2)


###################################################
### code chunk number 28: crates
###################################################
matshade(res$time,
         cbind(res$Crisk[,"Dead",],
               res$Crisk[,"Ins" ,],
               res$Crisk[,"OAD" ,]), plot = TRUE,
         xlim = c(0,10), xaxs = "i", yaxs = "i", las = 1,
         xlab = "Time since DM diagnosis (years)",
         ylab = "Cumulative probability",
          col = c("black","red","blue"))
text(8, 0.3 + c(1, 0, 2) / 25,
     c("Dead", "Ins", "OAD"),
     col = c("black", "red", "blue"), adj = 0, font = 2)


###################################################
### code chunk number 29: comp.rnw:809-811
###################################################
str(res$Crisk)
str(res$Srisk)


###################################################
### code chunk number 30: stack-ci
###################################################
zz <- mat2pol(res$Crisk[,c("Dead", "Ins", "OAD", "Surv"),1],
              x = res$time,
           xlim = c(0, 10), xaxs = "i", yaxs = "i", las = 1,
           xlab = "Time since DM diagnosis (years)",
           ylab = "Probability",
            col =  c("black","red","blue","forestgreen") )
text(9, mp(zz["9",]), c("Dead", "Ins", "OAD", "DM"), col = "white" )
matshade(res$time,
         cbind(res$Srisk[, 1, ],
               res$Srisk[, 2, ],
               res$Srisk[, 3, ]),
         col = 'transparent', col.shade = "white", alpha = 0.4)


###################################################
### code chunk number 31: comp.rnw:851-854
###################################################
s510 <- res$Stime[c("5", "10"),,]
dimnames(s510)[[1]] <- c(" 5 yr","10 yr")
round(ftable(s510, row.vars=1:2), 2)


###################################################
### code chunk number 32: comp.rnw:879-893
###################################################
data(DMlate)
dl <- mutate(DMlate, dofin = pmin(dodth, dooad, doins, dox, na.rm = TRUE),
                     xstat = factor(case_when(dofin == dodth ~ "Dth",
                                              dofin == dooad ~ "OAD",
                                              dofin == doins ~ "Ins",
                                                        TRUE ~ "DM"),
                                    levels = c("DM", "OAD", "Ins", "Dth")))
dmL <- Lexis(entry = list(per = dodm,
                          age = dodm - dobth,
                          tfD = dodm - dodm),
              exit = list(per = dofin),
       exit.status = xstat,
              data = dl)
summary(dmL, t = T)


###################################################
### code chunk number 33: boxes
###################################################
boxes(dmL, boxpos = TRUE)


###################################################
### code chunk number 34: comp.rnw:908-910
###################################################
sL <- splitLexis(dmL, time.scale = "age", breaks = seq(0, 120, 1/2))
summary(sL)


###################################################
### code chunk number 35: comp.rnw:912-915
###################################################
mOAD <- gamLexis(sL, ~ s(tfD, by=sex), to = "OAD")
mIns <- gamLexis(sL, ~ s(tfD, by=sex), to = "Ins")
mDth <- gamLexis(sL, ~ s(tfD, by=sex), to = "Dth")


###################################################
### code chunk number 36: comp.rnw:924-925
###################################################
nm <- data.frame(tfD = seq(0, 10, 1/20), sex = "M")


###################################################
### code chunk number 37: comp.rnw:929-936
###################################################
system.time(
cR <- ci.Crisk(list(OAD = mOAD,
                    Ins = mIns,
                    Dth = mDth),
               nB = 5000,
               nd = nm))
str(cR)


###################################################
### code chunk number 38: cR
###################################################
clr <- c("limegreen", "orange", "black")
matshade(cR$time, cbind(cR$Crisk[, "OAD", ],
                        cR$Crisk[, "Ins", ],
                        cR$Crisk[, "Dth", ]),
         col = clr, lty = 1, lwd = 2,
         plot = TRUE, ylim = c(0, 0.5), yaxs = "i",
         xlab = "Time since DM (years)", ylab = "Probability of state")
text(0, 1/3 - c(1,3,2)/30, c("OAD", "Ins", "Dth"),
     col = clr, adj = 0, font = 2)


###################################################
### code chunk number 39: Sr1
###################################################
matshade(cR$time, cbind(cR$Srisk[,1,],
                        cR$Srisk[,2,],
                        cR$Srisk[,3,]),
         col = "black", lty = 1, lwd = 2,
         plot = TRUE, ylim = c(0,1), xaxs = "i", yaxs = "i")
text(9, mp(c(0, cR$Srisk["9", , 1], 1)),
     rev(c(dimnames(cR$Crisk)[[2]])))
box(bty = "o")


###################################################
### code chunk number 40: Sr2
###################################################
zz <- mat2pol(cR$Crisk[, c("Dth", "Ins", "OAD", "Surv"), "50%"],
              x = cR$time,
           xlim = c(0,10), xaxs = "i", yaxs = "i", las = 1,
           xlab = "Time since DM diagnosis (years)",
           ylab = "Probability",
            col =  c(gray(0.3), "red", "blue", "limegreen"))
matshade(cR$time, cbind(cR$Srisk[,1,],
                        cR$Srisk[,2,],
                        cR$Srisk[,3,]),
         col = "transparent", col.shade = "white", alpha = 0.4)
text(9, mp(c(0, cR$Srisk["9", , 1], 1)),
     rev(c(dimnames(cR$Crisk)[[2]])), col = "white")
box(bty = "o", col ="white")


###################################################
### code chunk number 41: comp.rnw:1023-1024
###################################################
ftable(round(cR$Stime[paste(c(2, 5, 10)), , ], 2), row.vars = 1)


###################################################
### code chunk number 42: comp.rnw:1051-1070
###################################################
nm <- data.frame(tfD = seq(0, 10, 1/20), sex = "M")
nw <- data.frame(tfD = seq(0, 10, 1/20), sex = "F")
# set the seed
set.seed(1952)
mR <- ci.Crisk(list(OAD = mOAD,
                    Ins = mIns,
                    Dth = mDth),
               nd = nm,
               nB = 500,
          sim.res = "crisk" )
# REset the seed to get the same sequence
set.seed(1952)
wR <- ci.Crisk(list(OAD = mOAD,
                    Ins = mIns,
                    Dth = mDth),
               nd = nw,
               nB = 500,
          sim.res = "crisk" )
str(wR)


###################################################
### code chunk number 43: comp.rnw:1075-1080
###################################################
dS <- mR[,"Surv",] - wR[,"Surv",]
dS <- apply(dS, 1, quantile, probs = c(.5, .025, .975)) * 100
str(dS)
rS <- mR[,"Surv",] / wR[,"Surv",]
rS <- apply(rS, 1, quantile, probs = c(.5, .025, .975))


###################################################
### code chunk number 44: difrat
###################################################
par(mfrow = c(1,2))
matshade(as.numeric(colnames(dS)), t(dS), plot = TRUE,
         lwd = 3, ylim = c(-7, 7),
         xlab = "Time since DM diagnosis (years)",
         ylab = "Men - Women survival difference (%)")
abline(h = 0)
matshade(as.numeric(colnames(rS)), t(rS), plot = TRUE,
         lwd = 3, ylim = c(1/1.2, 1.2), log ="y",
         xlab = "Time since DM diagnosis (years)",
         ylab = "Men - Women survival ratio")
abline(h = 1)


###################################################
### code chunk number 45: comp.rnw:1106-1116
###################################################
fR <- ci.Crisk(list(OAD = mOAD,
                    Ins = mIns,
                    Dth = mDth),
               nd = nw,
               nB = 500,
          sim.res = "crisk" )
dxS <- mR[,"Surv",] - fR[,"Surv",]
dxS <- apply(dxS, 1, quantile, probs = c(.5, .025, .975)) * 100
rxS <- mR[,"Surv",] / fR[,"Surv",]
rxS <- apply(rxS, 1, quantile, probs = c(.5, .025, .975))


###################################################
### code chunk number 46: difratx
###################################################
par(mfrow = c(1,2))
matshade(as.numeric(colnames(dS)), t(dS), plot = TRUE,
         lwd = 3, ylim = c(-7, 7),
         xlab = "Time since DM diagnosis (years)",
         ylab = "Men - Women drug-free survival difference (%)")
matshade(as.numeric(colnames(dxS)), t(dxS), 
         lty = "21", lend = "butt", lwd = 3, col = "forestgreen")
abline(h = 0)

matshade(as.numeric(colnames(rS)), t(rS), plot = TRUE,
         lwd = 3, ylim = c(1/1.2, 1.2), log ="y",
         xlab = "Time since DM diagnosis (years)",
         ylab = "Men - Women drug-free survival ratio")
matshade(as.numeric(colnames(rxS)), t(rxS),
         lty = "21", lend = "butt", lwd = 3, col = "forestgreen")
abline(h = 1)


###################################################
### code chunk number 47: comp.rnw:1153-1157
###################################################
ende <- Sys.time()
cat("  Start time:", format(anfang, "%F, %T"),
  "\n    End time:", format(  ende, "%F, %T"),
  "\nElapsed time:", round(difftime(ende, anfang, units = "mins"), 2), "minutes\n")