Cause-Specific Hazard vs Multi-State Model in R

Question

I'm following the examples set out in https://cran.r-project.org/web/packages/survival/vignettes/compete.pdf for this, although this problem is on a real dataset.

I have a competing risks scenario where patients can die of one of two causes. An example dataset can be set up using:

library(survival)
data(cancer)

mgus2$etime <- with(mgus2, ifelse(pstat == 0, futime, ptime))
event <- with(mgus2, ifelse(pstat == 0, 2 * death, 1))
mgus2$event <- factor(event, 0:2, labels = c("censor", "pcm", "death"))

In this example, patients can have either pcm or death, but not both - therefore a traditional competing risks situation. A multi-state model can be fitted using:

coxph(Surv(etime, event) ~ sex, data = mgus2, id = id)

which gives:

Call:
coxph(formula = Surv(etime, event) ~ sex, data = mgus2, id = id)

      
1:2        coef exp(coef) se(coef) robust se      z     p
  sexM -0.05938   0.94235  0.18723   0.18733 -0.317 0.751

      
1:3       coef exp(coef) se(coef) robust se     z        p
  sexM 0.22800   1.25608  0.06900   0.06869 3.319 0.000902

 States:  1= (s0), 2= pcm, 3= death 

Likelihood ratio test=11.11  on 2 df, p=0.003866
n= 1384, number of events= 975 

My understanding is that this fits two separate cause-specific cox models to the data. If this were true, I would be able to equivalently model the data using:

coxph(Surv(etime, event == "pcm") ~ sex, data = mgus2, id = id)
coxph(Surv(etime, event == "death") ~ sex, data = mgus2, id = id)

Which gives the very similar, but not identical, result:

> coxph(Surv(etime, event == "pcm") ~ sex, data = mgus2, id = id)
Call:
coxph(formula = Surv(etime, event == "pcm") ~ sex, data = mgus2, 
    id = id)

         coef exp(coef) se(coef)      z     p
sexM -0.05938   0.94235  0.18723 -0.317 0.751

Likelihood ratio test=0.1  on 1 df, p=0.7511
n= 1384, number of events= 115 
> coxph(Surv(etime, event == "death") ~ sex, data = mgus2, id = id)
Call:
coxph(formula = Surv(etime, event == "death") ~ sex, data = mgus2, 
    id = id)

      coef exp(coef) se(coef)     z        p
sexM 0.228     1.256    0.069 3.304 0.000952

Likelihood ratio test=11.01  on 1 df, p=0.0009061
n= 1384, number of events= 860

My question is, why do these two results differ? Am I not actually fitting a cause-specific model in the first instance, or have I mis-understood what is required to fit these models? The reason I want to model them separately is that I will be doing subgroup analyses, followed by post-estimation using emmeans. So far, I can't make this work with the multi-state method, only single failure cox models.

Answer 1

The difference is that in your competing risks model, robust variances are being calculated. From the coxph docs, the robust argument is described as

robust
should a robust variance be computed? The default is TRUE if: there is a cluster argument, there are case weights that are not 0 or 1, or there are id values with more than one event.

If robust is set to TRUE, the calculated robust variances will also be used to calculate the z scores and p values of the coefficient table. There are no robust variances calculated in the individual endpoint models because there is only a single event per ID. This results in using the standard variances to calculate the z score and p value, hence the differences in the coefficient table output.

You can set robust = FALSE in your multi-endpoint model to make the coefficient tables exactly match the individual endpoint models, or robust = TRUE in your individual endpoint models to get them to match the multi-endpoint version.

coxph(Surv(etime, event) ~ sex, data = mgus2, id = id, robust = FALSE)
#> Call:
#> coxph(formula = Surv(etime, event) ~ sex, data = mgus2, robust = FALSE, 
#>     id = id)
#> 
#>       
#> 1:2        coef exp(coef) se(coef)      z     p
#>   sexM -0.05938   0.94235  0.18723 -0.317 0.751
#> 
#>       
#> 1:3     coef exp(coef) se(coef)     z        p
#>   sexM 0.228     1.256    0.069 3.304 0.000952
#> 
#>  States:  1= (s0), 2= pcm, 3= death 
#> 
#> Likelihood ratio test=11.11  on 2 df, p=0.003866
#> n= 1384, number of events= 975

Compared to

coxph(Surv(etime, event == "pcm") ~ sex, data = mgus2, id = id)
#> Call:
#> coxph(formula = Surv(etime, event == "pcm") ~ sex, data = mgus2, 
#>     id = id)
#> 
#>          coef exp(coef) se(coef)      z     p
#> sexM -0.05938   0.94235  0.18723 -0.317 0.751
#> 
#> Likelihood ratio test=0.1  on 1 df, p=0.7511
#> n= 1384, number of events= 115

coxph(Surv(etime, event == "death") ~ sex, data = mgus2, id = id)
#> Call:
#> coxph(formula = Surv(etime, event == "death") ~ sex, data = mgus2, 
#>     id = id)
#> 
#>       coef exp(coef) se(coef)     z        p
#> sexM 0.228     1.256    0.069 3.304 0.000952
#> 
#> Likelihood ratio test=11.01  on 1 df, p=0.0009061
#> n= 1384, number of events= 860