Intervention Analysis with Mirrors

`$\def\Dir{\text{Dir}}$`
`$\renewcommand{\exp}[1]{\operatorname{exp}\left(#1\right)}$`
`$\def\Mult{\text{Mult}}$`
`$\def\*#1{\mathbf{#1}}$`
`$\def\m#1{\boldsymbol{#1}}$`
`$\def\Unif{\text{Unif}}$`
`$\def\win{\tilde{w}_{\text{in}}}$`
`$\def\reals{\mathbb{R}}$`
`$\def\simplex{\Delta}$`
`$\def\*#1{\mathbf{#1}}$`
`$\def\m#1{\boldsymbol{#1}}$`
`$\def\PD{\mathrm{PD}}$`
`$\def\FDP{\mathrm{FDR}}$`
`$\newcommand\mbb[1]{\mathbb{#1}}$`
`$\newcommand\mbf[1]{\mathbf{#1}}$`
`$\def\mc#1{\mathcal{#1}}$`
`$\def\mrm#1{\mathrm{#1}}$`
`$\def\absarg#1{\left|#1\right|}$`
<div style = "position:fixed; visibility: hidden">
`$$\require{color}\definecolor{myred}{rgb}{0.705882352941177, 0.341176470588235, 0.36078431372549}$$`
`$$\require{color}\definecolor{mygreen}{rgb}{0.352941176470588, 0.541176470588235, 0.501960784313725}$$`
</div>

## Microbiome Intervention Analysis with `mbtransfer`

<div id="subtitle">
Kris Sankaran <br/>
(w/ Pratheepa Jeganathan)
06 | August | 2024 <br/>
Lab: <a href="https://go.wisc.edu/pgb8nl">go.wisc.edu/pgb8nl</a> <br/>
</div>

<div id="subtitle_right">
Advanced Methods for Microbiome Data Analysis<br/>
Joint Statistical Meetings</br>
Slides: <a href="https://go.wisc.edu/h5ow0o">go.wisc.edu/h5ow0o</a><br/>
Paper: <a href=" https://go.wisc.edu/5x8k89">go.wisc.edu/5x8k89</a>
</div>

---

### Motivation: Interaction in Dynamic Communities

Microbiome communities are dynamic [1; 2]. More precise knowledge of how it changes 
can guide microbiome design [3; 4].

<img src="figure/curtis-opening.png"/>
<span style="font-size: 18px;">
John Curtis was a UW - Madison Botany
professor who studied how forest ecosystems change in response to environmental
pressures (and created the Bray-Curtis distance). Figure from [5].
</span>

---

### Interventions

To disentangle these dynamics, longitudinal experimental designs with natural or
induced interventions have gained popularity. The questions are:

1. **Who**: Which taxa are affected by an intervention?
1. **When**: Are the effects immediate? Lagged? Do they persist?
1. **Why**: Are there host factors that mediate the effect?

<img src="figure/motivation.png"/>
<span style="font-size: 18px;">
Examples of microbiomes responding to environmental shifts, discussed in our supplement S2 [6]. 
</span>

---

---

### Ingredients

1. **Transfer Functions** [7]: Model
intervention and interaction effects across time horizons using a
high-dimensional linear model.

1. **Mirror Statistics** [8]: Use agreement across data
splits to allow false discovery rate-controlled inference of intervention
effects.

---

### Transfer Functions

`$$y_{t} = \sum_{p = 1}^{P} A_{p} y_{t - p} + \sum_{q = 0}^{Q - 1} B_{q}w_{t - q} + \epsilon_{t}$$`

where `$y_{t}$` is the series of interest and `$w_{t}$` encodes the intervention.

Unlike gLV models [22; 23], they can model lagged intervention effects.
]

.pull-right[
<img src="figure/pulse.png" width="340" style="display: block; margin: auto;" />
<span style="font-size: 18px;">
In Box and Tiao [7], `$y_{t}$` are ozone concentrations and `$w_{t}$` are new regulations.
 </span>

]

---

### Transfer Functions for Microbial Communities

We estimate a separate transfer function model `$f_{j}$` for each taxon:

`$$\mathbf{y}_{t}^{(i)} = \mathbf{f}\left(\mathbf{Y}^{(i)}_{t - 1}, \mathbf{W}^{(i)}_{t}, \mathbf{z}^{(i)} \right) + \mathbf{\epsilon}_{t}^{(i)}$$`

For each coordinate, we use a regularized linear boosting model [9] with candidate interactions found in an initial prescreen [10].

|Variable| Interpretation |
|--|---|
| `$\mathbf{y}_{t}^{(i)} \in \reals^{J}$` | The microbiome community profiles in subject `$i$` at time `$t$`. |
| `$\mathbf{Y}^{(i)}_{t - 1} \in \reals^{J \times P}$` | The length `$P$` memory of the `$i^{th}$` subject's community. |
| `$\mathbf{W}^{(i)}_{t - 1} \in \reals^{D \times Q}$` | The length `$Q$` memory of the `$i^{th}$` subject's interventions. |
| `$\mathbf{z}^{(i)} \in \reals^{S}$` | Subject-level covariates that can interact with `$\mathbf{Y}_{t - 1}^{(i)}$` and `$\mathbf{W}_{t}^{(i)}$`.|

---

### Mirror Statistics

1. Which taxa `$j$` are affected by the external perturbation, and how long do the
perturbations last?

2. We can use the approach of Dai et al. [8] to gauge statistical significance of effects that are not available in model
coefficients.
.center[
<img src="figure/trajectories_significance.png" width=700/>
]

---

### General Approach: Mirror Statistics

Split the data into `$\mathcal{D}^{(1)} = \left(\mathbf{X}^{(1)}, \mathbf{y}^{(1)}\right)$` and `$\mathcal{D}^{(2)} = \left(\mathbf{X}^{(2)}, \mathbf{y}^{(2)}\right)$` and check for agreement in the estimates `$\hat{\mathbf{\beta}}^{(1)}$` and `$\hat{\mathbf{\beta}}^{(2)}$` across splits.

---

### Test Statistic for Linear Model

To detect whether `$\beta_{j} \neq 0$`, we can use the statistic:

.center[
`$M_{j} = \text{sign}\left(\hat{\beta}^{(1)}_{j}\hat{\beta}^{(2)}_{j}\right)\left[\left|\hat{\beta}_{j}^{(1)}\right| + \left|\hat{\beta}_{j}^{(2)}\right|\right]$`
]

---

### FDR Estimation

To choose a threshold `$\tau$`, estimate the false discovery rate:

.center[
`$$\widehat{\FDP}\left(\tau\right) = \frac{\left|\{j : M_{j} < -\tau\}\right|}{\left|\left\{j : M_{j} > \tau\right\}\right|}$$`
]

---

### FDR Estimation

To choose a threshold `$\tau$`, estimate the false discovery rate:

.center[
`$$\widehat{\FDP}\left(\tau\right) = \frac{\left|\{j : M_{j} < -\tau\}\right|}{\left|\left\{j : M_{j} > \tau\right\}\right|}$$`
]

.pull-left[
<img src="figure/mirror-explanation2.png" width=500/>
]
.pull-right[
Assumptions:
* Symmetry of `$M_{j}$` under the null
* Independence of `$M_{j}$`
]

---

### Partial Dependence

For transfer function models, we need an alternative to `$\beta_{j}$`. We estimate the partial dependence [9; 11] of the intervention `$\mathbf{W}_{t}$` across time horizons.

.center[
`$$\PD_{j}^{(s)} = \frac{1}{\absarg{\mathcal{D}^{(s)}}}\sum_{d_{t}^{(i)} \in \mathcal{D}^{(s)}}\left[\myred{\hat{f}_{j}^{(s)}\left(\*Y^{(i)}_{t}, \*1_{Q}, \*z^{(i)}\right)} - \mygreen{\hat{f}_{j}^{(s)}\left(\*Y^{(i)}_{t}, \*0_{Q}, \*z^{(i)}\right)}\right].$$`
<img src="figure/pd-1.png" width=650/>
]

---

### Partial Dependence

For transfer function models, we need an alternative to `$\beta_{j}$`. We estimate the partial dependence [9; 11] of the intervention `$\mathbf{W}_{t}$` across time horizons.

---

### Partial Dependence

For transfer function models, we need an alternative to `$\beta_{j}$`. We estimate the partial dependence [9; 11] of the intervention `$\mathbf{W}_{t}$` across time horizons.

---

### Partial Dependence Mirrors

Our definition uses `$\mathrm{PD}_{j}^{(s)}$` in place of `$\hat{\beta}_{j}^{(s)}$` from the
usual definition.

`\begin{align*}
	M_{j} = \text{sign}\left(\PD_{j}^{(1)}\PD_{j}^{(2)}\right)\left[\absarg{\PD_{j}^{(1)}} + \absarg{\PD_{j}^{(2)}}\right]
\end{align*}`

---

---

### Evaluation Criteria

1. **Forecasting**: How do methods compare in their ability to predict microbial
community response to interventions?

1. **Selection**: How do methods compare in their FDR control and power in
detecting affected taxa?

---

### Simulation Parameters

<img src="figure/dimensionality.png" width=50/> Number of (nonnull) taxa: How do dimensionality and proportion of nulls influence performance?

<img src="figure/normal_curve.png" width=80/> Signal strength: How are inference and forecasting are affected by intervention strength?

<img src="figure/transform.png" width=40/> Normalization strategy: We can apply methods to either raw or transformed data.

<img src="figure/phylo.png" width=50/> Phylogenetic correlation: We include correlated taxa that reflect shared evolutionary ancestry.<br/>

<img src="figure/reads.png" width=90/> Sequencing depth: In a typical experiment, not all samples are sequenced to the same read depth.

---

### Generative Mechanism

Negative binomial models often fit microbiome data well [12; 13; 14]. We use an
autoregressive version:

`$$\mathbf{y}_{t}^{(i)}\vert\theta_{t}^{(i)},\mathbf{\varphi},d^{(i)}\sim\text{NB}\left(d^{(i)}\exp{\theta_{t}^{(i)}},\mathbf{\varphi}\right)\\\theta_{t}^{(i)}=\sum_{p=1}^{P}A_{p}\theta_{t-p}^{(i)}+\sum_{q=1}^{Q}\left(B_{q}+C_{q}\odot \mathbf{z}^{(i)}\right)\mathbf{w}_{t-q}^{i}+\mathbf{\epsilon}_{t}^{(i)}\\d^{(i)}\sim\Gamma\left(10,\lambda\right)\\\epsilon_{t}^{(i)}\sim\mathcal{N}\left(0,\Sigma\left(\alpha\right)\right)$$`

---

### Comparison

Forecasting

1. `MDSINE2` [15; 16; 17]: A Bayesian nonparametric negative binomial autoregressive model for discovering taxonomic interactions.
1. `fido` [18; 19]: An extended multinomial regression model that we adapted to forecasting.

Selection

1. `DESeq2` [20] with main and interaction effects for `$\mathbf{w}_{t}$` and `$\mathbf{z}$`.
1. Two-Sample `$t$`-test comparing the four `$\mathbf{y}_{t}$` for the first four `$t$` after `$\mathbf{w}_{t} = 1$` and with all other `$\mathbf{y}_{t}$`

---

### Prediction Performance: `$\alpha = 0.1, \lambda = 0.1$`

---

### Prediction Performance: `$\alpha = 10, \lambda = 0.1$`

---

### All Settings

---

### Power and FDR

---

---

### Data Description

David et al. [21] investigated short-term change in
microbiome composition and function in response to diet perturbations.

.pull-three-quarters-left[
* 20 participants observed for 15 days.
* 15 days = 4 pre + 5 intervention + 6 post.
* Specially prepared animal- (10 subjects) and plant-based (10 subjects) diets.
]

---

### In- and Out-of-Sample Prediction

<img src="figure/Fig7.png" width=750/><br/>
<span style="font-size: 18px;">
Forecasting performance across subsets of taxa and time lags. Out-of-sample
performance is strongest among highly abundant taxa (top row) and short-term
time lags (left column).
</span>

---

### Mirror Statistics

<span style="font-size: 18px;">
Boxplots of `$M_{j}$` that are larger than zero give evidence for significant intervention effects.
</span>

---

### Estimates and Data

For each taxon, we can compare the partial dependence profiles when the diet
intervention is or is not applied.

<img src="figure/Fig8.png" width=900/>
<span style="font-size: 18px;"><br/>
Estimates `$\mathrm{PD}_{j}$` (top panel) and underlying data (bottom panel) for taxa significant effects.
</span>

---

### Takeaways

`mbtransfer` supports FDR-controlled inference across multiple time horizons for
models of dynamic microbiome data.

.pull-three-quarters-right[
[Paper](https://go.wisc.edu/5x8k89)<br/>
<img src="figure/qr-paper.png" width=80/>
<br/>
<br/>
<br/>
[Package](https://krisrs1128.github.io/mbtransfer)<br/>
<img src="figure/qr-package.png" width=80/>
]

---

* Collaborators: Pratheepa Jeganathan (McMaster University)
* Lab Members: Margaret Thairu, Hanying Jiang, Shuchen Yan, Yuliang Peng, Kai Cui, Sam Merten, and Kobe Uko
* Funding: NIGMS R01GM152744

---

### References

[1] L. Dethlefsen et al. "Incomplete recovery and individualized
responses of the human distal gut microbiota to repeated antibiotic
perturbation". In: _Proceedings of the National Academy of Sciences_
108.supplement_1 (Sep. 2010), p. 4554–4561. ISSN: 1091-6490. DOI:
[10.1073/pnas.1000087107](https://doi.org/10.1073%2Fpnas.1000087107).
URL:
[http://dx.doi.org/10.1073/pnas.1000087107](http://dx.doi.org/10.1073/pnas.1000087107).

[2] G. K. Gerber. "The dynamic microbiome". In: _FEBS Letters_ 588.22
(Feb. 2014), p. 4131–4139. ISSN: 1873-3468. DOI:
[10.1016/j.febslet.2014.02.037](https://doi.org/10.1016%2Fj.febslet.2014.02.037).
URL:
[http://dx.doi.org/10.1016/j.febslet.2014.02.037](http://dx.doi.org/10.1016/j.febslet.2014.02.037).

[3] In: _Nature Microbiology_ 9.7 (Jul. 2024), p. 1633–1633. ISSN:
2058-5276. DOI:
[10.1038/s41564-024-01763-1](https://doi.org/10.1038%2Fs41564-024-01763-1).
URL:
[http://dx.doi.org/10.1038/s41564-024-01763-1](http://dx.doi.org/10.1038/s41564-024-01763-1).

[4] S. Tian et al. "Making gut microbiomes from scratch". In: _Cell
Host &amp; Microbe_ 30.11 (Nov. 2022), p. 1508–1509. ISSN: 1931-3128.
DOI:
[10.1016/j.chom.2022.10.005](https://doi.org/10.1016%2Fj.chom.2022.10.005).
URL:
[http://dx.doi.org/10.1016/j.chom.2022.10.005](http://dx.doi.org/10.1016/j.chom.2022.10.005).

[5] J. T. Curtis. _The vegetation of Wisconsin-ordination of plant
communities_. Chapter 7. Madison, WI: University of Wisconsin Press,
Dec. 1959.

[6] K. Sankaran et al. "mbtransfer: Microbiome intervention analysis
using transfer functions and mirror statistics". In: _PLOS
Computational Biology_ 20.6 (Jun. 2024). Ed. by N. Nagarajan, p.
e1012196. ISSN: 1553-7358. DOI:
[10.1371/journal.pcbi.1012196](https://doi.org/10.1371%2Fjournal.pcbi.1012196).
URL:
[http://dx.doi.org/10.1371/journal.pcbi.1012196](http://dx.doi.org/10.1371/journal.pcbi.1012196).

[7] G. E. Box et al. "Intervention analysis with applications to
economic and environmental problems". In: _Journal of the American
Statistical association_ 70.349 (1975), pp. 70-79.

[8] C. Dai et al. "False Discovery Rate Control via Data Splitting".
In: _Journal of the American Statistical Association_ 118.544 (May.
2022), p. 2503–2520. ISSN: 1537-274X. DOI:
[10.1080/01621459.2022.2060113](https://doi.org/10.1080%2F01621459.2022.2060113).
URL:
[http://dx.doi.org/10.1080/01621459.2022.2060113](http://dx.doi.org/10.1080/01621459.2022.2060113).

[9] J. H. Friedman. "Greedy function approximation: A gradient
boosting machine." In: _The Annals of Statistics_ 29.5 (Oct. 2001).
ISSN: 0090-5364. DOI:
[10.1214/aos/1013203451](https://doi.org/10.1214%2Faos%2F1013203451).
URL:
[http://dx.doi.org/10.1214/aos/1013203451](http://dx.doi.org/10.1214/aos/1013203451).

[10] G. Thanei et al. "The xyz algorithm for fast interaction search
in high-dimensional data". In: _Journal of Machine Learning Research_
19.37 (2018), pp. 1-42. URL:
[http://jmlr.org/papers/v19/16-515.html](http://jmlr.org/papers/v19/16-515.html).

[11] P. Biecek. "DALEX: Explainers for Complex Predictive Models in
R". In: _Journal of Machine Learning Research_ 19.84 (2018), pp. 1-5.
URL:
[http://jmlr.org/papers/v19/18-416.html](http://jmlr.org/papers/v19/18-416.html).

[12] X. Zhang et al. "Negative binomial mixed models for analyzing
microbiome count data". In: _BMC Bioinformatics_ 18.1 (Jan. 2017).
ISSN: 1471-2105. DOI:
[10.1186/s12859-016-1441-7](https://doi.org/10.1186%2Fs12859-016-1441-7).
URL:
[http://dx.doi.org/10.1186/s12859-016-1441-7](http://dx.doi.org/10.1186/s12859-016-1441-7).

[13] S. Jiang et al. "A Bayesian zero-inflated negative binomial
regression model for the integrative analysis of microbiome data".
In: _Biostatistics_ 22.3 (Dec. 2019), p. 522–540. ISSN: 1468-4357.
DOI:
[10.1093/biostatistics/kxz050](https://doi.org/10.1093%2Fbiostatistics%2Fkxz050).
URL:
[http://dx.doi.org/10.1093/biostatistics/kxz050](http://dx.doi.org/10.1093/biostatistics/kxz050).

---

### References

[14] A. K. Mishra et al. "Negative binomial factor regression with
application to microbiome data analysis". In: _Statistics in
Medicine_ 41.15 (Apr. 2022), p. 2786–2803. ISSN: 1097-0258. DOI:
[10.1002/sim.9384](https://doi.org/10.1002%2Fsim.9384). URL:
[http://dx.doi.org/10.1002/sim.9384](http://dx.doi.org/10.1002/sim.9384).

[15] V. Bucci et al. "MDSINE: Microbial Dynamical Systems INference
Engine for microbiome time-series analyses". In: _Genome Biology_
17.1 (Jun. 2016). ISSN: 1474-760X. DOI:
[10.1186/s13059-016-0980-6](https://doi.org/10.1186%2Fs13059-016-0980-6).
URL:
[http://dx.doi.org/10.1186/s13059-016-0980-6](http://dx.doi.org/10.1186/s13059-016-0980-6).

[16] T. Gibson et al. "Robust and Scalable Models of Microbiome
Dynamics". In: _Proceedings of the 35th International Conference on
Machine Learning_. Ed. by J. Dy and A. Krause. Vol. 80. Proceedings
of Machine Learning Research. PMLR, 2018, pp. 1763-1772. URL:
[https://proceedings.mlr.press/v80/gibson18a.html](https://proceedings.mlr.press/v80/gibson18a.html).

[17] T. E. Gibson et al. "Microbial dynamics inference at
ecosystem-scale".  (Dec. 2021). DOI:
[10.1101/2021.12.14.469105](https://doi.org/10.1101%2F2021.12.14.469105).
URL:
[http://dx.doi.org/10.1101/2021.12.14.469105](http://dx.doi.org/10.1101/2021.12.14.469105).

[18] J. D. Silverman et al. "Dynamic linear models guide design and
analysis of microbiota studies within artificial human guts".  (Apr.
2018). DOI: [10.1101/306597](https://doi.org/10.1101%2F306597). URL:
[http://dx.doi.org/10.1101/306597](http://dx.doi.org/10.1101/306597).

[19] J. D. Silverman et al. "Bayesian Multinomial Logistic Normal
Models through Marginally Latent Matrix-T Processes". In: _Journal of
Machine Learning Research_ 23.7 (2022), pp. 1-42. URL:
[http://jmlr.org/papers/v23/19-882.html](http://jmlr.org/papers/v23/19-882.html).

[20] M. I. Love et al. "Moderated estimation of fold change and
dispersion for RNA-seq data with DESeq2". In: _Genome Biology_ 15.12
(Dec. 2014). ISSN: 1474-760X. DOI:
[10.1186/s13059-014-0550-8](https://doi.org/10.1186%2Fs13059-014-0550-8).
URL:
[http://dx.doi.org/10.1186/s13059-014-0550-8](http://dx.doi.org/10.1186/s13059-014-0550-8).

[21] L. A. David et al. "Diet rapidly and reproducibly alters the
human gut microbiome". In: _Nature_ 505.7484 (Dec. 2013), p. 559–563.
ISSN: 1476-4687. DOI:
[10.1038/nature12820](https://doi.org/10.1038%2Fnature12820). URL:
[http://dx.doi.org/10.1038/nature12820](http://dx.doi.org/10.1038/nature12820).

[22] D. Gonze et al. "Microbial communities as dynamical systems".
In: _Current Opinion in Microbiology_ 44 (Aug. 2018), p. 41–49. ISSN:
1369-5274. DOI:
[10.1016/j.mib.2018.07.004](https://doi.org/10.1016%2Fj.mib.2018.07.004).
URL:
[http://dx.doi.org/10.1016/j.mib.2018.07.004](http://dx.doi.org/10.1016/j.mib.2018.07.004).

[23] T. A. Joseph et al. "Compositional Lotka-Volterra describes
microbial dynamics in the simplex". In: _PLOS Computational Biology_
16.5 (May. 2020). Ed. by V. Dakos, p. e1007917. ISSN: 1553-7358. DOI:
[10.1371/journal.pcbi.1007917](https://doi.org/10.1371%2Fjournal.pcbi.1007917).
URL:
[http://dx.doi.org/10.1371/journal.pcbi.1007917](http://dx.doi.org/10.1371/journal.pcbi.1007917).

---

### Phylogenetic Covariance

---

### Example Trajectories

---

### Computation Time

---

### Transfer Functions

.pull-left[
1. The pattern of autoregressive coefficients determines `$y_{t}$`'s response to
 "pulse" and "step" interventions.

1. Note that, unlike standard generalized Lotka-Volterra models [22; 23], they can model lagged intervention effects.
]

.pull-right[
<img src="figure/intervention-analysis.png" width="400" style="display: block; margin: auto;" />
]

---

### Generative Mechanism

.center[
<img src="figure/simulated_data.png" width=800/><br/>
]
<span style="font-size: 18px">
Taxonomic abundances in the simulated negative binomial data, transformed to
taxon-specific quantiles. The intervention windows are surrounded in red. Only
the top three rows have true intervention effects.
</span>

---

### Figure Attributions

math function by Ralf Schmitzer from <a href="https://thenounproject.com/browse/icons/term/math-function/" target="_blank" title="math function Icons">Noun Project</a> (CC BY 3.0)

biological tree by Leslie Coonrod from <a href="https://thenounproject.com/browse/icons/term/biological-tree/" target="_blank" title="biological tree Icons">Noun Project</a> (CC BY 3.0)

standard normal distribution by Nick Taras from <a href="https://thenounproject.com/browse/icons/term/standard-normal-distribution/" target="_blank" title="standard normal distribution Icons">Noun Project</a> (CC BY 3.0)

microbiome by Stefanie Peschel from <a href="https://thenounproject.com/browse/icons/term/microbiome/" target="_blank" title="microbiome Icons">Noun Project</a> (CC BY 3.0)