What does the output of flexsurvreg, for a Gompertz proportional hazards model, mean?

Question

I am estimating a Gompertz proportional hazards model in R using the package "flexsurvreg", but I'm having a hard time understanding the output of this function. My dataset is collected from skeletal remains, so my event is estimated age at death and my explanatory variable in community dwelling- either in monastic or urban communities during the medieval period in the UK. In the past, I've used the Cox proportional hazards model from the "survival" package, but the standard in my field is a Gompertz model. The problem is that the output of flexsurvreg is different than the coxph function so I don't know how to interpret the results.

Here's the output from flexsurvreg:

Call:
flexsurvreg(formula = Surv(age) ~ monastery, data = dat, dist = "gompertz")

Estimates: 
                data mean  est       L95%      U95%      se        exp(est)  L95%    
shape                 NA   0.072851  0.069204  0.076498  0.001861        NA        NA
rate                  NA   0.002167  0.001817  0.002584  0.000194        NA        NA
monasteryUrban  0.322643   0.338239  0.205303  0.471174  0.067826  1.402475  1.227897
                U95%    
shape                 NA
rate                  NA
monasteryUrban  1.601874

N = 1029,  Events: 1029,  Censored: 0
Total time at risk: 42770.5
Log-likelihood = -4113.269, df = 3
AIC = 8232.537

And here's the output from coxph on the same data

Call:
coxph(formula = Surv(age) ~ monastery, data = dat)

  n= 1029, number of events= 1029 

                  coef exp(coef) se(coef)    z Pr(>|z|)   
monasteryUrban 0.21742   1.24286  0.06753 3.22  0.00128 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

               exp(coef) exp(-coef) lower .95 upper .95
monasteryUrban     1.243     0.8046     1.089     1.419

Concordance= 0.531  (se = 0.011 )
Likelihood ratio test= 10.11  on 1 df,   p=0.001
Wald test            = 10.37  on 1 df,   p=0.001
Score (logrank) test = 10.41  on 1 df,   p=0.001

For the flexsurvreg output, how do I interpret the overall performance of the model- is log-likelihood the same thing as a likelihood ratio test? How can I tell whether the hazards between my two groups (monastic vs. urban dwellers) are significantly different?

Answer 1

Unlike the Cox model, the model based on the Gompertz distribution estimates parameter values describing the shape of the baseline hazard as well as the association of your covariate with survival. So the first two rows of output represent the estimated shape and rate of the underlying Gompertz distribution, which (unlike some other distributions) is compatible with proportional hazards; see the flexsurvreg() and dgompertz() documentation in the flexsurv package.

The third row represents results for the regression coefficient and corresponding hazard ratio (HR), exp(est), associated with your covariate, monastery. That's for the difference (resp. HR) between the Urban category and the reference category (presumably monastic-dwelling). That's analogous to the coefficient and HR for the Cox model. The 95% confidence intervals for that coefficient show that it's significantly different from 0 (and that the associated HR is different from 1), the usual criterion for significance as in the Cox model.

With a parametric model you should show that the data fit the modeled shape of the curve adequately. As you have expectations particular to your field of study you should follow its standard procedures for documenting the adequacy of the shape of the fit.