Is it possible to perform lasso regression (glmnet with "cox") for variable selection and then conduct Cox regression using selected variables?

What is the difference between analyzing with lasso regression only AND Cox regression with selected variables? I want to use Cox regression which has more functions in post-prediction.

Can anyone provide an example of the predict of the glmnet Cox model? In the glmnet vignette on CRAN, there are linear and logistic examples, but there is no example of Cox.

  • $\begingroup$ Please say more about the scale of your study (numbers of cases and events, number of potential predictors) and the nature of the "post-prediction" you wish to perform. The best answer to your question might be different depending on those details. $\endgroup$
    – EdM
    Jul 5 '19 at 16:38
  • $\begingroup$ Dear Mihai Chelaru, Thank you for your comments. I have 208 cases and 35 events (6months to 10years survival data). The number of potential predictors is about 40. I just want to get Harrel's C index with confidence interval. $\endgroup$
    – Tae Lee
    Jul 6 '19 at 6:09

If you simply take the predictors returned by LASSO and use them in a fresh Cox model, you have ignored the fact that you used the data to select the predictors. The model thus may be overfit; p-values and such will not be reliable. You will have removed LASSO's penalization of regression coefficients that minimizes overfitting. See this answer for some more details.

The vignette for Cox modeling with the glmnet package claims that "the Cox Model is rarely used for actual prediction" so the authors haven't provided prediction functionality for Cox models beyond the linear-predictor values and the related "relative-risks." There seems to be no way to extract a full Cox model object or its baseline cumulative hazard from glmnet.

You have a few ways to proceed.

First, if all you care about is the C-index as suggested by your comment, you can get the information you need to calculate it from the linear-predictor values returned by predict.glmnet. The C-index is simply the fraction of pairs of comparable cases for which the event order matches the order of linear-predictor values. Although censored cases can't be compared against each other, all cases with events can be compared against each other (in your data, 35*17 = 595 pairs) and any case with an event can be compared against censored cases having later censoring times. This page shows some ways to do this calculation yourself.

Second, you can force coxph to provide a Cox model with the coefficients returned by glmnet. This link outlines how. Instead of using the default coxph settings, you provide a vector of the glmnet coefficient values with the init= argument, and then prevent any additional fitting with an iter.max=0 argument, which is passed on to the coxph.control function. The output will be a coxph.object suitable for many purposes but it will not have a coefficient covariance matrix, so you may get error messages in some applications and functionality for things like confidence intervals (questionable in any event for a model fit by LASSO) is limited.

Third, if you don't need anything beyond simple LASSO you could try the glmpath package. It might be slower than glmnet but with a small study like this that shouldn't make much difference. Its predict function for coxpath objects has a type="coxph" argument option that returns a coxph.object for which "the components of a coxph object such as coefficients, variance, and the test statistics are adjusted to the shrinkage." So that might be the most generally useful approach for simple LASSO. I haven't used this myself, and I don't know how those adjustments for shrinkage are made. You will need to make sure that the shrinkage value specified by the s= argument matches the mode= argument for specifying the single point along the regularization path for which you want the Cox model.

With your original question answered, however, do think more about whether LASSO is the best approach to answering your scientific question. With only 35 events LASSO is probably only returning (at the optimal cross-validated partial likelihood) about 3 non-zero coefficients among your 40 possible predictors. In clinical data there are often highly correlated predictors, so the particular predictors chosen can depend heavily on the particular data sample at hand (which you can examine by repeating your modeling on multiple bootstrap sample of your data). And if you have a multi-category predictor simple LASSO might just choose one category, which might not always make sense.

You might be better off using your knowledge of the subject matter to combine similar predictors and/or select predictors that are most likely related to outcome. If you have a specific hypothesis about one predictor that you would like to test but need to correct for other covariates, look at the recommendations in this paper; in that case, using ridge-regression penalties for all predictors except for that one might be better.

  • $\begingroup$ I really want to express my gratitude for your kind and profound answer!! $\endgroup$
    – Tae Lee
    Jul 6 '19 at 19:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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