3
$\begingroup$

I am trying to fit a survival model in a loop between a set of ~3000 proteins and a disease outcome, with technical covariates of those proteins that are important to their measurement/estimation as covariates. The model, generally, looks like this:

survreg(Surv(time_to_event, diagnosis) ~ Age + Sex + Location + Ethnicity + Leukocyte_Count + Platelet_Count + Monocyte_Count + Lymphocyte_Count + Neutrophil_Count + Basophil_Count + Monocyte_Count + protein, data=data, na.action = na.exclude)

These data are coming from the UK Biobank, and while there is not perfect overlap between the outcome data and proteomics data, we still have an n>10k individuals with several hundred cases for the disease outcome. Importantly, I am using exactly the same set of participants for each run of the model - the only thing that changes is the protein measured. For some proteins in the set, everything works perfectly and the model completes without issue; for others, it seems that regardless of the number of Newton-Raphson Iterations it goes through (I have tested up to 1000 iterations), the Log(scale) term never converges (-Inf Z-score); and for others still, I get the following error:

Error in if (any(singular)) fit$coefficients[singular] <- NA : missing value where TRUE/FALSE needed

I have looked at the survreg source code and identified that these lines of code are throwing this error:

# Do this before attaching the na.action, so that residuals() won't
    #   reinsert missing values under na.exclude
    if (robust) {
        fit$naive.var <- fit$var
        if (!model) fit$model <- m  #temporary addition, so resid doesn't
                                    # have to reconstruct
        if (length(cluster))
             fit$var <- crossprod(rowsum(residuals.survreg(fit, 'dfbeta'), 
                                         cluster))
        else fit$var <- crossprod(residuals.survreg(fit, 'dfbeta'))
        if (!model) fit$model <- NULL  # take it back out
        }

    # set singular coefficients to NA
    #  this is purposely not done until the residuals, etc. are computed
    singular <- (diag(fit$var)==0)[1:length(fit$coefficients)]
    if (any(singular)) fit$coefficients[singular] <- NA

My relatively naive interpretation of this is that fit$var is NULL for certain proteins, but I'm not sure why this would be the case. I have tried fitting a coxph() model instead of survreg on the same outcome and predictor data and the coxph model completes entirely without issue. However, for my uses I must use survreg.

Does anyone have any insight into why I could be getting these non-convergence warnings and the error in question? I understand that it can be difficult or impossible to diagnose errors like this without the data itself, but I don't know how to even begin to diagnose the cause of this error because I can't look at how the model is being fit (no model object is created with the error in question), so any common sense advice on what to look for would be extremely helpful.

Edit: apologies for the original lack of clarity - my goal here it to perform mediation analysis with the mediate package. For time to event data, the mediate package can only take survreg models, not cox models. Therefore, for the analysis I am currently doing, I must use the survreg function to fit the model.

$\endgroup$
2
  • $\begingroup$ This does look weird. 1st step would be to make sure everything, i.e. packages and R itself, are up to date. Then check out. Then we can check the issues on the github I find this: github.com/therneau/survival/issues/237, so does your model also not work when you fit outside of a function in a fresh environment? $\endgroup$ Mar 4 at 15:57
  • $\begingroup$ @LukasLohse Thank you for the response! I have tried it both locally and on a remote cluster with a fresh installation of the relevant packages and most up to date version of R, the results are consistent - i.e. the same proteins cause the error $\endgroup$
    – womy
    Mar 4 at 17:24

1 Answer 1

5
$\begingroup$

This is probably not a good way to accomplish what you want, for several reasons.

First, you might not even have enough disease cases to do this properly. You have 11 separate predictors beyond the protein value, of which 7 are continuous values (blood counts). To avoid overfitting in a clinical survival model, you typically need about 15 events per coefficient that you are estimating. A single model as you have written it thus needs at least 180 events, more depending on the nature of the Location and Ethnicity predictors. So a lot depends on just how many "several hundred cases" represent.

Second, your model as written assumes a strictly linear association between the time-acceleration (for accelerated failure time models; or log-hazard for proportional hazard models) of developing the disease and each of the continuous predictors separately, including the protein value. That's highly unlikely to be the case. There are ways to model continuous predictors flexibly; see Section 2.4 (among others) of Frank Harrell's Regression Modeling Strategies. Those approaches require estimating more coefficients, however, and thus more events unless combined with some penalization (e.g., ridge regression). Or you might consider a way to combine the various blood counts into a smaller set of predictors; see Harrell's discussion of data reduction.

Third, it's not at all clear why you "must use survreg." What you are fitting with survreg() under its default settings is a Weibull model, which is a specific form of a proportional hazards model, one that also happens to have an interpretation in terms of accelerated failure time. A Cox model allows for a more flexible baseline hazard without that particular (and possibly unrealistic) restriction of the shape of the baseline hazard. As Cox models seem to work OK for you, part of the problem here might be your trying to force data to fit a baseline hazard function that is inappropriate.

Fourth, with 3000 proteins evaluated separately you have a massive problem with multiple comparisons. Even with false discovery rate control you are going to need some of your 3000 proteins to be extraordinarily "statistically significant." That will be particularly hard to accomplish if you stick with your current implicit assumption that each protein has a strictly linear association with the log-hazard/time-acceleration of outcome.

Finally, to address your specific question: coxph()internally centers and scales continuous predictors to avoid convergence problems arising from the exponentiations involved in fitting. I'm not sure whether survreg() also does that. Try centering and scaling protein values yourself. If that doesn't work, then look at the distributions of protein values for those that don't converge, and consider whether you really "must use survreg."

$\endgroup$
5
  • $\begingroup$ 1) I have over 300 cases, but I have also tried this model reducing the number of predicting covariates - this has resolved the error in some cases, but not in others; 2) Thank you for the suggestion of combining some of the covariates, I will attempt this! 3) I have tested a few underlying distributions, including exponential and weibull. Is there some way to check which distribution would be best? I must use "survreg" because I am doing mediation analysis with the mediate package, it only takes survival models fit with survreg 4) I am doing FDR correction, many of the associations remain $\endgroup$
    – womy
    Mar 4 at 17:33
  • $\begingroup$ Thank you for your response, I was running into character limits in my first comment. As a follow up question, if the fundamental issue is related to having too many covariates or that it is improper to force my data to fit the weibull distribution, why would the model complete for some proteins but not others? If my approach is flawed in these ways, shouldn't those issues appear in every model since the time to event data is being forced to the same distribution for all models? Only the protein is changing. $\endgroup$
    – womy
    Mar 4 at 17:58
  • 1
    $\begingroup$ The fundamental statistical issue, as I see it, is that you are probably not using the correct functional forms for relating continuous variables to outcomes.That's separate from numerical convergence problems. Even if an incorrect model converges the results won't be useful, particularly if you are attempting something so difficult as mediation analysis. Chapters 18 and 19 of Harrell's Regression Modeling Strategies discuss choice of parametric model families; plots of "survival curves" of censored residuals are helpful. $\endgroup$
    – EdM
    Mar 4 at 19:37
  • 2
    $\begingroup$ For mediation analysis, you maybe shouldn't blindly rely on the mediate package. Please read the extensive coverage of evaluating mediation and how complicated things get in survival studies, in Causal Inference: What If by Hernán and Robins. Consider getting some highly experienced local statistical consultation to keep you from going astray. $\endgroup$
    – EdM
    Mar 4 at 19:41
  • 1
    $\begingroup$ Thank you very much for the useful resources and valuable suggestions! $\endgroup$
    – womy
    Mar 4 at 21:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.