### R code from vignette source 'mort.rnw' ################################################### ### code chunk number 1: mort.rnw:2-9 (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: mort.rnw:11-19 ################################################### 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: mort.rnw:23-27 ################################################### library(Epi) library(popEpi) library(survival) library(tidyverse) ################################################### ### code chunk number 4: mort.rnw:29-31 ################################################### clear() setwd("C:/Bendix/teach/AdvCoh/courses/IDEG2025/pracs") ################################################### ### code chunk number 5: mort.rnw:33-38 ################################################### 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: mort.rnw:79-85 ################################################### data(DMlate) set.seed(1952) DMlate <- DMlate[sample(1:nrow(DMlate), 2000),] rownames(DMlate) <- 1:2000 str(DMlate) head(DMlate) ################################################### ### code chunk number 7: mort.rnw:106-107 ################################################### (y <- with(DMlate, sum(dox - dodm))) ################################################### ### code chunk number 8: mort.rnw:115-116 ################################################### (d <- with(DMlate, sum(!is.na(dodth)))) ################################################### ### code chunk number 9: mort.rnw:122-123 ################################################### d / y ################################################### ### code chunk number 10: mort.rnw:126-127 ################################################### d / y * 100 ################################################### ### code chunk number 11: mort.rnw:136-145 ################################################### tt <- xtabs(cbind(D = !is.na(dodth), Y = dox - dodm) ~ agr, data = mutate(DMlate, agr = cut(dodm - dobth, seq(0, 100, 10), right = FALSE))) tt cbind(mort = tt[, "D"] / tt[, "Y"] * 100) ################################################### ### code chunk number 12: mort.rnw:160-161 ################################################### cbind(tt[, 1] / tt[, 2] * 100) ################################################### ### code chunk number 13: mort.rnw:194-201 ################################################### Lx <- Lexis(entry = list(age = dodm - dobth), exit = list(age = dox - dobth), exit.status = factor(!is.na(dodth), labels = c("Alive","Dead")), data = DMlate) subset(DMlate, near(dodm, dox)) summary(Lx) head(Lx) ################################################### ### code chunk number 14: mort.rnw:224-227 ################################################### Sx <- splitLexis(Lx, breaks = seq(0, 100, 5)) summary(Lx) summary(Sx) ################################################### ### code chunk number 15: mort.rnw:240-246 ################################################### tt <- xtabs(cbind(D = lex.Xst == "Dead", Y = lex.dur) ~ I(floor(age / 5) * 5), data = Sx) tt (rt <- cbind(mort = tt[, "D"] / tt[, "Y"] * 100)) ################################################### ### code chunk number 16: morta ################################################### plot(rownames(tt), rt, log ="y", type = "o", xlab = "Age", pch = 16) ################################################### ### code chunk number 17: mort.rnw:264-264 ################################################### ################################################### ### code chunk number 18: mortas ################################################### plot(rownames(tt), rt, log ="y", type = "s", xlab = "Age") points(rownames(tt), rt, pch = 16) ################################################### ### code chunk number 19: mort.rnw:298-301 ################################################### Sx <- splitLexis(Lx, breaks = 0:100, time.scale = "age") summary(Lx) summary(Sx) ################################################### ### code chunk number 20: mort.rnw:308-310 ################################################### subset(Lx, lex.id %in% 2:3) subset(Sx, lex.id %in% 2:3) ################################################### ### code chunk number 21: mort.rnw:314-315 ################################################### mL <- glmLexis(Sx, ~ Ns(age, knots = seq(40, 80, 10))) ################################################### ### code chunk number 22: mort.rnw:324-330 ################################################### mP <- glm((lex.Xst == "Dead") ~ Ns(age, knots = seq(40, 80, 10)), offset = log(lex.dur), family = poisson, data = Sx) round(cbind(ci.exp(mL), ci.exp(mP)), 4) ################################################### ### code chunk number 23: mort.rnw:346-349 ################################################### nd <- data.frame(age = 30:95) pr.rates <- ci.pred(mL, nd) * 100 head(pr.rates) ################################################### ### code chunk number 24: mortaspl ################################################### matshade(nd$age, pr.rates, plot = TRUE, log = "y", lwd = 3, xlab = "Attained age", ylab = "Mortality rate per 100 PY") ################################################### ### code chunk number 25: mort.rnw:370-376 ################################################### Lx <- Lexis(entry = list(tfd = dodm - dodm), exit = list(tfd = dox - dodm), exit.status = factor(!is.na(dodth), labels = c("Alive", "Dead")), data = DMlate) summary(Lx) head(Lx) ################################################### ### code chunk number 26: mort.rnw:386-389 ################################################### Sx <- splitLexis(Lx, breaks = seq(0, 20, 0.5)) summary(Lx) summary(Sx) ################################################### ### code chunk number 27: mort.rnw:395-397 ################################################### subset(Lx, lex.id %in% 2:3) subset(Sx, lex.id %in% 2:3) ################################################### ### code chunk number 28: mort.rnw:401-402 ################################################### tL <- glmLexis(Sx, ~ Ns(tfd, knots = c(0, 1, 3, 6, 10))) ################################################### ### code chunk number 29: mort.rnw:408-411 ################################################### nd <- data.frame(tfd = seq(0, 12, 0.2)) pr.rates <- ci.pred(tL, nd) head(pr.rates) ################################################### ### code chunk number 30: mortdspl ################################################### matshade(nd$tfd, pr.rates * 100, plot = TRUE, log = "y", lwd = 3, xlab = "Time since DM diagnosis", ylab = "Mortality rate per 100 PY") ################################################### ### code chunk number 31: mort.rnw:451-454 ################################################### head(cbind(nd, pr.rates)) surv <- exp(-cumsum(pr.rates[, 1]) * 0.2) plot(nd$tfd, surv, type = "l", ylim = c(0,1)) ################################################### ### code chunk number 32: surv ################################################### matshade(nd$tfd, ci.surv(tL, nd), plot = TRUE, lwd = 2, ylim = c(0.5,1), yaxs = "i", xlab = "Time since diagnosis (years)", ylab = "Survival probability") lines(survfit(Surv(dox - dodm, !is.na(dodth)) ~ 1, data = Lx), col = "red")