49. Bayes factors

A p-value can only ever produce evidence against a null. It has no mechanism for supporting one hypothesis over another, and p = 0.6 tells you nothing about whether H_0 is true (§5).

The Bayesian answer to model comparison fixes both problems.

The definition

Compare two hypotheses by how well each predicted the data:

BF_{10} = \frac{p(\text{data} \mid H_1)}{p(\text{data} \mid H_0)}

The Bayes factor is the ratio of marginal likelihoods — the probability each hypothesis assigned to the data actually observed.

Its role becomes clear in the odds form of Bayes' theorem:

\underbrace{\frac{P(H_1 \mid \text{data})}{P(H_0 \mid \text{data})}}_{\text{posterior odds}} = \underbrace{\frac{p(\text{data} \mid H_1)}{p(\text{data} \mid H_0)}}_{\text{Bayes factor}} \times \underbrace{\frac{P(H_1)}{P(H_0)}}_{\text{prior odds}}

The Bayes factor is exactly the factor by which the data updates your odds. It is the evidence, cleanly separated from your prior beliefs — which is why it can be reported to an audience that disagrees with you about the prior.

Reading the number

Jeffreys' conventional scale:

BF_{10}	Evidence for H_1
1–3	barely worth mentioning
3–10	substantial
10–30	strong
30–100	very strong
> 100	decisive

And symmetrically: BF_{10} = 0.1 (i.e. BF_{01} = 10) is strong evidence for H_0.

That symmetry is the key advantage over p-values. A Bayes factor can come out clearly below 1, which is positive support for the null — something no p-value can express. "We found evidence that the effect is absent" becomes a statable conclusion rather than a fallacy.

The boundaries are conventions, not laws, and inherit all the arbitrariness of \alpha = 0.05. Report the number.

The marginal likelihood

For a composite hypothesis with parameters, p(\text{data} \mid H) is not a single likelihood — it's an average over the prior:

p(\text{data} \mid H) = \int p(\text{data} \mid \theta, H)\,p(\theta \mid H)\,d\theta

This is where the automatic Occam's razor comes from.

A complex model spreads its prior over many parameter values. It can fit a wide range of possible datasets, so it must assign only a little probability to each — including the one you actually saw.

A simple model concentrates its prior. If the data lands where it predicted, it assigns high probability there and wins the comparison.

Complexity is penalised automatically, with no explicit penalty term. A model that can explain anything explains nothing in particular.

The critical caveat: because the marginal likelihood integrates over the prior, the Bayes factor depends on the prior even with infinite data. This is genuinely different from posterior estimation, where the prior washes out. A vague prior on \theta under H_1 spreads the mass thinly and drives BF_{10} toward 0 — the Jeffreys–Lindley paradox. You cannot use an improper or arbitrarily diffuse prior here and expect a meaningful answer.

Worked example

A coin: 60 heads in 100 flips. Compare H_0: p = 0.5 against H_1: p is unknown, with a uniform prior.

Under H_0, the parameter is fixed, so the marginal likelihood is just the Binomial probability:

p(\text{data} \mid H_0) = \binom{100}{60}(0.5)^{100} \approx 0.0108

Under H_1, average over the uniform prior. The integral is a Beta function:

p(\text{data} \mid H_1) = \int_0^1 \binom{100}{60}p^{60}(1-p)^{40}\,dp = \binom{100}{60}B(61, 41) = \frac{1}{101} \approx 0.0099

(A pleasant fact: with a uniform prior the marginal likelihood is 1/(n+1) regardless of k — every possible outcome count is equally likely a priori.)

BF_{10} = \frac{0.0099}{0.0108} \approx 0.91

Slightly favouring H_0 — barely worth mentioning either way.

Now compare with the p-value. A two-sided test gives z = 2 and p \approx 0.046: "significant", reject the null.

The two disagree, and the disagreement is instructive. The p-value says the data is somewhat unusual under fairness. The Bayes factor asks a different question — did the alternative predict this data better? — and the answer is no, because H_1 had to spread its bets across all values of p from 0 to 1, while H_0 put everything on 0.5.

This is the Jeffreys–Lindley effect in miniature. It's also why p \approx 0.05 results replicate so poorly: they are much weaker evidence than the threshold suggests.

A fairer alternative. If we use a more focused H_1 — say p \sim \text{Beta}(20, 20), expressing "if biased, probably mildly so" — the comparison changes, because the alternative is no longer wasting prior mass on implausible values like p = 0.99. The code below runs this.

Doing it in Python

The coin comparison, exactly:

import numpy as np
from scipy.stats import binom, beta as beta_dist
from scipy.special import betaln

k, n = 60, 100

# H0: p = 0.5 exactly
m0 = binom.pmf(k, n, 0.5)

# H1: p ~ Beta(a, b), integrated out analytically
def marginal(k, n, a, b):
    log_m = (betaln(a + k, b + n - k) - betaln(a, b)
             + np.log(binom.pmf(k, n, 0.5)) - k*np.log(0.5) - (n-k)*np.log(0.5))
    return np.exp(log_m)

m1_flat = marginal(k, n, 1, 1)
print(f"p(data | H0)          = {m0:.6f}")
print(f"p(data | H1, uniform) = {m1_flat:.6f}   (= 1/(n+1) = {1/(n+1):.6f})")
print(f"\nBF10 = {m1_flat/m0:.4f}   -> mild support for H0")

How the choice of alternative changes the verdict:

import numpy as np
from scipy.stats import binom
from scipy.special import betaln

k, n = 60, 100
m0 = binom.pmf(k, n, 0.5)

def bayes_factor(k, n, a, b):
    log_bf = (betaln(a + k, b + n - k) - betaln(a, b)
              - (k * np.log(0.5) + (n - k) * np.log(0.5)))
    return np.exp(log_bf)

print(f"{'alternative prior':>26} {'BF10':>10} {'favours':>10}")
for name, (a, b) in [
    ("Beta(1,1)   very vague", (1, 1)),
    ("Beta(5,5)   broad", (5, 5)),
    ("Beta(20,20) focused", (20, 20)),
    ("Beta(50,50) tight", (50, 50)),
    ("Beta(0.5,0.5) Jeffreys", (0.5, 0.5)),
]:
    bf = bayes_factor(k, n, a, b)
    print(f"{name:>26} {bf:>10.4f} {'H1' if bf > 1 else 'H0':>10}")

print("\nThe vaguer the alternative, the more it is penalised -- Occam's razor,")
print("and the Jeffreys-Lindley effect, in one table.")

Bayes factor against p-value across sample sizes, holding the proportion fixed:

import numpy as np
from scipy.stats import binom, norm
from scipy.special import betaln

def bayes_factor(k, n, a=1, b=1):
    log_bf = (betaln(a + k, b + n - k) - betaln(a, b)
              - n * np.log(0.5))
    return np.exp(log_bf)

print(f"{'n':>8} {'k':>8} {'p-value':>10} {'BF10':>12} {'agree?':>10}")
for n in (100, 400, 1_600, 10_000):
    k = int(0.6 * n)
    z = (k - n/2) / np.sqrt(n * 0.25)
    pval = 2 * (1 - norm.cdf(abs(z)))
    bf = bayes_factor(k, n)
    agree = ("yes" if (pval < 0.05) == (bf > 3) else "NO")
    print(f"{n:>8} {k:>8} {pval:>10.2e} {bf:>12.3e} {agree:>10}")

print("\nWith a real 10-point effect, both eventually agree overwhelmingly.")

And the case where they disagree — a tiny effect with a huge sample:

import numpy as np
from scipy.stats import norm
from scipy.special import betaln

def bayes_factor(k, n, a=1, b=1):
    return np.exp(betaln(a + k, b + n - k) - betaln(a, b) - n * np.log(0.5))

print("A barely-there effect (p = 0.51) at increasing sample size:\n")
print(f"{'n':>10} {'p-value':>10} {'BF10':>10} {'verdict':>28}")
for n in (1_000, 10_000, 100_000, 1_000_000):
    k = int(round(0.51 * n))
    z = (k - n/2) / np.sqrt(n * 0.25)
    pval = 2 * (1 - norm.cdf(abs(z)))
    bf = bayes_factor(k, n)
    verdict = f"p {'sig' if pval < 0.05 else 'ns':>3}, BF favours {'H1' if bf > 1 else 'H0'}"
    print(f"{n:>10} {pval:>10.4f} {bf:>10.3f} {verdict:>28}")

print("\nThe p-value chases significance as n grows. The Bayes factor asks")
print("whether the alternative actually predicted better -- a different question.")

Your turn

1. BF_{10} = 20. What does that mean, and what do you need to add to get the posterior odds?

2. Why can't a p-value give evidence for H_0?

3. Why does a very vague prior under H_1 push BF_{10} toward zero?

Solutions

1. The data is 20 times more likely under H_1 than under H_0 — strong evidence on Jeffreys' scale.

To get posterior odds you need the prior odds:

\text{posterior odds} = 20 \times \frac{P(H_1)}{P(H_0)}

If you started at even odds (1:1), you end at 20:1, i.e. P(H_1 \mid \text{data}) = 20/21 \approx 0.952.

But if H_1 was a wild claim you'd have given 1:1000 beforehand, the posterior odds are 20/1000 = 0.02 — still 50:1 against. Strong evidence isn't enough to rescue an implausible hypothesis, which is exactly the base-rate lesson from §0 in a new setting.

2. Because of what a p-value measures: $P(\text{data at least this extreme} \mid H_0)$. It's computed entirely under H_0 and never consults the alternative.

A large p-value means the data is unsurprising under H_0. But it's very often equally unsurprising under a small-effect alternative — a study with low power produces large p-values whether or not an effect exists. The statistic cannot distinguish "no effect" from "not enough data to see one".

A Bayes factor compares both hypotheses explicitly, so BF_{10} = 0.05 is a meaningful statement that H_0 predicted the data 20 times better.

The frequentist tools for this purpose are equivalence testing (TOST) or a confidence interval narrow enough to exclude every effect you'd care about — both of which require you to state what "no effect" means quantitatively.

3. Because the marginal likelihood is an average over the prior:

p(\text{data} \mid H_1) = \int p(\text{data} \mid \theta)\,p(\theta)\,d\theta

A very diffuse prior spreads its mass over a huge range of \theta, the vast majority of which predict the observed data terribly. Those near-zero likelihoods are included in the average and drag it down. The good fit near \hat\theta gets weighted by the small prior mass sitting there.

Push this to the limit — a prior spread over (-\infty, \infty) — and p(\text{data} \mid H_1) \to 0, so BF_{10} \to 0 and the null wins no matter what the data says. That's the Jeffreys–Lindley paradox.

Two practical consequences:

• Improper priors cannot be used for Bayes factors, even though they're often fine for estimation. There's no normalising constant, so the ratio is undefined.
• The prior under H_1 must represent a real belief about the effect's plausible size. That's a feature — it forces you to say what effect you're actually testing for, rather than "anything other than zero".

Check yourself in code

Compute the Bayes factor for 60 heads in 100 flips under two different alternatives and compare with the p-value.

Print exactly this:

BF10 vague 0.913
BF10 focused 2.2546
p-value 0.0455
vague favours H0: True

Use a Beta(1,1) prior for the vague alternative and Beta(20,20) for the focused one, a two-sided Normal-approximation p-value, and round every value to 4 decimal places.

import numpy as np
from scipy.stats import norm
from scipy.special import betaln

k, n = 60, 100

def bayes_factor(k, n, a, b):
    return np.exp(betaln(a + k, b + n - k) - betaln(a, b) - n * np.log(0.5))

print("BF10 vague", round(float(bayes_factor(k, n, 1, 1)), 4))

# Print the focused Beta(20,20) Bayes factor, the two-sided p-value from the
# normal approximation, and whether the vague alternative favours H0 (BF < 1).

import numpy as np
from scipy.stats import norm
from scipy.special import betaln

k, n = 60, 100

def bayes_factor(k, n, a, b):
    return np.exp(betaln(a + k, b + n - k) - betaln(a, b) - n * np.log(0.5))

vague = float(bayes_factor(k, n, 1, 1))
print("BF10 vague", round(vague, 4))
print("BF10 focused", round(float(bayes_factor(k, n, 20, 20)), 4))

z = (k - n / 2) / np.sqrt(n * 0.25)
print("p-value", round(float(2 * (1 - norm.cdf(abs(z)))), 4))
print("vague favours H0:", bool(vague < 1))

---

A Bayes factor is the ratio of marginal likelihoods — how much the data shifts your odds between two hypotheses. Unlike a p-value it is symmetric, so it can support a null as readily as reject it, and it penalises complexity automatically because a model that spreads its predictions thinly assigns little probability to any one outcome. The price is that the prior under H_1 matters permanently, so it must be chosen to mean something.

That closes §7. Next: processes that evolve in time.