Institut für Medizinische Biometrie
und Statistik (IMBI)

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

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

• library(mvna) my.nelaal <- mvna(my.data, c("0", "1", "2"), tra, "cens")

• 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() }

• my.fit.surv <- survfit(Surv(time, to != "cens") ~ 1, data = my.data, conf.type = "log-log")

• plot(my.fit.surv, xlab = "Time", ylab = "Survival Probability", mark.time= FALSE, lwd = 2) curve(exp(-0.9*x), add = TRUE, lwd = 2)  

• tra.km <- matrix(FALSE, ncol=2, nrow=2) dimnames(tra.km) <- list(c("0", "1"), c("0", "1")) tra.km[1, 2] <- TRUE

• my.data.km <- my.data my.data.km$to <- ifelse(my.data.km$to == "cens", "cens", 1)

• library(etm) km.etm <- etm(my.data.km, c("0", "1"), tra.km, "cens", s = 0)

• lines(x = km.etm$time, y = km.etm$est[1, 1, ], type = "s", col = "red")

• require(cmprsk) my.cif <- cuminc(my.data$time, my.data$to, cencode = "cens")

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

• cif.etm <- etm(my.data, c("0", "1", "2"), tra, "cens", s = 0)

• cif.etm$est[, , 1]

• lt.times <- rgamma(100, shape= 0.5, rate = 2)

• 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, ]

• my.nelaal2 <- mvna(my.data2, c("0", "1", "2"), tra, "cens")

• 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() }

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

• my.fit.surv2 <- survfit(Surv(entry,exit, to != "cens") ~ 1, data = my.data2, conf.type = "log-log")

• cif.etm2 <- etm(my.data2, c("0", "1", "2"), tra, "cens", s = 0)

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

• 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

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

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

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

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

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

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

• data(abortion)

• 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

• tra <- matrix(FALSE, nrow = 4, ncol = 4) tra[1, 2:4] <- TRUE

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

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