Zu den Inhalten springen

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