Sie sind hier: R code for chapter 10

# R code for chapter 10

Section 10.2

## Section 10.2

### Section 10.2.1

Cox proportional hazards model: first method

```cox.icu.pneu01 <- coxph(Surv(entry, exit, to == 1) ~ age +
sex, my.icu.pneu, subset = from == 0)
cox.icu.pneu02 <- coxph(Surv(entry, exit, to == 2) ~ age +
sex, my.icu.pneu, subset = from == 0)
cox.icu.pneu12 <- coxph(Surv(entry, exit, to == 2) ~ age +
sex, my.icu.pneu, subset = from == 1)

summary(cox.icu.pneu01)
summary(cox.icu.pneu02)
summary(cox.icu.pneu12)
```

For the second method, we need an extended data set and transition specific covariates.

```from0 <- subset(my.icu.pneu, from == 0)
my.icu.pneu.ext <- rbind(from0, from0)

my.icu.pneu.ext\$new.to <- c(rep(1, nrow(from0)), rep(2, nrow(from0)))
my.icu.pneu.ext\$new.status <- as.numeric(c(from0\$to == 1, from0\$to == 2))

from1 <- subset(my.icu.pneu, from == 1)
from1\$new.to <- 2
from1\$new.status <- as.numeric(from1\$to == 2)

my.icu.pneu.ext <- rbind(my.icu.pneu.ext, from1)
my.icu.pneu.ext\$trans <- c(rep(1, nrow(from0)), rep(2, nrow(from0)),
rep(3, nrow(from1)))
my.icu.pneu.ext <- my.icu.pneu.ext[, -c(5, 6)]

my.icu.pneu.ext\$age.1 <- with(my.icu.pneu.ext, age * (trans == 1))
my.icu.pneu.ext\$age.2 <- with(my.icu.pneu.ext, age * (trans == 2))
my.icu.pneu.ext\$age.3 <- with(my.icu.pneu.ext, age * (trans == 3))
my.icu.pneu.ext\$male.1 <- with(my.icu.pneu.ext, (sex == "M") * (trans == 1))
my.icu.pneu.ext\$male.2 <- with(my.icu.pneu.ext, (sex == "M") * (trans == 2))
my.icu.pneu.ext\$male.3 <- with(my.icu.pneu.ext, (sex == "M") * (trans == 3))
```

Cox model:

```fit.cox.icu.ext <- coxph(Surv(entry, exit, new.status) ~
age.1 + age.2 + age.3 +
male.1 + male.2 + male.3 +
strata(trans), my.icu.pneu.ext)
summary(fit.cox.icu.ext)
```

The next analyses will use the mstate package.

```## possible transitions
mat <- trans.illdeath()
```

Data sets for two hypothetical individuals:

```woman.medianage <- data.frame(age = rep(median.age, 3), male
= rep(0, 3), trans = 1:3)
attr(woman.medianage, "trans") <- mat
class(woman.medianage) <- c("msdata", "data.frame")
woman.medianage <- expand.covs(woman.medianage,
c("age", "male"))
woman.medianage\$strata <- 1:3

man.medianage <- data.frame(age = rep(median.age, 3), male
= rep(1, 3), trans = 1:3)
attr(man.medianage, "trans") <- mat
class(man.medianage) <- c("msdata", "data.frame")
man.medianage <- expand.covs(man.medianage, c("age", "male"))
man.medianage\$strata <- 1:3
```

Cox model with breslow's method for handling ties:

```fit.cox.icu.ext.bres <- coxph(Surv(entry, exit, new.status)  ~
age.1 + age.2 + age.3 +
male.1 + male.2 + male.3 +
strata(trans), my.icu.pneu.ext,
method = "breslow")
```

Predicted cumulative hazards using `msfit`

```msfit.woman <- msfit(fit.cox.icu.ext.bres, woman.medianage,
trans = mat)
msfit.man <- msfit(fit.cox.icu.ext.bres, man.medianage,
trans = mat)
```

Predicted transition probabilities:

```pt.woman <- probtrans(msfit.woman, 0)
pt.man <- probtrans(msfit.man, 0)
```

Figure 10.1

Computation of the confidence intervals of the predicted transition probabilities

```cis <- function(p, se) {
alpha <- qnorm(0.975)
lower <- 1 - (1 - p)^(exp(alpha * (se/((1 - p) * log(1 - p)))))
upper <- 1 - (1 - p)^(exp(-alpha * (se/((1 - p) * log(1 - p)))))
data.frame(lower, upper)
}

cis.woman <- cis(pt.woman[[1]]\$pstate2, pt.woman[[1]]\$se2)
cis.man <- cis(pt.man[[1]]\$pstate2, pt.man[[1]]\$se2)
```

Code for the figure:

```olpar <- par(no.readonly = TRUE)
nf <- layout(matrix(c(1, 2), ncol = 2), height = c(1, 1), width = c(1, 1))

par(mar = c(5.1, 5, 4, 0))
plot(pt.woman[[1]]\$time, pt.woman[[1]]\$pstate2, type = "s", lwd = 2,
ylim = c(0, 0.1), xlim = c(0, 100), xlab = "Days",
ylab = expression(hat(P)["01"](0, t, z)), main = "Woman")
lines(pt.woman[[1]]\$time, cis.woman\$upper, type = "s", lwd = 2, lty = 3)
lines(pt.woman[[1]]\$time, cis.woman\$lower, type = "s", lwd = 2, lty = 3)

par(mar = c(5.1, 3, 4, 2))
plot(pt.man[[1]]\$time, pt.man[[1]]\$pstate2, type = "s", lwd = 2,
ylim = c(0, 0.1), xlim = c(0, 100), xlab = "Days",
ylab = "", main = "Man", yaxt = "n")
lines(pt.man[[1]]\$time, cis.man\$upper, type = "s", lwd = 2, lty = 3)
lines(pt.man[[1]]\$time, cis.man\$lower, type = "s", lwd = 2, lty = 3)

par(oldpar)
```

### Section 10.2.2

Cox models for 0 → 1 and 1 → 0 transitions with robust variance:

```cox.ventil.01 <- coxph(Surv(time, to == 1) ~ age + sex +
cluster(id), sir.cont,
subset = from == 0)
cox.ventil.10 <- coxph(Surv(time, to == 0) ~ age + sex +
cluster(id), sir.cont,
subset = from == 1)
summary(cox.ventil.01)
summary(cox.ventil.10)
```