# Bayesian modeling and FDR correction

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.

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)

• Thanks for the answer, I will read the mentioned article and get back to you. In the meantime could you please explain what you don't understand in my question? Oct 5, 2015 at 13:50
• @yannick. I understand your question. I just think that my answer shifts a bit your intial main question (e.g. seek a combination) but to my knowledge (which is limited) what my answer answers your intial motivation Oct 5, 2015 at 14:47
• I aligned your equations for better readability, but if you don't like it feel free to rollback or reject the edit. Oct 6, 2015 at 13:48
• Thanks, it was indeed a very good read! Could you simply add one more sentence to explain what you mean with "so that the margin posterior for the λi (and the resulting bayes factor) are no longer independant each other" ? Specifically, Gelman does not mention Bayes factors in this paper, and does not like them, so I guess you must be having something in mind. Oct 12, 2015 at 14:24
• @yannick Happy if it helps. In substance, I think that there is no limitation in using such a pooled model for Bayes factor or hpd testing (as in Gelman paper). The two strategies can be applied and are to my knowledge, two different answers to a same problematic. To my knowledge, Bayes factor is prefered by most authors. In any cases, the two strategies rely on writing the posterior which is the element I suggest to modify by including a hierarchical model for the lambda_i. I updated my answer. Feel free to ask if not clear enough Oct 13, 2015 at 13:36