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)

Load the example data

rm(list = ls())

data(simdat3) # see the end of the script for simulation code 

Restructuring data

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")])

Some descriptives

descriptives <- survcheck(Surv(start, end, type) ~ 1, 
                   data = simdat3, id = id, istate = from)

descriptives

Multistate estimation

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)

Dataset simulation

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


timonelmer/dena documentation built on April 15, 2023, 11:51 p.m.