Zu den Inhalten springen

Institute of Medical Biometry and Statistics (IMBI)

HIER SPENDEN Hier spenden
Home › Institut für Medizinische Biometrie und Statistik Neu › englisch › Institute of Medical Biometry and Statistics › Educational Books › Competing Risks and Multistate Models with R › R code › R code for chapter 9

R code for chapter 9

Section 9.2

Section 9.2.1

Subsection: Combined endpoint analysis

For the next analysis, we have to transform the icu.pneu data included in the kmi package into a format which allows us to use the mvna and etm package, as the data set is primarily tailored for a cox analysis with a time-dependent covariate.

Load the data.

  • data(icu.pneu)
Then, we perfom the transformation.

  • my.icu.pneu <- icu.pneu my.icu.pneu <- my.icu.pneu[order(my.icu.pneu$id, my.icu.pneu$start), ] masque <- diff(my.icu.pneu$id) my.icu.pneu$from <- 0 my.icu.pneu$from[c(1, masque) == 0] <- 1 my.icu.pneu$to2 <- my.icu.pneu$event my.icu.pneu$to2[my.icu.pneu$status == 0] <- "cens" my.icu.pneu$to2[c(masque, 1) == 0] <- 1
Event indicator variable to for a first event analysis:

  • my.icu.pneu$to <- ifelse(my.icu.pneu$to2 %in% c(2, 3), 2, my.icu.pneu$to2)
Keep the variables of interest

  • my.icu.pneu <- my.icu.pneu[, c("id", "start", "stop", "from", "to", "to2", "age", "sex")] names(my.icu.pneu)[c(2, 3)] <- c("entry", "exit")
Creation of a matrix with logical entries defining the possible transitions:

  • tra.idm <- matrix(FALSE, 3, 3, dimnames = list(c(0, 1, 2), c(0, 1, 2))) tra.idm[1, 2:3] <- TRUE tra.idm[2, 3] <- TRUE
Computation of the cumulative transition hazards using mvna:

  • mvna.idm <- mvna(my.icu.pneu, c("0", "1", "2"), tra.idm, "cens")
Figure 9.2

  • xyplot(mvna.idm, tr.choice = c("0 1", "0 2", "1 2"), lwd = 3, strip=strip.custom(bg="white", factor.levels = c(expression(0 %->% 1), expression(0 %->% 2), expression(1 %->% 2)), par.strip.text = list(cex = 1.2)), ylim = c(0, 9), xlab = "Days", layout = c(3, 1), xlim = c(0, 100))
Estimation of the transition probabilities using the etm package:

  • etm.idm <- etm(my.icu.pneu, c("0", "1", "2"), tra.idm, "cens", s = 0)
Figure 9.3

  • plot(etm.idm, tr.choice = "0 1", conf.int = TRUE, lwd = 2, legend = FALSE, ylim = c(0, 0.1), xlim = c(0, 100), xlab = "Days", ci.fun = "cloglog")
Landmark analysis
Creation of the landmark time points:

  • time.points <- c(seq(3, 10, 1), 15)
Computation of transition probabilities with s = time.points

  • landmark.etm <- lapply(time.points, function(start.time) { etm(my.icu.pneu, c("0", "1", "2"), tra.idm, "cens", start.time) })
Figure 9.4

  • nf <- layout(matrix(1:9, ncol = 3, byrow = TRUE), width = rep(1, 9), height = rep(1, 9)) for (i in seq_along(time.points)) { plot(landmark.etm[[i]], tr.choice = "0 2", lwd = 2, legend = FALSE, main = paste(expression("s ="), time.points[i], sep = " "), xlim = c(0, 100), lty = 2, ylab = "Probability of end-of-stay", xlab = "Days") lines(landmark.etm[[i]], tr.choice = "1 2", lwd = 2, col = 1, lty = 1, xlim = c(0, 100)) if (i == 1) { legend(45, 0.4, c(expression(P["02"](s, t)), expression(P["12"](s, t))), col = 1, lty = c(2, 1), lwd = 2, bty = "n", cex = 1.3) } }

Subsection: Analysis of competing endpoints in a progressive model.

For the analysis of the progressive model, we have to transform the data again into a progressive model.

  • my.icu.pneu.prog <- my.icu.pneu masque <- c(1, diff(my.icu.pneu.prog$id)) to <- my.icu.pneu.prog$to2[masque != 0 & my.icu.pneu.prog$to2 != 1] to <- ifelse(to == "cens", "cens", as.numeric(to) + 2) my.icu.pneu.prog$to <- my.icu.pneu.prog$to2 my.icu.pneu.prog$to[masque != 0 & my.icu.pneu.prog$to2 != 1] <- to
Then, we can create the matrix of logical values defining the possible transitions.

  • tra.prog <- matrix(FALSE, 6, 6, dimnames = list(as.character(0:5), as.character(0:5))) tra.prog[1, c(2, 5:6)] <- TRUE tra.prog[2, 3:4] <- TRUE
Cumulative transitions hazards:

  • mvna.prog <- mvna(my.icu.pneu.prog, as.character(0:5), tra.prog, "cens")
Figure 9.6

  • nf <- layout(matrix(c(1, 2), nrow = 1), width = c(1, 1), height = c(1, 1)) oldpar <- par(no.readonly = TRUE) ## Neslon Aalen estimator of the cumulative transition hazard of discharge par(mar = c(5, 4, 4, 0)) plot(mvna.prog, tr.choice = "0 5", conf.int = FALSE, legend = FALSE, lty = 2, lwd = 2, xlab = "Days", ylab = "Cumulative transition hazard", main = "Discharge", xlim = c(0, 100), ylim = c(0, 5)) lines(mvna.prog, tr.choice = "1 3", conf.int = FALSE, lty = 1, lwd = 2) legend("bottomright", c("Infected", "Non-infected"), col = 1, lty = c(1, 2), bty = "n", cex = 1.3, lwd = 2) ## Neslon Aalen estimator of the cumulative transition hazard of death par(mar = c(5, 2, 4, 2)) plot(mvna.prog, tr.choice = "0 4", conf.int = FALSE, legend = FALSE, lty = 2, lwd = 2, xlab = "Days", ylab = "", main = "Death", xlim = c(0, 100), ylim = c(0, 5), yaxt = "n") lines(mvna.prog, tr.choice = "1 2", conf.int = FALSE, lty = 1, lwd = 2) par(oldpar)
Computation of matrix of transition probabilities:

  • etm.prog <- etm(my.icu.pneu.prog, as.character(0:5), tra.prog, "cens", 0)
Figure 9.7

  • xyplot(etm.prog, tr.choice = c("0 2", "0 3", "0 4", "0 5"), strip = strip.custom(bg = "white", factor.levels = c(expression(P["02"](0, t)), expression(P["03"](0, t)), expression(P["04"](0, t)), expression(P["05"](0, t)))), xlim = c(-5, 100), par.strip.text = list(cex = 1.1, fontface = 2), xlab = "Days")

Section 9.2.2

A little modification of the data to avoid entry times equal to exit times, which throws an error

  • data(sir.cont) sir.cont <- sir.cont[order(sir.cont$id, sir.cont$time), ] for (i in 2:nrow(sir.cont)) { if (sir.cont$id[i]==sir.cont$id[i-1]) { if (sir.cont$time[i]==sir.cont$time[i-1]) { sir.cont$time[i-1] <- sir.cont$time[i-1] - 0.5 } } }
Definition of matrix of logical values specifying possible transitions:

  • tra.ventil <- matrix(FALSE, 3, 3, dimnames = list(c("0", "1", "2"), c("0", "1", "2"))) tra.ventil[1, c(2, 3)] <- TRUE tra.ventil[2, c(1, 3)] <- TRUE
Estimation of cumulative transition hazards:

  • mvna.ventil <- mvna(sir.cont, c("0", "1", "2"), tra.ventil, "cens")
Figure 9.8

  • xyplot(mvna.ventil, tr.choice = c("0 2", "1 2"), aspect = 1, strip = strip.custom(bg = "white", factor.levels = c("No ventilation -- Discharge/Death", "Ventilation -- Discharge/Death"), par.strip.text = list(cex = 1.1)), scales = list(alternating = 1, xlab = "Days", ylab = "Nelson-Aalen estimates"), xlim = c(-5, 100))