1 Load packages

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")

2 Compute distances in the data

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")

3 Explore relationships within the mediation triangle variables

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