Bayesian Hypothesis Testing
Under the hierarchical model, the null hypothesis of equal genetic proportions between two populations, say \(k = 1,2\), is \(H_0: \boldsymbol \rho_1 = \boldsymbol \rho_2 = \boldsymbol \rho_{12}\). Our Bayesian hypothesis-testing strategy is implemented in two steps:
- Sample from the posterior distribution \(p(\boldsymbol \alpha\mid \boldsymbol Y, H_0)\).
- Sample from the posterior distribution of either classical test statistic, \(T = \mathcal X\) or \(\Lambda\), thus obtaining a posterior p-value \[ \textrm{p}_\textrm{v}^{\textrm{post}}= \textrm{Pr}(T > T_\textrm{obs} \mid \boldsymbol Y, H_0). \]
MCMC Sampling from the Posterior Distribution
In order to estimate \(\boldsymbol
\alpha\) under \(H_0\), we apply
the Dirichlet-Multinomial distribution to \(K-1\) lakes, with the first two being
collapsed into a common count vector \(\boldsymbol Y_{12} = \boldsymbol Y_1 + \boldsymbol
Y_2\) with sample size \(N_{12} = N_1 +
N_2\). MCMC sampling is implemented via a Hybrid Monte Carlo
(HMC) variant provided by the R package
rstan. In MADPop, sampling from \(p(\boldsymbol \alpha\mid \boldsymbol Y,
H_0)\) is accomplished with the function hUM.post(),
where “hUM” stands for hierarchical Unconstrained Multinomial
model. We start by merging the two samples from equal distributions
under \(H_0\):
## [1] "Dickey" "Simcoe"eqId0 <- paste0(popId, collapse = ".") # merged lake name
popId0 <- as.character(fish215$Lake) # all lake names
popId0[popId0 %in% popId] <- eqId0
table(popId0, dnn = NULL) # merged lake counts## Crystal Dickey.Simcoe Hogan Kingscote Macdonald
## 20 40 21 18 19
## Manitou Michipicoten Opeongo Seneca Slate
## 20 20 20 18 19Next, we sample from \(p(\boldsymbol \alpha\mid \boldsymbol Y, H_0)\) using rstan:
X0 <- cbind(Id = popId0, fish215[-1]) # allele data with merged lakes
nsamples <- 1e4
fit0 <- hUM.post(nsamples = nsamples, X = X0,
rhoId = eqId0, # output rho only for merged lakes
chains = 1, # next two arguments are passed to rstan
warmup = min(1e4, floor(nsamples/10)),
full.stan.out = FALSE)##
## SAMPLING FOR MODEL 'DirichletMultinomial' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 5.9e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 1001 / 10000 [ 10%] (Sampling)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Sampling)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Sampling)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Sampling)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.83 seconds (Warm-up)
## Chain 1: 7.557 seconds (Sampling)
## Chain 1: 8.387 seconds (Total)
## Chain 1:rstan typically outputs a number of messages and
warnings, many of which are harmless. The MADPop
package lists rstan as a dependency, such that its
sophisticated tuning mechanism is exposed. Optionally,
hUM.post() returns the entire rstan output
(full.stan.out = TRUE), which can be run through
rstan’s MCMC convergence diagnostics.
By setting the rhoId argument of
hUM.post(), the fit0 object contains samples
from \(p(\boldsymbol \alpha, \boldsymbol
\rho_{12} \mid \boldsymbol Y, H_0)\). Boxplots of the 50 highest
posterior probability genotypes in the common distribution \(\boldsymbol \rho_{12}\) are displayed
below:
rho.post <- fit0$rho[,1,] # p(rho12 | Y)
# sort genotype counts by decreasing order
rho.count <- colSums(Xsuff$tab[popId,])
rho.ord <- order(colMeans(rho.post), decreasing = TRUE)
# plot
nbxp <- 50 # number of genotypes for which to make boxplots
clrs <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2",
"#D55E00", "#CC79A7")
rho.ct <- rho.count[rho.ord[1:nbxp]]
rho.lgd <- unique(sort(rho.ct))
rho.box <- rho.post[,rho.ord[1:nbxp]]
par(mar = c(5,4.5,.5,.5)+.1, cex.axis = .7)
boxplot(x = rho.box,
las = 2, col = clrs[rho.ct+1], pch = 16, cex = .2)
title(xlab = "Genotype", line = 4)
title(ylab = expression(p(bold(rho)[(12)]*" | "*bold(Y))))
legend("topright", legend = rev(rho.lgd), fill = rev(clrs[rho.lgd+1]),
title = "Counts", ncol = 2, cex = .8)
Calculating the Posterior P-Value
This is calculated analogously to the bootstrapped p-value, by generating \(M\) contingency tables with \(\boldsymbol Y_k^{(m)} \stackrel{\textrm{ind}}{\sim}\textrm{Multinomial}(N_k, \boldsymbol \rho_{12}^{(m)}\), \(m = 1, \ldots, M\). However, for each contingency table, the common probability vector \(\boldsymbol \rho_{12}^{(m)}\) is different, corresponding to a random draw from \(p(\boldsymbol \rho_{12} \mid \boldsymbol Y, H_0)\). The test statistic \(T^{(m)}\) is then calculated and we have \[ \textrm{p}_\textrm{v}^{\textrm{post}}= \frac 1 M \sum_{m=1}^M \delta(T^{(m)} \ge T_{\textrm{obs}}). \] \(\textrm{p}_\textrm{v}^{\textrm{post}}\) is calculated with MADPop as follows:
## user system elapsed
## 0.387 0.032 0.422We can now compare the three types of p-values (asymptotic, bootstrapped, posterior) for each test statistic (\(\mathcal X\) and \(\Lambda\)):
| Asymptotic | Bootstrap | Posterior | |
|---|---|---|---|
| \(\mathcal X\) | 33.0 | 0.68 | 13 |
| \(\Lambda\) | 4.1 | 0.59 | 12 |
We see that \(\textrm{p}_\textrm{v}^{\textrm{post}}\) is much more conservative than \(\textrm{p}_\textrm{v}^{\textrm{boot}}\).