Bayesian Testing of Equal Genotype Proportions between Multiple Populations

Multinomial Model

In the two-population problem we have K = 2 lakes with N₁ and N₂ fish sampled from each. Let Y_k = (Y_k1, …Y_kC) denote the counts for each genotype observed in lake k, such that $\sum_{i=1}^C Y_{ki} = N_k$. The sampling model for these data is $$ \boldsymbol Y_k \stackrel{\textrm{ind}}{\sim}\textrm{Multinomial}(N_k, \boldsymbol \rho_k), \quad k = 1,2, $$ where ρ_k = (ρ_k1, ..., ρ_kC) are the population proportions of each genotype, and $\sum_{i=1}^C \rho_{ki} = 1$.

Hypothesis Testing

Our objective is to test $$ \begin{split} H_0 &: \textrm{The two populations have the same genotype proportions} \\ & \phantom{: } \iff \boldsymbol \rho_1 = \boldsymbol \rho_2. \end{split} $$ The classical test statistics for assessing H₀ are Pearson’s Chi-Square statistic 𝒳 and the Likelihood Ratio statistic Λ, $$ \mathcal X = \sum_{k=1}^2 \sum_{i=1}^C \frac{(N_k\hat\rho_i - Y_{ki})^2}{N_k\hat\rho_i}, \qquad \Lambda = 2 \sum_{k=1}^2 \sum_{i=1}^C Y_{ki} \log\left(\frac{Y_{ki}}{N_k\hat\rho_i}\right), \qquad \hat \rho_i = \frac{Y_{1i} + Y_{2i}}{N_1 + N_2}. $$ Under H₀, the asymptotic distribution of either of these test statistics T = 𝒳 or Λ is χ_{(C − 1)}², such that the p-value p_v = Pr(T > T_obs ∣ H₀) for an observed value of T_obs can be estimated as follows:

# observed values of the test statistics
chi2.obs <- chi2.stat(ctab) # Pearson's chi^2
LRT.obs <- LRT.stat(ctab) # LR test statistic
T.obs <- c(chi2 = chi2.obs, LRT = LRT.obs)
# p-value with asymptotic calculation
C <- ncol(ctab)
pv.asy <- pchisq(q = T.obs, df = C-1, lower.tail = FALSE)
signif(pv.asy, 2)

##  chi2   LRT 
## 0.330 0.041

The Chi-Square and LR tests are asymptotically equivalent and so should give roughly the same p-values. The huge discrepancy observed here indicates that the sample sizes are too small for asymptotics to kick in. A more reliable p-value estimate can be obtained by the Bootstrap method, which in this case consists of simulating M contigency tables with $\boldsymbol Y_k^{(m)} \stackrel{\textrm{ind}}{\sim}\textrm{Multinomial}(N_k, \hat{\boldsymbol \rho})$, m = 1, …, M, where $\hat{\boldsymbol \rho}$ is the estimate of the common probability vector ρ₁ = ρ₂ under H₀. For each contingency table (Y₁^(m), Y₂^(m)), we calculate the test statistic T^(m), and the bootstrapped p-value is defined as $$ \textrm{p}_\textrm{v}^{\textrm{boot}}= \frac 1 M \sum_{m=1}^M \delta(T^{(m)} \ge T_\textrm{obs}). $$ p_v^boot is calculated with MADPop as follows:

N1 <- sum(ctab[1,])                     # size of first sample
N2 <- sum(ctab[2,])                     # size of second sample
rho.hat <- colSums(ctab)/(N1+N2)        # common probability vector
# bootstrap distribution of the test statistics
# set verbose = TRUE for progress output
system.time({
  T.boot <- UM.eqtest(N1 = N1, N2 = N2, p0 = rho.hat, nreps = 1e4,
                      verbose = FALSE)
})

##    user  system elapsed 
##   0.459   0.008   0.467

# bootstrap p-value
pv.boot <- rowMeans(t(T.boot) >= T.obs)
signif(pv.boot, 2)

##   chi2    LRT 
## 0.0068 0.0060

Note that the bootstrap p-values for both tests are roughly the same and decisively reject H₀, whereas the less reliable asymptotic p-values both failed to reject (at quite different significance levels).

Parameter Estimation

The likelihood function for this model corresponds to a Dirichlet-Multinomial distribution, $$ \begin{split} \mathcal L(\boldsymbol \alpha\mid \boldsymbol Y) = \prod_{k=1}^K p(\boldsymbol Y_k \mid \boldsymbol \alpha) & = \prod_{k=1}^K \int p(\boldsymbol Y_k \mid \boldsymbol \rho_k) p(\boldsymbol \rho_k \mid \boldsymbol \alpha) \,\mathrm{d}\boldsymbol \rho_k \\ & = \prod_{k=1}^K \left[\frac{N_k!\cdot\Gamma(\alpha_0)}{\Gamma(N_k+\alpha_0)} \prod_{i=1}^C \frac{\Gamma(Y_{ki} + \alpha_i)}{Y_{ki}!\cdot\Gamma(\alpha_i)}\right], \end{split} $$ from which α can be estimated by maximum likelihood (Minka 2000). Alternatively we estimate α using a Bayesian approach, which is more readily extensible to more complex genotype models, for which ℒ(α ∣ Y) has no closed form (see Future Work).

First, we specify a prior distribution $$ \pi(\alpha_0, \bar{\boldsymbol \alpha}) = \frac{1}{(1+ \alpha_0)^2}, $$ which is a uniform distribution on the prior probability vector $\bar{\boldsymbol \alpha}$, and a uniform distribution on the variance-like parameter A = (1 + α₀)⁻¹, in the sense that the prior mean and variance of a given genotype population probability are E[ρ_kj] = α̃_j,   var(ρ_kj) = A ⋅ α̃_j(1 − α̃_j). Then, we sample from the posterior distribution p(α ∣ Y) ∝ ℒ(α ∣ Y)π(α) using a Markov chain Monte Carlo (MCMC) algorithm provided by the Stan programming language (Stan Development Team 2016), as detailed below.

MCMC Sampling from the Posterior Distribution

In order to estimate α under H₀, we apply the Dirichlet-Multinomial distribution to K − 1 lakes, with the first two being collapsed into a common count vector Y₁₂ = Y₁ + Y₂ with sample size N₁₂ = N₁ + N₂. MCMC sampling is implemented via a Hybrid Monte Carlo (HMC) variant provided by the R package rstan. In MADPop, sampling from p(α ∣ Y, H₀) 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₀:

popId                                   # equal lakes under H0

## [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            19

Next, we sample from p(α ∣ Y, H₀) 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 0.000173 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.73 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: 1.732 seconds (Warm-up)
## Chain 1:                15.58 seconds (Sampling)
## Chain 1:                17.312 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(α, ρ₁₂ ∣ Y, H₀). Boxplots of the 50 highest posterior probability genotypes in the common distribution ρ₁₂ 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, …, M. However, for each contingency table, the common probability vector ρ₁₂^(m) is different, corresponding to a random draw from p(ρ₁₂ ∣ Y, H₀). 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}}). $$ p_v^post is calculated with MADPop as follows:

system.time({
  T.post <- UM.eqtest(N1 = N1, N2 = N2, p0 = rho.post, nreps = 1e4,
                      verbose = FALSE)
})

##    user  system elapsed 
##   0.889   0.012   0.902

# posterior p-value
pv.post <- rowMeans(t(T.post) >= T.obs)

We can now compare the three types of p-values (asymptotic, bootstrapped, posterior) for each test statistic (𝒳 and Λ):

We see that p_v^post is much more conservative than p_v^boot.