R code for chapter 4

Quicklinks for Chapter 4

Section 4.2

We start by computing the cumulative hazards using the mvna package.

First step: Creation of a matrix of logical values indicating the possible transition types.

tra <- matrix(FALSE, ncol = 3, nrow = 3) dimnames(tra) <- list(c("0", "1", "2"), c("0", "1", "2")) tra[1, 2:3] <- TRUE tra

We have to do a slight modification of the data set to comply with the mvna function.

id <- seq_along(obs.cause) from <- rep(0, length(obs.cause)) to <- as.factor(ifelse(obs.cause == 0, "cens", obs.cause)) my.data <- data.frame(id, from, to, time = obs.times)

First load the library. Then, you can compute the cumulative hazards.

library(mvna) my.nelaal <- mvna(my.data, c("0", "1", "2"), tra, "cens")

Figure 4.2

xyplot(my.nelaal, strip = strip.custom(bg = "white"), ylab="Cumulative Hazard", lwd = 2) ### Add theorical quantities to the plot ## computing the true cumulative hazards x <- list(seq(0, max(my.nelaal$"0 1"$time), 0.01), seq(0, max(my.nelaal$"0 2"$time), 0.01)) y <- list(0.3 * x[[1]], 0.6 * x[[2]]) ## and plot them tcl <- trellis.currentLayout() for (i in 1:ncol(tcl)) { trellis.focus("panel", i, 1, highlight=FALSE) panel.lines(x[[i]], y[[i]], col="black", lty=1, lwd=1) trellis.unfocus() }

Subsection: Survival function

Estimation of P(T > t) using the survfit function of the survival package.

my.fit.surv <- survfit(Surv(time, to != "cens") ~ 1, data = my.data, conf.type = "log-log")

Figure 4.3

plot(my.fit.surv, xlab = "Time", ylab = "Survival Probability", mark.time= FALSE, lwd = 2) curve(exp(-0.9*x), add = TRUE, lwd = 2)

Subsection: Estimating the survival function from a multistate perspective

Estimation of the survival function using etm. and creation of a matrix with logical components defining the possible transitions.

tra.km <- matrix(FALSE, ncol=2, nrow=2) dimnames(tra.km) <- list(c("0", "1"), c("0", "1")) tra.km[1, 2] <- TRUE

Transformation of the data. New event indicator to takes 1 if an event (of any type), "cens" otherwise.

my.data.km <- my.data my.data.km$to <- ifelse(my.data.km$to == "cens", "cens", 1)

Call to etm().

library(etm) km.etm <- etm(my.data.km, c("0", "1"), tra.km, "cens", s = 0)

You can now check the results and compare them to Figure 4.3.

lines(x = km.etm$time, y = km.etm$est[1, 1, ], type = "s", col = "red")

Subsection: Cumulative incidence functions

Estimation of the CIFs using the cuminc function in the cmprsk package.

require(cmprsk) my.cif <- cuminc(my.data$time, my.data$to, cencode = "cens")

Figure 4.4

plot(my.cif, xlab="Time", wh=c(-6,-6), lwd=2, ylab="Cumulative event probability") ## True quantities curve((1 - exp(-0.9 * x))/3, add = TRUE) curve((1 - exp(-0.9 * x)) * 2/3, add = TRUE) text(2.5, 0.7, c("Event type 2"), cex=1.2) text(2.5, 0.2, c("Event type 1"), cex=1.2) box()

Subsection: Estimating the cumulative incidence functions from a multistate perspective

Estimation of the CIFs using etm.

cif.etm <- etm(my.data, c("0", "1", "2"), tra, "cens", s = 0)

Estimates at t₁

cif.etm$est[, , 1]

Subsection: Left-truncation

We add left-truncation (Gamma distributed) to the simulated data set.

lt.times <- rgamma(100, shape= 0.5, rate = 2)

The entry times will be defined as the left-truncation times lt.times. People for whom entry > exit will not be included (formally, they are never oberved in the 'study').

my.data2 <- my.data[,1:3] ##do not keep my.data$time my.data2$entry <- lt.times my.data2$exit <- my.data$time my.data2 <- my.data2[my.data2$entry < my.data2$exit, ]

Next, we estimate the cumulative hazards.

my.nelaal2 <- mvna(my.data2, c("0", "1", "2"), tra, "cens")

Figure 4.5

print(xyplot(my.nelaal2, strip=strip.custom(bg="white"), ylab="Cumulative Hazard", lwd=2)) ## theoretical quantities x <- list(seq(0, max(my.nelaal2$"0 1"$time), 0.01), seq(0, max(my.nelaal2$"0 2"$time), 0.01)) y <- list(0.3 * x[[1]], 0.6 * x[[2]]) tcl <- trellis.currentLayout() for (i in 1:ncol(tcl)) { trellis.focus("panel", i, 1, highlight=FALSE) panel.lines(x[[i]], y[[i]], col="black", lty=1, lwd=1) trellis.unfocus() }

Figure 4.6

plot(x = c(0, ceiling(max(obs.times))), y = c(1, 100), type = "n", xlab = "Time", ylab = "Individuals", axes = FALSE) axis(side = 1, at = seq(0, ceiling(max(obs.times)), by = 1)) axis(side = 2, at = c(1, seq(10, 100,by = 10))) box() for(i in 1:100){ if(lt.times[i] < obs.times[i]) { curve(0 * x + i, from = lt.times[i], to = obs.times[i], add = TRUE, lwd = 1.5) if (obs.cause[i] == 0) { points(obs.times[i], i, pch = 21) } } }

Estimation of the survival function using survfit:

my.fit.surv2 <- survfit(Surv(entry,exit, to != "cens") ~ 1, data = my.data2, conf.type = "log-log")

and using etm.

cif.etm2 <- etm(my.data2, c("0", "1", "2"), tra, "cens", s = 0)

Section 4.3

First, load the data.

data(sir.adm)

Cause-specific hazards

We transform sir.adm into a multistate-type data set.

to <- ifelse(sir.adm$status == 0, "cens", ifelse(sir.adm$status == 1, 2, 1)) my.sir.data <- data.frame(id = sir.adm$id, from = 0, to, time = sir.adm$time, pneu = sir.adm$pneu)

Estimation of the cumulative CSH, stratified pneumonia status on admission:

## no pneumonia my.nelaal.nop <- mvna(my.sir.data[my.sir.data$pneu == 0, ], c("0", "1", "2"), tra, "cens") ## with pneumonia my.nelaal.p <- mvna(my.sir.data[my.sir.data$pneu == 1, ], c("0", "1", "2"), tra, "cens")

Figure 4.7

ltheme <- canonical.theme(color = FALSE) ## in-built B&W theme ltheme$strip.background$col <- "white" ## change strip bg lattice.options(default.theme = ltheme) ## set as default

Plot for "no pneumonia":

dessin.nop <- xyplot(my.nelaal.nop, tr.choice = c("0 2", "0 1"), lwd = 2, layout = c(1, 2), strip = strip.custom( factor.levels = c("No pneumomia: Discharge", "No pneumomia: Death"), par.strip.text = list(font = 2)), ylim = c(0, 9), xlim = c(0, 190), xlab = "Days", scales = list(alternating = 1, x = list(at = seq(0, 150, 50)), y = list(at = seq(0, 8, 2))))

Plot for those having pneumomia:

dessin.p <- xyplot(my.nelaal.p, tr.choice = c("0 2", "0 1"), lwd = 2, layout = c(1, 2), strip = strip.custom( factor.levels = c("Pneumomia: Discharge", "Pneumomia: Death"), par.strip.text = list(font = 2)), ylab = "", ylim = c(0, 9), xlim = c(0, 190), xlab = "Days", scales = list(alternating = 1, x = list(at = seq(0, 150, 50)), y = list(at = seq(0, 8, 2))))

Display the 2 plots together.

print(dessin.nop, split = c(1, 1, 2, 1), more = TRUE, position = c(0, 0, 1.07, 1)) print(dessin.p, split = c(2, 1, 2, 1), position = c(-0.07, 0, 1, 1))

Cumulative incidence functions

Estimation of the cumulative incidence function using the cuminc function.

my.sir.cif <- cuminc(my.sir.data$time, my.sir.data$to, group = my.sir.data$pneu, cencode = "cens")

Using etm, stratified on pneumomia status.

my.sir.etm.nop <- etm(my.sir.data[my.sir.data$pneu == 0, ], c("0", "1", "2"), tra, "cens", s = 0) my.sir.etm.p <- etm(my.sir.data[my.sir.data$pneu == 1, ], c("0", "1", "2"), tra, "cens", s = 0)

Figure 4.9

op <- par(mfrow = c(1, 2)) ## Death plot(my.sir.etm.nop, tr.choice = "0 1", conf.int = FALSE, lwd = 2, lty = 1, xlab = "Days", ylab = "Probability", bty = "n", legend = FALSE) lines(my.sir.etm.p, tr.choice = "0 1", conf.int = FALSE, lwd = 2, lty = 2) legend(0, 0.6, c("No pneumonia", "Pneumonia"), col = 1, lty = c(1, 2), bty = "n", lwd = 2) mtext("Death", cex = 1.1, font = 2, line = 2) axis(1, at = seq(0, 200, 50)) ##Discharge plot(my.sir.etm.nop, tr.choice = "0 2", conf.int = FALSE, lwd = 2, lty = 1, xlab = "Days", ylab = "Probability", bty = "n", legend = FALSE) lines(my.sir.etm.p, tr.choice = "0 2", conf.int = FALSE, lwd = 2, lty = 2) axis(1, at = seq(0, 200, 50)) mtext("Discharge", cex = 1.1, font = 2, line = 2) par(op)

Figure 4.10

plot(my.sir.etm.p, tr.choice = '0 1', col = 1, lwd = 2, conf.int = TRUE, ci.fun = "cloglog", legend = FALSE, ylab="Probability", xlim=c(0,190)) lines(my.sir.etm.nop, tr.choice = '0 1', col = "gray", lwd = 2, conf.int = TRUE, ci.fun = "cloglog")

Section 4.4

Load the data.

data(abortion)

Figure 4.11

yo.demp <- ecdf(abortion$entry[abortion$group==1]) yo.demp2 <- ecdf(abortion$entry[abortion$group==0]) split.screen(figs = matrix(c(c(0, 0.495), c(0.505, 1), c(0,0), c(1,1)), ncol = 4)) screen(1) plot(x = c(0, 40), y = c(0, 1), xlab = "Week of pregnancy", ylab="Empirical distribution of study entry times", type = "n", axes = FALSE, main = "Exposed") axis(1, at = seq(0, 40, 5)) axis(2, at = seq(0, 1, 0.2)) box() lines(x = c(0, sort(unique(abortion$entry[abortion$group==1]))), y = c(0, yo.demp(sort(unique(abortion$entry[abortion$group==1])))), type="s") screen(2) plot(x = c(0, 40), y = c(0, 1), xlab = "Week of pregnancy", ylab="Empirical distribution of study entry times", type = "n", axes = FALSE, main = "Control") axis(1, at = seq(0, 40, 5)) axis(2, at = seq(0, 1, 0.2)) box() lines(x = c(0,sort(unique(abortion$entry[abortion$group==0]))), y = c(0,yo.demp2(sort(unique(abortion$entry[abortion$group==0])))), type="s") close.screen(all.screens = TRUE)

CIFs

First, you have to modify the data set.

my.abort <- abortion my.abort$from <- rep(0, nrow(my.abort)) names(my.abort)[5] <- c("to") # rename cause

Matrix of logicals defining the possible transitions.

tra <- matrix(FALSE, nrow = 4, ncol = 4) tra[1, 2:4] <- TRUE

Call to etm for computing the CIFs, stratified on treatment.

## Control my.abort.etm.nocd <- etm(my.abort[my.abort$group == 0,], c("0", "1", "2", "3"), tra, NULL, s = 0) ## exposed my.abort.etm.cd <- etm(my.abort[my.abort$group == 1,], c("0", "1","2","3"), tra, NULL, s = 0)

Figure 4.12

split.screen(figs = matrix(c(c(0, 0.495), c(0.505, 1), c(0, 0), c(1, 1)), ncol=4)) ## Exposed screen(1) plot(x = c(0, 40), y = c(0, 0.4), xlab="Week of pregnancy", ylab="Cumulative incidence function", type="n", axes=FALSE, main="Exposed") axis(1, at=seq(0, 40, 5)) axis(2, at=seq(0, 0.4, 0.05)) box() lines(x = my.abort.etm.cd$time, y = my.abort.etm.cd$est[1, 2,], type="s", lwd = 2, lty = 2)#induced lines(x = my.abort.etm.cd$time, y = my.abort.etm.cd$est[1, 4,], type="s", lwd = 2)#spont ## Control screen(2) plot(x = c(0, 40), y = c(0, 0.4), xlab = "Week of pregnancy", ylab="Cumulative incidence function", type="n", axes=FALSE, main="Control") axis(1, at=seq(0, 40, 5)) axis(2, at=seq(0, 0.4, 0.05)) box() lines(x = my.abort.etm.nocd$time, y = my.abort.etm.nocd$est[1, 2,], type="s", lwd=2, lty=2)#induced lines(x = my.abort.etm.nocd$time, y = my.abort.etm.nocd$est[1, 4,], type="s", lwd=2)#spont legend(10, 0.3, lty=1:2, legend=c("spontaneous", "induced"), bty = "n") close.screen(all.screens = TRUE)

Institute of Medical Biometry and Statistics (IMBI)