Bayesian Statistics
5. Bayes Factors (BIC approximation)
Final nail in the coffin for the p-value
- Another thing often heard about the p-value is that it can only be used to reject the null hypothesis (effect = 0), not as an index of evidence in favour of the null
- This is because the p-value, by definition, is uniformly distributed under the null hypothesis
- This means that, if the null hypothesis is true, all values of the p-value are equally likely; p = .001 is just as likely as p = 1
p <- c()
for(i in 1:10000) {
x <- rnorm(100)
y <- rnorm(100)
test <- cor.test(x, y)
p <- c(p, test$p.value)
}
data.frame(p = p) |>
ggplot(aes(x=p)) +
geom_histogram(bins=40, color="white", fill="grey") +
labs(x="p-values", y="Count") +
theme_bw()
Evidence against the null hypothesis?
- Altough Frequentists recently came up with some workarounds (region of practical equivalence testing), it is still seen as a limitation of the Frequentist framework
- This is problematic in science because we are often interested in disproving a hypothesis, or proving that there is no effect
- Does Uncle Bayes have a trick in his hat?
The “Likelihood” as a known model
- Bayes’ Theorem: \(P(H|data) =P(H)~P(data|H)\)
- Bayes’ theorem relates the probability of a hypothesis (H) given some observed evidence (data X) to the probability of data given the hypothesis, as well as the prior probability of the hypothesis
- \(P(data|H)\) is the likelihood of the data given a hypothesis
- It typically corresponds to a “model of the world”
- In the previous example, we had “the model” of coin tosses. We know the likelihood of outcomes given various coin types because it corresponds to a physical reality (a fair coin gives 50% heads). So the “model” was obvious and easy to define.
- Conceptually, in this experiment:
- The likelihood was a statistical model linking a parameter x (the probability by which a coin is biased towards heads: 0.3, 0.4, 0.5, 0.6, 0.7) to the outcome (heads).
- The prior corresponded to the distribution of our credibility across the different possible x values.
- The posterior corresponded to the updated distribution of our credibility across the different possible x values after observing the data (coin tosses).
The “Likelihood” as a model with probabilistic parameters
- We often have general beliefs about probabilistic models.
- What are the odds that more type of coints exist?
- What are the odds that my model of the coin is True
- …
- In practice, the model’s likelihood depends on uncertain parameter over which we have other priors, or a combination of priors (e.g., Intercept + Slope)
- The “likelihood”, as well as the “hypothesis”, can represent different things (a statistical model, beliefs about whether a general hypothesis is True, etc.)
- Problems can be specified in different ways (it is flexible framework)
- In typical real-life statistical scenarios, the likelihood term \(P(data|H)\) is very complex to establish mathematically (i.e., it is hard to compute/approximate it)
- But it is not impossible. And accessing the likelihood is very useful…
Model 0: qsec ~ mpg
head(mtcars[c("mpg", "qsec", "wt")]) mpg qsec wt
Mazda RX4 21.0 16.46 2.620
Mazda RX4 Wag 21.0 17.02 2.875
Datsun 710 22.8 18.61 2.320
Hornet 4 Drive 21.4 19.44 3.215
Hornet Sportabout 18.7 17.02 3.440
Valiant 18.1 20.22 3.460mtcars |>
ggplot(aes(x=mpg, y=qsec)) +
geom_point() +
geom_smooth(method="lm", formula = "y ~ x") +
theme_bw()
model0 <- lm(qsec ~ mpg, data = mtcars)
summary(model0)
Call:
lm(formula = qsec ~ mpg, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-2.8161 -1.0287 0.0954 0.8623 4.7149
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.35477 1.02978 14.911 2.05e-15 ***
mpg 0.12414 0.04916 2.525 0.0171 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.65 on 30 degrees of freedom
Multiple R-squared: 0.1753, Adjusted R-squared: 0.1478
F-statistic: 6.377 on 1 and 30 DF, p-value: 0.01708mpgis a significant predictor ofqsec(\(\beta\) = 0.12, p < .05)- Is it worth it to add another predictor? Does it improve the model?
Model 1: qsec ~ mpg + wt
mtcars |>
ggplot(aes(x=wt, y=qsec)) +
geom_point() +
geom_smooth(method="lm", formula = "y ~ x") +
theme_bw()
- Do you think adding
wtwill improve the model?
model1 <- lm(qsec ~ mpg + wt, data = mtcars)
summary(model1)
Call:
lm(formula = qsec ~ mpg + wt, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-2.1092 -0.9617 -0.2343 0.8399 4.2769
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.92949 3.53431 1.961 0.0596 .
mpg 0.32039 0.09138 3.506 0.0015 **
wt 1.39324 0.56286 2.475 0.0194 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.524 on 29 degrees of freedom
Multiple R-squared: 0.3191, Adjusted R-squared: 0.2722
F-statistic: 6.797 on 2 and 29 DF, p-value: 0.003796- Be careful! A model with multiple predictors is different from individual models taken separately (because variables interact with each other)
- Is model1 better than model0?
- Higher R2, but more parameters (more complex model). Comparing BICs is an alternative that takes into account the number of parameters
Pizza plot
plot(bayestestR::bayesfactor_models(model1, denominator=model0)) +
scale_fill_manual(values = c("red", "yellow")) 
The End (for now)
Thank you!
