Sie sind hier: Startseite Competing Risks and … R Code R code for chapter 10

R code for chapter 10

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