knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) load(file = "../m_cmm_vignette.RData") # and other packages that we need for visualization library(survival)
# load the dena package with the data and functions library(dena)
rm(list = ls()) data(simdat3) # see the end of the script for simulation code
first we need to restructure the data a bit (i.e., add a variable that indicates the previous state (in the tutorial example we assumed that the previous state is always the alone state))
# for(id in unique(simdat3$id)){ tmp <- simdat3[simdat3$id %in% id,] tmp$from <- c(NA,as.character(tmp$type[-nrow(tmp)])) simdat3[simdat3$id %in% id,"from"] <- tmp$from # and some restructuring of the time variable for(row in 1:nrow(tmp)){ if(row == 1){ tmp$start[1] <- 0 tmp$end[1] <- tmp$time[1] }else{ tmp$start[row] <- tmp$end[row-1] tmp$end[row] <- tmp$start[row]+tmp$time[row] } } simdat3[simdat3$id %in% id,"start"] <- tmp$start simdat3[simdat3$id %in% id,"end"] <- tmp$end }
Now we have all the necessary variables in the correct format for mutlsitate analysis:
head(simdat3[,c("id","start","end","from","type","Covariate1","Covariate2")])
descriptives <- survcheck(Surv(start, end, type) ~ 1, data = simdat3, id = id, istate = from) descriptives
simdat3$from <- factor(simdat3$from, levels = c("censored", "alone","family","friend","partner")) simdat3$type <- factor(simdat3$type, levels = c("censored", "alone","family","friend","partner")) # multistate method from the survival package (package manual, Therneau, 2020, p 61) m1 <-coxph(Surv(start, end, type) ~ Covariate1 + Covariate2, data = simdat3, id = id, istate = from) m1
Lets look at the model output:
m1
The numbers 1 to 5 correspond to the states (1) Alone, (2) Family, (3) Friend, and (4) Partner. The model output shows the hazard ratios for the covariates separate for each transition between the three states. For example, in the transition from (4) Partner to (2) Family, only the Covariate 1 shows a significant effect (beta = -0.195, SE = 0.078, p = 0.005).
With the following function, one can also transform the multistate output model into a more reader-friendly table
to.paper.table <- function(m1, digits = 3){ tmp <- as.data.frame(summary(m1)$coefficients) tmp$names <- rownames(summary(m1)$cmap) transitions <- do.call(rbind.data.frame, strsplit(rownames(tmp),"_"))[2] tmp$from <- summary(m1)$states[as.numeric(substr(as.vector(transitions[,1]),1,1))] tmp$to <- summary(m1)$states[as.numeric(substr(as.vector(transitions[,1]),3,3))] tmp$sig <- ifelse(tmp$`Pr(>|z|)`> 0.05, "",ifelse(tmp$`Pr(>|z|)` > 0.01,"*",ifelse(tmp$`Pr(>|z|)` > 0.001,"**","***"))) tmp$coef.LB <- tmp$coef-1.96*tmp$`robust se` tmp$coef.UB <- tmp$coef+1.96*tmp$`robust se` round_df(tmp[,c("from", "to","names","coef","sig","robust se","Pr(>|z|)","coef.LB","coef.UB")],digits = digits) } to.paper.table(m1)
set.seed(5689) simdat3 <- frailtySurv::genfrail(N = 50, K = 40, beta = c(log(0.5), -log(1)), frailty = "gamma", theta = 3, censor.rate = 0, lambda_0=function(t, tau=4.6, C=0.01) (tau*(C*t)^tau)/t) colnames(simdat3) <- c("id","rep","time","event","Covariate1","Covariate2") simdat3$type <- factor(sample(c("partner", "family","friend", "alone"), nrow(simdat3), replace = T, prob = c(0.2,0.2,0.2,0.3))) #save(simdat3, file = "../data/simdat3.RData")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.