Answers.Body: <p>Final answer. Keep in mind I work in Python so I'm trying to translate here.</p>
  <ol>
  <li>The output of <code>coef(cv.lassoModel)</code> ... those are your genes of interest thats your <strong>feature selection</strong>. <code>head</code> I assume gives you the top values ... I think thats right (I need to check its not <code>coef(lassoModel)</code>) ... Basically thats where your results are.</li>
  </ol>
  <p>This might be over simplistic but <code>co.se</code> without the <code>head</code> is your answer. Not every gene will carry a weight BTW. If the entirety of that is (700) <code>.</code> ... try point 2 below and then point 4. I don't think that will happen (see <strong>Whats happened</strong>)</p>
  <ol start="2">
  <li>I would redo this with <code>alpha=0</code> (ridge regression) and <code>alpha=0.5</code> elasticnet, but you don't need to</li>
  <li>Generate an accuracy score</li>
  </ol>
  <pre><code>    y_predicted &lt;- predict(best_model, s = best_lambda, newx = x)
      sst &lt;- sum((y - mean(y))^2)
      sse &lt;- sum((y_predicted - y)^2)
      rsq &lt;- 1 - sse/sst
      rsq
  </code></pre>
  <p>Hopefully this is large e.g. around 0.8 (in Python this automatic BTW). You could ignore that, but its what I'd do</p>
  <ol start="4">
  <li>The intercept is simply the weight when your gene expression score is 0.</li>
  <li>Lasso regression handles sparse data - so no worries about <code>&quot;.&quot;</code>. I don't whether that Python's version of &quot;nan&quot; (not a number), or if its zero.  ... @haci knows what this means and his answer is authoritative for the second part of your question.</li>
  <li>In Python transformations are common, but that doesn't seem an option here, so find providing your accuracy score is reasonable, personally I'm happy.</li>
  </ol>
  <p>In Python what you want is feature selection
  <code>print(pd.DataFrame({'Features': X_train.columns, 'Coefficients': lasso.coef_}))</code></p>
  <p><strong>Whats happened?</strong> The coefficients for <code>head(co.se)</code> have been shrunk to zero, thats only the top 10? results, you've 700 weights. The full list <em>should</em> (famous last words) contain a minority of really big coefficients ... those are really important and those are your answer. So by printing only the top few results ... not the whole lot the results are still there you've just not seen them. Thats it, basically you were almost there and you could happily settle on that as the final answer. Basically, alot of genes ain't doing much a few genes are doing all the work - thats the answer. Below I stress repeatedly I don't think anything is wrong here and I personally like the result.</p>
  <hr />
  <p>Background</p>
  <p><strong>Lasso regression</strong>, it's a cousin of ridge regression and key parameter is lambda (who'd guessed). What your plots are doing is tuning the parameter and finding the associated weight individually or collectively. Fairly straight forward. So like whats the value of lambda for your data set ('answered' below)?  It's worth doing some reading on how lambda is calculated to understand the plot, if you're not sure how it works.</p>
  <p><strong>Imbalance</strong>, well to be honest I'm not convinced there's anything wrong with your analysis. Okay there's a lack of stationarity for some of your weights ... but some are fine and gist is pretty clear. Yeah okay there's a problem with stabilisation of lambda - thats not necessarily imbalance.</p>
  <p><strong>Imbalance</strong> is when it's 1 in 10 and over. You've got 1 in 6, 1 in 7, not huge. However, - <em>you might be right</em> - and if you're worried simply repeat the analysis excluding <code>RP</code> (value 7). Then check whether the parameterisation plot (lambda against weight) is more stable or more &quot;homogeneous&quot; (multinomial plot), i.e. the genes behaviour in a more uniform way. If that is does work - it gives you nicer results AND <code>RP</code> is essential to the analysis thats a separate question about augmentation. That is complicated. For example, the multinomial plot with a 3 target categorical model might give a clear value for lambda.</p>
  <p><strong>Transformation</strong> Here, I would suggest if you want to stabilise the multinomial plot (if it's possible) I would look at a transformation. The data you presented doesn't look like it needs a transformation, but we can't see the whole data. In Python there's a <em>max-min</em> transformation - thats where I might start.</p>
  <p><strong>Overall good results</strong> Personally, I'd say these were good results, except the weights ain't homogeneous for some genes (most are fine BTW), but the patterns clear. The <em>red</em> (top), <em>blue</em> (bottom) and cyan (negative y-axis) are real results showing that within the 700 genes <strong>3 genes</strong> are having a massive impact. Sure you have positive and negative and thats how it works one gene is positively contributing to the classification, others are negatively contributing to it. Then you've the black and 'deep pink' lines which have stabilised around 5E4. <strong>From 700 genes down to 5 genes</strong>. Whats wrong with that?</p>
  <p><strong>Whats the answer to the data set?</strong></p>
  <p>Its ...</p>
  <blockquote>
  <p>NRT
  RT
  NRP
  RP</p>
  </blockquote>
  <p>... is due to 'red', 'blue', 'cyan', 'black' and 'deep pink' genes. Thats it.</p>
  <p><strong>Frank assessment</strong> To be perfectly honest, looks good. No-one is going to say &quot;Wooaaa lambda ain't stable, destroy! ;-)&quot;, they'll say &quot;Wow thats cool and clear&quot;.</p>
  <hr />
  <p><strong>How I'd do it..</strong> I would look at basic ML robustness statistics and use this to assess the model rather than the parameterisation of lambda. Its horses for courses as they say. Upvotes.</p>
  <p><strong>Note</strong> I use Python not R for ML, both work well however. <code>glmnet</code> within ML is new to me BTW.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-22T09:02:50.540Z
Answers.Text: Thank you for your reply and detailed explanation. What I do not understand is why I am getting output as Intercept instead of gene names.
Answers.UserDisplayName: Angelo
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-22T09:43:27.533Z
Answers.Text: Another puzzling aspect is when I do co.se$RT!=0, I get a gene list with "." next to it :(
Answers.UserDisplayName: Angelo
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-22T11:14:58.627Z
Answers.Text: Technically you are getting a `sparse Matrix` with `co.se$RT!=0` and the "dots" are representing zeros or in this specific case `FALSE` values. You can see yourself by `as.matrix(co.se$RT!=0)`.
Answers.UserDisplayName: haci
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-22T11:54:27.233Z
Answers.Text: All of them :( I tried to normalise the values by getting zscores. HES4 is a gene name and its reading in different samples.
Answers.UserDisplayName: Angelo
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-22T13:59:23.700Z
Answers.Text: @Angelo I've addressed the comments in the first part of the answer.
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-22T04:39:07.903Z
Answers.LastActivityDate: 2023-07-22T16:23:48.990Z
Answers.LastEditDate: 2023-07-22T16:23:48.990Z
Answers.UserDisplayName: M__
Body: <p>I am new to glmnet and other ML techniques, so apologise in advance if this sounds a trivial question. However, your guidance is very much appreciated.</p>
  <p>So, I have nearly 700 genes with their expression profile across 4 conditions. I want to identify signature genes for each of the 4 conditions. I am using glmnet but I am not getting anywhere near the results.</p>
  <p>so my data looks like this</p>
  <pre class="lang-r prettyprint-override"><code>           Subtype     HES4   FNDC10    PRKCZ PRKCZ.AS1
  Sample1     NRT 5.452955 6.241387 8.546059  6.260245
  Sample2     RT 5.125350 7.038719 8.900024  6.267325
  Sample3     NRT 5.689560 5.903811 8.029434  4.213996
  Sample4     NRP 5.834468 5.722589 8.061339  5.072274
  Sample5     RP 5.401623 7.272111 9.108342  6.599629
  </code></pre>
  <p>purely in terms of numbers, the subtypes have following number of samples.</p>
  <pre><code>NRT 52
  RT 42
  NRP 30
  RP 7
  </code></pre>
  <p>Q1. May be the groups are imbalanced and caution should be taken?</p>
  <p>I have applied glmnet as following</p>
  <pre class="lang-r prettyprint-override"><code>modellinglasso=read.table(&quot;LRMatrix.txt&quot;, header=T, row.names = 1, sep=&quot;\t&quot; )
  modellinglasso$Subtype=as.factor(modellinglasso$Subtype)
  modellinglasso$Subtype
  #build a single pass (single fold) lasso-penalised model
  lassoModel &lt;- glmnet(
    x=data.matrix(modellinglasso[,-1]),
    y=modellinglasso$Subtype,
    standardize=TRUE,
    alpha=1.0,
    family=&quot;multinomial&quot;)
  plot(lassoModel, xvar=&quot;lambda&quot;)
  Warning message:
  In lognet(xd, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  :
    one multinomial or binomial class has fewer than 8  observations; dangerous ground
  </code></pre>
  <p><a href="https://i.sstatic.net/0lUP0.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/0lUP0.png" alt="enter image description here" /></a></p>
  <pre class="lang-r prettyprint-override"><code>#perform 10-fold cross validation
  cv.lassoModel &lt;- cv.glmnet(
    x=data.matrix(modellinglasso[,-1]),
    y=modellinglasso$Subtype,
    standardize=TRUE,
    alpha=1.0,
    nfolds=10,
    family=&quot;multinomial&quot;,
    parallel=TRUE)
  # plot variable deviances vs. shrinkage parameter, λ (lambda)
  plot(cv.lassoModel)
  </code></pre>
  <p><a href="https://i.sstatic.net/E02wY.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/E02wY.png" alt="enter image description here" /></a></p>
  <pre class="lang-r prettyprint-override"><code>idealLambda &lt;- cv.lassoModel<span class="math-container">$lambda.min
  idealLambda1se &lt;- cv.lassoModel$</span>lambda.1se
  print(idealLambda); print(idealLambda1se)
  
  [1] 0.1186125
  [1] 0.1186125
  
  # derive coefficients for each gene
  co &lt;- coef(cv.lassoModel, s=idealLambda, exact=TRUE)
  head(co)
  
  co.se &lt;- coef(cv.lassoModel, s=idealLambda1se, exact=TRUE)
  head(co.se)
  
  # identify predictors for each sub-type
  rownames(co.se<span class="math-container">$PR)[which(co.se$</span>NRT != 0)]
  
  [1] &quot;(Intercept)&quot;
  </code></pre>
  <p>output for</p>
  <pre><code>co.se$RT!=0
  6 x 1 sparse Matrix of class &quot;lgCMatrix&quot;
              1
  (Intercept) |
  HES4        .
  FNDC10      .
  PRKCZ       .
  PRKCZ.AS1   .
  MMEL1       .
  </code></pre>
  <p>It doesn't matter which subtype I put, I get the same output &quot;(Intercept)&quot; instead of gene names.</p>
  <p>Any help to get the gene signatures for the subtypes is sincerely appreciated.</p>
  <p>Thank  you</p>
  
Tags: <r><rna-seq><machine-learning>
Title: Applying glmnet to identify predictors for subtypes
ThreadText: # Applying glmnet to identify predictors for subtypes
  
  (Asked by: username_0 on 2023-07-21)
  
  I am new to glmnet and other ML techniques, so apologise in advance if this sounds a trivial question. However, your guidance is very much appreciated.
  
  So, I have nearly 700 genes with their expression profile across 4 conditions. I want to identify signature genes for each of the 4 conditions. I am using glmnet but I am not getting anywhere near the results.
  
  so my data looks like this
  
  ```r
             Subtype     HES4   FNDC10    PRKCZ PRKCZ.AS1
  Sample1     NRT 5.452955 6.241387 8.546059  6.260245
  Sample2     RT 5.125350 7.038719 8.900024  6.267325
  Sample3     NRT 5.689560 5.903811 8.029434  4.213996
  Sample4     NRP 5.834468 5.722589 8.061339  5.072274
  Sample5     RP 5.401623 7.272111 9.108342  6.599629
  
  ```
  
  purely in terms of numbers, the subtypes have following number of samples.
  
  ```
  NRT 52
  RT 42
  NRP 30
  RP 7
  
  ```
  
  Q1. May be the groups are imbalanced and caution should be taken?
  
  I have applied glmnet as following
  
  ```r
  modellinglasso=read.table("LRMatrix.txt", header=T, row.names = 1, sep="\t" )
  modellinglasso$Subtype=as.factor(modellinglasso$Subtype)
  modellinglasso$Subtype
  #build a single pass (single fold) lasso-penalised model
  lassoModel <- glmnet(
    x=data.matrix(modellinglasso[,-1]),
    y=modellinglasso$Subtype,
    standardize=TRUE,
    alpha=1.0,
    family="multinomial")
  plot(lassoModel, xvar="lambda")
  Warning message:
  In lognet(xd, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  :
    one multinomial or binomial class has fewer than 8  observations; dangerous ground
  
  ```
  
  [![enter image description here](https://i.sstatic.net/0lUP0.png)](https://i.sstatic.net/0lUP0.png)
  
  ```r
  #perform 10-fold cross validation
  cv.lassoModel <- cv.glmnet(
    x=data.matrix(modellinglasso[,-1]),
    y=modellinglasso$Subtype,
    standardize=TRUE,
    alpha=1.0,
    nfolds=10,
    family="multinomial",
    parallel=TRUE)
  # plot variable deviances vs. shrinkage parameter, λ (lambda)
  plot(cv.lassoModel)
  
  ```
  
  [![enter image description here](https://i.sstatic.net/E02wY.png)](https://i.sstatic.net/E02wY.png)
  
  ```r
  idealLambda <- cv.lassoModel$lambda.min
  idealLambda1se <- cv.lassoModel$lambda.1se
  print(idealLambda); print(idealLambda1se)
  
  [1] 0.1186125
  [1] 0.1186125
  
  # derive coefficients for each gene
  co <- coef(cv.lassoModel, s=idealLambda, exact=TRUE)
  head(co)
  
  co.se <- coef(cv.lassoModel, s=idealLambda1se, exact=TRUE)
  head(co.se)
  
  # identify predictors for each sub-type
  rownames(co.se$PR)[which(co.se$NRT != 0)]
  
  [1] "(Intercept)"
  
  ```
  
  output for
  
  ```
  co.se$RT!=0
  6 x 1 sparse Matrix of class "lgCMatrix"
              1
  (Intercept) |
  HES4        .
  FNDC10      .
  PRKCZ       .
  PRKCZ.AS1   .
  MMEL1       .
  
  ```
  
  It doesn't matter which subtype I put, I get the same output "(Intercept)" instead of gene names.
  
  Any help to get the gene signatures for the subtypes is sincerely appreciated.
  
  Thank you
  
  **Tags:** r, rna-seq, machine-learning
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2023-07-22. Score: 2 | Accepted answer
  
  Final answer. Keep in mind I work in Python so I'm trying to translate here.
  
  1. The output of `coef(cv.lassoModel)` ... those are your genes of interest thats your **feature selection**. `head` I assume gives you the top values ... I think thats right (I need to check its not `coef(lassoModel)`) ... Basically thats where your results are.
  
  This might be over simplistic but `co.se` without the `head` is your answer. Not every gene will carry a weight BTW. If the entirety of that is (700) `.` ... try point 2 below and then point 4. I don't think that will happen (see **Whats happened**)
  
  2. I would redo this with `alpha=0` (ridge regression) and `alpha=0.5` elasticnet, but you don't need to
  3. Generate an accuracy score
  
  ```
      y_predicted <- predict(best_model, s = best_lambda, newx = x)
      sst <- sum((y - mean(y))^2)
      sse <- sum((y_predicted - y)^2)
      rsq <- 1 - sse/sst
      rsq
  
  ```
  
  Hopefully this is large e.g. around 0.8 (in Python this automatic BTW). You could ignore that, but its what I'd do
  
  4. The intercept is simply the weight when your gene expression score is 0.
  5. Lasso regression handles sparse data - so no worries about `"."`. I don't whether that Python's version of "nan" (not a number), or if its zero. ... @username_2 knows what this means and his answer is authoritative for the second part of your question.
  6. In Python transformations are common, but that doesn't seem an option here, so find providing your accuracy score is reasonable, personally I'm happy.
  
  In Python what you want is feature selection
  `print(pd.DataFrame({'Features': X_train.columns, 'Coefficients': lasso.coef_}))`
  
  **Whats happened?** The coefficients for `head(co.se)` have been shrunk to zero, thats only the top 10? results, you've 700 weights. The full list *should* (famous last words) contain a minority of really big coefficients ... those are really important and those are your answer. So by printing only the top few results ... not the whole lot the results are still there you've just not seen them. Thats it, basically you were almost there and you could happily settle on that as the final answer. Basically, alot of genes ain't doing much a few genes are doing all the work - thats the answer. Below I stress repeatedly I don't think anything is wrong here and I personally like the result.
  
  ---
  
  Background
  
  **Lasso regression**, it's a cousin of ridge regression and key parameter is lambda (who'd guessed). What your plots are doing is tuning the parameter and finding the associated weight individually or collectively. Fairly straight forward. So like whats the value of lambda for your data set ('answered' below)? It's worth doing some reading on how lambda is calculated to understand the plot, if you're not sure how it works.
  
  **Imbalance**, well to be honest I'm not convinced there's anything wrong with your analysis. Okay there's a lack of stationarity for some of your weights ... but some are fine and gist is pretty clear. Yeah okay there's a problem with stabilisation of lambda - thats not necessarily imbalance.
  
  **Imbalance** is when it's 1 in 10 and over. You've got 1 in 6, 1 in 7, not huge. However, - *you might be right* - and if you're worried simply repeat the analysis excluding `RP` (value 7). Then check whether the parameterisation plot (lambda against weight) is more stable or more "homogeneous" (multinomial plot), i.e. the genes behaviour in a more uniform way. If that is does work - it gives you nicer results AND `RP` is essential to the analysis thats a separate question about augmentation. That is complicated. For example, the multinomial plot with a 3 target categorical model might give a clear value for lambda.
  
  **Transformation** Here, I would suggest if you want to stabilise the multinomial plot (if it's possible) I would look at a transformation. The data you presented doesn't look like it needs a transformation, but we can't see the whole data. In Python there's a *max-min* transformation - thats where I might start.
  
  **Overall good results** Personally, I'd say these were good results, except the weights ain't homogeneous for some genes (most are fine BTW), but the patterns clear. The *red* (top), *blue* (bottom) and cyan (negative y-axis) are real results showing that within the 700 genes **3 genes** are having a massive impact. Sure you have positive and negative and thats how it works one gene is positively contributing to the classification, others are negatively contributing to it. Then you've the black and 'deep pink' lines which have stabilised around 5E4. **From 700 genes down to 5 genes**. Whats wrong with that?
  
  **Whats the answer to the data set?**
  
  Its ...
  
  > NRT
  > RT
  > NRP
  > RP
  
  ... is due to 'red', 'blue', 'cyan', 'black' and 'deep pink' genes. Thats it.
  
  **Frank assessment** To be perfectly honest, looks good. No-one is going to say "Wooaaa lambda ain't stable, destroy! ;-)", they'll say "Wow thats cool and clear".
  
  ---
  
  **How I'd do it..** I would look at basic ML robustness statistics and use this to assess the model rather than the parameterisation of lambda. Its horses for courses as they say. Upvotes.
  
  **Note** I use Python not R for ML, both work well however. `glmnet` within ML is new to me BTW.
  
  > **Comments:**
  > **username_0**: Thank you for your reply and detailed explanation. What I do not understand is why I am getting output as Intercept instead of gene names.  
  > **username_0**: Another puzzling aspect is when I do co.se$RT!=0, I get a gene list with "." next to it :(  
  > **username_2**: Technically you are getting a `sparse Matrix` with `co.se$RT!=0` and the "dots" are representing zeros or in this specific case `FALSE` values. You can see yourself by `as.matrix(co.se$RT!=0)`.  
  > **username_0**: All of them :( I tried to normalise the values by getting zscores. HES4 is a gene name and its reading in different samples.  
  > **username_1**: @username_0 I've addressed the comments in the first part of the answer.  
  
  ---
  
  