Bayesian modeling and FDR correction

Question

Question

Assume that we have 2 experiments, each of which yields $N$ observations. We label them $(c_1^i,\ldots,c_N^i)$ with $i=1,2$. Each observation is a positive integer, and can be modeled as a Poisson variate. We wish to know where the two experiments disagree.

Assume further that in the classical context, the appropriate analysis for this experiment would be to

1. Perform a rate ratio test on all pairs of observations, asking whether $c_i^1$ and $c_i^2$ have the same Poisson rate ($\mathcal{H}_0$) or not ($\mathcal{H}_1$)
2. Perform a FDR (Benjamin-Hochberg) correction on the $p$-values to account for multiple testing
3. Select those pairs with $q<0.05$

I would like to know what is the Bayesian take on this, in particular at step 2.

Ideas

• I have pursued a model comparison approach, but am stuck at step 2

First, I write the probability for two observations and one rate, model $M_1$: $$p(M_1|c_i^1,c_i^2) = \int \text{d}\lambda p(c_i^1,c_i^2|\lambda)p(\lambda)$$ Or two rates, model $M_2$: $$p(M_2|c_i^1,c_i^2) = \int \text{d}\lambda_1\text{d}\lambda_2 p(c_i^1|\lambda_1)p(c_i^2|\lambda_2)p(\lambda_1)p(\lambda_2)$$ I use a proper prior for $p(\lambda)$, which I don't want to discuss here.

Then, I compute the Bayes factor for these two models: $K_i = p(M_2|c_i^j)/p(M_1|c_i^j)$ assuming equal prior odds. The Bayes factor for all observations then factorizes into the individual pairs: $K=K_1\cdots K_N$

Now, the big question, is how do we assign most probable models to each pair?

a) Maximize K, e.g. give 1 rate to points for which $K_i<1$ and 2 rates to the others. To me, this sounds like no $p$-value correction for multiple testing, i.e. a bad idea.

b) Seek a combination to reach, say, $K\simeq 3$. But how do I do that? That sounds much closer to FDR to me.

• Write down the posterior probability of one rate being larger than the other, and use this probability in a setting similar to the definition of a $p$-value. I did not try this yet but thought I could write it here.

Answer 1

I am not sure to fully answer your question but think that the following can help.

What you can try is to a design multilevel model including a partial pooling of the $\lambda_i$ by considering for example: $$ \begin{align} c_i &\sim P(\lambda_i) \\ λ_i &\sim Gamma(a,b) \\ a &\sim \text{Find a good prior here} \\ b &\sim findTheGoodPriorHere \end{align}$$

Or another possibility: $$\begin{align} c_i &\sim P(\lambda_i). \\ log(\lambda_i) &=\alpha + \theta_i \\ \mbox{with } \theta_i &\sim N(0,\sigma^2) \\ \alpha &\sim u_{R} \\ \sigma &\sim gamma(\epsilon,\epsilon) \end{align}$$

or any other suitable pooling, so that the margin posterior for the $\lambda_i$: $$ p(\lambda_i|(c_j)_{j=1:N}) \propto \int_{R}[ p(c_i|\lambda_i) p(\lambda_i|\beta) d\beta ] \cdot \prod_{j\ne i} [\int_{R^+} \int_{R} p(c_j|\lambda_j) p(\lambda_j|\beta) d\beta d\lambda_j] $$ (calling for generality $\beta$ the hyperparameter of $p(\lambda_j|\beta)$)

are no longer independant each other (because the hyperparameter $\beta$ is common to all the $\lambda_i$) and does no more result in independent inferences/comparisons (while this needs a dedicated discussion e.g. Why don't Bayesian methods require multiple testing corrections?). The important question is weither or not the partial pooling model is suitable for your design.

Here (http://www.stat.columbia.edu/~gelman/research/unpublished/multiple2.pdf) is a reference for such a consideration by Gelman (it does not use Bayes factor but to my knowledge, there is no limitation in using such a pooling model with Bayes factor)