library(survival)
library(tidyverse)Estimating and Testing Cause Specific Hazard Ratio Using R
Objective
In this document we present how to estimate and test cause specific hazard ratio for the probability of experiencing a certain event at a given time in a competing risks model in R. We focus on the basic model where each subject experiences only one out of k possible events as depicted in the figure below.
As this document aims to provide syntax for estimating and testing cause-specific hazard ratios using Cox’s PH model for competing risks, we assume that readers have working knowledge of a competing risks framework. The Reference below list a few literature for a quick refresher on this topic.
The syntax given here produce results match that produced by the default settings of SAS PROC PHREG (see the companion SAS document). This is usually necessary if validating results from the two software is the objective.
R Package
We use the survival package in this document.
Data used
The bone marrow transplant (BTM) dataset as presented by Guo & So (2018) is used. The dataset has the following variables:
Grouphas three levels, indicating three disease groups.Tis the disease-free survival time in days. A derived variableTYears = T/365.25is used in the analysis.Statushas value 0 ifTis censored; 1 ifTis time to relapse; 2 ifTis time to death.WaitTimeis the waiting time to transplant in days.For illustration, a categorical variable
waitCatis created fromwaitTimeaswaitCat = TRUEifwaitTime > 200, andFALSEotherwise.
bmt <- haven::read_sas(file.path("../data/bmt.sas7bdat")) %>%
mutate(
Group = factor(
Group,
levels = c(1, 2, 3),
labels = c('ALL', 'AML-Low Risk', 'AML-High Risk')
),
Status = factor(
Status,
levels = c(0, 1, 2),
labels = c('Censored', 'Relapse', 'Death')
),
TYears = T / 365.25,
waitCat = (WaitTime > 200),
ID = row_number()
)Estimating and testing the cause specific hazard ratio
Syntax-wise there are two ways to generate the estimates and related outputs using survival::coxph(). They produce essentially the same results except that the global null hypotheses are different.
Syntax 1: All competing events in one go
csh.1 <- survival::coxph(
survival::Surv(TYears, Status) ~ Group + survival::strata(waitCat),
data = bmt,
id = ID,
ties = 'breslow', ## default is 'efron'
robust = FALSE ## default is TRUE
)
summary(csh.1)Call:
survival::coxph(formula = survival::Surv(TYears, Status) ~ Group +
survival::strata(waitCat), data = bmt, ties = "breslow",
robust = FALSE, id = ID)
n= 137, number of events= 83
coef exp(coef) se(coef) z
GroupAML-Low Risk_1:2 -0.9271 0.3957 0.4493 -2.063
GroupAML-High Risk_1:2 0.6227 1.8640 0.3635 1.713
survival::strata(waitCat)waitCat=TRUE_1:2 -0.1310 0.8772 0.3240 -0.404
GroupAML-Low Risk_1:3 -0.3982 0.6715 0.3936 -1.012
GroupAML-High Risk_1:3 0.1177 1.1250 0.4034 0.292
survival::strata(waitCat)waitCat=TRUE_1:3 -0.2487 0.7798 0.3456 -0.720
Pr(>|z|)
GroupAML-Low Risk_1:2 0.0391 *
GroupAML-High Risk_1:2 0.0867 .
survival::strata(waitCat)waitCat=TRUE_1:2 0.6860
GroupAML-Low Risk_1:3 0.3117
GroupAML-High Risk_1:3 0.7704
survival::strata(waitCat)waitCat=TRUE_1:3 0.4717
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95
GroupAML-Low Risk_1:2 0.3957 2.5272 0.1640
GroupAML-High Risk_1:2 1.8640 0.5365 0.9142
survival::strata(waitCat)waitCat=TRUE_1:2 0.8772 1.1400 0.4649
GroupAML-Low Risk_1:3 0.6715 1.4891 0.3105
GroupAML-High Risk_1:3 1.1250 0.8889 0.5103
survival::strata(waitCat)waitCat=TRUE_1:3 0.7798 1.2824 0.3961
upper .95
GroupAML-Low Risk_1:2 0.9546
GroupAML-High Risk_1:2 3.8005
survival::strata(waitCat)waitCat=TRUE_1:2 1.6553
GroupAML-Low Risk_1:3 1.4525
GroupAML-High Risk_1:3 2.4801
survival::strata(waitCat)waitCat=TRUE_1:3 1.5351
Concordance= 0.631 (se = 0.031 )
Likelihood ratio test= 18.08 on 6 df, p=0.006
Wald test = 16.52 on 6 df, p=0.01
Score (logrank) test = 18.71 on 6 df, p=0.005In the output, rows with suffix 1:2 are for Status = 2, or Relapse; and 1:3 are for Status = 3, or Death. As usual, censoring must be the lowest level in Status, which in this example is coded 0.
Since both events (Relapse and Death) are model together, the global tests have 4 degrees of freedom, with the null hypothesis that there is no difference among different levels of Group in either Relapse or Death.
Syntax 2: One event at a time
csh.2 <- survival::coxph(
survival::Surv(TYears, Status == 'Relapse') ~
Group + survival::strata(waitCat),
data = bmt,
id = ID,
ties = 'breslow', ## default is 'efron'
robust = FALSE ## default is TRUE
)
summary(csh.2)Call:
survival::coxph(formula = survival::Surv(TYears, Status == "Relapse") ~
Group + survival::strata(waitCat), data = bmt, ties = "breslow",
robust = FALSE, id = ID)
n= 137, number of events= 42
coef exp(coef) se(coef) z
GroupAML-Low Risk -0.9271 0.3957 0.4493 -2.063
GroupAML-High Risk 0.6227 1.8640 0.3635 1.713
survival::strata(waitCat)waitCat=TRUE -0.1310 0.8772 0.3240 -0.404
Pr(>|z|)
GroupAML-Low Risk 0.0391 *
GroupAML-High Risk 0.0867 .
survival::strata(waitCat)waitCat=TRUE 0.6860
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
GroupAML-Low Risk 0.3957 2.5272 0.1640 0.9546
GroupAML-High Risk 1.8640 0.5365 0.9142 3.8005
survival::strata(waitCat)waitCat=TRUE 0.8772 1.1400 0.4649 1.6553
Concordance= 0.69 (se = 0.037 )
Likelihood ratio test= 16.04 on 3 df, p=0.001
Wald test = 14.49 on 3 df, p=0.002
Score (logrank) test = 16.66 on 3 df, p=8e-04The results are identical to those labeled with 1:2 in the earlier outputs under Syntax 1 above. However, since only Relapse is modeled, the global tests have only 2 degrees of freedom with the null hypothesis that there is no difference among different levels of Group for Relapse.
Summary
In
survival::coxph()the default method for handling ties isties = 'efron'. To match results with SAS, this needs to be changed toties =breslow`.For multi-state models such as a competing risk analysis,
survival::coxph()by default estimate the standard errors of parameter estimates with a robust sandwich estimator. To match default results with SAS, this needs to be set torobust = FALSE.Due to differences in internal numerical estimation methods of R and SAS, results only match up to the 4th decimal places. However, overall consistency can be established between the two for estimating and testing cause-specific hazard ratio using Cox’s PH model.
Session Info
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.4.2 (2024-10-31)
os macOS Sequoia 15.6.1
system aarch64, darwin20
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Europe/London
date 2025-08-29
pandoc 3.4 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
─ Packages ───────────────────────────────────────────────────────────────────
! package * version date (UTC) lib source
cmprsk 2.2-12 2024-05-19 [1] CRAN (R 4.4.0)
P survival * 3.7-0 2024-06-05 [?] CRAN (R 4.4.0)
tidycmprsk 1.1.0 2024-08-17 [1] CRAN (R 4.4.0)
[1] /Users/christinafillmore/Documents/GitHub/CAMIS/renv/library/macos/R-4.4/aarch64-apple-darwin20
[2] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
P ── Loaded and on-disk path mismatch.
──────────────────────────────────────────────────────────────────────────────Reference
Guo C and So Y. (2018). “Cause-specific analysis of competing risks using the PHREG procedure.” In Proceedings of the SAS Global Forum 2018 Conference. Cary, NC: SAS Institute Inc. https://support.sas.com/resources/papers/proceedings18/2159-2018.pdf.
Pintilie M. (2006). Competing Risks: A Practical Perspective. Wiley. http://dx.doi.org/10.1002/9780470870709
Therneau T, Crowson C, and Atkinson E. (2024). “Multi-state models and competing risks.” https://cran.r-project.org/web/packages/survival/vignettes/compete.pdf