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 5

R code for chapter 5

Section 5.2.2

Subsection Model specification, motivation and analysis of the model specification

Figure 5.1

  • split.screen(figs = c(1, 2)) screen(1) plot(x = c(0, 70), y = c(0, 0.06), xlab=expression(paste(Time, " ", italic(t))), ylab="Cause-specific infection hazard", type="n", axes=FALSE, main="", cex.main=1.7, cex.lab=1.5) axis(1, at=seq(0, 70, 10), cex.axis=1.5) axis(2, at=seq(0, 0.06, 0.01), cex.axis=1.25) box() curve(0.825 * 0.09/(x+1), from=0, to=80, lwd=2,lty=1, add=TRUE,n=8000) curve(0.09/(x+1), from=0, to=80, lwd=2, lty=2, add=TRUE,n=8000) legend(20, 0.02, c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty=2:1,bty="n", cex=1.2, lwd=2) screen(2) plot(x=c(0, 70), y=c(0, 0.06), xlab=expression(paste(Time, " ", italic(t))), ylab="Cause-specific end-of-neutropenia hazard", type="n", axes=FALSE, main="", cex.main=1.7, cex.lab=1.5) axis(1, at=seq(0, 70, 10), cex.axis=1.5) axis(2, at=seq(0, 0.06, 0.01), cex.axis=1.25) box() curve(0.024 * x, from=0, to=45, type="l", lwd=2, lty=2, add=TRUE,n=4000) curve(0.2 * 0.024 * x, from=0, to=45, type="l", lwd=2, add=TRUE,n=4000) legend(20,0.02,c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty=2:1,bty="n", cex=1.2, lwd=2) close.screen(all.screens=TRUE)
Figure 5.2

  • split.screen(figs = c(1,2)) screen(1) plot(x=c(0, 70), y=c(0, 1), xlab=expression(paste(Time, " ", italic(t))), ylab="One minus waiting time distribution", type="n", axes=FALSE, main="", cex.main=1.7, cex.lab=1.5) axis(1, at=seq(0, 70, 10), cex.axis=1.5) axis(2, at=seq(0, 1, 0.1), cex.axis=1.25) box() legend(30,0.7,c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty=2:1,bty="n", cex=1.2, lwd=2) curve(exp(-0.09 * log(x+1) - 0.5 * 0.024 * x^2),from=0, to=80,lwd=2,lty=2,add=TRUE,n=8000) curve(exp(-0.825 * 0.09 * log(x+1) - 0.2 * 0.5 * 0.024 * x^2), from=0,to=80,lwd=2,lty=1,add=TRUE,n=8000) ## Binomial probability screen(2) plot(x=c(0, 70), y=c(0, 1), xlab=expression(paste(Time, " ", italic(t))), ylab="Binomial probability", type="n", axes=FALSE, main="", cex.main=1.7, cex.lab=1.5) axis(1, at=seq(0, 70, 10), cex.axis=1.5) axis(2, at=seq(0, 1, 0.1), cex.axis=1.25) box() curve(0.825 * 0.09/(x+1)/(0.825 * 0.09/(x+1) + 0.2 * 0.024 * x), from=0, to=80, lwd=2,lty=1, add=TRUE,n=8000) curve(0.09/(x+1)/(0.09/(x+1) + 0.024 * x), from=0, to=80, lwd=2,lty=2, add=TRUE,n=8000) legend(30,0.7,c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty=2:1,bty="n", cex=1.2, lwd=2) close.screen(all.screens=TRUE
Figure 5.3

  • split.screen(figs = c(1,2)) screen(1) plot(c(0, 70), c(0, 1), xlab = expression(paste(Time, " ", italic(t))), ylab = "", type = "n", axes = FALSE, cex.lab=1.5) axis(1, at = seq(0, 70, 10), cex.axis = 1.5) axis(2, at = seq(0, 1, 0.1), cex.axis = 1.25) legend(0, 0.4, c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty = 2:1, bty = "n", cex = 1.2, lwd = 2) box() mtext(text="CIF for infection", side = 2, line = 3, cex = 1.25) curve(0.01 * cumsum(exp(-(0.09 * log(x+1) + 0.5 * 0.024 * x^2))* 0.09/(x+1)), from = 0,to = 80, lwd = 2, lty = 2, add = TRUE, n = 8000) curve(0.01 * cumsum(exp(-(0.825 * 0.09 * log(x+1) + 0.2 * 0.5 * 0.024 * x^2))* 0.825 * 0.09/(x+1)), from = 0, to = 80, lwd = 2, lty = 1, add = TRUE, n = 8000) ## Plot 2 screen(2) plot(c(0, 70), c(0, 1), xlab = expression(paste(Time, " ", italic(t))), ylab = "", type = "n", axes = FALSE, cex.lab=1.5) axis(1, at = seq(0, 70, 10), cex.axis = 1.5) axis(2, at = seq(0, 1, 0.1), cex.axis = 1.25) box() mtext(text="CIF for end-of-neutropenia", side = 2, line = 3, cex = 1.25) box() legend(0, 1, c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty = 2:1, bty = "n", cex = 1.2, lwd = 2) curve(0.01 * cumsum(exp(-(0.09 * log(x+1) + 0.5 * 0.024 * x^2))* 0.024 * x), from = 0,to = 80, lwd = 2, lty = 2, add = TRUE, n = 8000) curve(0.01 * cumsum(exp(-(0.825 * 0.09 * log(x+1) + 0.2 * 0.5 * 0.024 * x^2))* 0.2 * 0.024 * x), from = 0, to = 80, lwd = 2, lty = 1, add = TRUE, n = 8000) close.screen(all.screens=TRUE)
Figure 5.4

  • plot(c(0, 30), c(0, 0.2), xlab = expression(paste(Time, " ", italic(t))), ylab = "", type = "n", axes = FALSE, cex.lab=1.5) axis(1, at = seq(0, 30, 5), cex.axis = 1.5) axis(2, at = seq(0, 0.2, 0.05), cex.axis = 1.25) legend(0, 0.4, c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty = 2:1, bty = "n", cex = 1.2, lwd = 2) box() mtext(text="CIF for infection", side = 2, line = 3, cex = 1.25) curve(0.01 * cumsum(exp(-(0.09 * log(x+1) + 0.5 * 0.024 * x^2))* 0.09/(x+1)), from = 0,to = 80, lwd = 2, lty = 2, add = TRUE, n = 8000) curve(0.01 * cumsum(exp(-(0.825 * 0.09 * log(x+1) + 0.2 * 0.5 * 0.024 * x^2))* 0.825 * 0.09/(x+1)), from = 0, to = 80, lwd = 2, lty = 1, add = TRUE, n = 8000)
Figure 5.5

  • split.screen(figs=c(1,2)) screen(1) plot(c(0, 70), c(0, 1), xlab = expression(paste(Time, " ", italic(t))), ylab = "", type = "n", axes = FALSE, cex.lab = 1.5) axis(1, at = seq(0, 70, 10), cex.axis = 1.5) axis(2, at = seq(0, 1, 0.1), cex.axis = 1.25) box() mtext(text="CIF for end-of-neutropenia", side = 2, line = 3, cex = 1.25) box() legend(0, 1, c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty = 2:1, bty = "n", cex = 1.2, lwd = 2) curve(0.01 * cumsum(exp(-(0.09 * log(x+1) + 0.5 * 0.024 * x^2))* 0.024 * x), from = 0, to = 80, lwd = 2, lty = 2, add = TRUE, n = 8000) curve(0.01 * cumsum(exp(-(0.825 * 0.09 * log(x+1) + 0.7 * 0.5 * 0.024 * x^2))* 0.7 * 0.024 * x), from = 0, to = 80, lwd = 2, lty = 1, add = TRUE, n = 8000) text(35, 0.4, c(expression(paste(plain(exp), "(",beta[paste(0,1)],")"," =0.825"))), cex = 1.2) text(33.0, 0.3, c(expression(paste(plain(exp), "(",beta[paste(0,2)],")"," =0.7"))), cex = 1.2) screen(2) plot(c(0, 70), c(0, 1), xlab = expression(paste(Time, " ", italic(t))), ylab = "", type = "n", axes = FALSE, cex.lab = 1.5) axis(1, at = seq(0, 70, 10), cex.axis = 1.5) axis(2, at = seq(0, 1, 0.1), cex.axis = 1.25) box() mtext(text="CIF for end-of-neutropenia", side = 2, line = 3, cex = 1.25) box() legend(0, 1, c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), lty = 2:1, bty = "n", cex = 1.2, lwd = 2) curve(0.01 * cumsum(exp(-(0.09 * log(x+1) + 0.5 * 0.024 * x^2)) * 0.024 * x), from = 0,to = 80, lwd = 2, lty = 2, add = TRUE, n = 8000) curve(0.01 * cumsum(exp(-(0.825 * 0.09 * log(x+1) + 0.825 * 0.5 * 0.024 * x^2)) * 0.825 * 0.024 * x), from = 0, to = 80, lwd = 2, lty = 1, add = TRUE, n = 8000) text(35, 0.4, c(expression(paste(plain(exp), "(",beta[paste(0,1)],")"," =0.825"))), cex = 1.2) text(35.0, 0.3, c(expression(paste(plain(exp), "(",beta[paste(0,2)],")"," =0.825"))), cex = 1.2) close.screen(all.screens=TRUE)

Subsection: Simulation and analysis of simulated data

What follows is the actual code for simulating the ONKO-KISS data. cum.haz.0* are functions defining the all-cause hazards for Z = 0 and Z = 1, respectively. They will be used for creating the event times via numerical inversion (hence the little y).

  • cum.haz.01 <- function(t, y) { 0.09 * log(t + 1)+ 0.024 * t * t/ 2 - y } cum.haz.02 <- function(t, y) { 0.825 * 0.09 * log(t + 1) + 0.2 * 0.024 * t * t / 2 - y }

Cause-specific hazards

Z=0:

  • alpha.0.01 <- function(t) { 0.09 / (t + 1) } alpha.0.02 <- function(t) { 0.024 * t }
Z=1 :

  • alpha.1.01 <- function(t) { 0.825 * alpha.0.01(t) } alpha.1.02 <- function(t) { 0.2*alpha.0.02(t) }
create.times is a function to generate the event times

  • create.times <- function(n, ch, sup.int = 200) { times <- numeric(n) i <- 1 while (i <= n) { u <- runif(1) if (ch(0, -log(u)) * ch(sup.int, -log(u)) < 0) { times[i] <- uniroot(ch, c(0, sup.int), tol = 0.0001, y= -log(u))$root i <- i + 1 } else { cat("pos") } } times }
Binomial experiments to decide on the event type.

  • binom.status <- function(ftime, n, a01, a02, size = 1) { prob <- a01(ftime) / (a01(ftime) + a02(ftime)) out <- rbinom(n, size, prob) out }
Next, a function for generating the data.

  • gen.data <- function(n, PROB) { n1 <- n * PROB n0 <- n - n1 t0 <- create.times(n0, cum.haz.01) t1 <- create.times(n1, cum.haz.02) fs0 <- binom.status(t0, n0, alpha.0.01, alpha.0.02) fs1 <- binom.status(t1, n1, alpha.1.01, alpha.1.02) cov <- c(rep(0, n0), rep(1, n1)) ft <- c(t0, t1) fs <- c(fs0, fs1) fs <- (-1) * fs + 2 matrix(c(ft, fs, cov), ncol=3) }
Create the data:

  • set.seed(261339) y <- gen.data(n = 1500, PROB = 0.5) x <- data.frame(id = 1:length(y[, 1]), T = y[, 1], X.T = y[, 2], Z = y[, 3])
and add light censoring:

  • cens <- runif(length(x$T), max = 100) x$C <- cens x$TandC <- pmin(x$T, x$C) x$status <- as.numeric(x$T <= x$C) * x$X.T

First-event analysis

Cox model for the all-cause hazard:

  • fit <- coxph(Surv(TandC, status != 0) ~ Z, data = x) sfit <- summary(fit) sfit
A quick check of the confidence intervals:

  • exp(fit$coefficients + qnorm(0.975) *sqrt(fit$var)) exp(fit$coefficients - qnorm(0.975) *sqrt(fit$var))
Analysis of the cause-specific hazards by fitting two separate Cox models

  • fit01 <- coxph(Surv(TandC, status == 1) ~ Z, data = x) fit02 <- coxph(Surv(TandC, status == 2) ~ Z, data = x

Analysis of the duplicated data set

First, we create the duplicated data set.

  • xl <- rbind(x,x) xl$eventtype <- c(rep("interest",1500), rep("competing", 1500)) xl$newstat <- as.numeric(c(x$status == 1, x$status == 2))
1st cox model: first event analysis:

  • summary(coxph(Surv(TandC, newstat != 0) ~ Z, data = xl))
2nd cox model: first event analysis using a strata term:

  • summary(coxph(Surv(TandC, newstat != 0) ~ Z + strata(eventtype), data = xl))
Creation of the cause-specific covariates:

  • xl$Z.01 <- xl$Z * (xl$eventtype == "interest") xl$Z.02 <- xl$Z * (xl$eventtype == "competing")
Cox cause-specific hazards analysis with the duplicated data set:

  • summary(coxph(Surv(TandC, newstat != 0) ~ Z.01 + Z.02 + strata(eventtype), data = xl))

Modelling of common effects

Now, we create a random variable which has nothing in common with the way the data were simulated.

  • cv <- round(0.5 * rnorm(length(x$id),1)) xl$cv <- rep(cv, 2)
cox model with a common effect on both CSHs for cv

  • ## event of interest summary(coxph(Surv(TandC, status == 1) ~ Z + cv, data = x)) ## competing event summary(coxph(Surv(TandC, status == 2) ~ Z + cv,data = x))
Duplicated data set (both α01 and α02)

  • summary(coxph(Surv(TandC, newstat != 0) ~ Z.01 + Z.02 + cv + strata(eventtype), data = xl))

Breslow estimators, Nelson-Aalen estimators and interpretation

Breslow estimates of the baseline hazards:

  • a01.0 <- basehaz(coxph(Surv(TandC, status == 1) ~ Z,data = x), centered = FALSE) # csh01 a02.0 <- basehaz(coxph(Surv(TandC, status == 2) ~ Z, data = x), centered = FALSE) # csh02
Figure 5.6

  • split.screen(figs=c(1,2)) screen(1) plot(c(0, 10), c(0, 1), xlab = expression(paste(Time, " ", italic(t))), ylab= "", type = "n", axes = FALSE, cex.lab=1.5) axis(1, at = seq(0, 10, 5), cex.axis = 1.5) axis(2, at=seq(0, 1, 0.2), cex.axis = 1.25) box() lines(x=a02.0$time, y=a02.0$hazard, type="s", lwd=2) lines(x=a01.0$time, y=a01.0$hazard, type="s", lwd=2, col="darkgrey") legend(0,1, c(expression(paste(A[0], phantom()[1], phantom()[scriptstyle(";")], phantom()[0],"(",italic(t),")")), expression(paste(A[0],phantom()[2],phantom()[scriptstyle(";")], phantom()[0],"(",italic(t),")"))), lty=c(1,1),col=c("darkgrey", "black"),bty="n", cex=1.2, lwd=2) screen(2) plot(c(0, 70), y=c(0, 35), xlab = expression(paste(Time, " ", italic(t))), ylab="", type="n", axes=FALSE, cex.lab=1.5) axis(1, at=seq(0, 70, 10), cex.axis=1.5) axis(2, at=seq(0, 35, 5), cex.axis=1.25) box() lines(x=a02.0$time, y=a02.0$hazard, type="s", lwd=2) lines(x=a01.0$time, y=a01.0$hazard, type="s", lwd=2, col="darkgrey") legend(0,35,c(expression(paste(A[0],phantom()[1],phantom()[scriptstyle(";")], phantom()[0],"(",italic(t),")")), expression(paste(A[0],phantom()[2],phantom()[scriptstyle(";")], phantom()[0],"(",italic(t),")"))), lty=c(1,1),col=c("darkgrey", "black"),bty="n", cex=1.2, lwd=2) close.screen(all.screens=TRUE) mtext("Breslow estimates of the cumulative cause-specific hazards", cex=1.5, line=1, font = 2)
Baseline hazards using the duplicated data set:

  • basehaz(coxph(Surv(TandC, newstat != 0) ~ Z.01 + Z.02 + strata(eventtype), data = xl), centered = FALSE)
Figure 5.7

  • ## matrix of logical describing the possible transitions (for using mvna) tra <- matrix(FALSE, ncol = 3, nrow = 3) dimnames(tra) <- list(c("0", "1", "2"), c("0", "1", "2")) tra[1, 2:3] <- TRUE ## Tranformation of the data set into mvna format x.mvna <- x x.mvna$from <- rep(0, nrow(x)) x.mvna$to <- ifelse(x.mvna$status == 0, "cens", x.mvna$status) x.mvna$time <- x$TandC ## cumulative CSHs na.z0 <- mvna(x.mvna[x.mvna$Z == 0, ], c("0", "1", "2"), tra, "cens") na.z1 <- mvna(x.mvna[x.mvna$Z == 1, ], c("0", "1", "2"), tra, "cens") plot(x=c(0, 10), y=c(0, 1), xlab=expression(paste(Time, " ", italic(t))), ylab="", type="n", axes=F, main="", cex.main=1.7, cex.lab=1.5) axis(1, at=seq(0, 10, 5), cex.axis=1.5) axis(2, at=seq(0, 1, 0.2), cex.axis=1.25) box() lines(na.z0[["0 1"]][,"time"], na.z0[["0 1"]][,"na"], type="s", col="darkgrey", lwd=2) lines(na.z0[["0 2"]][,"time"], na.z0[["0 2"]][,"na"], type="s", lwd=2) lines(na.z1[["0 1"]][,"time"], na.z1[["0 1"]][,"na"], type="s", col="darkgrey", lwd=2, lty=2) lines(na.z1[["0 2"]][,"time"], na.z1[["0 2"]][,"na"], type="s", lwd=2, lty=2) legend(0,0.989,lty=c(1,1), rep(expression(paste(phantom("Event of interest"))),2), bty="n", col=c("darkgrey","black"),lwd=2) legend(1,1,lty=c(2,2), c("Event of interest","Competing event"), bty="n", col=c("darkgrey","black"),lwd=2, cex=1.2) legend(0,0.789,lty=c(2,1), rep(expression(paste(phantom("Event of interest"))),2), bty="n", col=c("darkgrey","darkgrey"),lwd=2) legend(1,0.8,lty=c(2,1), c(expression(paste(italic(Z)[i]," =0")), expression(paste(italic(Z)[i]," =1"))), bty="n", lwd=2, cex=1.2)

Prediction of the cumulative incidence functions

Prediction of the CIFs is performed using the mstate package, using the extended data set xl.

  • require(mstate) ## matrix indicating the possible transitions tmat <- trans.comprisk(2)
Add a variable trans that indicates the transition type:

  • xl$trans <- c(rep(1, 1500), rep(2, 1500))
Next, we fit the cox csh model with the "Breslow" method for handling ties.

  • fit <- coxph(Surv(TandC, newstat) ~ Z.01 + Z.02 + strata(trans), data = xl, method = "breslow")

New data frames for making predictions:

An individual with Z=0:

  • newdat.z0 <- data.frame(Z.01 = c(0, 0), Z.02 = c(0, 0), strata = c(1, 2))
An individual with Z=1:

  • newdat.z1 <- data.frame(Z.01 = c(1, 0), Z.02 = c(0, 1), strata = c(1, 2))
Predicted cumulative hazards:

  • msf.z0 <- msfit(fit, newdat.z0, trans = tmat) msf.z1 <- msfit(fit, newdat.z1, trans = tmat)
Predicted CIFs:

  • pt.z0 <- probtrans(msf.z0, 0)[[1]] pt.z1 <- probtrans(msf.z1, 0)[[1]]
Figure 5.8

  • ## modification of the data for using etm() x.etm <- data.frame(id = x$id, from = rep(11, nrow(x)), to = x$status, time = x$TandC, Z = x$Z) aa0 <- etm(x.etm[x.etm$Z == 0, ], c('11', '1', '2'), matrix(c(FALSE, TRUE, TRUE, rep(FALSE, 6)), ncol = 3, nrow = 3, byrow = TRUE), '0', 0) aa1 <- etm(x.etm[x.etm$Z == 1, ], c('11', '1', '2'), matrix(c(FALSE, TRUE, TRUE, rep(FALSE, 6)), ncol = 3, nrow = 3, byrow = TRUE), '0', 0) ## plot oldpar <- par(mfrow = c(1, 2)) plot(aa0, tr.choice = "11 1", col = "gray", lty = 1, legend = FALSE, lwd = 3, xlim = c(0, 30), ylab = "Cumulative Incidence") lines(aa1, tr.choice = "11 1", col = 1, lty = 1, lwd = 3) lines(pt.z0$time, pt.z0$pstate2, col = "gray", lty = 2, lwd = 3) lines(pt.z1$time, pt.z1$pstate2, col = 1, lty = 2, lwd = 3) legend(x = 0, y = 0.83, c("Z = 0", "Z = 1"), bty = "n", col = c("gray", 1), lty = 1, lwd = 3) legend(x = 0, y = 0.7, c("Aalen-Johansen", "Predicted"), col = 1, lty = c(1, 2), bty = "n", lwd = 3) plot(aa0, tr.choice = "11 2", col = "gray", lty = 1, legend = FALSE, lwd = 3, xlim = c(0, 30), ylab = "Cumulative Incidence") lines(aa1, tr.choice = "11 2", col = 1, lty = 1, lwd = 3) lines(pt.z0$time, pt.z0$pstate3, col = "gray", lty = 2, lwd = 3) lines(pt.z1$time, pt.z1$pstate3, col = 1, lty = 2, lwd = 3) par(oldpar)
Figure 5.9

  • plot(aa0, tr.choice = "11 1", col = "gray", lty = 1, legend = FALSE, lwd = 3, xlim = c(0, 30), ylab = "Cumulative Incidence", ylim = c(0, 0.2)) lines(aa1, tr.choice = "11 1", col = 1, lty = 1, lwd = 3) lines(pt.z0$time, pt.z0$pstate2, col = "gray", lty = 2, lwd = 3) lines(pt.z1$time, pt.z1$pstate2, col = 1, lty = 2, lwd = 3) legend(x = 15, y = 0.1, c("Z = 0", "Z = 1"), bty = "n", col = c("gray", 1), lty = 1, lwd = 3) legend(x = 15, y = 0.05, c("Aalen-Johansen", "Predicted"), col = 1, lty = c(1, 2), bty = "n", lwd = 3)

Subsection: Analysis of hospital data: Impact of pneumonia status on admission on intensive care unit mortality

Cox models

  • fit.pneu.01 <- coxph(Surv(time, to == 1) ~ pneu, my.sir.data) fit.pneu.02 <- coxph(Surv(time, to == 2) ~ pneu, my.sir.data) summary(fit.pneu.02) summary(fit.pneu.02)

Figure 5.10

Baseline hazards:

  • a01.0 <- basehaz(fit.pneu.01, centered=FALSE) a02.0 <- basehaz(fit.pneu.02, centered=FALSE
Plot:

  • split.screen(figs=c(1,2)) screen(1) plot(c(0, 50), c(0, 5), xlab = expression(paste(Time, " ", italic(t))), ylab = "Cumulative cause-specific hazard", type = "n", axes = FALSE, main = "No pneumonia", cex.main = 1.5, cex.lab = 1.5) axis(1, at=seq(0, 50, 10), cex.axis=1.5) axis(2, at=seq(0, 5, 1), cex.axis=1.25) box() lines(a02.0$time, a02.0$hazard, type="s", lwd=2, lty=2) lines(a01.0$time, a01.0$hazard, type="s", lwd=2) lines(my.nelaal.nop, conf.int = FALSE, col = rep("darkgray", 2), lty = c(1, 2), lwd = 2) legend(0,5,c("Death", "Discharge"), lty=1:2,bty="n", cex=1.2, lwd=2) screen(2) plot(x=c(0, 50), y=c(0, 5), xlab=expression(paste(Time, " ", italic(t))), ylab="Cumulative cause-specific hazard", type="n", axes=F, main="Pneumonia", cex.main=1.5, cex.lab=1.5) axis(1, at=seq(0, 50, 10), cex.axis=1.5) axis(2, at=seq(0, 5, 1), cex.axis=1.25) box() lines(a02.0$time, hr2 * a02.0$hazard, type="s", lwd=2, lty=2) lines(a01.0$time, hr1 * a01.0$hazard, type="s", lwd=2) lines(my.nelaal.p, conf.int = FALSE, col = rep("darkgray", 2), lty = c(1, 2), lwd = 2) legend(0,5,c("Death", "Discharge"), lty=1:2,bty="n", cex=1.2, lwd=2) close.screen(all.screens=TRUE)

Subsection: Analysis of pregnancy outcome data: Impact of coumarin derivatives on spontaneous and induced abortion

Proportional cause-specific hazards models

  • coxph(Surv(entry,exit, cause == 1) ~ group, data = abortion)# induced coxph(Surv(entry,exit, cause == 3) ~ group, data = abortion)# spontaneous coxph(Surv(entry,exit, cause == 2) ~ group, data = abortion) # life birth

Section 5.3.3

Subsection: Simulated data

Using standard Cox software for administratively censored data

Coding of the censored subdistribution failure time:

  • x$thetaandC <- ifelse(x$status==2, x$C, x$TandC)
With coxph and the administrative censoring times

  • summary(coxph(Surv(thetaandC,status==1)~Z,data=x))
Using the cmprsk package

  • fit.sh <- crr(ftime = x$TandC, fstatus = x$status, cov1 = x$Z, failcode = 1, cencode = 0)

Figure 5.13

Predicted CIFs from the proportional subdistribution hazards model:

  • daddeln <- predict.crr(fit.sh, cov1 = matrix(c(0, 1), nrow = 2))
Plot:

  • split.screen(figs=c(1,2)) screen(1) plot(c(0, 25), c(0, 0.2), xlab = expression(paste(Time, " ", italic(t))), ylab = "", type = "n", axes=FALSE, main = expression(paste(italic(Z)[i]," =0")), cex.main = 1.5, cex.lab = 1.5) axis(1, at = seq(0, 25, 5), cex.axis = 1.5) axis(2, at = seq(0, 0.2, 0.05), cex.axis = 1.25) box() mtext(text="CIF for infection", side = 2, line = 3, cex = 1.25) lines(cif[["0 1"]]$time, cif[["0 1"]]$est, type = "s", lwd = 2, lty = 1) ## predicted curves lines(daddeln[,1], daddeln[,2], type = "s", lwd = 2, lty = 1, col = "darkgrey") screen(2) plot(c(0, 25), c(0, 0.2), xlab = expression(paste(Time, " ", italic(t))), ylab = "", type = "n", axes = FALSE, main = expression(paste(italic(Z)[i]," =1")), cex.main = 1.5, cex.lab = 1.5) axis(1, at = seq(0, 25, 5), cex.axis = 1.5) axis(2, at = seq(0, 0.2, 0.05), cex.axis = 1.25) box() mtext(text="CIF for infection", side = 2, line = 3, cex = 1.25) lines(cif[["1 1"]]$time, cif[["1 1"]]$est, type = "s", lwd = 2,lty = 1) ## predicted curves lines(daddeln[,1], daddeln[,3], type = "s", lwd = 2, col = "darkgrey") close.screen(all.screens=TRUE)
Using standard Cox software and multiple imputation through the kmi package

  • require(kmi) ## imputed data sets imputed.data <- kmi(Surv(TandC, status != 0) ~ 1, data = x, etype = status, failcode = 1, nimp = 10) ## Cox models on the imputed data sets fit.impu <- cox.kmi(Surv(TandC, status == 1) ~ Z, imputed.data) summary(fit.impu)
With bootstrap:

  • imputed.data.boot <- kmi(Surv(TandC, status != 0) ~ 1, data = x, etype = status, failcode = 1, nimp = 10, bootstrap = TRUE, nboot = 100) summary(cox.kmi(Surv(TandC, status == 1) ~ Z, imputed.data.boot))

Subsection: Analysis of hospital data: Impact of pneumonia status on admission on intensive care unit mortality

Proportional subdistribution hazards model:

  • crr(ftime = my.sir.data$time, fstatus = my.sir.data$to, cov1 = my.sir.data$pneu, failcode = "1", cencode = "cens")

Section 5.4

Censoring of the simulated data at time 10:

  • x$TandC.art <- ifelse(x$TandC <= 10, x$TandC, 10) x$status.art <- ifelse(x$TandC <= 10, x$status, 0)
Subdistribution hazards analysis:

  • fit.sh.art <- crr(ftime = x$TandC.art, fstatus = x$status.art, cov1 = x$Z, failcode = 1, cencode = 0)
cause-specific hazards analysis

  • summary(coxph(Surv(TandC.art, status.art == 1) ~ Z,data = x)) summary(coxph(Surv(TandC.art, status.art == 2) ~ Z, data = x))

Section 5.5

Figure 5.14

In the following, we need to be careful about the event times, hence the use of findInterval.

  • ## all-cause hazard fit00 <- coxph(Surv(TandC, status != 0) ~ Z, data = x) surv.o <- survfit(Surv(TandC, status != 0) ~ Z, data = x) ### Plot: old.par <- par(no.readonly = TRUE) nf <- layout(t(matrix(c(1, 2, 3))), width = c(1, 1, 1)) ## All cause hazard times <- sort(surv.o$time) all.na.z0 <- cumsum(surv.o[1]$n.event / surv.o[1]$n.risk) all.na.z1 <- cumsum(surv.o[2]$n.event / surv.o[2]$n.risk) ind.z0 <- findInterval(times, surv.o[1]$time) ind.z0[ind.z0 == 0] <- NA all.na.z0 <- all.na.z0[ind.z0] ind.z1 <- findInterval(times, surv.o[2]$time) ind.z1[ind.z1 == 0] <- NA all.na.z1 <- all.na.z1[ind.z1] plot(all.na.z0, all.na.z1, type = "s", lwd = 2, xlim = c(0, 6), ylim = c(0, 2), xlab = expression(hat(A)["0.;0"](t)), ylab = expression(hat(A)["0."](t, "Z = 1")), main = "All-cause hazard") abline(a = 0, b = exp(fit00$coef), col = "darkgray", lwd = 2) ## alpha01 times <- sort(c(na.z0$"0 1"$time, na.z1$"0 1"$time)) ind.z0 <- findInterval(times, na.z0$"0 1"$time) ind.z0[ind.z0 == 0] <- NA na01.z0 <- na.z0$"0 1"$na[ind.z0] ind.z1 <- findInterval(times, na.z1$"0 1"$time) ind.z1[ind.z1 == 0] <- NA na01.z1 <- na.z1$"0 1"$na[ind.z1] plot(na01.z0, na01.z1, type = "s", lwd = 2, xlab = expression(hat(A)["01;0"](t)), ylab = expression(hat(A)["01"](t, "Z = 1")), main = "Event of interest") abline(a = 0, b = exp(fit01$coef), col = "darkgray", lwd = 2) ## alpha 02 times <- sort(c(na.z0$"0 2"$time, na.z1$"0 2"$time)) ind.z0 <- findInterval(times, na.z0$"0 2"$time) ind.z0[ind.z0 == 0] <- NA na02.z0 <- na.z0$"0 2"$na[ind.z0] ind.z1 <- findInterval(times, na.z1$"0 2"$time) ind.z1[ind.z1 == 0] <- NA na02.z1 <- na.z1$"0 2"$na[ind.z1] plot(na02.z0, na02.z1, type = "s", lwd = 2, xlim = c(0, 5), ylim = c(0, 2), xlab = expression(hat(A)["02;0"](t)), ylab = expression(hat(A)["02"](t, "Z = 1")), main = "Competing Event") abline(a = 0, b = exp(fit02$coef), col = "darkgray", lwd = 2)

Figure 5.15

Computation of the cumulative subdistribution hazard:

  • cif.etm.z0 <- etm(x.mvna[x.mvna$Z == 0, ], c("0", "1", "2"), tra, "cens", 0) cif.etm.z1 <- etm(x.mvna[x.mvna$Z == 1, ], c("0", "1", "2"), tra, "cens", 0) times <- sort(c(cif.etm.z0$time, cif.etm.z1$time)) cif.z0 <- trprob(cif.etm.z0, tr.choice = "0 1", timepoints = times) cif.z1 <- trprob(cif.etm.z1, tr.choice = "0 1", timepoints = times) sub.haz.z0 <- cumsum(1 -((1 - cif.z0) / (1 - c(0, cif.z0[-length(cif.z0)])))) sub.haz.z1 <- cumsum(1 -((1 - cif.z1) / (1 - c(0, cif.z1[-length(cif.z1)]))))
Plot:

  • plot(sub.haz.z0, sub.haz.z1, lwd = 2, type = "s", xlab = expression(hat(Lambda)(t, "Z = 0")), ylab = expression(hat(Lambda)(t, "Z = 1"))) abline(a = 0, b = exp(fit.sh$coef), col = "darkgray", lwd = 2)