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

R code for chapter 4

Quicklinks for Chapter 4

Section 4.2

 

Section 4.3

Section 4.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 t1

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)