# How can I use Bayesian statistics to test this particular hypothesis?

There is a set $$R=\{r_1, r_2, ..., r_K\}$$ of $$N$$ ranks (where $$N>> K$$). I test the hypothesis that the ranks in $$R$$ are not homogeneously distributed in $$\{1, 2, ..., N\}$$. As I am interested in particular in the case where ranks are low (that is, the ranks from $$R$$ tend to be at the beginning of the numbers), I calculate the following statistic:

$$c_R = -2 \sum_{i=1}^{K}\ln( \frac{r_i}{N})$$

It can be easily shown that this statistic has a $$\chi^2_{2\cdot K}$$ distribution (same as proof of Fisher's method).

Now, whatever you think about the theory behind that, this test works wonders in my practical application when testing for gene set enrichments. However, similar to many other approaches in the field, it has a relatively high FP rate; clearly, the assumption of independence of ranks in $$R$$ is not fulfilled for genes. Also, gene sets vary in their number; some are very small (5-15 genes compared to 10-20 thousand of other genes), some are quite large (200-500 genes). This results in a kind of Lindley's paradox, in that even small deviations from uniform distribution for the large gene sets get very low value (it is possible to calculate an effect size and in fact this is an important number to look at when analysing the results; but that is beyond the scope of this question).

I would like to try to use a bayesian approach to calculate a posterior probability after observing $$c_R$$, using a uniform prior. Having little experience with bayesian statistics, how would I go about it?

H_0 is that the ranks in $$R$$ are uniformly distributed (better would be: uniformly distributed, but partly correlated). $$H_1$$ is that they aren't.

Say, I would like to compute

$$P(H_0|c_k) = \frac{P(c_k|H_0)\cdot P(H_0)}{P(c_k|H_0)\cdot P(H_0) + P(c_k|H_1)\cdot P(H_1)}$$

$$P(c_k|H_0)$$ is simply the p-value from the $$\chi^2_{2N}$$ distribution. I am somewhat hazy about the remaining terms. I read "Statistical Rethinking" a while ago, and I think that I roughly know what to do, but I am unable to pin the details. What would be an alternative to the uniform distribution? I guess a Poisson distribution with $$\lambda$$ varying from 1 to $$N$$. How should $$\lambda$$ be distributed? Heck if I know. How should I calculate $$P(H_0)$$?

I think that once I can sample from $$H_1$$, I can estimate the $$P(c_t|H_1)$$ similarly using random sampling.

Actually, it would be best if I could introduce a measure for how tightly the genes are correlated and sample from that, but I'm not sure how to do it.

Rather than using a software package like (R)stan, I would prefer to be able to get at a solution myself in order to understand the process.

So, my question: how do I in this setup calculate $$P(H_0)$$? How should I sample the $$H_1$$?

• I don't know if I can help you. I started by simulating $c_R$ using $N=100$ and $K=10$. I made $10^6$ draws. The mean was about 19.4 and the variance was about 34.1. This is not a chi-square distribution where the variance is twice the mean and the mean equals the degrees of freedom. Am I doing something wrong?
– mef
Oct 21 '20 at 12:05
• Mean of about 20 is right. As for the variance, this is caused by drawing ranks without replacement. So the statement that ranks follow a uniform distribution is not true if you draw without replacement (no tied ranks), but is asymptotically true when k/n -> 0. Which is, fortunately, the case for me. Try k=50 and n=20000. Oct 23 '20 at 7:21
• So the degrees of freedom are more like $2K$ than $2N$.
– mef
Oct 23 '20 at 16:04
• Damn, my mistake. Of course 2K. Oct 26 '20 at 14:54

## 1 Answer

I think it might be helpful to recast the problem (and change the notation slightly) in order to focus on the Bayesian framework and avoid confusion with the frequentist framework.

Suppose there are two possible models for how an observation is generated: $$M_1$$ and $$M_2$$. (There is no notion of null and alternative: there are just two models.) Given an observation $$y$$, what is the probability it was generated by $$M_1$$? The answer is provided by Bayes rule: $$$$p(M_1|y) = \frac{p(y|M_1)\,p(M_1)}{p(y|M_1)\,p(M_1) + p(y|M_2)\,p(M_2)} .$$$$ Note that $$p(y|M_i)$$ is the value of a density; it is not a tail probability such a a $$p$$-value. So if $$y \sim \chi_{2K}^2$$ according to $$M_1$$, then $$$$p(y|M_1) = \textsf{Chi-square}(y|2K) = \frac{1}{2^K\,\Gamma(K)}\,y^{K-1}\,e^{-y/2} .$$$$ The analyst must specify the two remaining components: the prior model probabilities $$p(M_i)$$ and the distribution for the observation given the other model $$p(y|M_2)$$.

The prior model probabilities $$p(M_i)$$ cannot (and must not) be computed from the current observation $$y$$. They are derived from "non-sample" information, which includes what has be learned from other experiments. You must say what they are.

You must also specify $$M_2$$. In the Bayesian framework, it takes a model to beat a model. It is not sufficient to say "$$M_2$$ is not $$M_1$$". It is possible that $$M_2$$ is quite complicated, being (for example) and average of a number of sub-models. That's fine; but in the end you must be able to say what the density $$p(y|M_2)$$ is. You might simulate $$y$$ according to $$M_2$$ and use a kernel density approach to approximate $$p(y|M_2)$$. One way or another you must come up with $$M_2$$ and $$p(y|M_2)$$.

One final comment: The posterior odds ratio can be expressed as the Bayes factor (BF) times the prior odds ratio: $$$$\frac{p(M_1|y)}{p(M_2|y)} = \underbrace{\frac{p(y|M_1)}{p(y|M_2)}}_{\text{BF}} \times \frac{p(M_1)}{p(M_2)} .$$$$ This expression shows the contributions of the sample information and non-sample information. The Bayes factor is the ratio of two densities. You can compute it without saying what the prior model probabilities are.