# Model evaluation for feature selection

I have a dataset of gene expression data and I'm trying to find genes related to particular diseases. My labels are dichotomous (sick - not sick) and I used a Logistic regression with LASSO regularization in order to extract meaningful features (basically taking all the genes with coefficient different from zero). Hyperparameter lambda has been found using leave-one-out cross-validation. In order to find my coefficient should I train a new model using my best lambda on all the data?

I think there are some drawbacks in my approach but since I'm relatively new to feature reduction I cannot formulate a query on google that tells me if or where I am wrong. Can you help me?

• How many genes $p$ do you have, and what's the number $n$ of experiments? – Edgar Oct 4 '19 at 12:26
• 4300 genes and 260 samples – Giuseppe Minardi Oct 4 '19 at 12:37

1. It's better to do the standard $$10$$-fold cross-validation in your case, since you have $$n=260$$ observations. Leave-one-out cross-validation is more appropriate if the $$n$$ is smaller.

2. You're completely right, after the cross-validation you need to train a new model on the full data with the hyperparameter $$\lambda$$ found by the cross-validation.

3. If you're using R for this analysis (which is highly recommended), everything should be easy to implement with the cv.glmnet() function from the glmnet package, with parameter family="binomial" for the logistic regression.

4. Take care that you normalized your gene expression values reasonably. Look up "RPKM" and "TPM" for a further understanding of how to model with gene expression values.

Below a minimum working example based on your comments, this should give you a smooth start with the cv.glmnet function (X is your gene expression matrix, y your 0-1-encoded disease status):

library(glmnet)

my_glmnet <- cv.glmnet(X, y, family="binomial") # default CV with 10 folds

plot(my_glmnet) # here you can see the CVMSE and the optimal lambdas, as well as the number of non-zero coefficients at the top

which_opti_lambda <- which(my_glmnet$$lambda==my_glmnet$$lambda.min) # this gives you the index of the optimal lambda (you can also try lambda.1se)

opti_coefficients <- my_glmnet$$glmnet.fit$$beta[, which_opti_lambda] # get the coefficients of the final run for the optimal lambda

which(opti_coefficients!=0) # these are the genes that were chosen by lasso cv logistic regression


• Yes, I'm using R and yes, I'm using glmnet. My advisor told me to use cross-validation also in order to test my model but I do not know how to use cv.glmnet in order to do cross-validation without estimating lambda. he function only works using different values of lambda. I'm currenntly using a simple train-test split but my advisor told me to use cross validation – Giuseppe Minardi Oct 4 '19 at 12:56
• cv.glmnet does all the lambda calculations for you, you can recover the lambda values from the glmnet object, as well as the lambda.min (or lambda.1se) value for which the optimal cross-validation error was achieved. – Edgar Oct 4 '19 at 12:58
• I know, but should I train my model again with my best lambda or can I use the mse given by cv.glmnet? – Giuseppe Minardi Oct 4 '19 at 16:27
• If you do 10-fold CV with cv.glmnet, it gives you the fitted object for the (11th) full model, that means, it trains the model again for you already (but with a full lambda sequence, too). The best lambda is either lambda.min or lambda.1se or any lambda of your choice -- you will have to find the appropriate index in the lambda sequence and choose the corresponding column of the coefficient matrix. – Edgar Oct 4 '19 at 16:40
• I updated my answer with a minimal working example. – Edgar Oct 4 '19 at 16:57