Sie sind hier: Startseite Competing Risks and … R Code R code for chapter 9

R code for chapter 9

Quicklinks for Chapter 9

Section 9.2

Section 9.2.1

Section 9.2.2

 

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))