library(brms)
f <- brms::brmsformula(qsec ~ 0 + Intercept + mpg)Bayesian Statistics
9. Linear Models (Posterior)
MCMC - NUTS
- By default,
brmsuses MCMC- In particular, the No-U-Turn Sampler (NUTS) algorithm, which is a self-adjusting variant of Hamiltonian Monte Carlo (HMC)
- The key idea behind HMC is that it treats sampling as a physics problem, where we imagine the posterior distribution space as a hilly landscape and use momentum to explore it efficiently
- Instead of setting a fixed number of steps, NUTS keeps “rolling the ball” until it detects a U-turn (where the trajectory starts going back toward its starting point)
MCMC vs. Variational Inference (VI)
- Stan (the language used by
brmsin the background) also implements variational inference algorithms that sample from an approximation of the posterior, which are faster but less accurate- You can use them by specifying:
brm(..., algorithm = "meanfield")(or"fullrank"or"pathfinder") - Unlike MCMC, which samples from the true posterior, Variational Inference (VI) approximates it using simpler distribution (e.g., Gaussian) and minimizes the difference between the two (e.g., optimizing the Kullback-Leibler KL divergence)
- You can use them by specifying:
- Pros:
- Can be useful for “quick” checks and adjustments but should not necessarily be used for final results
- If the uncertainty estimation is not the main goal but speed is
- Cons:
- Less accurate
- Fragile (algorithm can easily fail)
- VI posteriors may be overconfident (narrower CI)
Variational Inference (VI) in the Brain
- The brain is thought to be a Bayesian inference machine (Friston et al.), constantly predicting sensory input and updating beliefs based on errors
- Prediction errors (mismatch between expectation and reality) drive learning, perception, and action
- Computing exact posteriors would be too computationally expensive for the brain
- Instead, it might approximates the posterior in a way similar to VI:
- The brain represents beliefs with simpler distributions (e.g., Gaussian-like approximations).
- Instead of exploring all possibilities, it minimizes a function (“Free Energy”) to find the best fit
- The brain minimizes Free Energy to improve its model could be mathematically identical to VI minimizing the Kullback-Leibler (KL) divergence between an approximate posterior and the true posterior[^In fact, Friston’s “Variational Free Energy” is mathematically equivalent to what computer scientists call the Evidence Lower Bound (ELBO) in machine learning.]
- Cognitive and neural processing could be optimized for efficiency rather than exactness, much like how VI trades accuracy for speed
Sampling Diagnostics - Trace Plots
- HAIRY CATERPILLARS
- The MCMC algorithm is supposed to draw independent samples from the posterior
- There should be no autocorrelation between successive samples (i.e., no patterns in the trace plot)
- They should look like hairy caterpillars
- We typically draw multiple independent chains (useful for parallel computing, e.g., one chain per core)
plot(model)
- The plots shows 1) the posterior distribution of the parameters and 2) the trace plot of the MCMC algorithm
- One trace per chain: there should not be any pattern
Posterior Predictive Check - Distribution
- We can see whether the distribution of predicted values are closer to the observed values (in particular, its mean when using linear models)
ggplot(mtcars, aes(x=qsec)) +
geom_line(data=pred_long, aes(x=iter_value, group=iter_group),
stat="density", alpha=0.3) +
geom_density(fill="skyblue", alpha=0.5) +
theme_minimal()
- Note the use of
geom_line(stat="density")for the lines instead ofgeom_density()to be able to modify the alpha of the lines (and not the filling) - Posterior predictive checks of the distribution is particularly useful for more complex models to see if the model manages to reproduce the
Visualizing the Effects - Plot (1)
- Step 1: plot the raw data points
ggplot(mtcars, aes(x=mpg)) +
geom_point(aes(y=qsec), size=3) 
- Note that we specify the “y-aesthetic” in the
geom_point()layer and not in theggplot()layer, because other layers will different variable names for their y-axis (e.g.,Predictedfor the predicted line)
Visualizing the Effects - Plot (2)
- Step 2: add the median predicted line with the CI bounds
- Using the information from the “wide” format
ggplot(mtcars, aes(x=mpg)) +
geom_point(aes(y=qsec), size=3) +
geom_ribbon(data=pred_wide, aes(ymin=CI_low, ymax=CI_high), fill="darkred", alpha=0.2) +
geom_line(data=pred_wide, aes(y=Predicted), linewidth=2, color="orange")
Visualizing the Effects - Plot (3)
- Step 3: add the individual lines
ggplot(mtcars, aes(x=mpg)) +
geom_point(aes(y=qsec), size=3) +
geom_ribbon(data=pred_wide, aes(ymin=CI_low, ymax=CI_high), fill="darkred", alpha=0.2) +
geom_line(data=pred_wide, aes(y=Predicted), linewidth=2, color="orange") +
geom_line(data=pred_long, aes(y=iter_value, group=iter_group), alpha=0.1)
- Each line corresponds to one “plausible” regression line (given our posterior beliefs about the model parameters)
- Most of them indicate a positive relationship between
mpgandqsec
What about… Model Assumptions
- Most model assumptions that applied for Frequentist models still apply for Bayesian models (e.g., linearity, homoscedasticity, normality of residuals…)
- They are not assumptions of the estimation process, but rather assumptions of the model itself
- For instance, the assumption of Homoskedasticity directly comes from the fact that there’s a single fixed \(\sigma\) parameter for the model (which corresponds to the general uncertainty in predictions).
- The Normality of the residuals comes from the distribution that is being fit to the data (a Gaussian one, see lecture on OLS/MLE).
- However, the Bayesian framework is much more flexible and allows to easily relax some of these assumptions by using different families of models.
check_model(model)
