Answers.Body: <p>confusion: Why when I am putting each of 23 genes individually in glm model, 5 of them returns p-value &lt; 0.05 but when summary of final model containing these 5 genes says only p-value of one of them is &lt; 0.05.</p>
  <p>Ans. Pr(&gt;|z|) simply suggests that whether z-value deviates beyond average z-value and if so,what is the extent of its probability. less than .05 indicates a 95% confidence in computed z value i.e. an acceptable significance level. The summary of 5 genes simply tells significance of over-all model.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-03-19T16:49:07.870Z
Answers.LastActivityDate: 2024-03-19T16:49:07.870Z
Answers.UserDisplayName: Subhash C. Davar
Body: <p>I want to find some good predictors (genes). I have checked the Spearman correlation of the expression of each of 23 genes with dependent variables (responders Vs non-responders and I saw only the 5 genes identified initially showed significant correlation with the dependent variable.</p>
  <p>This is my data, log transformed RNA-seq:</p>
  <pre><code>              TRG    CDK6 EGFR  KIF2C CDC20
      Sample 1  TRG12  11.39 10.62  9.75 10.34
      Sample 2  TRG12  10.16  8.63  8.68  9.08
      Sample 3  TRG12   9.29 10.24  9.89 10.11
      Sample 4  TRG45  11.53  9.22  9.35  9.13
      Sample 5  TRG45   8.35 10.62 10.25 10.01
      Sample 6  TRG45  11.71 10.43  8.87  9.44
      ...
  </code></pre>
  <p>23 differentially expressed genes between responder patients (TRG12) and non-responder patients (TRG45).</p>
  <p>1- I tested each of 23 genes individually in this code and each of them gives p-value &lt; 0.05 remained as a good predictor; For example for CDK6 I have done</p>
  <pre><code>&gt; glm=glm(TRG ~ CDK6, data = df, family = binomial(link = 'logit'))
  &gt; 
  &gt; summary(glm)
  
  Call:
  glm(formula = TRG ~ CDK6, family = binomial(link = &quot;logit&quot;), 
      data = df)
  
  Deviance Residuals: 
      Min       1Q   Median       3Q      Max  
  -1.8286  -0.9317   0.6382   0.9500   1.5966  
  
  Coefficients:
              Estimate Std. Error z value Pr(&gt;|z|)   
  (Intercept)  26.8950     9.6627   2.783  0.00538 **
  CDK6         -2.6932     0.9765  -2.758  0.00582 **
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 77.347  on 55  degrees of freedom
  Residual deviance: 64.062  on 54  degrees of freedom
  AIC: 68.062
  
  Number of Fisher Scoring iterations: 5
  </code></pre>
  <p>Finally I obtained five genes and I put them in this model:</p>
  <pre><code>final &lt;- glm(TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2 , data = df, family = binomial(link = 'logit'))
  </code></pre>
  <p>but</p>
  <pre><code>&gt; summary(final)
  
  Call:
  glm(formula = TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2, family = binomial(link = &quot;logit&quot;), 
      data = df)
  
  Deviance Residuals: 
      Min       1Q   Median       3Q      Max  
  -1.5696  -0.9090   0.2233   0.7125   2.0131  
  
  Coefficients:
              Estimate Std. Error z value Pr(&gt;|z|)  
  (Intercept)   3.2258    13.3877   0.241   0.8096  
  CDK6         -2.6101     1.0860  -2.403   0.0162 *
  CXCL8         0.8702     0.9783   0.889   0.3737  
  IL6           0.9591     0.8883   1.080   0.2803  
  ISG15         0.9574     0.5101   1.877   0.0605 .
  PTGS2        -0.4623     1.1005  -0.420   0.6744  
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 77.347  on 55  degrees of freedom
  Residual deviance: 53.917  on 50  degrees of freedom
  AIC: 65.917
  
  Number of Fisher Scoring iterations: 6
  </code></pre>
  <p>I am just seeing CDK6 as significant. I done cross validation</p>
  <pre><code>&gt; cv.glm(df, final, K=nrow(df))$delta
  [1] 0.2130050 0.2125748
  </code></pre>
  <p>Likely the delta values are not greatly differ (Model is right ???)</p>
  <p>Then I tried Sensitivity and Specificity</p>
  <pre><code>&gt; threshold=0.5
  &gt; predicted_values&lt;-ifelse(predict(final,type=&quot;response&quot;)&gt;threshold,1,0)
  &gt; actual_values&lt;-final$y
  &gt; conf_matrix&lt;-table(predicted_values,actual_values)
  &gt; conf_matrix
                  actual_values
  predicted_values  0  1
                 0 16  7
                 1 10 23
  &gt; sensitivity(conf_matrix)
  [1] 0.6153846
  &gt; specificity(conf_matrix)
  [1] 0.7666667
  </code></pre>
  <p>2- I went this way for step-wise method</p>
  <pre><code>full.model &lt;- glm(TRG ~., data = df, family = binomial)
  
  library(MASS)
  step.model &lt;- full.model %&gt;% stepAIC(trace = FALSE)
  coef(step.model)
  
  &gt; coef(step.model)
  (Intercept)        CDK6        CHGA        CSF3       CXCL8       ERBB2         IL6       ISG15       KRT14       KRT17 
  -23770.6261   -469.3660    376.3528   -238.5697   1404.5905    590.2586    789.5295    842.2058    423.0402   -651.3403 
       MAGEA1         OSM       PTGS2 
     496.6933   -637.0583   -500.1040 
  </code></pre>
  <p>In addition to 5 genes another genes remains as predictors but likely none of them are significant by Pr(&gt;|z|) &gt; 0.05</p>
  <pre><code>&gt; summary(step.model)
  
  Call:
  glm(formula = TRG ~ CDK6 + CHGA + CSF3 + CXCL8 + ERBB2 + IL6 + 
      ISG15 + KRT14 + KRT17 + MAGEA1 + OSM + PTGS2, family = binomial, 
      data = df)
  
  Deviance Residuals: 
         Min          1Q      Median          3Q         Max  
  -1.317e-04  -2.100e-08   2.100e-08   2.100e-08   1.396e-04  
  
  Coefficients:
               Estimate Std. Error z value Pr(&gt;|z|)
  (Intercept)  -23770.6  2922188.9  -0.008    0.994
  CDK6           -469.4    64920.6  -0.007    0.994
  CHGA            376.4    46191.2   0.008    0.993
  CSF3           -238.6    30857.1  -0.008    0.994
  CXCL8          1404.6   172558.7   0.008    0.994
  ERBB2           590.3    71989.5   0.008    0.993
  IL6             789.5   100130.7   0.008    0.994
  ISG15           842.2   102934.0   0.008    0.993
  KRT14           423.0    52565.5   0.008    0.994
  KRT17          -651.3    79489.6  -0.008    0.993
  MAGEA1          496.7    60319.6   0.008    0.993
  OSM            -637.1    82121.8  -0.008    0.994
  PTGS2          -500.1    73541.7  -0.007    0.995
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 7.7347e+01  on 55  degrees of freedom
  Residual deviance: 1.4741e-07  on 43  degrees of freedom
  AIC: 26
  
  Number of Fisher Scoring iterations: 25
  </code></pre>
  <p>3- I went to lasso regression</p>
  <pre><code>&gt; library(glmnet)
  &gt; lassoModel &lt;- glmnet(
  +   x=data.matrix(df[,-1]),
  +   y=df<span class="math-container">$TRG,
  +   standardize=TRUE,
  +   alpha=1.0,
  +   family="multinomial")
  &gt; #Derive coefficients for each gene to each subtype
  &gt; co &lt;- coef(lassoModel, s=idealLambda, exact=TRUE)
  &gt; co
  $</span>`TRG12`
  24 x 1 sparse Matrix of class &quot;dgCMatrix&quot;
                         1
  (Intercept)  4.441322326
  ALB         -0.010867331
  AQP9        -0.051908908
  CALML5       0.126359364
  CCL4        -0.059869284
  CDK6         0.283781437
  CHGA        -0.088283426
  CREB3L3      0.109988758
  CSF3        -0.012311176
  CXCL5       -0.030532485
  CXCL6       -0.063610112
  CXCL8       -0.084385732
  ERBB2       -0.102303095
  FGFRL1       0.148759702
  IL1B        -0.072225606
  IL6         -0.122671026
  ISG15       -0.297386587
  KLK5         0.039891341
  KRT14       -0.004536929
  KRT17        0.041693665
  MAGEA1      -0.154075458
  OSM          0.028987648
  PTGS2       -0.034960405
  S100A7A     -0.008190105
  
  $TRG45
  24 x 1 sparse Matrix of class &quot;dgCMatrix&quot;
                         1
  (Intercept) -4.441322326
  ALB          0.010867331
  AQP9         0.051908908
  CALML5      -0.126359364
  CCL4         0.059869284
  CDK6        -0.283781437
  CHGA         0.088283426
  CREB3L3     -0.109988758
  CSF3         0.012311176
  CXCL5        0.030532485
  CXCL6        0.063610112
  CXCL8        0.084385732
  ERBB2        0.102303095
  FGFRL1      -0.148759702
  IL1B         0.072225606
  IL6          0.122671026
  ISG15        0.297386587
  KLK5        -0.039891341
  KRT14        0.004536929
  KRT17       -0.041693665
  MAGEA1       0.154075458
  OSM         -0.028987648
  PTGS2        0.034960405
  S100A7A      0.008190105
  
  &gt;
  </code></pre>
  <p>Likely none of genes is significant so I just tried 5 genes from first glm</p>
  <pre><code>&gt; finalLasso &lt;- glm(df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5 , data=df, family=binomial(link=&quot;logit&quot;))
  &gt; summary(finalLasso)
  
  Call:
  glm(formula = df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5, family = binomial(link = &quot;logit&quot;), 
      data = df)
  
  Deviance Residuals: 
      Min       1Q   Median       3Q      Max  
  -1.6780  -0.6733   0.1643   0.6925   1.9182  
  
  Coefficients:
              Estimate Std. Error z value Pr(&gt;|z|)  
  (Intercept)  -0.4178    14.7456  -0.028   0.9774  
  CDK6         -2.2672     1.1537  -1.965   0.0494 *
  IL6           1.2524     0.9197   1.362   0.1733  
  ISG15         1.1763     0.5479   2.147   0.0318 *
  CXCL8         0.5413     0.6443   0.840   0.4008  
  CALML5       -0.7890     0.4219  -1.870   0.0615 .
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 77.347  on 55  degrees of freedom
  Residual deviance: 49.839  on 50  degrees of freedom
  AIC: 61.839
  
  Number of Fisher Scoring iterations: 6
  
  &gt; 
  </code></pre>
  <p>That 2 genes are significant; really confusing :(</p>
  <p>Please somebody save me from this horrible confusion, finally which genes are significant remaining in model and which not</p>
  <blockquote>
  <p>First confusion: Why when I am putting each of 23 genes individually
  in glm model, 5 of them returns p-value &lt; 0.05 but when summary of
  final model containing these 5 genes says only p-value of one of them
  is &lt; 0.05</p>
  <p>Second confusion: When I am using stepwise regression, summary of
  model says none of my genes is significant while some of them had been
  remained as predictors by <code>coef(step.model)</code></p>
  <p>Third confusion: CDK6 as my significant genes in two of model shows
  negative coefficient while I know that show be correlated with
  dependent variable positively</p>
  </blockquote>
  
Tags: <r><rna-seq><linear-regression>
Title: How do I resolve disagreements in the determination of significant genes?
ThreadText: # How do I resolve disagreements in the determination of significant genes?
  
  (Asked by: username_0 on 2019-01-15)
  
  I want to find some good predictors (genes). I have checked the Spearman correlation of the expression of each of 23 genes with dependent variables (responders Vs non-responders and I saw only the 5 genes identified initially showed significant correlation with the dependent variable.
  
  This is my data, log transformed RNA-seq:
  
  ```
                TRG    CDK6 EGFR  KIF2C CDC20
      Sample 1  TRG12  11.39 10.62  9.75 10.34
      Sample 2  TRG12  10.16  8.63  8.68  9.08
      Sample 3  TRG12   9.29 10.24  9.89 10.11
      Sample 4  TRG45  11.53  9.22  9.35  9.13
      Sample 5  TRG45   8.35 10.62 10.25 10.01
      Sample 6  TRG45  11.71 10.43  8.87  9.44
      ...
  
  ```
  
  23 differentially expressed genes between responder patients (TRG12) and non-responder patients (TRG45).
  
  1- I tested each of 23 genes individually in this code and each of them gives p-value < 0.05 remained as a good predictor; For example for CDK6 I have done
  
  ```
  > glm=glm(TRG ~ CDK6, data = df, family = binomial(link = 'logit'))
  > 
  > summary(glm)
  
  Call:
  glm(formula = TRG ~ CDK6, family = binomial(link = "logit"), 
      data = df)
  
  Deviance Residuals: 
      Min       1Q   Median       3Q      Max  
  -1.8286  -0.9317   0.6382   0.9500   1.5966  
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)   
  (Intercept)  26.8950     9.6627   2.783  0.00538 **
  CDK6         -2.6932     0.9765  -2.758  0.00582 **
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 77.347  on 55  degrees of freedom
  Residual deviance: 64.062  on 54  degrees of freedom
  AIC: 68.062
  
  Number of Fisher Scoring iterations: 5
  
  ```
  
  Finally I obtained five genes and I put them in this model:
  
  ```
  final <- glm(TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2 , data = df, family = binomial(link = 'logit'))
  
  ```
  
  but
  
  ```
  > summary(final)
  
  Call:
  glm(formula = TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2, family = binomial(link = "logit"), 
      data = df)
  
  Deviance Residuals: 
      Min       1Q   Median       3Q      Max  
  -1.5696  -0.9090   0.2233   0.7125   2.0131  
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
  (Intercept)   3.2258    13.3877   0.241   0.8096  
  CDK6         -2.6101     1.0860  -2.403   0.0162 *
  CXCL8         0.8702     0.9783   0.889   0.3737  
  IL6           0.9591     0.8883   1.080   0.2803  
  ISG15         0.9574     0.5101   1.877   0.0605 .
  PTGS2        -0.4623     1.1005  -0.420   0.6744  
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 77.347  on 55  degrees of freedom
  Residual deviance: 53.917  on 50  degrees of freedom
  AIC: 65.917
  
  Number of Fisher Scoring iterations: 6
  
  ```
  
  I am just seeing CDK6 as significant. I done cross validation
  
  ```
  > cv.glm(df, final, K=nrow(df))$delta
  [1] 0.2130050 0.2125748
  
  ```
  
  Likely the delta values are not greatly differ (Model is right ???)
  
  Then I tried Sensitivity and Specificity
  
  ```
  > threshold=0.5
  > predicted_values<-ifelse(predict(final,type="response")>threshold,1,0)
  > actual_values<-final$y
  > conf_matrix<-table(predicted_values,actual_values)
  > conf_matrix
                  actual_values
  predicted_values  0  1
                 0 16  7
                 1 10 23
  > sensitivity(conf_matrix)
  [1] 0.6153846
  > specificity(conf_matrix)
  [1] 0.7666667
  
  ```
  
  2- I went this way for step-wise method
  
  ```
  full.model <- glm(TRG ~., data = df, family = binomial)
  
  library(MASS)
  step.model <- full.model %>% stepAIC(trace = FALSE)
  coef(step.model)
  
  > coef(step.model)
  (Intercept)        CDK6        CHGA        CSF3       CXCL8       ERBB2         IL6       ISG15       KRT14       KRT17 
  -23770.6261   -469.3660    376.3528   -238.5697   1404.5905    590.2586    789.5295    842.2058    423.0402   -651.3403 
       MAGEA1         OSM       PTGS2 
     496.6933   -637.0583   -500.1040 
  
  ```
  
  In addition to 5 genes another genes remains as predictors but likely none of them are significant by Pr(>|z|) > 0.05
  
  ```
  > summary(step.model)
  
  Call:
  glm(formula = TRG ~ CDK6 + CHGA + CSF3 + CXCL8 + ERBB2 + IL6 + 
      ISG15 + KRT14 + KRT17 + MAGEA1 + OSM + PTGS2, family = binomial, 
      data = df)
  
  Deviance Residuals: 
         Min          1Q      Median          3Q         Max  
  -1.317e-04  -2.100e-08   2.100e-08   2.100e-08   1.396e-04  
  
  Coefficients:
               Estimate Std. Error z value Pr(>|z|)
  (Intercept)  -23770.6  2922188.9  -0.008    0.994
  CDK6           -469.4    64920.6  -0.007    0.994
  CHGA            376.4    46191.2   0.008    0.993
  CSF3           -238.6    30857.1  -0.008    0.994
  CXCL8          1404.6   172558.7   0.008    0.994
  ERBB2           590.3    71989.5   0.008    0.993
  IL6             789.5   100130.7   0.008    0.994
  ISG15           842.2   102934.0   0.008    0.993
  KRT14           423.0    52565.5   0.008    0.994
  KRT17          -651.3    79489.6  -0.008    0.993
  MAGEA1          496.7    60319.6   0.008    0.993
  OSM            -637.1    82121.8  -0.008    0.994
  PTGS2          -500.1    73541.7  -0.007    0.995
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 7.7347e+01  on 55  degrees of freedom
  Residual deviance: 1.4741e-07  on 43  degrees of freedom
  AIC: 26
  
  Number of Fisher Scoring iterations: 25
  
  ```
  
  3- I went to lasso regression
  
  ```
  > library(glmnet)
  > lassoModel <- glmnet(
  +   x=data.matrix(df[,-1]),
  +   y=df$TRG,
  +   standardize=TRUE,
  +   alpha=1.0,
  +   family="multinomial")
  > #Derive coefficients for each gene to each subtype
  > co <- coef(lassoModel, s=idealLambda, exact=TRUE)
  > co
  $`TRG12`
  24 x 1 sparse Matrix of class "dgCMatrix"
                         1
  (Intercept)  4.441322326
  ALB         -0.010867331
  AQP9        -0.051908908
  CALML5       0.126359364
  CCL4        -0.059869284
  CDK6         0.283781437
  CHGA        -0.088283426
  CREB3L3      0.109988758
  CSF3        -0.012311176
  CXCL5       -0.030532485
  CXCL6       -0.063610112
  CXCL8       -0.084385732
  ERBB2       -0.102303095
  FGFRL1       0.148759702
  IL1B        -0.072225606
  IL6         -0.122671026
  ISG15       -0.297386587
  KLK5         0.039891341
  KRT14       -0.004536929
  KRT17        0.041693665
  MAGEA1      -0.154075458
  OSM          0.028987648
  PTGS2       -0.034960405
  S100A7A     -0.008190105
  
  $TRG45
  24 x 1 sparse Matrix of class "dgCMatrix"
                         1
  (Intercept) -4.441322326
  ALB          0.010867331
  AQP9         0.051908908
  CALML5      -0.126359364
  CCL4         0.059869284
  CDK6        -0.283781437
  CHGA         0.088283426
  CREB3L3     -0.109988758
  CSF3         0.012311176
  CXCL5        0.030532485
  CXCL6        0.063610112
  CXCL8        0.084385732
  ERBB2        0.102303095
  FGFRL1      -0.148759702
  IL1B         0.072225606
  IL6          0.122671026
  ISG15        0.297386587
  KLK5        -0.039891341
  KRT14        0.004536929
  KRT17       -0.041693665
  MAGEA1       0.154075458
  OSM         -0.028987648
  PTGS2        0.034960405
  S100A7A      0.008190105
  
  >
  
  ```
  
  Likely none of genes is significant so I just tried 5 genes from first glm
  
  ```
  > finalLasso <- glm(df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5 , data=df, family=binomial(link="logit"))
  > summary(finalLasso)
  
  Call:
  glm(formula = df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5, family = binomial(link = "logit"), 
      data = df)
  
  Deviance Residuals: 
      Min       1Q   Median       3Q      Max  
  -1.6780  -0.6733   0.1643   0.6925   1.9182  
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
  (Intercept)  -0.4178    14.7456  -0.028   0.9774  
  CDK6         -2.2672     1.1537  -1.965   0.0494 *
  IL6           1.2524     0.9197   1.362   0.1733  
  ISG15         1.1763     0.5479   2.147   0.0318 *
  CXCL8         0.5413     0.6443   0.840   0.4008  
  CALML5       -0.7890     0.4219  -1.870   0.0615 .
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 77.347  on 55  degrees of freedom
  Residual deviance: 49.839  on 50  degrees of freedom
  AIC: 61.839
  
  Number of Fisher Scoring iterations: 6
  
  > 
  
  ```
  
  That 2 genes are significant; really confusing :(
  
  Please somebody save me from this horrible confusion, finally which genes are significant remaining in model and which not
  
  > First confusion: Why when I am putting each of 23 genes individually
  > in glm model, 5 of them returns p-value < 0.05 but when summary of
  > final model containing these 5 genes says only p-value of one of them
  > is < 0.05
  >
  > Second confusion: When I am using stepwise regression, summary of
  > model says none of my genes is significant while some of them had been
  > remained as predictors by `coef(step.model)`
  >
  > Third confusion: CDK6 as my significant genes in two of model shows
  > negative coefficient while I know that show be correlated with
  > dependent variable positively
  
  **Tags:** r, rna-seq, linear-regression
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2024-03-19. Score: 0
  
  confusion: Why when I am putting each of 23 genes individually in glm model, 5 of them returns p-value < 0.05 but when summary of final model containing these 5 genes says only p-value of one of them is < 0.05.
  
  Ans. Pr(>|z|) simply suggests that whether z-value deviates beyond average z-value and if so,what is the extent of its probability. less than .05 indicates a 95% confidence in computed z value i.e. an acceptable significance level. The summary of 5 genes simply tells significance of over-all model.
  
  
  ---
  
  