Allow flexible specification of how to center the linear predictor

[therneau]

[Edit] Allow a flexible specification of the reference observation, most likely by providing a data frame, and re-evaluate the choice of default for it.

---

This really belongs as an email to the maintainer, but finding an email address has been elusive. Hopefully I can get an email back and start a conversation.

I saw a reference to censored tidymodels when browsing the rstudio conference, which led me to https://www.tidyverse.org/blog/2021/11/survival-analysis-parsnip-adjacent/

There are some very interesting ideas there. One very bad one, however, is the following sentence. " While we think through a more general approach to specifying such a reference observation or baseline hazard, we ensure that all engines in censored use the same approach: a reference observation of 0". I've been down that road, and it often ends in overlow/underflow of the exp() function. On the other side of the coin, using predict() to get a survival curve, rather than a survfit method as coxph does, is the design I would also use (if I could go back in time.)

In general, I think it would be useful for the package authors to engage in a proactive conversation with the developers of the main packages you are building off of, don't just hope we all will notice your blogs and jump on. I've made more than a few mistakes over the years and there is no reason for you to re-invent them. An interactive forum may well be a better format for this than email, but I don't naturally live in that space so need targeted advice on which and how.

Terry T., author of rpart and survival

[hfrick]

Hi Terry, thank you for the feedback! We are all out of office this week but will respond once we're back.

[therneau]

One of the few things I did right in my long career was that "out of the office" = "no email". Neither read it or respond. Enjoy your time away.

[hfrick]

Thanks for the understanding! Regarding the centering of the linear predictor: I've read through the Details section of ?predict.coxph which also brings up the numeric over/underflow

Second, downstream calculations depend on the risk score exp(linear predictor), which will fall prey to numeric overflow for a linear predictor greater than .Machine$double.max.exp.

Could you elaborate more on which downstream calculations you mean? I'm trying to understand better in which use cases this would be particularly biting and how common those are in a prediction context.

An aspect that is particularly important to us for {censored} is consistency. We are currently not taking advantage of predict.coxph()'s reference argument in favor of consistency with predictions from the other model types and implementations we wrap. As mentioned in the blog post, we'd eventually like to offer a general approach to specifying a reference observation or baseline hazard. Do you have any additional thoughts on different use cases for baseline hazard parametrizations? I've seen that predict.coxph() offers overall means and strata means.

[therneau]

I'm still trying to find out how to have a more general conversation. A bug list just doesn't seem correct. Comments a. There is no vignette here, nor any reference to an external overview document.
b. I just saw the R-conference talk. I saw a couple of good things, and at least one that made me cringe. What is the forum to communicate that? c. What is the forum, if any, to talk about high level issues? (I avoid mailing lists, in general, due to bandwidth, and by default not following Rstudio ones). d. What group is providing the expert input on survival analysis? I think that such input is needed, at the design phase.

I will try to gather a reproducable example for you wrt centering. Detailed discussion is found in the noweb source of the survival package. Some of it is pretty nitpicky, written to remind me of nasty data sets I've seen.

[hfrick]

Given that censored is an extension package for parnsip, there is a lot more documentation on parsnip than on censored. The parnsip pkgdown site at https://parsnip.tidymodels.org contains a Get Started section as well as links to further articles (via the navbar on top). That said, the pkgdown site for censored contains the article Fitting and Predicting with censored which illustrates basic functionality for survival models.

We track any issues with the package here so this is a good place. We typically use RStudio Community as a secondary place for discussions, the blog post you originally read links out to a discussion thread there. That would be a good venue as well!

Thank you for the pointer, I'll take a look at the noweb folder of survival then!

[therneau]

The issue of centering is one that resurfaces every dozen years in a new guise. As an example look at the longley data (type "longley" in R); James Longley, An appraisal of least squares programs from the electronic computer from the point of view of the user, J Amer Stat Assoc, 1967. In short, Longley found that for this data set, which was being used in numerous papers, different computing systems at the universiy gave quite different answers, some of which had 0 digits in common. You can see the problem by a simple graph: fit <- lm(Employed ~ Population, longley) plot(Employed ~ Population, longley, xlim=c(0, 131), ylim=c(8, 70)) abline(fit)

The intercept at x=0 is very far from the data, leading to large correlations and consequent problems with the fit. Standard advice has long been to recenter covariates like this, one one hand, while on the other this particular data set motivated major updates in the computational medthods for linear models. Using the QR method, as lm does, tames the issue (for this data set at least).

The same issue arises in a Cox model, but made worse by the presence of exp(). Internally, the coxph routine recenters covariates, and in fact sometimes needs to re-recenter them during the course of the compuation. On github.com/therneau/survival is a directory 'fail' that documents some of the data sets that have driven these internal changes. But none of this can fix a basic issue with survival curves, namely the well known equation below shows that if you have the predicted survival curve for some covariate set "x0", you can easily obtain the curve for any other covariate choice "x". log(S(t;x)) = log(S(t; x0)) + exp(beta *(x - x0)) This leads a lot of people to the idea of saving just one curve and reusing it. This is exacerbated by the tendency of textbooks to call x0= 0 the "baseline" survival curve.

If you have any covariates in the model that are far from 0 in a similar fashion to the longley variables (small sd relative to the mean), and use x0=0, then the exp() part of the above equation can overflow. Worse, the compuation of S(x0) itself is likely to extremely inaccurate or even fail outright. My strong advice has been to let the survfit() code produce predicted survival curves. It is fast enough that the efficiency argument has no legs (and if you also want std deviations there isn't a shortcut formula for that part anyway), and more importantly it is reliably accurate. As a package author, accuracy is more important than convenience --- I don't want to ever be behind a false medical 'result' that kills people.

If you absolutely must use the shortcut, then pick a rational x0. The coxph result includes a "mean" component and uses that to recenter. It is not at all critical to use a mean: any centering point within the range of the covariate is fine; and yes/no or categorical variables don't need to be recentered at all.

[hfrick]

The choice of the reference observation here was made with the aim of standardization and as a temporary choice on the way to a more flexible specification of that reference observation. I'll edit the issue slightly to reflect that. Thank you for providing additional context!