You are here: Home R code for chapter 4

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