Cumulative baseline hazard for relative survival models

[heor-robyoung]

With the publication of flexsurvcure, this seems more important than ever. For specification of relative survival models, only the baseline hazard is passed through to the likelihood function, and is summed for those patients who have events to the hazard from the parametric distribution.

Does this not ignore the contribution of the cumulative baseline hazard to the likelihood of the censored observations?

[chjackson]

Sorry for leaving this issue open for so long!

The likelihood function that flexsurv uses for relative survival models is from the paper by Nelson et al..

The likelihood for someone with an event is

log(h^*(t_i) + \lambda(t_i)) + log(R(t_i))

where h* is the baseline hazard, lambda is the excess hazard, and R the relative survival.

Or for a right-censored observation, the likelihood is log(R(t_i))

In parametric relative survival models, the parametric model is placed on the relative survival R or excess hazard lambda , rather than the absolute hazard or survival. Note that -log(R(t)) equals cumulative total hazard - cumulative baseline hazard = cumulative excess hazard.

So the likelihood contribution for a censored person doesn't depend on the baseline hazard.

I've now put a reference to this paper in the help page for flexsurvreg.

[heor-robyoung]

Hi Chris,

The term relating to the accumulation of hazard from the baseline function is dropped because it is a constant term, and so isn't influential on the location of the maximum likelihood. However, the constant contribution should still be there when considering the absolute value of the likelihood function as when computing AIC; i.e. eqn (5) within the reference is fine for the search, but is misnamed as it does not describe the full lnLi; equation (4) is the full evaluation.

Best,

Rob.

[chjackson]

Right I see. I'm wondering when this matters for model comparison though. Would we ever want to compare models with two different baseline hazard offsets, and could this be done validly with likelihood (4)?

In any case I should document where the "likelihood" value comes from, at least for comparison with other software. I don't suppose you know what the Stata packages for this do here?

I haven't ever used relative survival models, but I see them getting used more and more, so I'd like this to be right!

[heor-robyoung]

Hi Chris,

I think we'd commonly see this in comparison of models with 0 baseline hazard versus relative survival models; generally I've just mentioned that the two are incomparable due to this issue, but provided that the full likelihood is evaluated I don't see that comparison by information criteria would be any more inappropriate than for simple distributions. I'm afraid I don't know about Stata - I'll ask around some of my Stata-based colleagues and see if they know.

Best,

Rob.

[chjackson]

Yes, agreed, I see it as like having a covariate with a known effect.

Anyway I've now implemented the corrected likelihood in https://github.com/chjackson/flexsurv-dev/commit/6b9445561714473a6903a35b6cfbb3f4542321a9 . Let me know if this looks sound.

[heor-robyoung]

Just scanning the changes, looks to me like you're assuming that the bhazard is applied constantly from baseline, which is only going to be degenerately true. A common case will be an integration of piecewise constant hazards from lifetables, or some other arbitrary baseline hazard function. I don't really see a way around it but to ask for the accumulated hazard directly; easy enough for the simple right-censored case, but with a further argument for the counting form.

[chjackson]

I don't understand what you mean by "applied constantly from baseline, which is only going to be degenerately true"

Doesn't this change implement the likelihood more accurately for the models that exist in the package already?

Allowing the baseline hazard to be supplied in a different form seems like an additional feature to me. It may be a worthwhile feature though if people often have data in that form.

[chjackson]

Ah I get you - it's assuming the hazard is constant when calculating the cumulative hazard...

[chjackson]

I think a piecewise constant baseline hazard would be possible in the current setup in the same way as models with time dependent covariates are implemented - by using the counting process form for the data: Surv(start,stop,status).

The papers for the Stata commands stgenreg and strcs suggest they ignore the constant term in the likelihood.

[heor-robyoung]

I can confirm that stgenreg fails to correctly report the log-likelihood; however, flexsurvreg doesn't match the maximum likelihood value from streg, so I'm less concerned about what Stata does than maintaining internal consistency.

` bhaz <- 0.1 mdl_0 <- flexsurvreg(Surv(time/365, status == 2) ~ 1, data = lung, dist = "exp")

mdl_1 <- flexsurvreg(Surv(time/365, status == 2) ~ 1, dist = "exp", data = lung, inits = mdl_0$res[, "est"] - bhaz, fixedpars = TRUE, bhazard = rep(bhaz, nrow(lung)))

mdl_0$loglik == mdl_1$loglik `

This clearly needs to be true - the same model is being evaluated at the same point.

For the piecewise approach, I still think the cumulative hazard just needs to be asked for. This is asking for a lot of work for the user in dataset manipulation for what is still a special case; the general case is any arbitrary integrable baseline hazard. We either need to ask for the result of that integration directly, or a function by which to integrate.

[chjackson]

Having considered this more, I think the main problem is clear documentation and specification of what the software is designed to do.

Currently flexsurvreg(...bhazard=) is designed to estimate the relative survival function - the fitted parameter values, and predicted survival and hazard curves (returned by summary.flexsurvreg) all refer to the relative survival. Beyond the data, all that is required are the baseline hazards at the event times for people who have had an event.

If, in addition, you want to compute the absolute survival function, you have to supply a baseline survival function, and combine it by hand with the fitted relative survival. Likewise if you want to compute the full likelihood, you have to supply something extra: the cumulative baseline hazard values for everyone.

There is an analogy with Cox regression, which estimates hazard ratios by maximising a partial likelihood, without needing to specify a model for the baseline hazard. The baseline hazard is only needed if you want the full survival curve for a particular covariate value. survival::coxph returns the partial likelihood, and users who need the full survival distribution have to combine the returned coxph object with an estimate of the baseline hazard, e.g. by using basehaz.

So in https://github.com/chjackson/flexsurv-dev/commit/b6e688b968467476fe0baccbda892cded1ca3f40 I've clarified the help page, print output, and added a new vignette section, to explain relative survival models in more detail. I've stated that the likelihood returned by flexsurvreg(...bhazard=) is a partial likelihood, and explained how that could be converted that to a full likelihood, or how the absolute survival function can be obtained. For the likelihood, I'd rather not clutter the interface of flexsurvreg with an extra argument just for the sake of doing such a simple addition.

Though I can believe that people would benefit from nicer tools to deal with these models. I'm not familiar enough with this kind of modelling to be confident proposing a design, so I'm happy for new issues to be submitted to discuss those.

[heor-robyoung]

Hi Chris,

I think I'm happy with this approach. If you could flag up the AIC in the print as well as the log-likelihood, I think we're good to close.

[chjackson]

Thanks, Rob - done in https://github.com/chjackson/flexsurv-dev/commit/0f817fcd49fac2b7e56ae35b6a1726cf1d3d1dbd .