Bayesian Statistics
5. Bayes Factors
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?
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 taking separately (because variable interact with each other)
- Is model1 better than model0?
- Higher R2, but more parameters (more complex model). Comparing BICs is an alterantive that takes into account the number of parameters
Pizza plot
plot(bayestestR::bayesfactor_models(model1, denominator=model0)) +
scale_fill_manual(values = c("red", "yellow")) 
t-Tests
- True Bayes factor (i.e., not approximations using Frequentist models) can be also computed for “simple” tests, such as t-tests
- t-tests test for the difference of means of a variable between two groups (or against a constant value)
- This is a frequentist t-test:
# `am` is a binary variable (0 and 1)
t.test(mtcars$mpg ~ mtcars$am)
Welch Two Sample t-test
data: mtcars$mpg by mtcars$am
t = -3.7671, df = 18.332, p-value = 0.001374
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-11.280194 -3.209684
sample estimates:
mean in group 0 mean in group 1
17.14737 24.39231 mtcars |>
ggplot(aes(x=as.factor(am), y=mpg)) +
geom_jitter(width=0.1, size=3) 
- We conclude that there is a significance difference between the two groups (am == 1 > am == 0, p < .01)
Prior Distribution
- In the previous examples, we had priors in the form of discrete categories (e.g., 1/3 that the coin is fair, 1/3 that it’s rigged, …)
- But how to specify a prior for a continuous variable like a correlation coefficient?
- Using a distribution!
- Let’s imagine taking a normal distribution as a prior for \(\rho\) (i.e., we think that \(\rho\) is normally distributed)
- Imagine that someone proposes the following parameters
- e.g., \(\mu = 0.42, \sigma = 0.7\)
xspace <- seq(-1.5, 1.5, by=0.01)
prob <- dnorm(xspace, mean=0.42, sd=0.333)
data <- data.frame(x = xspace, y = prob)
data |>
ggplot(aes(x=x, y=y)) +
geom_line(color="blue", linewidth=3) 
Exercice
- Do you think this is a good prior?
- Why is this not a very good prior?
- Try to come up with a suitable prior for the correlation coefficient \(\rho\)
Solution
- Priors represent your beliefs before seeing the data
- In can be from previous studies, from your own experience, from your intuition, …
- Do we have a hypothesis about the direction (positive vs. negative) of the correlation between
dratandqsec- No. The correlation could be positive or negative.
- Do we have reasons to expect a non-null correlation?
- No.
- Our prior could be symmetric around 0, with a higher probability on 0 (i.e., we think that the correlation is more likely to be 0 than any other value)
- What about \(Normal(0, 0.5)\)
data$y <- dnorm(xspace, mean=0, sd=0.5)
data |>
ggplot(aes(x=x, y=y)) +
geom_line(color="blue", linewidth=3) 
- What’s the problem?
What about Uniform?
- \(Uniform(-1, 1)\) means that all values between -1 and 1 are equally likely
# Compute the uniform distribution
data$y <- dunif(xspace, -1, 1)
data |>
ggplot(aes(y=y)) +
geom_line(aes(x=x), color="blue", linewidth=3) +
theme_bw()
- This is a non-informative prior
- i.e., we don’t introduce any prior knowledge about the correlation coefficient (other than domain knowledge about what a correlation coefficient is)
- As we will see later, uniform priors tend to give results that are similar to frequentist maximum likelihood estimates (i.e., Bayesian models with uniform priors tend to have very similar results to frequentist models)
- But they are non-recommended.
- They are often non realistic: we tend to have some expectations about our parameters, and priors should reflect these beliefs
What other prior?
- What is another distribution that has bounds?
- The beta distribution!
- But what’s the problem?
- It is defined on the interval [0, 1]
- Solution?
- Distributions can often be scaled and shifted
Shifted and Scaled Beta
- We can have a shifted & scaled beta distribution that is defined on the interval [-1, 1] with a peak at 0
- E.g., How about \(Beta(1.5, 1.5)\) rescaled and shifted to [-1, 1]?
# Compute the beta distribution
data$y <- dbeta(xspace, 1.5, 1.5)
# Rescale its x-axis
data$x2 <- (data$x - 0.5) * 2
data |>
ggplot(aes(y=y)) +
geom_line(aes(x=x), color="blue", linewidth=3) +
geom_line(aes(x=x2), color="orange", linewidth=2, alpha=0.9) +
coord_cartesian(xlim=c(-1.2, 1.2)) +
theme_bw()
- Original distribution in blue and shifted and scaled in orange
Meaning
- This means that we think that the correlation coefficient is more likely to be 0, and then “gradually” less likely as we move away from 0, with very low likelihood for values close to -1 or 1
Visualize BayesFactor’s defaults
# List of rscales
rscales <- list(
"medium.narrow" = 1/sqrt(27),
"medium" = 1/3,
"wide" = 1/sqrt(3),
"ultrawide" = 1
)# Compute the beta distribution following the formula
# And store in the same dataframe
xspace <- seq(0, 1, by=0.01) # Original beta is ]0, 1[
data <- data.frame() # Inititalize empty dataframe
for(scale in names(rscales)) {
rscale <- rscales[[scale]]
dat <- data.frame(
y = dbeta(xspace, 1/rscale, 1/rscale),
x = (xspace - 0.5) * 2, # Rescale its x-axis
scale = scale
)
data <- rbind(data, dat)
}- Exercice: plot that!
data |>
ggplot(aes(x=x, y=y)) +
geom_line(aes(color=scale), linewidth=2) +
coord_cartesian(xlim=c(-1.2, 1.2)) 
The End (for now)
Thank you!
