In this article we illustrate how to fit cause specific hazards models to competing risks data. The standard way to estimate cause specific hazards is to create one data set for each event type and fit a separate model. However, it is also possible to create one combined data set and enter the event type as a covariate (with interactions), such that it is possible to estimate shared effects (i.e., effects that contribute equally to the hazard of multiple event types).
For illustration we use the fourD
data set from the etm
package. The data set contains time-constant covariates like age
and sex
as well as time-to-event (time
) and event type indicator status
(0 = censored, 1 = death from cardiovascular events, 2 = death from other causes).
## id sex age medication status time treated
## 1 5002 Male 60 Placebo 0 5.8480493 0
## 4 5006 Female 68 Placebo 0 5.2539357 0
## 7 5011 Female 70 Placebo 1 2.9541410 0
## 9 5014 Male 69 Placebo 1 0.9856263 0
## 10 5017 Female 58 Placebo 1 0.2902122 0
## 11 5018 Male 63 Placebo 1 3.9452430 0
The data transformation required to fit PAMMs to competing risks data is similar to the transformation in the single event case (see the data transformation vignette for details). In fact, internally the standard transformation is applied to each event type using as_ped
, however, some choices have to be made
For cause specific hazards without shared effects the combination of cause specific interval split points and list output is usually sufficient. For models with shared effects we need to stack the individual data sets and use split points common for all event types.
Below we transform the data set for the case without shared effects. By specifying cobmine = FALSE
, the individual data sets are not stacked but rather returned in a list.
cut <- sample(fourD$time, 100)
ped <- fourD %>%
select(-medication, - treated) %>%
as_ped(Surv(time, status)~., id = "id", cut = cut, combine = FALSE)
str(ped, 1)
## List of 2
## $ cause = 1:Classes 'ped' and 'data.frame': 32917 obs. of 9 variables:
## ..- attr(*, "breaks")= num [1:100] 0.203 0.208 0.26 0.271 0.323 ...
## ..- attr(*, "id_var")= chr "id"
## ..- attr(*, "intvars")= chr [1:6] "id" "tstart" "tend" "interval" ...
## ..- attr(*, "trafo_args")=List of 3
## ..- attr(*, "time_var")= chr "time"
## ..- attr(*, "status_var")= chr "status"
## $ cause = 2:Classes 'ped' and 'data.frame': 32917 obs. of 9 variables:
## ..- attr(*, "breaks")= num [1:100] 0.203 0.208 0.26 0.271 0.323 ...
## ..- attr(*, "id_var")= chr "id"
## ..- attr(*, "intvars")= chr [1:6] "id" "tstart" "tend" "interval" ...
## ..- attr(*, "trafo_args")=List of 3
## ..- attr(*, "time_var")= chr "time"
## ..- attr(*, "status_var")= chr "status"
## - attr(*, "class")= chr [1:4] "ped_cr_list" "ped_cr" "ped" "list"
## - attr(*, "trafo_args")=List of 5
## - attr(*, "risks")= int [1:2] 1 2
# data set for event type 1 (death from cardiovascular events)
head(ped[[1]])
## id tstart tend interval offset ped_status
## 1 5002 0.0000000 0.2026010 (0,0.2026009582] -1.596517 0
## 2 5002 0.2026010 0.2080767 (0.2026009582,0.2080766598] -5.207435 0
## 3 5002 0.2080767 0.2600958 (0.2080766598,0.2600958248] -2.956143 0
## 4 5002 0.2600958 0.2710472 (0.2600958248,0.2710472279] -4.514288 0
## 5 5002 0.2710472 0.3230664 (0.2710472279,0.3230663929] -2.956143 0
## 6 5002 0.3230664 0.5475702 (0.3230663929,0.5475701574] -1.493863 0
## sex age cause
## 1 Male 60 1
## 2 Male 60 1
## 3 Male 60 1
## 4 Male 60 1
## 5 Male 60 1
## 6 Male 60 1
# data set for event type 2 (death from other causes)
head(ped[[2]])
## id tstart tend interval offset ped_status
## 1 5002 0.0000000 0.2026010 (0,0.2026009582] -1.596517 0
## 2 5002 0.2026010 0.2080767 (0.2026009582,0.2080766598] -5.207435 0
## 3 5002 0.2080767 0.2600958 (0.2080766598,0.2600958248] -2.956143 0
## 4 5002 0.2600958 0.2710472 (0.2600958248,0.2710472279] -4.514288 0
## 5 5002 0.2710472 0.3230664 (0.2710472279,0.3230663929] -2.956143 0
## 6 5002 0.3230664 0.5475702 (0.3230663929,0.5475701574] -1.493863 0
## sex age cause
## 1 Male 60 2
## 2 Male 60 2
## 3 Male 60 2
## 4 Male 60 2
## 5 Male 60 2
## 6 Male 60 2
To fit the model, we could loop through the list entries and fit the model of interest, however, there is also a convenience function, that recognizes the data type and fits the models accordingly:
library(mgcv)
pam_csh <- map(ped, ~ pamm(ped_status ~ s(tend) + sex + age, data = .x))
map(pam_csh, summary)
## $`cause = 1`
##
## Family: poisson
## Link function: log
##
## Formula:
## ped_status ~ s(tend) + sex + age
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.192363 0.566010 -5.640 1.7e-08 ***
## sexMale -0.144349 0.127828 -1.129 0.2588
## age 0.020247 0.008188 2.473 0.0134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(tend) 1.004 1.008 7.97 0.00483 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -0.00834 Deviance explained = 0.559%
## UBRE = -0.91428 Scale est. = 1 n = 32917
##
## $`cause = 2`
##
## Family: poisson
## Link function: log
##
## Formula:
## ped_status ~ s(tend) + sex + age
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.91075 0.84671 -8.162 3.30e-16 ***
## sexMale 0.21104 0.17409 1.212 0.225
## age 0.06311 0.01193 5.290 1.22e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(tend) 1.026 1.052 11.46 0.000791 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -0.00535 Deviance explained = 2.25%
## UBRE = -0.94944 Scale est. = 1 n = 32917
The data transformation is performed as before, but setting combine=TRUE
(the default), the interval cut points are created based on all event times (event times of all event types, here) and stacked:
ped_stacked <- fourD %>%
select(-medication, - treated) %>%
as_ped(Surv(time, status)~., id = "id", cut = cut) %>%
mutate(cause = as.factor(cause))
head(ped_stacked)
## id tstart tend interval offset ped_status
## 1 5002 0.0000000 0.2026010 (0,0.2026009582] -1.596517 0
## 2 5002 0.2026010 0.2080767 (0.2026009582,0.2080766598] -5.207435 0
## 3 5002 0.2080767 0.2600958 (0.2080766598,0.2600958248] -2.956143 0
## 4 5002 0.2600958 0.2710472 (0.2600958248,0.2710472279] -4.514288 0
## 5 5002 0.2710472 0.3230664 (0.2710472279,0.3230663929] -2.956143 0
## 6 5002 0.3230664 0.5475702 (0.3230663929,0.5475701574] -1.493863 0
## sex age cause
## 1 Male 60 1
## 2 Male 60 1
## 3 Male 60 1
## 4 Male 60 1
## 5 Male 60 1
## 6 Male 60 1
Model for cause specific hazards with shared effects is performed by inclusion of interaction effects:
pam_csh_shared <- pamm(
formula = ped_status ~ s(tend, by = cause) + sex + sex:cause + age + age:cause,
data = ped_stacked)
summary(pam_csh_shared)
##
## Family: poisson
## Link function: log
##
## Formula:
## ped_status ~ s(tend, by = cause) + sex + sex:cause + age + age:cause
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.178455 0.585401 -5.430 5.65e-08 ***
## sexMale -0.136946 0.132264 -1.035 0.30048
## age 0.020086 0.008468 2.372 0.01770 *
## sexFemale:cause2 -3.796266 1.059605 -3.583 0.00034 ***
## sexMale:cause2 -3.490538 1.013118 -3.445 0.00057 ***
## age:cause2 0.044233 0.015049 2.939 0.00329 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(tend):cause1 1.001 1.003 7.402 0.006550 **
## s(tend):cause2 1.025 1.050 11.377 0.000824 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -0.00686 Deviance explained = 2.01%
## UBRE = -0.93355 Scale est. = 1 n = 63338
The results indicate that cause specific terms (interactions) are necessary in this case and the two models largely agree. For example, the age effect for the two causes are very similar for both models:
pamm_csh1
) vs. 0.02 (pamm_csh_shared
)pamm_csh2
) vs. 0.02 + 0.044= 0.064 (pamm_csh_shared
)Finally, in many cases we will want to calculate and visualize the cumulative incidence functions for different covariate combinations. In pammtools
this can be again achieved using make_newdata
and using the appropriate add_*
function, here add_cif
:
ndf <- ped_stacked %>%
make_newdata(tend = unique(tend), cause = unique(cause)) %>%
group_by(cause) %>%
add_cif(pam_csh_shared)
ndf %>%
select(tend, cause, cif, cif_lower, cif_upper) %>%
group_by(cause) %>%
slice(1:3)
## # A tibble: 6 × 5
## # Groups: cause [2]
## tend cause cif cif_lower cif_upper
## <dbl> <fct> <dbl> <dbl> <dbl>
## 1 0.203 1 0.0213 0.0168 0.0266
## 2 0.208 1 0.0219 0.0173 0.0274
## 3 0.260 1 0.0274 0.0216 0.0342
## 4 0.203 2 0.0100 0.00710 0.0139
## 5 0.208 2 0.0103 0.00730 0.0142
## 6 0.260 2 0.0129 0.00914 0.0178
ggplot(ndf, aes(x = tend, y = cif)) +
geom_line(aes(col = cause)) +
geom_ribbon(
aes(ymin = cif_lower, ymax = cif_upper, fill = cause),
alpha = .3)
Similar to other applications of add_*
functions, we can additionally group by other covariate values:
ndf <- ped_stacked %>%
make_newdata(tend = unique(tend), cause = unique(cause), sex = unique(sex))
ndf <- ndf %>%
group_by(cause, sex) %>%
add_cif(pam_csh_shared)
The estimated CIFs can then be compared w.r.t. to cause
for each category of sex
:
ggplot(ndf, aes(x = tend, y = cif)) +
geom_line(aes(col = cause)) +
geom_ribbon(
aes(ymin = cif_lower, ymax = cif_upper, fill = cause),
alpha = .3) +
facet_wrap(~sex, labeller = label_both)
or w.r.t. to sex
for each cause:
ggplot(ndf, aes(x = tend, y = cif)) +
geom_line(aes(col = sex)) +
geom_ribbon(
aes(ymin = cif_lower, ymax = cif_upper, fill = sex),
alpha = .3) +
facet_wrap(~cause, labeller = label_both)