### R code from vignette source 'rmst.rnw' ################################################### ### code chunk number 1: rmst.rnw:3-10 (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: rmst.rnw:12-20 ################################################### 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: rmst.rnw:23-28 ################################################### library(Epi) library(popEpi) library(survival) library(tidyverse) clear() ################################################### ### code chunk number 4: rmst.rnw:31-33 ################################################### anfang <<- now() cat("Start time:", format(anfang, "%F, %T"), "\n") ################################################### ### code chunk number 5: rmst.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: rmst.rnw:60-72 ################################################### data(DMlate) set.seed(19540803) DMlate <- DMlate[sample(1:nrow(DMlate), 1000), ] Lx <- Lexis(entry = list(age = dodm - dobth, tfd = 0), exit = list(tfd = dox - dodm), exit.status = factor(!is.na(dodth), labels = c("DM", "Dead")), data = DMlate) sL <- splitLexis(Lx, seq(0, 15, 0.5), "tfd") summary(Lx) summary(sL) ################################################### ### code chunk number 7: rmst.rnw:77-81 ################################################### m1 <- glmLexis(sL, ~ Ns(age, knots = c(30, 50, 70)) + Ns(tfd, knots = c(0, 1, 4, 10)) + sex) round(ci.exp(m1, subset = "sex"), 3) ################################################### ### code chunk number 8: rmst.rnw:115-120 ################################################### surv.arr <- NArray(list(adx = c(50, 60, 70), sex = c("M", "F"), tfd = tfd <- seq(0, 15, .1), surv = c("surv", "lo", "up"))) str(surv.arr) ################################################### ### code chunk number 9: rmst.rnw:142-150 ################################################### for(adx in c(50, 60, 70)) for( sx in c("M", "F")) { nd <- data.frame(tfd = tfd, age = adx + tfd, sex = sx) surv.arr[paste(adx), sx, , ] <- ci.surv(m1, nd) } ################################################### ### code chunk number 10: rmst.rnw:171-179 ################################################### ftable(surv.arr[,,c("5","10"),], row.vars = c(1,3,2)) # difference in survival probabilities (%) round(100 * ftable(surv.arr[,"F",c("5","10"),1] - surv.arr[,"M",c("5","10"),1]), 2) # remember the brackets when using the pipe, |> ((surv.arr[,"F",c("5","10"),1] - surv.arr[,"M",c("5","10"),1]) * 100) |> ftable() |> round(2) ################################################### ### code chunk number 11: rmst.rnw:191-195 ################################################### head(nd) msM <- ci.Crisk(list(Mort = m1), mutate(nd, sex = "M"))$Stime msF <- ci.Crisk(list(Mort = m1), mutate(nd, sex = "F"))$Stime str(msF) ################################################### ### code chunk number 12: rmst.rnw:202-205 ################################################### msM[c("5","10"), "Surv", ] msF[c("5","10"), "Surv", ] msF[c("5","10"), "Surv", 1] - msM[c("5","10"), "Surv", 1] ################################################### ### code chunk number 13: rmst.rnw:225-235 ################################################### head(nd) nB <- 10000 # no of bootstrap samples ain <- 5:7 * 10 sex <- c("M", "F") simres <- NArray(list(adx = ain, sex = sex, tfd = nd$tfd, sim = 1:nB)) str(simres) length(simres) ################################################### ### code chunk number 14: rmst.rnw:244-255 ################################################### for (adx in ain) for ( sx in sex) { set.seed(20250503) simres[paste(adx), sx, , ] <- ci.Crisk(list(Mort = m1), nd = mutate(nd, sex = sx, age = adx + tfd), nB = nB, sim.res = "crisk")[, "Surv", ] cat(adx, sx, format(Sys.time(), format = "%T"), "\n") } ################################################### ### code chunk number 15: rmst.rnw:287-294 ################################################### mid <- function(z) z[-1] - diff(z) / 2 get5.10 <- function(x) { c("5" = sum(mid(x[1:51])), "10"= sum(mid(x[1:101]))) / 10 } ################################################### ### code chunk number 16: rmst.rnw:309-312 ################################################### str(simres) zz <- apply(simres, c(1,2,4), get5.10) str(zz) ################################################### ### code chunk number 17: rmst.rnw:322-334 ################################################### round(ftable(apply(zz, 1:3, quantile, c(50, 5, 95)/100), col.vars = 2:1), 2) # the tidyverse version would be apply(zz, 1:3, quantile, c(50, 5, 95)/100) |> ftable(col.vars = 2:1) |> round(2) # the usual tidyverse version takes up a bit more space zz |> apply(1:3, quantile, c(50, 5, 95) / 100) |> ftable(col.vars = 2:1) |> round(2) ################################################### ### code chunk number 18: rmst.rnw:350-354 ################################################### m1 <- glmLexis(sL, ~ Ns(age, knots = c(50, 60, 70, 80)) + Ns(tfd, knots = c(0, 1, 5)) + sex) round(ci.exp(m1, subset = "sex"), 3) ################################################### ### code chunk number 19: m1 ################################################### pd <- expand.grid(tfd = c(NA, seq(0, 35, .5)), ain = c(50, 55, 60, 65, 70, 75)) |> mutate(age = ain + tfd, age = ifelse(age < 86, age, NA)) head(pd) matshade(pd$age, ci.pred(m1, mutate(pd, sex = "M")) * 1000, xlab = "Attained age", ylab = "Mortality per 1000 PY", plot = TRUE, col = "blue", lwd = 3, log = "y") matshade(pd$age, ci.pred(m1, mutate(pd, sex = "F")) * 1000, col = "red", lwd = 3, log = "y") ################################################### ### code chunk number 20: rmst.rnw:377-381 ################################################### mi <- glmLexis(sL, ~ (Ns(age, knots = c(50, 60, 70, 80)) + Ns(tfd, knots = c(0, 1, 5))) * sex) anova(m1, mi, test = "Chisq") ################################################### ### code chunk number 21: mi ################################################### round(ci.exp(mi)[,1, drop = FALSE], 3) par(mfrow = c(2,2)) # main effects model matshade(pd$age, ci.pred(m1, mutate(pd, sex = "M")) * 1000, plot = TRUE, ylim = c(2, 200), xlab = "Attained age", ylab = "Mortality per 1000 PY", col = "blue", lwd = 2, log = "y", alpha = .1) matshade(pd$age, ci.pred(m1, mutate(pd, sex = "F")) * 1000, col = "red", lwd = 2, alpha = .1) matshade(pd$age, ci.pred(mi, mutate(pd, sex = "M")) * 1000, plot = TRUE, ylim = c(2, 200), xlab = "Attained age", ylab = "Mortality per 1000 PY", col = "blue", lwd = 2, log = "y", alpha = .1) matshade(pd$age, ci.pred(mi, mutate(pd, sex = "F")) * 1000, col = "red", lwd = 2, alpha = .1) # interaction model matshade(pd$tfd, ci.pred(m1, mutate(pd, sex = "M")) * 1000, plot = TRUE, ylim = c(2, 200), xlab = "Diabetes duration", ylab = "Mortality per 1000 PY", col = "blue", lwd = 2, log = "y", alpha = .1) matshade(pd$tfd, ci.pred(m1, mutate(pd, sex = "F")) * 1000, col = "red", lwd = 2, alpha = .1) matshade(pd$tfd, ci.pred(mi, mutate(pd, sex = "M")) * 1000, plot = TRUE, ylim = c(2, 200), xlab = "Diabetes duration", ylab = "Mortality per 1000 PY", col = "blue", lwd = 2, log = "y", alpha = .1) matshade(pd$tfd, ci.pred(mi, mutate(pd, sex = "F")) * 1000, col = "red", lwd = 2, alpha = .1) ################################################### ### code chunk number 22: rmst.rnw:427-451 ################################################### simresi <- simres * NA # a quick way to get an array with same dimensions for (adx in ain) for ( sx in sex) { set.seed(20250503) cat(adx, sx, format(Sys.time(), format = "%T"), "\n") simresi[paste(adx), sx, , ] <- ci.Crisk(list(Mort = mi), nd = mutate(nd, sex = sx, age = adx + tfd), nB = nB, sim.res = "crisk")[, "Surv", ] } # # tidyverse apply(simresi, c(1, 2, 4), get5.10) |> # bootstrap sample of 5 and 10 y RMST apply(1:3, quantile, c(50, 5, 95)/100) |> # quantiles of these ftable(col.vars = 2:1) |> # nice formatting of resulting array round(2) # # old fashioned: one at a time aa <- apply(simresi, c(1,2,4), get5.10) bb <- apply(aa, 1:3, quantile, c(50, 5, 95)/100) cc <- ftable(bb, col.vars = 2:1) round(cc, 2) ################################################### ### code chunk number 23: rmst.rnw:466-473 ################################################### str(simres) simdif <- simres[,"F",,] - simres[,"M",,] str(simdif) dd <- apply(simdif, c(1,3), get5.10) str(dd) dd <- apply(dd, 1:2, quantile, c(50, 5, 95)/100) round(ftable(dd, col.vars = 2:1), 2) ################################################### ### code chunk number 24: rmst.rnw:476-485 ################################################### simdifi <- simresi[,"F",,] - simresi[,"M",,] ii <- apply(simdifi, c(1,3), get5.10) ii <- apply(ii, 1:2, quantile, c(50, 5, 95)/100) # # main effects model (proportional hazards) round(ftable(dd, col.vars = 2:1), 2) # # interaction model round(ftable(ii, col.vars = 2:1), 2) ################################################### ### code chunk number 25: rmst.rnw:500-508 ################################################### str(simresi) simadif <- simresi[2:3,,,] - simresi[1:2,,,] dimnames(simadif)[[1]] <- paste0(dimnames(simresi)[[1]][2:3], " - ", dimnames(simresi)[[1]][1:2]) str(simadif) aa <- apply(simadif, c(1,2,4), get5.10) bb <- apply(aa, 1:3, quantile, c(50, 5, 95)/100) round(ftable(bb, col.vars = 2:1), 2) ################################################### ### code chunk number 26: rmst.rnw:536-579 ################################################### nB <- 10000 # no of bootstrap samples nd <- data.frame(tfd = seq(0, 30, .1)) (ain <- seq(50, 75, 5)) sex <- c("M", "F") simres <- NArray(list(adx = ain, sex = sex, tfd = nd$tfd, sim = 1:nB)) str(simres) # # simulated cumulative risks for (adx in ain) for ( sx in sex) { set.seed(20250503) simres[paste(adx), sx, , ] <- ci.Crisk(list(Mort = m1), nd = mutate(nd, sex = sx, age = adx + tfd), nB = nB, sim.res = "crisk")[, "Surv", ] cat(adx, sx, format(Sys.time(), format = "%T"), "\n") } # # function to derive integrals from function values 0.1 apart getRMST <- function(x) { c("10" = sum(mid(x[1:101])), "20" = sum(mid(x[1:201])), "30" = sum(mid(x[1:301]))) / 10 } str(simres) # # compute integrals to chosen horizons aa <- apply(simres, c(1,2,4), getRMST) str(aa) # # confidence intervals for RMST bb <- apply(aa, 1:3, quantile, c(50, 5, 95) / 100) str(bb) round(ftable(bb, col.vars = 2:1), 2) ################################################### ### code chunk number 27: forest ################################################### par(xaxt = "n") plotEst(t(bb[,"10",,"M"]), y=6:1, col = "blue", xlim = c(0,30), xaxt = "n", grid = seq(5, 30, 5), xlab = "RMST the next 10, 20 or 30 years") linesEst(t(bb[,"20",,"M"]), y=6:1 - 0.1, col = "blue") linesEst(t(bb[,"30",,"M"]), y=6:1 - 0.2, col = "blue") linesEst(t(bb[,"10",,"F"]), y=6:1 - 0.3, col = "red") linesEst(t(bb[,"20",,"F"]), y=6:1 - 0.4, col = "red") linesEst(t(bb[,"30",,"F"]), y=6:1 - 0.5, col = "red") par(xaxt = "s") axis(side = 1, at = seq(5, 25, 5)) axis(side = 1, at = 0:30, labels = NA, tcl = -0.3) ################################################### ### code chunk number 28: rmst.rnw:615-619 ################################################### ende <- now() 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")