# Songbird on paired samples
**URL:** https://forum.qiime2.org/t/songbird-on-paired-samples/14379
**Category:** Library Support
**Tags:** best-of-the-forum, q2-songbird
**Created:** [April 6, 2020, 8:03pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379 "2020-04-06T20:03:37Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@jwdebelius](https://forum.qiime2.org/u/jwdebelius)
#### Post date: [April 6, 2020, 8:03pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/1 "2020-04-06T20:03:37Z")
Hi all, especially @mortonjt ,
I'm trying to figure out Songbird on paired samples. Ive been through the tutorial and looked at the differential code from the original paper. Im slightly confused how to use this for paired samples, since as far as I can tell by the formula, the code is looking at a bulk comparison between before and after tooth brushing without conditioning on the individual. (The code is `C(brushing_event)` and I might have expected something closer to the `(C(brushing_event) | host_subject_id)` in an LME... except that I m not totally sure that works in this model.) I might also be missing something with pre-processing where the pairwise comparison is already set up.
Alternatively, if there is another kosher way to do paired sample tests for relative abundance deltas, that would be super useful. I found one paper, but it was a UniFrac-like approach and I'd like differential abundance for individual features or clades.
Thanks,
Justine
---
### Author:  [@mortonjt](https://forum.qiime2.org/u/mortonjt)
#### Post date: [April 6, 2020, 11:23pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/2 "2020-04-06T23:23:48Z")
Hey @jwdebelius, very good question. One important thing to know about Songbird is that it doesn't actually perform null hypothesis testing. The main point was to show that ranks are a useful concept to embrace. The reason why we showed that paired testing works is because the mean calculation in paired t-test is identical to a standard t-test - the only difference here is the hypothesis test, which we tackled using [paired\_ttest](https://github.com/knightlab-analyses/reference-frames/blob/master/ipynb/ttest.py) method (which was adapted from scipy) after identifying taxa and performing the appropriate log-ratio transformation.
For your case, there are 3 possibilities:
1. You could do this approach of first identifying microbes, followed by either a paired t-test or [LME](https://www.statsmodels.org/stable/mixed_linear.html).
2. You could try to adapt this [code](https://mortonjt.github.io/probable-bug-bytes/probable-bug-bytes/differential-abundance/) to design your own statistical test. Note that this can be tricky (I don't have the answer to this at the moment).
3. There are other tools that better handle longitudinal data such as [MIMIX](https://github.com/nsgrantham/mimix) and [stray](https://github.com/jsilve24/stray), but you may run into similar problems as discussed in (2).
---
### Author:  [@jwdebelius](https://forum.qiime2.org/u/jwdebelius)
#### Post date: [April 6, 2020, 11:50pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/3 "2020-04-06T23:50:19Z")
Thank you! I will check those out and maybe be back with questions.
Best,
Justine
---
### Author:  [@vrbana](https://forum.qiime2.org/u/vrbana)
#### Post date: [November 15, 2020, 9:17pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/4 "2020-11-15T21:17:58Z")
@jwdebelius just curious, what did you end up trying for this? Anything work well? I had the exact same question and was glad to see you had already asked it 😆
I've been messing around with edgeR and other RNASeq-related programs which handle paired samples well, but I'm not sure I fully trust those for microbiome data and would like to use ranks/ratios. Also can't really use those for non-count based metabolome data.
---
### Author:  [@jwdebelius](https://forum.qiime2.org/u/jwdebelius)
#### Post date: [November 16, 2020, 8:23am UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/5 "2020-11-16T08:23:46Z")
Hi @vrbana,
I hope you're well! @mortonjt helped me out with some code and ideas! He shared a baysian approach. Maybe ask him directly?
Best,
Justine
---
### Author:  [@mortonjt](https://forum.qiime2.org/u/mortonjt)
#### Post date: [November 16, 2020, 8:35pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/6 "2020-11-16T20:35:32Z")
Hi @vrbana, I gave @jwdebelius some Stan scripts to build a linear mixed effects negative binomial model : [Show how to fit a random effects models with negative binomial · GitHub](https://gist.github.com/mortonjt/9d336566c83541e9c80f0f78c3a97750)
And see my blog post here to see how to build your own model : [How to build your own differential abundance tool in 15 MINUTES! – Jamie Morton – Genomics, ecology and statistics](https://mortonjt.github.io/probable-bug-bytes/probable-bug-bytes/differential-abundance/)
So this will allow you to build something like DESeq2, but will give you a little more flexibility in statistical testing (DESeq2 only tests against the mean, whereas this will allow you to do ratio tests).
---
### Author:  [@vrbana](https://forum.qiime2.org/u/vrbana)
#### Post date: [November 16, 2020, 11:29pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/7 "2020-11-16T23:29:52Z")
Thanks Jamie! What model do you typically use for continuous data (i.e. metabolome or metagenome rel abund)? It's my understanding that the negative binomial is only appropriate for count data.
---
### Author:  [@mortonjt](https://forum.qiime2.org/u/mortonjt)
#### Post date: [November 16, 2020, 11:55pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/8 "2020-11-16T23:55:30Z")
It's a good question, one that no one has a great grasp on.
Metabolomics data is generated from discrete abundances -- molecular abundances are inherently discrete quantities. It's the fragmentation that converts these discrete quantities into continuous qualities. Because the underlying quantities of interest is inherently discrete, I typically stick with discrete models. Is it a hack? Most definitely. But there are very few tools that can deal with discrete latent variables in this fashion ( I don't believe it is currently tractable in this case). So your options are to either model everything in a discrete world (with some hacky approximation) or model everything continuously.
But you can certainly experiment with [different distributions](https://mc-stan.org/docs/2_18/functions-reference/positive-continuous-distributions.html) in Stan -- Stan gives you quite a bit of freedom in this regard. You could even swap out the Negative
Binomial with a Gaussian distribution or even better a [Truncated Gaussian distribution](https://mc-stan.org/docs/2_18/stan-users-guide/truncated-data-section.html) to model the abundances. As long as the confounding sequencing depth is accounted for (and your predicted abundances are strictly positive), these modeling choices should be ok.
---
### Author:  [@vrbana](https://forum.qiime2.org/u/vrbana)
#### Post date: [November 17, 2020, 12:40am UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/9 "2020-11-17T00:40:15Z")
Interesting, even though the abundance is MS1 peak area? Seems like people treat it both ways in other papers; I wish we had discussed/argued about this in a joint lab meeting with the metabolomics folks! But thanks, I will try out some different distributions
---
### Author:  [@mortonjt](https://forum.qiime2.org/u/mortonjt)
#### Post date: [November 17, 2020, 2:22am UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/10 "2020-11-17T02:22:49Z")
If you have access to replicates / dilution experiments - that can also help in deciding what an appropriate error distribution would be.
But I think you get an idea about the conceptual challenges - definitely curious to hear what you find out!
---
### Author:  [@vrbana](https://forum.qiime2.org/u/vrbana)
#### Post date: [November 19, 2020, 11:42pm UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/12 "2020-11-19T23:42:32Z")
Sorry I can't for the life of me figure this out and can't find any good examples, how do you incorporate subject into your DE abundance tool example? For example, in an R-style formula I would just add (1 | subject) and for statsmodels it is just groups=dataframe['subject'].
---
### Author:  [@mortonjt](https://forum.qiime2.org/u/mortonjt)
#### Post date: [November 20, 2020, 12:04am UTC](https://forum.qiime2.org/t/songbird-on-paired-samples/14379/13 "2020-11-20T00:04:57Z")
In my example I have subject ids directly coded : [Show how to fit a random effects models with negative binomial · GitHub](https://gist.github.com/mortonjt/9d336566c83541e9c80f0f78c3a97750#file-negative-binomial-random-effects-stan-L9)
This would require a small update to the input to the stan model highlighted in the above tutorial via
```python
dat = {
'N' : X.shape[0],
'D' : table.shape[0],
'p' : X.shape[1],
'depth' : np.log(table.sum(axis='sample').values),
'x' : X.values,
'y' : table.values.astype(np.int64)
'subj_ids' :
}
```
where `subj_ids` is some numeric conversion of your actual host subject ids.
There are also quite a bit of documentation on stan that I recommend checking out : [Stan User’s Guide](https://mc-stan.org/docs/2_25/stan-users-guide/index.html).
Also stan has its own channel ; for instance this page highlights an example fitting mixed effects models : [Linear mixed effects model with varying intercept and slope in matrix notation - #3 by Wade - Modeling - The Stan Forums](https://discourse.mc-stan.org/t/linear-mixed-effects-model-with-varying-intercept-and-slope-in-matrix-notation/7308/3)