# Logistic Regression to identify significant covariates(variables)

I am working on logistic regression to identify the most cucial genes that can predict the response of stimuli.

I have a workflow that looks fine to me but I want inputs from the people who have more experience on similar problems.

I have a data set that looks like below:

          Resp    CDK6 EGFR  KIF2C CDC20
Sample 1  pos  11.39 10.62  9.75 10.34
Sample 2  pos  10.16  8.63  8.68  9.08
Sample 3  pos   9.29 10.24  9.89 10.11
Sample 4  neg  11.53  9.22  9.35  9.13
Sample 5  neg   8.35 10.62 10.25 10.01
Sample 6  pos  11.71 10.43  8.87  9.44
...


This a table of dimension 130 * 82.

What I have done is calculated each genes individually. Below is the code for that and selected the ones with that gave p-value <= 0.05.

 glm=glm(Resp ~ CDK6, data = dataframe, family = binomial(link = 'logit'))


Once I have x number of genes (say 35 genes) that have p-value <=0.05, I use the model/code below to extract the most significant signatures of response.

glm2 <- glm(Resp ~ Gene1 + Gene2 + Gene3 + Gene4 + Gene5 , data = dataframe, family = binomial(link = 'logit'))


The above model gives me a 9 genes that are significant with p-value <=0.01.

Then I perform cross validation, like below

cv.glm(dataframe, glm2, K=nrow(dataframe))$delta  This was followed by confusion matrix and sensitivity & specificity scores  threshold=0.5 predicted_values<-ifelse(predict(glm2,type="response")>threshold,1,0) actual_values<-glm2$y
conf_matrix<-table(predicted_values,actual_values)

> sensitivity(conf_matrix)
[1] 0.9153846
> specificity(conf_matrix)
[1] 0.9666667


My Question is, I started with 82 genes/variables and I got it down to 9 genes with the above approach. Is this approach ok? or am I doing something wrong or this can be further tuned up?

• Using p-values to select features tends to be problematic. Are you looking to make accurate predictions, or just check which genes are most strongly associated with the outcome, or perform causal inference? In the first two cases, you might want to consider using LASSO (as in the glmnet package) to select relevant features, instead of p-values. Commented Aug 10, 2023 at 20:46
• Hi, Thank you for your comment I tried glmnet it did not work well with my dataset. the SE output came as "." for some strange reason. I am looking for genes that can predict the outcome of the dosage. In addition, I googled and stumbled upon this post bioinformatics.stackexchange.com/questions/6813/… I have a kind of similar problem but not the same (My dataset does not work with glmnet) and if choose a model to put all genes together in first step istead of 1 by 1 the model doesn't converge (it does with gaussian) Commented Aug 10, 2023 at 20:57
• LASSO ought to work well with sparse and "wide" data that I've seen from genetics problems before. You aren't usually very interested in the standard errors with LASSO, you just look at the regression coefficients. As long as they're all on the same scale, the coefficients of the least important inputs will be reduced to 0, and the most coefficients of the most important inputs will be left with the largest coefficients. Commented Aug 10, 2023 at 21:29
• I'm also a little confused about what's happening in the question you linked to. You are not supposed to use LASSO to choose features for a separate model, it doesn't work like that. Your dataset probably works fine with LASSO. This is starting to seem like a completely separate question on a separate topic, however. Commented Aug 10, 2023 at 21:31
• Thank you @shadowtalker I have sucessfully implemented Lasso and I got the results that seemed elusive. Commented Aug 11, 2023 at 14:28

Another approach may be to use a Bayesian model. A laplace prior, centered at 0, in a Bayesian model has similar properties to a lasso penalized model, The scale parameter of the laplace prior is related to the penalty in the lasso (you will probably want to standardize or otherwise scale your inputs). Unlike glmnet, most Bayesian tools will not fit the multiple lambda values in one step, you would need to fit each penalty/scale individually. You could also look at the horseshoe prior which pulls very small coefficients towards 0 while not making much change to the larger ones. The brms package has an implementation of the horseshoe prior (other packages do as well, but brms is the one that I remember specifically).