library(phyloseq); packageVersion("phyloseq")
## [1] '1.28.0'
library(ggplot2); packageVersion("ggplot2")
## [1] '3.2.0'
Source the MODIMA: Multivariate Omnibus Distance Mediation Analysis code.
library(energy)
source("modima.R")
Load STAT data with metabolic phenotype variables.
load("STAT2.RData")
We are going to be testing if microbiota mediates effect of antibiotic exposure on metabolic phenotypes.
Exposures:
phy_pheno = subset_samples(phy2, Location == "cecal")
sample_data(phy_pheno)$Abx = factor(sample_data(phy_pheno)$Treatment != "C")
abx_dist = dist(cbind(sample_data(phy_pheno)$Abx))
tx_dist = as.dist(outer(sample_data(phy_pheno)$Treatment,
sample_data(phy_pheno)$Treatment,
"!="))
Response: percent fat
fat_dist = dist(sample_data(phy_pheno)$pFat)
Mediator: microbiota with bray curtis distance
bray_dist = phyloseq::distance(phy_pheno, method="jsd")
Exposure vs. response
ggplot(data=sample_data(phy_pheno), aes(x=Abx,
y = pFat)) +
geom_boxplot() +
geom_jitter(width = .05, size=1,
aes(color=Abx))
with(sample_data(phy_pheno), t.test(pFat~Abx))
##
## Welch Two Sample t-test
##
## data: pFat by Abx
## t = -4.1999, df = 30.851, p-value = 0.0002105
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.98912 -1.38088
## sample estimates:
## mean in group FALSE mean in group TRUE
## 20.170 22.855
Exposure vs. mediator
phy.o = ordinate(phy_pheno, method="MDS", distance=bray_dist)
plot_ordination(physeq = phy_pheno, ordination = phy.o,
color="Treatment")
# Recall
source("https://raw.githubusercontent.com/alekseyenko/WdStar/master/Wd.R")
WdS.test(dm = bray_dist,
f = sample_data(phy_pheno)$Treatment,
nrep = 999)
## $p.value
## [1] 0.001
##
## $statistic
## [1] 4.535805
##
## $nrep
## [1] 999
Tw2.posthoc.1vsAll.tests(dm = bray_dist,
f = sample_data(phy_pheno)$Treatment,
nrep = 999)
## N1 N2 p.value tw2.stat nrep
## C 40 10 0.001 4.717367 999
## P 40 10 0.001 4.788102 999
## T 40 10 0.009 2.529902 999
## V 40 10 0.001 4.119393 999
## VP 40 10 0.001 5.91846 999
Mediator vs. Response given exposure
#Mediator given response
library(ade4)
MgivenE = wca(dudi.pco(bray_dist, scannf=F, nf=3),
sample_data(phy_pheno)$Treatment,
scannf=F, nf=3)
## Warning in dudi.pco(bray_dist, scannf = F, nf = 3): Non euclidean distance
RgivenE = with(sample_data(phy_pheno), resid(lm(pFat ~ Treatment)))
plot(RgivenE, MgivenE$li[,1])
cor.test(RgivenE, MgivenE$li[,1])
##
## Pearson's product-moment correlation
##
## data: RgivenE and MgivenE$li[, 1]
## t = 0.65893, df = 48, p-value = 0.5131
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1886376 0.3634508
## sample estimates:
## cor
## 0.09468158
plot(RgivenE, MgivenE$li[,2])
cor.test(RgivenE, MgivenE$li[,1])
##
## Pearson's product-moment correlation
##
## data: RgivenE and MgivenE$li[, 1]
## t = 0.65893, df = 48, p-value = 0.5131
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1886376 0.3634508
## sample estimates:
## cor
## 0.09468158
Given that the correlations of the mediator (principal axes) and response given exposure is pretty small and not statistically significant, we expect that the mediation effect will be small too. We do not expect a significant test result.
modima(exposure = abx_dist,
mediator = bray_dist,
response = fat_dist, nrep=9999)
##
## MODIMA: a Method for Multivariate Omnibus Distance Mediation
## Analysis
##
## data: exposure = abx_dist
## mediator = bray_dist
## response = fat_dist
## number of permutations= 10000
## sample estimates are
## -bias-corrected distance correlation (energy::bcdcor) of indicated pairs and
## -partial distance correlation (energy::pdcor) of indicated triple
## MODIMA stat = 0.0022072, p-value = 0.3072
## sample estimates:
## exposure–mediator bcdcor exposure–response bcdcor
## 0.21410985 0.14728477
## mediator–response bcdcor mediator–response–exposure pdcor
## 0.04149499 0.01030876