I am trying to find a variable signature associated with a characteristic. Particularly I am looking to get a prognostic model from multi-variable data for gene expression. I have the "Time (survival years)" and "Status (dead/alive)" data for individuals and I also have expression data for those individuals.

Based on my reading about statistical models and posts 1 , 2 and 3, I decided to use multivariate Cox propotional-hazard model. I was reading how to select few best markers, and I found "backward" selection and "LASSO" method to achieve that. Further reading suggested "LASSO" could be a good choice and can be implemented using glmnet package in R.

I divided my data into 60-40 (TRAIN-TEST) proportions, and using the 60% to find those variables and then use them on 40% to find there success as prognostic markers.

The training data has 457 subjects (patients) and 612 covariates (continuous log transformed values). The survivability data is right censored with ~15% Dead cases (Status=1). I don't have any idea about what would be a good number of predictors.

For the TRAIN data, I used the code below.

md #expression data
mdsurv #surviablity data( years) and death (0 or 1)
#creating table
obsLen<-5 # right censored 5 years

keep<-mdsurv$Dead==0 & mdsurv$Status==1

cvfit<-cv.glmnet(mdt,mdsurv, family = "cox",alpha=1)
c<-coef(cvfit ,s='lambda.min')

I am getting 16 genes(covariates) with lambda.min cutoff for coef() and zero with lambda.1se. When I do univariate cox analysis. HR value for some of those genes is less than one. Shall I still use those genes for my test subject or just make a model using the one with HR>=1.

Thanks. P.S. Edited with new information.More edits

  • $\begingroup$ My previous post had a mistake because I had improperly formatted the data and hence was getting incorrect solution. Now I have updated the post. I am getting significant covariates in the LASSO and univariate cox analysis. I am sorry for the confusion. $\endgroup$ – Piyush Joshi Feb 27 '19 at 5:04

First, you should not be worried by the predictors for which hazard ratios (HR) are < 1. That simply means that increased expression of those predictors is associated with better outcome rather than worse outcome. There is no reason to remove them from further use.

Second, with about 70 events in your training set you typically would have power to fit a standard Cox model with 5 unpenalized predictors without overfitting. (Usual rule of thumb is about 15 events per predictor being considered.) LASSO has identified 16 predictors, but LASSO also has penalized their regression coefficients to lower magnitudes than they would have in a standard Cox model, to avoid overfitting. If you wish to use a training/test setup (see below for why that might not be wise) then you should use the penalized coefficients provided by LASSO to evaluate performance. Do not simply use those predictors as the basis of a standard Cox model, as you then lose the protection against overfitting.

Third, with so few events you typically lose too much information by setting aside separate training and test sets. A more powerful approach can be to develop the model on all of your data to obtain a LASSO model. You then demonstrate the quality of your model-building process on multiple bootstrap samples of the data.

See this page and its links for an overview of how you could proceed. Essentially, after you have your model, you repeat the entire model-building process (cross-validation to choose lambda, model with the LASSO predictor selections and coefficients based on that value of lambda) on multiple bootstrapped samples from all the original data. You evaluate the performance of each model on all the original data. Averaging your estimates of performance for multiple models developed on bootstrapped samples provides a useful measure of the quality of the model-building process. Note that with LASSO you will typically get different predictors selected for each bootstrap sample but this procedure will still document the quality of your modeling approach.

  • $\begingroup$ Thanks @EdM for your explanation. However, I am confused with your first explanation. How can they be best predictor but not significant? My readings also suggested that LASSO based cox analysis tends to be less conservative in terms of coefficient, but I didn't find information that running cox with and without LASSO could have different answers. Could you also guide me how can I do the bootstrap modeling or an alternative way of selecting variables for cox analysis. Some post suggested backward stepwise cox model isn't appropriate and suggested LASSO based variable selection. $\endgroup$ – Piyush Joshi Feb 26 '19 at 23:45
  • $\begingroup$ @PiyushJoshi although finding with LASSO a "best predictor but not significant" is no longer an issue in your data, it is a real possibility in general. Say that with the number of events and the variability of outcomes that you need to have an HR of 1.5 to pass a standard test of statistical significance in a Cox model. Say that the best predictor in your data has an HR of 1.2. LASSO could identify that predictor yet it would not pass a standard significance test. I've edited the answer to address your revised question. $\endgroup$ – EdM Feb 27 '19 at 16:23
  • $\begingroup$ Thanks @EdM. I really appreciate your help. I just need two more clarifications. 1) how can generate multiple bootstrapped models. Shall I randomly choose some % (say 60%) of sample, enlarge it original size with repetition, get LASSO based predictor. But then how can I apply it to original data. 2) After multiple bootstrapping, how can I choose final predictors? The ones that are common in all or some other criteria? And what about their coefficient, would that be average of coef obtain from application of all models to the original data? $\endgroup$ – Piyush Joshi Feb 28 '19 at 18:05
  • $\begingroup$ Can I use bootLasso[rdocumentation.org/packages/HDCI/versions/1.0-2/topics/… package and/or BootValidation [cran.r-project.org/web/packages/BootValidation/… package for the algorithm you are suggesting? $\endgroup$ – Piyush Joshi Feb 28 '19 at 19:12
  • $\begingroup$ @PiyushJoshi I don't have experience with BootValidation but it seems at first glance to do what I suggested. (HDCI doesn't seem to handle Cox models.) The predictors that you report would be those from your original model, with bootstrapping to document the quality of your model-building process. Note that a bootstrapped sample is the same size as your full data set but taken with replacement from your data so that about 1/3 of cases are duplicates. If that package doesn't do what you need, the boot package in R provides a general way to proceed. $\endgroup$ – EdM Feb 28 '19 at 19:31

I am doing Cox regression analyses using glmnet myself, and ran into this question. The penalization method used in Cox models built by glmnet is more generally elastic net. Elastic net could be defined as a mix between LASSO and ridge-regression, and regularizes the data to avoid overfitting as well. I think it would be helpful to ascertain this information, in case you are writing a scientific paper, or something like that.

The reference is: https://www.jstatsoft.org/article/view/v039i05

The alpha value setting determines the mix between LASSO and ridge-regression used by glmnet. The default is alpha=1, which is pure LASSO as used in this question. Other values can be chosen: alpha=0 is pure ridge regression, while intermediate values combine LASSO and ridge.

  • $\begingroup$ You are correct that the glmnet() function does allow for elastic net. That's a strength of the package. The tradeoff between LASSO and ridge regression is determined by the alpha parameter setting in the call to the function. The default value for that setting, however, is alpha=1, which is all LASSO, no ridge, according to the manual. That setting was used explicitly in the OP. Unless you provide a different value for alpha, the function will do a LASSO fit by default. $\endgroup$ – EdM Feb 14 at 18:14
  • 1
    $\begingroup$ Thank you very much for the clarification, @EdM! $\endgroup$ – Fla28 Feb 14 at 18:26
  • $\begingroup$ There's no need to apologize. Your answer made an important point--glmnet() does do elastic net, it's just that the LASSO extreme is what it chooses by default. And your link to the paper is helpful. Editing an answer when you realize that a prior version was incomplete or misleading is an important way that contributors to this site learn more about statistics, one that I have frequently engaged in and for which no apology is necessary. I'll edit your question to suggest a way to proceed, without apology. $\endgroup$ – EdM Feb 14 at 21:54
  • $\begingroup$ The edit is great, thank you! $\endgroup$ – Fla28 Feb 15 at 11:08

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.