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

###################################################
### code chunk number 1: surv.rnw:4-11 (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: surv.rnw:13-21
###################################################
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: surv.rnw:23-28
###################################################
library(Epi)
library(popEpi)
library(survival)
library(tidyverse)
clear()


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


###################################################
### code chunk number 5: surv.rnw:35-40
###################################################
vrs <-
data.frame(R = do.call(paste, c(R.Version()[c("major","minor")], sep =".")),
         Epi = packageVersion("Epi"),
      popEpi = packageVersion("popEpi"))
print(vrs, row.names = FALSE)


###################################################
### code chunk number 6: surv.rnw:65-71
###################################################
data(DMlate)
par(mfrow = 2:1)
hist(DMlate$dodm, breaks = seq(1995, 2010, 1/4), col = 1)
abline(v = 1995:2010, col = "red")
hist(DMlate$dodm, breaks = seq(1995, 2010, 1/12), col = 1)
abline(v = 1995:2010, col = "red")


###################################################
### code chunk number 7: surv.rnw:76-79
###################################################
set.seed(19540803)
DMlate <- DMlate[sample(1:nrow(DMlate), 1000), ]
str(DMlate)


###################################################
### code chunk number 8: KM
###################################################
sf <- survfit(Surv(dox - dodm, !is.na(dodth)) ~ 1, data = DMlate)
sf
plot(sf, yaxs ="i")
abline(h = 0.5, lty = 3)


###################################################
### code chunk number 9: surv.rnw:153-154
###################################################
(mort <- with(DMlate, sum(!is.na(dodth)) / sum(dox - dodm)))


###################################################
### code chunk number 10: KMexp
###################################################
plot(sf, yaxs = "i")
abline(h = 0.5, col = "gray")
t = seq(0, 15, 0.1)
lines(t, exp(-mort * t), col = "red", lwd =2)


###################################################
### code chunk number 11: surv.rnw:201-205
###################################################
Lx <- Lexis(exit = list(tfd = dox - dodm),
     exit.status = factor(!is.na(dodth),
                          labels = c("DM", "Dead")),
            data = DMlate)


###################################################
### code chunk number 12: surv.rnw:230-231
###################################################
summary(Lx)


###################################################
### code chunk number 13: surv.rnw:237-238
###################################################
Lx[1:5, ]


###################################################
### code chunk number 14: surv.rnw:246-249
###################################################
sL <- splitLexis(Lx, seq(0, 15, 0.5))
subset(Lx, lex.id %in% c(2,4))
subset(sL, lex.id %in% c(2,4))


###################################################
### code chunk number 15: surv.rnw:277-278
###################################################
m0 <- glmLexis(sL, ~ Ns(tfd, knots = c(0,1,4,10)))


###################################################
### code chunk number 16: surv.rnw:286-290 (eval = FALSE)
###################################################
## mr <- glm(cbind(lex.Xst == "Dead", lex.dur)
##           ~ Ns(tfd, knots = c(0,1,4,10)),
##           family = poisreg,
##             data = sL)


###################################################
### code chunk number 17: surv.rnw:295-300 (eval = FALSE)
###################################################
## mp <- glm(lex.Xst == "Dead"
##           ~ Ns(tfd, knots = c(0,1,4,10)),
##           offset = log(lex.dur),
##           family = poisson,
##             data = sL)


###################################################
### code chunk number 18: surv.rnw:317-320
###################################################
nd <- data.frame(tfd = seq(0, 15, 0.1))
pm <- ci.pred(m0, nd)
head(cbind(tfd = nd$tfd, pm))


###################################################
### code chunk number 19: mort
###################################################
matshade(nd$tfd, pm * 100, plot = TRUE, lwd = 2,
         ylim = c(0, 10), yaxs = "i",
         xlab = "Time since diabetes (years)",
         ylab = "Mortality rate (% / year)")
abline(v = c(0,1,4,10), lty = 3)


###################################################
### code chunk number 20: surv.rnw:351-353
###################################################
am <- with(sL, tapply((dodm - dobth) + tfd, floor(tfd), mean))
cbind(am, c(NA, diff(am)))


###################################################
### code chunk number 21: KMmod
###################################################
plot(sf, yaxs = "i", ylim = c(0.5, 1),
         xlab = "Time since diabetes (years)",
         ylab = "Survival probability")
abline(h = 0.5)
ps <- ci.surv(m0, nd)
matshade(nd$tfd, ps, col = "red", lwd = 2)


###################################################
### code chunk number 22: surv.rnw:413-417
###################################################
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")