The Mindfulness Study
In most microbiome mediation analysis studies, microbiome composition is thought to mediate some outcome of interest. Sometimes, though, we might want to treat the microbiome as an outcome. Our motivating example here looks at the relationship between mindfulness and the microbiome. In this setting, we might imagine that mindfulness training influences behaviors (diet, sleep, exercise) that could then change the microbiome. Clinical surveys that gauge those behaviors now play the role of mediator. Alternatively, mindfulness could trigger neurophysiological changes that influence the microbiome in more subtle ways. The point of mediation analysis is to provide evidence for one or the other of these possibilities.
library(ggdist)
library(ggrepel)
library(phyloseq)
library(tidyverse)
library(multimedia)
data(mindfulness)
set.seed(20231228)This package needs a mediation_data object to organize
the pretreatment, treatment, mediator, and outcome variables.
taxa_names(mindfulness) specifies that we want all the taxa
names to be outcome variables. We’ve already appended all the mediator
variable names with the prefix “mediator,” so we can specify them using
the starts_with selector. We’ll treat the
subject column as pretreatment variables. The study tracked
participants across a few timepoints, and we want to control for
different baseline behavioral and microbiome profiles.
exper <- mediation_data(
mindfulness,
taxa_names(mindfulness),
"treatment",
starts_with("mediator"),
"subject"
)
exper#> [Mediation Data]
#> 307 samples with measurements for,
#> 1 treatment: treatment
#> 4 mediators: mediator.sleepDistRaw, mediator.fatigueRaw, ...
#> 55 outcomes: Defluviitalea, Gordonibacter, ...
#> 1 pretreatment: subject
This block applies a logistic-normal multinomial model to the microbiome outcome. This is an example where the response variables are all modeled jointly. We’ll use a ridge regression model for each of the mediators. This is a better choice than ordinary linear models because we have to estimate many subject-level parameters, and our sample size is not that large.
# use Flavonifractor as the reference for the log-ratio transform
outcomes(exper) <- outcomes(exper) |>
select(Flavonifractor, everything())
model <- multimedia(
exper,
lnm_model(seed = .Random.seed[1]),
glmnet_model(lambda = 0.5, alpha = 0)
) |>
estimate(exper)#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.012067 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 120.67 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -243307.612 1.000 1.000
#> Chain 1: 200 -99580.409 1.222 1.443
#> Chain 1: 300 -90008.408 0.850 1.000
#> Chain 1: 400 -87376.289 0.645 1.000
#> Chain 1: 500 -85962.756 0.519 0.106
#> Chain 1: 600 -85132.538 0.434 0.106
#> Chain 1: 700 -84674.011 0.373 0.030
#> Chain 1: 800 -84249.244 0.327 0.030
#> Chain 1: 900 -84035.521 0.291 0.016
#> Chain 1: 1000 -83782.160 0.262 0.016
#> Chain 1: 1100 -83487.854 0.163 0.010 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
Evaluating Effects
With this model, we can estimate direct and indirect effects. We’re averaging over different settings of the mediator treatment variables. Since the LNM model outputs relative abundance profiles, all the estimated effects should be interpreted on a probability scale.
direct <- direct_effect(model, exper) |>
group_by(outcome) |>
summarise(direct_effect = mean(direct_effect)) |>
arrange(-abs(direct_effect))These are the indirect effects when we set mediators to the treatment setting one at a time.
indirect <- indirect_pathwise(model, exper) |>
group_by(mediator, outcome) |>
summarise(indirect_effect = mean(indirect_effect)) |>
arrange(-abs(indirect_effect))
indirect#> # A tibble: 220 × 3
#> # Groups: mediator [4]
#> mediator outcome indirect_effect
#> plot_mediators visualizes the relationship between
mediator and outcome variables. The title for each subplot is the
estimated indirect effect (note that it’s Control -
Treatment in the current definition). For example, in the
panel relating sleep disturbances and Blautia, there is a rough positive
association between these variables. The control group has a higher
number of sleep disturbances, which is translated into a positive
indirect effect for this pair. The model accounts for more than these
pairwise relationships, though – it controls for subject-level baselines
and includes direct effects of treatment – so these figures should be
interpreted within the larger model context.
exper_rela <- exper
outcomes(exper_rela) <- normalize(outcomes(exper))
plot_mediators(indirect, exper_rela)
Model Alterations
Hypothesis testing is built on the idea that a simple submodel might explain the data just as well as a full, more complicated one. In mediation analysis, we can zero out sets of edges to define submodels. For example, we can re-estimate a model that does not allow any relationship between the treatment and the mediators.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.008222 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 82.22 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -243307.612 1.000 1.000
#> Chain 1: 200 -99580.409 1.222 1.443
#> Chain 1: 300 -90008.408 0.850 1.000
#> Chain 1: 400 -87376.289 0.645 1.000
#> Chain 1: 500 -85962.756 0.519 0.106
#> Chain 1: 600 -85132.538 0.434 0.106
#> Chain 1: 700 -84674.011 0.373 0.030
#> Chain 1: 800 -84249.244 0.327 0.030
#> Chain 1: 900 -84035.521 0.291 0.016
#> Chain 1: 1000 -83782.160 0.262 0.016
#> Chain 1: 1100 -83487.854 0.163 0.010 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
These modified models can be used to simulate synthetic null data that help contextualize the effects seen in real data. For example, we can compare imaginary cohort drawn from the original and the altered models. Notice that I re-assigned all the treatments. This removes potential confounding between the treatment and subject level effects – this effect doesn’t exist on average across the population because the study was randomized, but we want to remove any associations present in our finite sample.
new_assign <- treatments(exper)[sample(nrow(treatments(exper))), , drop = FALSE]
profile <- setup_profile(model, new_assign, new_assign)
m1 <- sample(model, profile = profile, pretreatment = pretreatments(exper))
m2 <- sample(altered, profile = profile, pretreatment = pretreatments(exper))The sampled objects above have the same class as exper.
We can extract the real and simulated mediator data to see the treatment
effects along those edges. As expected, the altered model has no
systematic difference between treatment groups. This lends some evidence
to the conclusion that mindfulness does influence the mediators, which
was implicit in our discussion of indirect effects above.
list(
real = bind_cols(treatments(exper), mediators(exper), pretreatments(exper)),
original = bind_cols(new_assign, mediators(m1), pretreatments(exper)),
altered = bind_cols(new_assign, mediators(m2), pretreatments(exper))
) |>
bind_rows(.id = "source") |>
pivot_longer(starts_with("mediator"), names_to = "mediator") |>
ggplot() +
geom_boxplot(
aes(value, reorder(mediator, value, median), fill = treatment)
) +
facet_grid(~source)
Here is the analogous alteration that removes all direct effects. The x-axis is on a log scale, and the LNM’s estimated treatment effects for rare taxa are not practically significant. It’s interesting that even though the original data don’t show much difference between treatment groups, the model suggests that there are real effects, even for the abundant taxa. It’s possible that, in the raw boxplots, heterogeneity across subjects and mediators masks true effects, and the LNM is able to detect these effects by appropriately accounting for sources of heterogeneity.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.008601 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 86.01 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -286918.019 1.000 1.000
#> Chain 1: 200 -104917.320 1.367 1.735
#> Chain 1: 300 -92595.460 0.956 1.000
#> Chain 1: 400 -88954.112 0.727 1.000
#> Chain 1: 500 -87383.372 0.585 0.133
#> Chain 1: 600 -86258.965 0.490 0.133
#> Chain 1: 700 -85569.685 0.421 0.041
#> Chain 1: 800 -84971.656 0.369 0.041
#> Chain 1: 900 -84708.412 0.329 0.018
#> Chain 1: 1000 -84309.170 0.296 0.018
#> Chain 1: 1100 -84071.410 0.197 0.013
#> Chain 1: 1200 -83880.348 0.023 0.008 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.009878 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 98.78 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -266185.924 1.000 1.000
#> Chain 1: 200 -103322.121 1.288 1.576
#> Chain 1: 300 -90720.935 0.905 1.000
#> Chain 1: 400 -87659.862 0.688 1.000
#> Chain 1: 500 -86393.749 0.553 0.139
#> Chain 1: 600 -85556.774 0.462 0.139
#> Chain 1: 700 -84583.951 0.398 0.035
#> Chain 1: 800 -84086.834 0.349 0.035
#> Chain 1: 900 -83825.808 0.311 0.015
#> Chain 1: 1000 -83560.389 0.280 0.015
#> Chain 1: 1100 -83420.529 0.180 0.012
#> Chain 1: 1200 -83137.414 0.023 0.010 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
profile <- setup_profile(model, treatments(exper), treatments(exper))
y1 <- sample(model, profile = profile, pretreatment = pretreatments(exper))
y2 <- sample(
altered_direct,
profile = profile,
pretreatment = pretreatments(exper)
)
y3 <- sample(
altered_indirect,
profile = profile,
pretreatment = pretreatments(exper)
)taxa_order <- normalize(outcomes(exper)) |>
log() |>
apply(2, \(x) {
x[is.infinite(x)] <- NA
median(x, na.rm = TRUE)
}) |>
order()
combined_samples <- list(
"Original Data" = normalize(outcomes(exper)),
"Full Model" = normalize(outcomes(y1)),
"No Direct Effect" = normalize(outcomes(y2)),
"No Indirect Effect" = normalize(outcomes(y3))
) |>
bind_rows(.id = "source") |>
mutate(treatment = rep(treatments(exper)$treatment, 4)) |>
pivot_longer(colnames(outcomes(exper))) |>
mutate(name = factor(name, levels = colnames(outcomes(exper))[taxa_order]))
ggplot(combined_samples) +
geom_boxplot(
aes(log(value), name, fill = treatment),
size = 0.4, outlier.size = 0.8, position = "identity",
alpha = 0.7
) +
labs(x = "log(Relative Abundance)", y = "Genus", fill = "", col = "") +
facet_grid(~source, scales = "free") +
theme(panel.border = element_rect(fill = NA))
combined_samples |>
group_by(name, source, treatment) |>
summarise(presence = mean(value > 0)) |>
ggplot() +
geom_col(
aes(presence, name, fill = treatment, col = treatment),
position = "identity", width = 1,
alpha = 0.7
) +
facet_grid(~source, scales = "free") +
scale_x_continuous(expand = c(0, 0)) +
theme(panel.border = element_rect(fill = NA))
Synthetic Null Testing
A standard approach to inference is to apply a bootstrap (this is done at the end of the vignette). An alternative, which is less common, but well-suited to high-dimensions is to calibrate variable selection sets using synthetic null data. For example, this is the intuition behind the knockoff, Clipper, and mirror statistic algorithms. The general recipe is:
- Simulate synthetic null data from a model/sampling mechanism where we know for certain that an effect does not exist.
- Use the full model to estimate effects on this synthetic null data.
- Define a selection rule so that the fraction of selected synthetic null effects is kept below an false discovery rate threshold.
The synthetic null effects essentially serve as negative controls to calibrate inference.
In multimedia, null_contrast estimates
effects on real and synthetic null data, while fdr_summary
ranks and selects features. Each point in the figure below is the
directed effect for a single feature in either real or synthetic data.
We rank features according to their absolute effect size, and as soon as
negative control features enter, we begin increasing the estimated false
discovery rate for the rule that selects features up to that point.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.00872 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 87.2 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -286918.019 1.000 1.000
#> Chain 1: 200 -104917.320 1.367 1.735
#> Chain 1: 300 -92595.460 0.956 1.000
#> Chain 1: 400 -88954.112 0.727 1.000
#> Chain 1: 500 -87383.372 0.585 0.133
#> Chain 1: 600 -86258.965 0.490 0.133
#> Chain 1: 700 -85569.685 0.421 0.041
#> Chain 1: 800 -84971.656 0.369 0.041
#> Chain 1: 900 -84708.412 0.329 0.018
#> Chain 1: 1000 -84309.170 0.296 0.018
#> Chain 1: 1100 -84071.410 0.197 0.013
#> Chain 1: 1200 -83880.348 0.023 0.008 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.009524 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 95.24 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -584473.269 1.000 1.000
#> Chain 1: 200 -163530.307 1.787 2.574
#> Chain 1: 300 -107654.318 1.364 1.000
#> Chain 1: 400 -99725.734 1.043 1.000
#> Chain 1: 500 -97692.416 0.839 0.519
#> Chain 1: 600 -96484.766 0.701 0.519
#> Chain 1: 700 -95723.900 0.602 0.080
#> Chain 1: 800 -95435.003 0.527 0.080
#> Chain 1: 900 -94649.591 0.469 0.021
#> Chain 1: 1000 -94481.985 0.423 0.021
#> Chain 1: 1100 -94105.353 0.323 0.013
#> Chain 1: 1200 -93857.267 0.066 0.008 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
The indirect effect analog is shown below. It seems we have more evidence for direct than indirect effects. This is consistent with literature on the difficulty of estimating indirect effects in mediation analysis.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.009584 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 95.84 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -266185.924 1.000 1.000
#> Chain 1: 200 -103322.121 1.288 1.576
#> Chain 1: 300 -90720.935 0.905 1.000
#> Chain 1: 400 -87659.862 0.688 1.000
#> Chain 1: 500 -86393.749 0.553 0.139
#> Chain 1: 600 -85556.774 0.462 0.139
#> Chain 1: 700 -84583.951 0.398 0.035
#> Chain 1: 800 -84086.834 0.349 0.035
#> Chain 1: 900 -83825.808 0.311 0.015
#> Chain 1: 1000 -83560.389 0.280 0.015
#> Chain 1: 1100 -83420.529 0.180 0.012
#> Chain 1: 1200 -83137.414 0.023 0.010 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.010941 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 109.41 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -428484.112 1.000 1.000
#> Chain 1: 200 -123978.876 1.728 2.456
#> Chain 1: 300 -103202.512 1.219 1.000
#> Chain 1: 400 -99097.492 0.925 1.000
#> Chain 1: 500 -97271.647 0.744 0.201
#> Chain 1: 600 -96487.144 0.621 0.201
#> Chain 1: 700 -95306.014 0.534 0.041
#> Chain 1: 800 -94767.004 0.468 0.041
#> Chain 1: 900 -94420.053 0.416 0.019
#> Chain 1: 1000 -94194.245 0.375 0.019
#> Chain 1: 1100 -93791.799 0.275 0.012
#> Chain 1: 1200 -93423.006 0.030 0.008 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
fdr_indirect <- fdr_summary(contrast, "indirect_overall") |>
dplyr::rename(effect = indirect_effect)fdr_data <- bind_rows(fdr_direct, fdr_indirect, .id = "effect_type") |>
mutate(effect_type = c("Direct", "Indirect")[as.integer(effect_type)])
ggplot(fdr_data, aes(effect, fdr_hat)) +
geom_text_repel(
data = filter(fdr_data, keep),
aes(label = outcome, col = source), size = 3, force = 10
) +
labs(x = "Effect", y = expression(widehat(FDR)), col = "Data Source") +
geom_point(aes(col = source, size = keep)) +
scale_size_discrete("Selected", range = c(.5, 3), guide = FALSE) +
scale_y_continuous(limits = c(-0.1, 0.55)) +
facet_wrap(~effect_type, scales = "free_x") +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
strip.text = element_text(size = 14),
legend.text = element_text(size = 12)
)
bind_mediation(exper_rela) |>
select(
subject, starts_with("mediator"), treatment, Roseburia, Blautia,
Faecalibacterium
) |>
pivot_longer(starts_with("mediator"), names_to = "mediator") |>
mutate(mediator = str_remove(mediator, "mediator.")) |>
pivot_longer(
Roseburia:Faecalibacterium,
names_to = "taxon", values_to = "abundance"
) |>
group_by(subject, mediator, treatment, taxon) |>
summarise(abundance = mean(abundance), value = mean(value)) |>
ggplot() +
geom_point(aes(value, abundance, col = treatment)) +
facet_grid(taxon ~ mediator, scales = "free") +
labs(y = "Relative Abundance", x = "Mediator Value", col = "") +
theme(
axis.title = element_text(size = 14),
strip.text = element_text(size = 12),
strip.text.y = element_text(size = 12, angle = 0),
legend.text = element_text(size = 12),
panel.border = element_rect(fill = NA)
)
Here is the analog for indirect pathwise effects. Since these effects distinguish between individual mediators, we can facet by these variables.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.00966 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 96.6 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -266185.924 1.000 1.000
#> Chain 1: 200 -103322.121 1.288 1.576
#> Chain 1: 300 -90720.935 0.905 1.000
#> Chain 1: 400 -87659.862 0.688 1.000
#> Chain 1: 500 -86393.749 0.553 0.139
#> Chain 1: 600 -85556.774 0.462 0.139
#> Chain 1: 700 -84583.951 0.398 0.035
#> Chain 1: 800 -84086.834 0.349 0.035
#> Chain 1: 900 -83825.808 0.311 0.015
#> Chain 1: 1000 -83560.389 0.280 0.015
#> Chain 1: 1100 -83420.529 0.180 0.012
#> Chain 1: 1200 -83137.414 0.023 0.010 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.010124 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 101.24 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -423667.039 1.000 1.000
#> Chain 1: 200 -128902.987 1.643 2.287
#> Chain 1: 300 -103520.469 1.177 1.000
#> Chain 1: 400 -99192.987 0.894 1.000
#> Chain 1: 500 -97231.301 0.719 0.245
#> Chain 1: 600 -95809.891 0.602 0.245
#> Chain 1: 700 -95124.335 0.517 0.044
#> Chain 1: 800 -94517.067 0.453 0.044
#> Chain 1: 900 -94169.813 0.403 0.020
#> Chain 1: 1000 -93783.792 0.363 0.020
#> Chain 1: 1100 -93318.660 0.264 0.015
#> Chain 1: 1200 -93157.243 0.035 0.007 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
fdr <- fdr_summary(contrast, "indirect_pathwise")
ggplot(fdr, aes(indirect_effect, fdr_hat)) +
geom_text_repel(
data = filter(fdr, keep),
aes(label = outcome, fill = source), size = 4
) +
labs(x = "Indirect Effect", y = "Estimated False Discovery Rate") +
geom_point(aes(col = source, size = keep)) +
scale_size_discrete(range = c(.2, 3)) +
facet_wrap(~mediator)
Random Forests
The LNM makes relatively strong assumptions on the outcome model: The effects must enter linearly through a multinomial likelihood. As a nonparametric alternative, we can fit separate random forest models to each microbe’s abundance. Since this model doesn’t directly account for compositions, it will help to transform the response to relative abundances.
model <- multimedia(
exper_rela,
rf_model(num.trees = 50),
glmnet_model(lambda = 0.5, alpha = 0)
) |>
estimate(exper_rela)We can compute direct and indirect effects from the estimated model.
direct <- direct_effect(model, exper_rela)
direct_sorted <- direct |>
group_by(outcome, contrast) |>
summarise(direct_effect = mean(direct_effect)) |>
arrange(-direct_effect)
indirect_sorted <- indirect_pathwise(model, exper_rela) |>
group_by(outcome, mediator, contrast) |>
summarise(indirect_effect = mean(indirect_effect)) |>
arrange(-abs(indirect_effect))
direct_sorted#> # A tibble: 55 × 3
#> # Groups: outcome [55]
#> outcome contrast direct_effect
#> #> # A tibble: 220 × 4
#> # Groups: outcome, mediator [220]
#> outcome mediator contrast indirect_effect
#> We can visualize some of the features with the largest direct and indirect effects. Nonlinearity seems to have helped uncover a few mediator to microbe relationships, for example, the relationships between sleep disturbance/Ruminococcus2 and cereal/Faecalibacterium.
Bootstrap
The bootstrap can help evaluate uncertainty for more complex
statistics, and it’s a good option for our random forest analysis. In
each iteration, we randomly subsample timepoints from the original
experiment and compute a list of estimates, specified by the list
funs below. The output concatenates output across all
bootstrap iterations.
funs <- list(direct = direct_effect, indirect = indirect_overall)
# inference <- bootstrap(model, exper_rela, funs, B = 1e3)
inference <- readRDS(url("https://go.wisc.edu/ba3e87"))We can now visualize the bootstrap confidence intervals associated with the estimated effects. The largest interval corresponds to a 90% bootstrap confidence interval. The inner interval covers the median +/- one-third of the bootstrap distribution’s mass. Some associations certainly seem worth following up on, though more data will be needed for validation.
Notice that the direct effects seem generally stronger than indirect effects. Even with a flexible random forest model, it seems that changes in diet and sleep alone are not enough to explain the changes in microbiome composition seen in this study. This suggests that mechanisms related to the gut-brain axis may be at play, lying behind the observed direct effects.
p1 <- ggplot(inference$indirect) +
geom_vline(xintercept = 0) +
stat_slabinterval(
aes(indirect_effect, reorder(outcome, indirect_effect, median)),
.width = c(.66, .95)
) +
labs(
title = "Indirect Effects", y = "Taxon", x = "Effect (Rel. Abund Scale)"
)
p2 <- ggplot(inference$direct) +
geom_vline(xintercept = 0) +
stat_slabinterval(
aes(direct_effect, reorder(outcome, direct_effect, median)),
.width = c(.66, .95)
) +
labs(title = "Direct Effects", y = "Taxon", x = "Effect (Rel. Abund Scale)")
p1 | p2
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] phyloseq_1.49.0 ggrepel_0.9.6 multimedia_0.2.0 tidyselect_1.2.1
#> [5] ranger_0.16.0 glmnetUtils_1.1.9 vroom_1.6.5 lubridate_1.9.3
#> [9] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
#> [13] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 tidyverse_2.0.0
#> [17] patchwork_1.3.0 glue_1.7.0 ggdist_3.3.2 compositions_2.0-8
#> [21] brms_2.21.0 Rcpp_1.0.13 ggplot2_3.5.1 rmarkdown_2.28
#>
#> loaded via a namespace (and not attached):
#> [1] sys_3.4.2 tensorA_0.36.2.1
#> [3] jsonlite_1.8.8 shape_1.4.6.1
#> [5] magrittr_2.0.3 farver_2.1.2
#> [7] zlibbioc_1.51.1 vctrs_0.6.5
#> [9] multtest_2.61.0 htmltools_0.5.8.1
#> [11] S4Arrays_1.5.7 progress_1.2.3
#> [13] distributional_0.5.0 curl_5.2.2
#> [15] Rhdf5lib_1.27.0 SparseArray_1.5.36
#> [17] rhdf5_2.49.0 sass_0.4.9
#> [19] StanHeaders_2.32.10 bslib_0.8.0
#> [21] plyr_1.8.9 cachem_1.1.0
#> [23] buildtools_1.0.0 igraph_2.0.3
#> [25] lifecycle_1.0.4 iterators_1.0.14
#> [27] pkgconfig_2.0.3 Matrix_1.7-0
#> [29] R6_2.5.1 fastmap_1.2.0
#> [31] GenomeInfoDbData_1.2.12 MatrixGenerics_1.17.0
#> [33] digest_0.6.37 colorspace_2.1-1
#> [35] S4Vectors_0.43.2 miniLNM_0.1.0
#> [37] GenomicRanges_1.57.1 vegan_2.6-8
#> [39] labeling_0.4.3 timechange_0.3.0
#> [41] fansi_1.0.6 httr_1.4.7
#> [43] abind_1.4-8 mgcv_1.9-1
#> [45] compiler_4.4.1 bit64_4.0.5
#> [47] withr_3.0.1 backports_1.5.0
#> [49] inline_0.3.19 highr_0.11
#> [51] QuickJSR_1.3.1 pkgbuild_1.4.4
#> [53] MASS_7.3-61 bayesm_3.1-6
#> [55] DelayedArray_0.31.11 biomformat_1.33.0
#> [57] loo_2.8.0 permute_0.9-7
#> [59] tools_4.4.1 ape_5.8
#> [61] nlme_3.1-166 rhdf5filters_1.17.0
#> [63] grid_4.4.1 checkmate_2.3.2
#> [65] cluster_2.1.6 reshape2_1.4.4
#> [67] ade4_1.7-22 generics_0.1.3
#> [69] operator.tools_1.6.3 gtable_0.3.5
#> [71] tzdb_0.4.0 formula.tools_1.7.1
#> [73] data.table_1.16.0 hms_1.1.3
#> [75] tidygraph_1.3.1 utf8_1.2.4
#> [77] XVector_0.45.0 BiocGenerics_0.51.1
#> [79] foreach_1.5.2 pillar_1.9.0
#> [81] posterior_1.6.0 robustbase_0.99-4
#> [83] splines_4.4.1 lattice_0.22-6
#> [85] bit_4.0.5 survival_3.7-0
#> [87] maketools_1.3.0 Biostrings_2.73.1
#> [89] knitr_1.48 gridExtra_2.3
#> [91] V8_5.0.0 IRanges_2.39.2
#> [93] SummarizedExperiment_1.35.1 stats4_4.4.1
#> [95] xfun_0.47 bridgesampling_1.1-2
#> [97] Biobase_2.65.1 matrixStats_1.4.1
#> [99] DEoptimR_1.1-3 rstan_2.32.6
#> [101] stringi_1.8.4 UCSC.utils_1.1.0
#> [103] yaml_2.3.10 evaluate_1.0.0
#> [105] codetools_0.2-20 cli_3.6.3
#> [107] RcppParallel_5.1.9 munsell_0.5.1
#> [109] jquerylib_0.1.4 GenomeInfoDb_1.41.1
#> [111] coda_0.19-4.1 parallel_4.4.1
#> [113] rstantools_2.4.0 prettyunits_1.2.0
#> [115] bayesplot_1.11.1 Brobdingnag_1.2-9
#> [117] glmnet_4.1-8 mvtnorm_1.3-1
#> [119] scales_1.3.0 crayon_1.5.3
#> [121] rlang_1.1.4