How to compare two models of binomial distributions?

Question

Assume we have two series of indepedent success-failure observations, e.g. coin tosses $$ \boldsymbol x_1 \in \{H,T\}^{n_1} \\ \boldsymbol x_2 \in \{H,T\}^{n_2} $$

Also, let $k_i$ be the number of heads (=successes) in $\boldsymbol x_i$.

Now, I want to assume the observations are generated by two random variables $X_i$ drawn from a binomial distribution:

$$ X_i \sim B(n_i, p_i) $$

Goal (intuitively)

Basically, what I want to do is compare two different scenarios:

1. The observations have been generated by two different coins
2. The observations have been generated by the same coin

For the comparison I want to have a value $\Lambda \in [0, 1]$ that expresses the difference of the two models.

Goal (formally)

I want to estimate the model parameter $p_i$ and compare two different models $\boldsymbol \theta^{(j)}, j=1,2$ of the underlying binomial distribution for the joint probability of the observations, i.e.

$$ L(\boldsymbol\theta^{(j)}) = P(X_1=k_1, X_2=k_2; \boldsymbol\theta^{(j)}) = P(X_1=k_1; \theta_1^{(j)}) \cdot P(X_2 = k_2; \theta_2^{(j)}) $$

So, for the first scenario $j=1$, i.e. when the observations have been generated by two different coins, I assume the MLE for each binomial distribution:

$$ \theta_i^{(1)} = k_i/n_i $$

For the second scenario $j=2$, I assume that the underlying binomial distributions have the same parameter, namely:

$$ \theta_1^{(2)} = \theta_2^{(2)} = \frac{k_1 + k_2}{n_1 + n_2} $$

Approach (1)

In order to compare the two models, I could just use the likelihood ratio

$$ \Lambda := \frac{L(\boldsymbol\theta^{(2)})}{L(\boldsymbol\theta^{(1)})} $$

We know for sure that the first model is more likely than the second, since it has been estimated with the MLEs. This means that $\Lambda \in [0, 1]$.

Approach (2)

A colleague of mine insists that we can actually calculate the probability of the models, by using Bayes:

$$P(\boldsymbol \theta^{(j)}|X_1=k_1, X_2=k_2) = \frac{P(X_1=k_1,X_2=k_2|\boldsymbol\theta^{(j)})\cdot P(\boldsymbol \theta^{(j)})}{\sum_{l=1}^2 P(X_1=k_1,X_2=k_2|\boldsymbol\theta^{(l)})\cdot P(\boldsymbol \theta^{(l)})}$$

Since we don't know the priors $P(\boldsymbol\theta^{(j)})$ he suggests

$$ P(\boldsymbol\theta^{(1)}) = P(\boldsymbol\theta^{(2)}) = 0.5 $$

He suggests to use

$$ \Lambda := P(\boldsymbol \theta^{(2)}|X_1=k_1, X_2=k_2) \in [0,1] $$ Question

Is the second approach correct? If not, what assumption would I have to make in order to legitimate it?

In other words. Is it okay to constraint the entire hypotheses space to our two hypotheses $\Theta^C = \{\boldsymbol\theta^{(1)}, \boldsymbol\theta^{(2)}\}$?

Let me know what you think. Thank you in advance!

Answer 1

I agree with your colleague. Fitting MLE's to each experiment will always yield a likelihood greater than or equal to the model with a common probability of heads, so you need to account for this somehow.

You are comparing models with different sized parameter spaces. You have 2 parameters to fit in the first one, and only 1 in the second. Therefore, you need to account for the inevitable increase in likelihood gained by the extra degrees of freedom. There are a few common approaches:

1. AIC
2. BIC
3. Penalized likelihood (in general)
4. Bayesian inference with a prior over the two options.

Of course, what if $p_1 = p_2 + \epsilon, \; |\epsilon|\ll p_2$

You are ruling out this possibility in your approach, but you may not care. Whether or not it is "correct" depends on how sensitive you want your test to be.

The nice thing about (4) is you automatically get $\Lambda$ due to the fact that Bayesian inference works with probabilities.