<p>By fitting the full model <code>~FAB + Sex + Age + TMB + WBC + BM_percentage</code> you can estimate the effect of FAB (or any of the other variables) on RNA expression after adjusting for differences caused by Sex, Age, TMB, WBC, and BM percent. See <a href="https://stats.stackexchange.com/questions/17336/how-exactly-does-one-control-for-other-variables">this post</a> for exactly what it means to adjust for other variables</p>
<p>When you fit a reduced model (in your case <code>~Sex + Age</code>) and then conduct a likelihood ratio test (LRT), you are comparing the fit of the full model (<code>~FAB + Sex + Age + TMB + WBC + BM_percentage</code>) to the fit of the reduced model (<code>~Sex + Age</code>). If the fit is substantially better, the likelihood of the full model will be larger than the reduced model and the p-value will be significant. This means that adding the extra variables that are in the full model (FAB, TMB, WBC, BM_percentage) improved the fit of the model to the data because together these variables have some type of relationship with RNA expression even after adjusting for the variables in the reduced model.</p>
<p>The challenge with using a LRT is that it doesn't tell you which of the variables added in the full model are responsible for the improved fit - it just indicates that the full model is better.  This type of test is useful when you are trying to make broad statements about several variables having a relationship with RNA expression. From the DESeq2 vignette:</p>
<blockquote>
<p>The LRT is therefore useful for testing multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. The LRT for count data is conceptually similar to an analysis of variance (ANOVA) calculation in linear regression, except that in the case of the Negative Binomial GLM, we use an analysis of deviance (ANODEV), where the deviance captures the difference in likelihood between a full and a reduced model.</p>
</blockquote>
<p>Definitely read the DESeq2 vignette for more details. In your case if the LR test is significant, it means that FAB, TMB, WBC, and BM percent combined are associated with RNA expression.</p>
<p>Most people are interested in specific comparisons when performing differential expression analyses with DESeq2. These comparisons are usually like do males have higher gene expression than females (comparing sex differences), do patients with higher TMB have greater expression than patients with lower TMB (comparing effect of TMB on RNA expression). If you want to do a specific comparison, you will likely want to use the Wald test, which calculates p-values for individual coefficients (male vs female coefficient, TMB coefficient etc). The Wald test determines whether a coefficient is equal to 0 or some other null hypothesis (see <code>lfcThreshold</code> parameter).</p>
<p>The Wald statistic and accompanying p-values are reported by default when using <code>results()</code> and not specifying <code>test=&quot;LRT&quot;</code></p>
<p>There are some additional differences between Wald tests vs LRT, including accuracy differences when operating on small sample sizes and computational efficiency issues when operating with large sample sizes. <a href="https://stats.stackexchange.com/questions/193643/likelihood-ratio-vs-wald-test">Discussion on Wald vs LRT</a></p>
<p><strong>TL;DR</strong> Use LRT to test if several variables or factor levels are associated with RNA expression. Use the Wald test to determine if specific comparisons are significant.</p>

CC BY-SA 4.0
2022-06-02T21:42:29.580Z
thank you for the detailed write up...this is more than enough for my data reporting which would be now be more statistical rather than using loose terms thanks a ton
kcm
CC BY-SA 4.0
2022-06-15T06:43:27.280Z
Little doubt "The design formula should have all of the factors in your metadata that account for major sources of variation in your data. The last factor entered in the formula should be the condition of interest." are these saying the same thing you described?  https://hbctraining.github.io/DGE_workshop_salmon_online/lessons/04a_design_formulas.html
kcm
CC BY-SA 4.0
2022-06-16T18:10:14.493Z
The last factor in the formula is the variable whose comparison (if sex is the last variable then male vs female) will show up by default when you call `results()`. You can change which comparison you get with the `contrast` argument.
Tomas Bencomo
CC BY-SA 4.0
2022-06-24T01:09:23.123Z
if I have a single factor in my design for example `age` and i ran deseq2 what would be my interpretation  would I be saying here Im testing the effect of `age` on the differences in gene expression?is that statistically correct statement?
kcm
CC BY-SA 4.0
2022-06-25T06:20:51.423Z
Yes that is correct. If the age variable is continuous, the age coefficient for a gene $g$ would represent the log2FoldChange of $g$ when age increases by 1 unit. If age has no effect on gene expression for $g$, the coefficient will be very small and not significant (assuming adequate power)
Tomas Bencomo
CC BY-SA 4.0
2022-06-28T13:00:05.927Z
in case of numerical such as age what becomes my intercept? for categorical Im now clear ,in case of age is it the lowest age becomes the intercept ?
kcm
CC BY-SA 4.0
2022-06-30T04:48:28.940Z
The intercept would be the value when age is 0. Based on your questions it would probably be helpful to learn more about linear models as they are very common in bioinformatics. DESeq2 uses a negative binomial linear model. I would start with linear regression and then checkout [this book](https://bookdown.org/roback/bookdown-BeyondMLR/ch-MLRreview.html#introduction-to-beyond-multiple-linear-regression)
Tomas Bencomo
CC BY-SA 4.0
2022-06-30T10:43:53.390Z
'it would probably be helpful to learn more about linear models as they are very common in bioinformatics' yes i been vigorously digging and reading trying to wrap my head on the concept and the biology going on
kcm
CC BY-SA 4.0
2022-06-02T17:54:20.540Z
2022-06-02T17:54:20.540Z
Tomas Bencomo
<p>I would like to know how to interpret the output of the formula although i been thorough loads of literature but I'm not confident yet.</p>
<p>So this is my full model I'm running</p>
<pre><code>dds &lt;- DESeqDataSetFromMatrix(countData=rpkm_ordered, colData=coldata, design= ~ FAB + Sex + Age + TMB + WBC + BM_percentage)
</code></pre>
<p>So if I have to interpret the output of this how it would be ?</p>
<p>What I have read is this</p>
<p>Here in the above formula we are trying or controlling the <code>FAB</code> variable while we look for variation that is coming from other variables?</p>
<p>Now my confusion is when I made this <code>dds</code> object the <code>FAB</code> was the variable which I fixed while I look for variation coming from other variable.</p>
<p>Next thing is I would like to run a reduced model.</p>
<p>I tried this</p>
<pre><code>test.Age &lt;- as_tibble(results(DESeq(dds, test=&quot;LRT&quot;, reduced= ~Sex + Age )))
</code></pre>
<p>Now my question is since I'm using the same <code>dds</code> object here where the <code>FAB</code> is the variable i fixed and  with new formula which is</p>
<p><code>reduced= ~Sex + Age</code> this so how would i interpret this model?</p>
<p>Any suggestion or help would be really appreciated</p>

CC BY-SA 4.0
2022-05-30T17:03:01.500Z
You know that RPKM is inappropriate for DESeq2?
swbarnes2
CC BY-SA 4.0
2022-05-30T18:15:30.847Z
yes..are you `rpkm_ordered` telling about this i suppose its just dataframe i named in my code which i didn't change..
kcm
CC BY-SA 4.0
2022-05-30T18:16:44.137Z
@swbarnes2 sorry for this  `rpkm_ordered`  confusion since once i had to reorder my dataframe...which I did quite long time back so its just in my code..
kcm
<r><deseq2>
deseq2 full and reduced model interpretation
# deseq2 full and reduced model interpretation

(Asked by: username_0 on 2022-05-30)

I would like to know how to interpret the output of the formula although i been thorough loads of literature but I'm not confident yet.
So this is my full model I'm running
dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered, colData=coldata, design= ~ FAB + Sex + Age + TMB + WBC + BM_percentage)

So if I have to interpret the output of this how it would be ?
What I have read is this
Here in the above formula we are trying or controlling the FAB variable while we look for variation that is coming from other variables?
Now my confusion is when I made this dds object the FAB was the variable which I fixed while I look for variation coming from other variable.
Next thing is I would like to run a reduced model.
I tried this
test.Age <- as_tibble(results(DESeq(dds, test="LRT", reduced= ~Sex + Age )))

Now my question is since I'm using the same dds object here where the FAB is the variable i fixed and  with new formula which is
reduced= ~Sex + Age this so how would i interpret this model?
Any suggestion or help would be really appreciated

**Tags:** r, deseq2
> **Comments:**
> **username_1**: You know that RPKM is inappropriate for DESeq2?  
> **username_0**: yes..are you `rpkm_ordered` telling about this i suppose its just dataframe i named in my code which i didn't change..  
> **username_0**: @swbarnes2 sorry for this  `rpkm_ordered`  confusion since once i had to reorder my dataframe...which I did quite long time back so its just in my code..  

---

## 1 answers

### Answer by username_2 on 2022-06-02. Score: 2 | Accepted answer

By fitting the full model `~FAB + Sex + Age + TMB + WBC + BM_percentage` you can estimate the effect of FAB (or any of the other variables) on RNA expression after adjusting for differences caused by Sex, Age, TMB, WBC, and BM percent. See [this post](https://stats.stackexchange.com/questions/17336/how-exactly-does-one-control-for-other-variables) for exactly what it means to adjust for other variables

When you fit a reduced model (in your case `~Sex + Age`) and then conduct a likelihood ratio test (LRT), you are comparing the fit of the full model (`~FAB + Sex + Age + TMB + WBC + BM_percentage`) to the fit of the reduced model (`~Sex + Age`). If the fit is substantially better, the likelihood of the full model will be larger than the reduced model and the p-value will be significant. This means that adding the extra variables that are in the full model (FAB, TMB, WBC, BM\_percentage) improved the fit of the model to the data because together these variables have some type of relationship with RNA expression even after adjusting for the variables in the reduced model.

The challenge with using a LRT is that it doesn't tell you which of the variables added in the full model are responsible for the improved fit - it just indicates that the full model is better. This type of test is useful when you are trying to make broad statements about several variables having a relationship with RNA expression. From the DESeq2 vignette:

> The LRT is therefore useful for testing multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. The LRT for count data is conceptually similar to an analysis of variance (ANOVA) calculation in linear regression, except that in the case of the Negative Binomial GLM, we use an analysis of deviance (ANODEV), where the deviance captures the difference in likelihood between a full and a reduced model.

Definitely read the DESeq2 vignette for more details. In your case if the LR test is significant, it means that FAB, TMB, WBC, and BM percent combined are associated with RNA expression.

Most people are interested in specific comparisons when performing differential expression analyses with DESeq2. These comparisons are usually like do males have higher gene expression than females (comparing sex differences), do patients with higher TMB have greater expression than patients with lower TMB (comparing effect of TMB on RNA expression). If you want to do a specific comparison, you will likely want to use the Wald test, which calculates p-values for individual coefficients (male vs female coefficient, TMB coefficient etc). The Wald test determines whether a coefficient is equal to 0 or some other null hypothesis (see `lfcThreshold` parameter).

The Wald statistic and accompanying p-values are reported by default when using `results()` and not specifying `test="LRT"`

There are some additional differences between Wald tests vs LRT, including accuracy differences when operating on small sample sizes and computational efficiency issues when operating with large sample sizes. [Discussion on Wald vs LRT](https://stats.stackexchange.com/questions/193643/likelihood-ratio-vs-wald-test)

**TL;DR** Use LRT to test if several variables or factor levels are associated with RNA expression. Use the Wald test to determine if specific comparisons are significant.

> **Comments:**
> **username_0**: thank you for the detailed write up...this is more than enough for my data reporting which would be now be more statistical rather than using loose terms thanks a ton  
> **username_0**: Little doubt "The design formula should have all of the factors in your metadata that account for major sources of variation in your data. The last factor entered in the formula should be the condition of interest." are these saying the same thing you described?  https://hbctraining.github.io/DGE_workshop_salmon_online/lessons/04a_design_formulas.html  
> **username_2**: The last factor in the formula is the variable whose comparison (if sex is the last variable then male vs female) will show up by default when you call `results()`. You can change which comparison you get with the `contrast` argument.  
> **username_0**: if I have a single factor in my design for example `age` and i ran deseq2 what would be my interpretation  would I be saying here Im testing the effect of `age` on the differences in gene expression?is that statistically correct statement?  
> **username_2**: Yes that is correct. If the age variable is continuous, the age coefficient for a gene $g$ would represent the log2FoldChange of $g$ when age increases by 1 unit. If age has no effect on gene expression for $g$, the coefficient will be very small and not significant (assuming adequate power)  
> **username_0**: in case of numerical such as age what becomes my intercept? for categorical Im now clear ,in case of age is it the lowest age becomes the intercept ?  
> **username_2**: The intercept would be the value when age is 0. Based on your questions it would probably be helpful to learn more about linear models as they are very common in bioinformatics. DESeq2 uses a negative binomial linear model. I would start with linear regression and then checkout [this book](https://bookdown.org/roback/bookdown-BeyondMLR/ch-MLRreview.html#introduction-to-beyond-multiple-linear-regression)  
> **username_0**: 'it would probably be helpful to learn more about linear models as they are very common in bioinformatics' yes i been vigorously digging and reading trying to wrap my head on the concept and the biology going on  

---

