1.3 Three tools for golem engineering

This book teaches three tools for building better statistical models: Bayesian inference, multilevel models, and information criteria. These tools are not mutually exclusive — they reinforce each other and form the foundation of the book's approach to statistical modeling.

Tool 1: Bayesian inference

Bayesian inference is a framework for updating beliefs under uncertainty. Given a prior distribution $P(\theta)$ and a likelihood $P(D \mid \theta)$, Bayes' theorem tells us how to compute the posterior distribution:

$\underbrace{p(\theta \mid D)}_{\text{posterior}} \propto \underbrace{p(D \mid \theta)}_{\text{likelihood}} \cdot \underbrace{p(\theta)}_{\text{prior}}$

The posterior distribution is the complete answer to the question "what do we know about $\theta$ after seeing the data?" It is a distribution, not a point estimate, so it naturally captures uncertainty.

Bayesian inference has several advantages over frequentist methods:

• Full uncertainty: The posterior gives us the entire distribution of plausible values, not just a point estimate and a standard error.
• Prior information: We can formally incorporate what we already know, rather than pretending to be ignorant.
• Coherence: The rules of probability ensure that our conclusions are internally consistent.
• Flexibility: The same framework works for simple models and complex hierarchical models.

# Bayesian inference from scratch: grid approximation
# Estimate the proportion of water on a globe
# Data: 6 waters in 9 tosses

# 1. Define the parameter grid

p_grid <- seq(from = 0, to = 1, length.out = 100)

# 2. Define the prior: uniform (uninformative)

prior <- rep(1, 100)

# 3. Compute the likelihood: binomial

likelihood <- dbinom(6, size = 9, prob = p_grid)

# 4. Compute the posterior: multiply and normalize

posterior <- (likelihood _ prior) / sum(likelihood _ prior)

# 5. Summarize the posterior

cat("Posterior mean:", round(sum(p*grid * posterior), 3), "\n")
cat("Posterior std dev:", round(sqrt(sum((p*grid - sum(p_grid * posterior))^2 * posterior)), 3), "\n")
cat("95% credible interval: [",
round(p_grid[which(cumsum(posterior) >= 0.025)[1]], 3), ",",
round(p_grid[which(cumsum(posterior) >= 0.975)[1]], 3), "]\n")

# Plot the posterior

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "Proportion of water (p)",
ylab = "Posterior density",
main = "Bayesian posterior: 6 waters in 9 tosses")
polygon(c(p_grid, rev(p_grid)),
c(posterior, rep(0, length(posterior))),
col = rgb(0.27, 0.51, 0.71, 0.3), border = NA)
abline(v = 6/9, lty = 2, col = "firebrick")
text(0.72, max(posterior) * 0.9, "MLE = 6/9", col = "firebrick")

The posterior distribution captures everything we know about $p$ after seeing the data. The mean, standard deviation, and credible interval are all summaries of this distribution. Unlike a $p$-value, the posterior tells us the full range of plausible values and their relative probabilities.

Tool 2: Multilevel models

Multilevel models (also called hierarchical or mixed-effects models) handle data with structure: students nested in schools, measurements nested in patients, trials nested in subjects. Instead of assuming that all groups are identical (complete pooling) or completely independent (no pooling), multilevel models use partial pooling: each group gets its own parameters, but the group parameters are drawn from a common distribution.

The key idea is that the population distribution links the groups together. Groups with little data are pulled toward the population mean (shrinkage), while groups with lots of data are left mostly alone. This produces better estimates for all groups, especially those with small sample sizes.

# Demonstrate partial pooling with a simple simulation
# 5 tanks of tadpoles, each with different sample sizes
set.seed(42)
n_tanks <- 5
true_survival <- c(0.8, 0.5, 0.7, 0.3, 0.9)  # true survival per tank
tank_sizes <- c(10, 20, 5, 30, 8)  # different sample sizes

# Simulate data

survived <- rbinom(n_tanks, tank_sizes, true_survival)

# No pooling: estimate each tank independently

no_pool <- survived / tank_sizes

# Complete pooling: ignore tank differences

complete_pool <- rep(sum(survived) / sum(tank_sizes), n_tanks)

# Partial pooling (simplified): shrink toward overall mean

overall_mean <- sum(survived) / sum(tank_sizes)

# Shrinkage factor depends on sample size

shrinkage <- tank*sizes / (tank_sizes + 10) # 10 is a tuning parameter
partial_pool <- shrinkage * no*pool + (1 - shrinkage) * overall_mean

# Compare

comparison <- data.frame(
tank = 1:n_tanks,
true = true_survival,
n = tank_sizes,
no_pooling = round(no_pool, 3),
partial_pooling = round(partial_pool, 3),
complete_pooling = round(complete_pool, 3)
)
print(comparison)

# Plot

plot(1:n_tanks, true_survival, pch = 16, col = "black", cex = 1.5,
ylim = c(0, 1), xlab = "Tank", ylab = "Survival probability",
main = "Partial pooling: shrinkage toward the mean")
points(1:n_tanks, no_pool, pch = 17, col = "firebrick")
points(1:n_tanks, partial_pool, pch = 18, col = "steelblue")
points(1:n_tanks, complete_pool, pch = 15, col = "darkgreen")
legend("bottomleft",
legend = c("True", "No pooling", "Partial pooling", "Complete pooling"),
col = c("black", "firebrick", "steelblue", "darkgreen"),
pch = c(16, 17, 18, 15), cex = 0.8)

Notice how partial pooling pulls the estimates for small tanks (e.g., tank 3 with $n=5$) toward the overall mean, while leaving large tanks (e.g., tank 4 with $n=30$) mostly unchanged. This is the key benefit of multilevel models: they borrow strength across groups.

Tool 3: Information criteria

Information criteria (like WAIC and PSIS) estimate out-of-sample predictive accuracy. They tell us which model will predict new data best, without needing to collect new data. This is the cure for overfitting: instead of choosing the model that fits the sample best (which always favors more complex models), we choose the model that will generalize best.

The core idea is the deviance, which measures how well the model fits the data:

$\text{Deviance} = -2 \sum_i \log P(y_i \mid \theta)$

But the deviance of the training data always improves as we add parameters. The trick is to estimate the deviance of new data, which penalizes complexity. WAIC (Widely Applicable Information Criterion) does this by computing the log-pointwise-predictive-density and adding a penalty for the effective number of parameters:

$\text{WAIC} = -2(\text{lppd} - p_\text{WAIC})$

where $\text{lppd}$ is the log-pointwise-predictive-density and $p_\text{WAIC}$ is the effective number of parameters.

# Demonstrate overfitting: more parameters always fit better
set.seed(42)
n <- 20
x <- runif(n, 0, 10)
y <- 2 * x + rnorm(n, sd = 2)

# Fit polynomials of increasing degree

degrees <- 1:8
train_error <- numeric(length(degrees))

for (i in seq_along(degrees)) {
d <- degrees[i]
fit <- lm(y ~ poly(x, degree = d, raw = TRUE))
train_error[i] <- mean(resid(fit)^2) # training MSE
}

# Plot training error vs model complexity

plot(degrees, train_error, type = "b", pch = 16, lwd = 2,
col = "firebrick", xlab = "Polynomial degree",
ylab = "Training MSE",
main = "Overfitting: training error always decreases")
grid()

# Add a simulated test error curve

test*error <- train_error + 0.3 * degrees^2 _ 0.1
lines(degrees, test_error, type = "b", pch = 17, lwd = 2, col = "steelblue")
legend("topright",
legend = c("Training error (in-sample)", "Test error (out-of-sample)"),
col = c("firebrick", "steelblue"), pch = c(16, 17), lty = 1, lwd = 2)

The training error always decreases as we add parameters — the golem happily fits the noise. But the test error (estimated out-of-sample) has a minimum: beyond a certain complexity, the model starts overfitting and generalizes worse. Information criteria estimate this out-of-sample error without needing a separate test set.

How the three tools work together

The three tools form a coherent workflow:

1. Specify: Write down a scientific model and translate it into a Bayesian model (likelihood + priors).
2. Fit: Use Bayesian inference (quadratic approximation, MCMC) to compute the posterior.
3. Critique: Use posterior predictive checks to see if the model can reproduce the data.
4. Compare: Use information criteria to compare models and estimate which will predict best.

If the model fails the critique, we go back to step 1 and specify a better model. If multiple models pass the critique, we compare them and possibly average their predictions.

# The full workflow: specify, fit, critique, compare
# Model 1: Linear relationship
# Model 2: Quadratic relationship

set.seed(42)
x <- runif(50, 0, 10)
y <- 1 + 2 _ x + 0.3 _ x^2 + rnorm(50, sd = 3)

# Step 1 & 2: Specify and fit both models

m1 <- lm(y ~ x) # linear
m2 <- lm(y ~ x + I(x^2)) # quadratic

# Step 3: Critique - compare R-squared

cat("Model 1 (linear): R-squared =", round(summary(m1)$r.squared, 4), "\n")
cat("Model 2 (quadratic): R-squared =", round(summary(m2)$r.squared, 4), "\n")

# Step 4: Compare using AIC (a simple information criterion)

cat("\nAIC (lower is better):\n")
cat("Model 1:", round(AIC(m1), 2), "\n")
cat("Model 2:", round(AIC(m2), 2), "\n")
cat("\nDifference:", round(AIC(m1) - AIC(m2), 2), "\n")
cat("Model 2 is preferred by", round(AIC(m1) - AIC(m2), 2), "AIC units\n")

Summary of Section 1.3

• Bayesian inference updates beliefs under uncertainty using Bayes' theorem: posterior $\propto$ likelihood $\times$ prior.
• Multilevel models handle structured data through partial pooling: groups borrow strength from each other.
• Information criteria estimate out-of-sample predictive accuracy to prevent overfitting.
• The three tools form a workflow: specify, fit, critique, compare.
• These are engineering practices, not magic spells — they make modeling transparent and iterative.