You are here: Home R code for chapter 5

# R code for chapter 5

Section 5.2.2

Section 5.3.3

Section 5.4

Section 5.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,
curve(0.09/(x+1), from=0, to=80, lwd=2, lty=2,

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,

curve(exp(-0.825 * 0.09 * log(x+1) - 0.2 * 0.5 * 0.024 * x^2),

## 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),
curve(0.09/(x+1)/(0.09/(x+1) + 0.024 * x),

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

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

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