Sie sind hier: R code for chapter 9

# R code for chapter 9

Section 9.2

Section 9.2.1

Section 9.2.2

## Section 9.2

### 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.

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

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