In joint frailty model with competing risks and recurrent events - how to interpret the coefficients?

Question

I took the output from the example described in manual to one of R's packags for frailty modelling. Apart from the statistical package, how to interpret the outcome for recurrent events? Is this HR for "having more recurrences"? "for having the event and recurrences" (somehow adjusted)? How to interpret the HRs?

> modJoint.gap
Call:
frailtyPenal(formula = Surv(time, event) ~ cluster(id) + sex + 
    dukes + charlson + terminal(death), formula.terminalEvent = ~sex + 
    dukes + charlson, data = readmission, recurrentAG = FALSE, 
    n.knots = 10, kappa = c(100, 100), hazard = "Splines")


  Joint gamma frailty model for recurrent and a terminal event processes 
  using a Penalized Likelihood on the hazard function 

Recurrences:
------------- 
                 coef exp(coef) SE coef (H) SE coef (HIH)        z          p
sexFemale   -0.527310   0.59019    0.140818      0.137168 -3.74462 1.8067e-04
dukesC       0.397348   1.48787    0.154980      0.172323  2.56387 1.0351e-02
dukesD       1.274035   3.57525    0.202234      0.180836  6.29981 2.9801e-10
charlson1-2  0.390411   1.47759    0.256722      0.324360  1.52075 1.2832e-01
charlson3    0.433207   1.54220    0.136762      0.143414  3.16760 1.5370e-03

           chisq df global p
dukes    46.2610  2 9.01e-11
charlson 12.3464  2 2.08e-03

Terminal event:
---------------- 
                 coef exp(coef) SE coef (H) SE coef (HIH)        z          p
sexFemale   -0.340342  0.711527    0.220174      0.248636 -1.54579 1.2216e-01
dukesC       0.903630  2.468547    0.337757      0.337335  2.67539 7.4643e-03
dukesD       2.724323 15.246081    0.382510      0.373507  7.12223 1.0619e-12
charlson1-2  0.714284  2.042723    0.624573      0.565360  1.14364 2.5277e-01
charlson3    1.112842  3.042996    0.246157      0.263593  4.52086 6.1589e-06

           chisq df global p
dukes    62.0806  2 3.31e-14
charlson 20.7119  2 3.18e-05

 Frailty parameters: 
   theta (variance of Frailties, w): 0.737739 (SE (H): 0.104455 ) p = 8.1612e-13 
   alpha (w^alpha for terminal event): 0.733697 (SE (H): 0.218835 ) p = 0.00080016 
 
Penalized marginal log-likelihood = -0.48
   Convergence criteria: 
   parameters = 4.16e-06 likelihood = 5.66e-05 gradient = 4.24e-08 

Likelihood Cross-Validation (LCV) criterion in the semi parametrical case:
   approximate LCV = 0.0351718 

   n observations= 861  n subjects= 403
   n recurrent events= 458
   n terminal events= 109
   n censored events= 403
   number of iterations:  8 
   Number of nodes for the Gauss-Laguerre quadrature:  20 

   Exact number of knots used:  10 
   Value of the smoothing parameters:  100 100

Answer 1

In this type of model there are two types of event, one Recurrent (indicated by event in these data) and one Terminal (called death). There are two separate sets of regression coefficients estimated, one set for the Recurrences and another set for the Terminal event.*

Remember that what you show here with exp(coef) are hazard ratios. They represent the relative risks of having an event. In this model, the hazard of an individual having a Recurrent event is assumed independent of the number of prior Recurrent events. So one good way to think about the HR values from the output labeled Recurrences is that they are the relative hazards for having a (new) Recurrent event at any time, versus the hazard at baseline covariate values (here: Male, dukes A-B, charlson 0).

The HRs for the Terminal event are separate from those for the Recurrences. They have the usual interpretation for survival models of death.

You can inspect the baseline hazards of the two events over time with plot(modJoint.gap).

The individual gamma frailties $\omega_i$ (w), with modeled parameter theta, act multiplicatively on the covariate-adjusted hazards for the Recurrences. Those same individual frailties are raised to the power of another estimated parameter,alpha, for the Terminal event hazards.

---

*Function calls with this frailty pack package are a bit different from the usual survival package, as the terminal event is included here as a term +terminal(death) in the formula rather than having a multi-level factor for the event indicator in the Surv() function. That might lead to some confusion.