Zu den Inhalten springen

Institute of Medical Biometry and Statistics (IMBI)

R code for chapter 10

Quicklinks for Chapter 10

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)