Maximum likelihood in the morning, Bayes in the afternoon

Morning · Why we care about likelihood

Every model you fit in this course is, secretly, a statement of the form $Y \sim f(\theta)$: data $Y$ are random draws from a distribution whose shape is controlled by parameters $\theta$. Once you write down $f$, the data become a function of $\theta$:

Same formula, same numbers — only the interpretation flips. As a function of $y$ it is a probability (or probability density). As a function of $\theta$ with $y$ held fixed at the value you actually observed, it is the likelihood. The likelihood is not a probability over $\theta$. It is a measure of how compatible different $\theta$ values are with the data you have.

Read the full plain-language summary for Day 4 → — the same ideas worked through slowly, with examples and no math.

Core concepts you will use all morning

Walk-through 1 · Bernoulli warmup: proportion of successes

Cleanest possible example. We watch 50 cars approach an intersection on a yellow-then-red light; 30 run the red, 20 stop. Each car is an independent Bernoulli trial with unknown probability $p$. The data $y = (y_1, \dots, y_{50})$ are 0s and 1s, with $\sum y_i = 30$.

The likelihood is

Take the log:

Differentiate, set to zero:

The proportion of successes is the MLE under Bernoulli/Binomial sampling. Same answer if we used a single Binomial$(n=50, p)$ likelihood — only the binomial coefficient differs, and it doesn't depend on $p$.

p_grid <- seq(0.01, 0.99, by = 0.001)
L <- p_grid^30 * (1 - p_grid)^20
nLL <- -30 * log(p_grid) - 20 * log(1 - p_grid)

par(mfrow = c(1, 2))
plot(p_grid, L, type = "l", lwd = 2, xlab = "p", ylab = "Likelihood")
abline(v = 0.6, col = "red")
plot(p_grid, nLL, type = "l", lwd = 2, xlab = "p", ylab = "-log L")
abline(v = 0.6, col = "red")

Walk-through 2 · Grid search: Von Bertalanffy growth

The Von Bertalanffy growth function (VBGF) describes asymptotic growth in fishes:

where $a$ is age, $L_a$ is length at age $a$, $L_\infty$ is asymptotic maximum length, and $\kappa$ is the growth coefficient. Given $n$ aged individuals, the likelihood is

The original course code holds $\sigma = 10$ and $L_\infty = 100$ fixed and grids over $\kappa \in \{0.1, 0.2, 0.3, 0.4\}$ to show how the likelihood changes with the parameter. Here is that grid search, complete and runnable.

# Load the Vonbert data (columns: sex, age, length). Keep sex==1.
TheData <- read.csv("data/vonbert.csv")
TheData <- TheData[TheData$sex == 1, ]
age <- TheData$age; len <- TheData$length

# negative log-likelihood as a function of (Linf, kappa, sigma)
nll_vbgf <- function(Linf, kappa, sigma) {
 pred <- Linf * (1 - exp(-kappa * age))
 -sum(dnorm(len, mean = pred, sd = sigma, log = TRUE))
}

# Grid search over kappa with Linf and sigma fixed (as in the original lab).
kappa_grid <- seq(0.05, 0.50, by = 0.005)
nll_vec <- sapply(kappa_grid, function(k) nll_vbgf(Linf = 100, kappa = k, sigma = 10))
plot(kappa_grid, nll_vec, type = "l", lwd = 2,
 xlab = expression(kappa), ylab = "Negative log-likelihood")
abline(v = kappa_grid[which.min(nll_vec)], col = "red")

Walk-through 3 · Numerical optimization: mle() and optim()

The stats4::mle() wrapper around optim() takes a function of named parameters and a list of starting values. It returns the MLE, the Hessian-based standard errors, and lets you compute Wald or profile confidence intervals.

library(stats4)

nll_vbgf3 <- function(Linf, kappa, sigma) {
 pred <- Linf * (1 - exp(-kappa * age))
 -sum(dnorm(len, mean = pred, sd = sigma, log = TRUE))
}

fit <- mle(nll_vbgf3,
 start = list(Linf = 100, kappa = 0.2, sigma = 10),
 method = "L-BFGS-B",
 lower = c(10, 0.001, 0.01),
 upper = c(500, 2, 100))
summary(fit)
confint(fit) # profile-likelihood CIs (the right kind for nonlinear models)
-logLik(fit)[1] # minimum NLL at the MLE

Walk-through 4 · Two-component Poisson mixture (whales)

Bowhead whale blow counts come from a mix of males and females, and the two sexes blow at different average rates. Treat each whale's sex as an unobserved (latent) label and model the counts as a two-component mixture: a male component with rate $\lambda_1$ and a female component with rate $\lambda_2$. A whale is male with probability $\pi$ — the proportion of males in the population — and female with probability $1-\pi$. The goal is to estimate all three quantities: the male blow rate $\lambda_1$, the female blow rate $\lambda_2$, and the sex proportions $\pi$ (males) and $1-\pi$ (females). The likelihood for one count $y$ is a mixture:

Three parameters, no closed form. We minimize the joint negative log-likelihood with mle().

TheData <- read.csv("data/whale.csv")   # columns: sex (1/2, unused here), blows
blows <- TheData$blows

nll_mix <- function(lam1, lam2, p) {
 comp <- p * dpois(blows, lam1) + (1 - p) * dpois(blows, lam2)
 -sum(log(comp))
}

fit_mix <- mle(nll_mix,
 start = list(lam1 = 10, lam2 = 30, p = 0.4),
 method = "L-BFGS-B",
 lower = c(0.01, 0.01, 0.001),
 upper = c(Inf, Inf, 0.999))
summary(fit_mix)

Walk-through 5 · Likelihood ratio tests & profile likelihood

LRT: do males and females share a growth curve?

Fit two nested Vonbert models to the full dataset (sex 1 and sex 2 combined):

• Model A (null, 3 parameters): one $L_\infty$, one $\kappa$, one $\sigma$, shared across sexes.
• Model B (alternative, 5 parameters): sex-specific $L_\infty$ and $\kappa$, shared $\sigma$.

Likelihood-ratio statistic: $\Lambda = 2(\ell_B - \ell_A)$. Under $H_0$: $\Lambda \sim \chi^2_{\Delta\,\text{df} = 2}$.

TheData <- read.csv("data/vonbert.csv")
Sex <- TheData$sex; Age <- TheData$age; Len <- TheData$length

nll_A <- function(Linf, kappa, sigma) {
 pred <- Linf * (1 - exp(-kappa * Age))
 -sum(dnorm(Len, pred, sigma, log = TRUE))
}
nll_B <- function(Linf1, kappa1, Linf2, kappa2, sigma) {
 pred <- ifelse(Sex == 1,
 Linf1 * (1 - exp(-kappa1 * Age)),
 Linf2 * (1 - exp(-kappa2 * Age)))
 -sum(dnorm(Len, pred, sigma, log = TRUE))
}

fitA <- mle(nll_A, start = list(Linf = 100, kappa = 0.2, sigma = 10),
 method = "L-BFGS-B", lower = c(10, 0.001, 0.01))
fitB <- mle(nll_B, start = list(Linf1 = 100, kappa1 = 0.2, Linf2 = 100, kappa2 = 0.2, sigma = 10),
 method = "L-BFGS-B", lower = c(10, 0.001, 10, 0.001, 0.01))

LRT <- 2 * (-logLik(fitA)[1] - (-logLik(fitB)[1]))
pval <- pchisq(LRT, df = 2, lower.tail = FALSE)
c(LRT = LRT, p_value = pval)

# AIC ranking on the same models
AIC(fitA, fitB)

Profile likelihood confidence intervals

Most likelihood surfaces in ecology are not symmetric around the MLE, so Wald CIs (MLE ± 1.96·SE) lie. The profile likelihood CI is the set of parameter values for which

i.e. all values whose deviance is within 1.92 NLL units of the minimum. confint() on an mle object gives exactly this. You can also draw it.

nll_p <- function(p) -(30 * log(p) + 20 * log(1 - p))
phat <- 0.6
best <- nll_p(phat)
pp <- seq(0.40, 0.80, length.out = 400)
plot(pp, nll_p(pp) - best, type = "l", lwd = 2,
 ylab = expression(-Delta * log~L), xlab = "p")
abline(h = 1.92, lty = 2, col = "red")
abline(v = phat, col = "grey")

Read the 95% CI off where the curve crosses the horizontal line at 1.92.

Bonus walk-through 6 · MLE for logistic regression, by hand vs glm()

This example shows that glm(y ~ x, family=binomial) is not magic — it is iteratively minimizing the same negative log-likelihood you can write down in five lines of R. Simulating data lets us check that glm() and your hand-rolled MLE agree to many decimal places.

Simulate $n=200$ presence/absence observations of a songbird detected as a function of forest cover $x$:

library(stats4)
set.seed(123)
n <- 200
x <- rnorm(n)
beta0_true <- -0.5; beta1_true <- 1.2
p <- 1 / (1 + exp(-(beta0_true + beta1_true * x)))
y <- rbinom(n, size = 1, prob = p)

# --- 1. Hand-rolled MLE ----------------------------------------------
nll_logit <- function(beta0, beta1) {
 eta <- beta0 + beta1 * x
 mu <- 1 / (1 + exp(-eta)) # inverse-logit
 -sum(dbinom(y, size = 1, prob = mu, log = TRUE))
}
fit_mle <- mle(nll_logit, start = list(beta0 = 0, beta1 = 0))

# --- 2. The R one-liner ----------------------------------------------
fit_glm <- glm(y ~ x, family = binomial(link = "logit"))

# --- 3. Side-by-side -------------------------------------------------
cbind(MLE = coef(fit_mle), GLM = coef(fit_glm), truth = c(beta0_true, beta1_true))
# Standard errors:
cbind(MLE = sqrt(diag(vcov(fit_mle))),
 GLM = sqrt(diag(vcov(fit_glm))))

The estimates agree precisely. The standard errors agree because both methods invert the same observed Fisher information.

Bonus walk-through 7 · Zero-inflated Poisson (a non-standard mixture)

You count fledglings per nest. Some nests fail entirely (true structural zeros: nest abandoned, depredated). The successful nests produce a Poisson count of fledglings — including occasional Poisson zeros. The mixture has two sources of zeros and is called a zero-inflated Poisson (ZIP):

Here $\pi$ is the probability of structural zero and $\lambda$ is the Poisson mean among non-structural nests. Base R's glm() has no zero-inflated family, so ZIP is not a one-line glm(...) call — but ready-made ZIP fitters do exist: pscl::zeroinfl() is the standard one, and the VGAM and glmmTMB packages fit ZIP models too. We write the negative log-likelihood by hand anyway, because it is a single line and doing so is the transferable skill for the non-standard models that have no canned fitter:

library(stats4)
set.seed(11)
n <- 300
pi_true <- 0.35; lam_true <- 2.5
z <- rbinom(n, 1, pi_true) # 1 = structural zero
y <- ifelse(z == 1, 0, rpois(n, lam_true))

# Plain Poisson MLE underestimates the mean and misses the zero inflation:
mean(y); var(y)
mean(y == 0); dpois(0, mean(y)) # observed zeros >> Poisson-predicted zeros

# ZIP negative log-likelihood:
nll_zip <- function(p, lambda) {
 if (p <= 0 || p >= 1 || lambda <= 0) return(1e10)
 prob_zero <- p + (1 - p) * exp(-lambda)
 ll_pos <- log(1 - p) + dpois(y, lambda, log = TRUE)
 ll <- ifelse(y == 0, log(prob_zero), ll_pos)
 -sum(ll)
}
fit_zip <- mle(nll_zip, start = list(p = 0.2, lambda = 1),
 method = "L-BFGS-B",
 lower = c(0.001, 0.001),
 upper = c(0.999, 50))
summary(fit_zip)
confint(fit_zip)

Afternoon · From likelihood to posterior

You spent the morning maximizing $L(\theta \mid y)$ — finding the single $\theta$ that makes the data most plausible. Bayesian inference goes further: it gives you a distribution over $\theta$ that combines the data with what you knew (or assumed) before. Same likelihood, plus a prior, equals a posterior:

The marginal $p(y) = \int p(y\mid\theta) p(\theta)\,d\theta$ does not depend on $\theta$, so it acts only as a normalizing constant. For most realistic models it is impossible to compute in closed form — and that is exactly why MCMC exists. MCMC samples from the unnormalized product on the right, dodging the integral entirely.

Bayes' rule, worked out for a single coin

To anchor the algebra, fit a binomial model to a single observation: $y = 7$ successes in $n = 10$ trials. Use a flat prior $p \sim \text{Uniform}(0,1)$ — which, note, is exactly a $\text{Beta}(1,1)$ distribution (a constant density on $[0,1]$):

Recognize this kernel? It is a Beta(8, 4) density (up to a constant). The normalizing constant exists in closed form via the Beta function:

This is the magic of conjugacy: the prior family is chosen so that the posterior stays in the same family. The flat $\text{Uniform}(0,1)$ prior is a $\text{Beta}(1,1)$, so this really is a Beta prior + Binomial likelihood → Beta posterior — and the posterior $\text{Beta}(8,4) = \text{Beta}(1+7,\, 1+3)$ simply adds the 7 successes and 3 failures to the prior's two parameters.

theta <- seq(0.001, 0.999, length.out = 500)
prior <- dbeta(theta, 1, 1)
likelihood <- dbinom(7, size = 10, prob = theta)
posterior <- dbeta(theta, 8, 4)
plot(theta, posterior, type = "l", lwd = 2, col = "darkred",
 ylab = "density (or scaled likelihood)", xlab = expression(theta))
lines(theta, prior, lty = 2)
lines(theta, likelihood / max(likelihood) * max(posterior), col = "darkblue")
legend("topleft", lty = c(2, 1, 1), col = c("black", "darkblue", "darkred"),
 legend = c("Prior Beta(1,1)", "Likelihood (scaled)", "Posterior Beta(8,4)"), bty = "n")

Conjugacy 1 · Beta–Binomial in full

For $y \mid p \sim \text{Bin}(n, p)$ and $p \sim \text{Beta}(\alpha, \beta)$,

The posterior just adds successes to $\alpha$ and failures to $\beta$.

Worked example: spotted-owl nest survival. You monitored 12 nests; 9 fledged. Choose three priors and see how each posterior moves.

y <- 9; n <- 12
theta <- seq(0.001, 0.999, length.out = 500)
priors <- list(
 flat = c(1, 1), # Uniform
 skeptical = c(2, 8), # prior centered at 0.2
 informative = c(30, 10) # tight, centered at 0.75
)
par(mfrow = c(1, 3))
for (nm in names(priors)) {
 ab <- priors[[nm]]; post <- c(ab[1] + y, ab[2] + n - y)
 plot(theta, dbeta(theta, post[1], post[2]), type = "l", lwd = 2,
 main = paste(nm, ": Beta(", post[1], ", ", post[2], ")", sep = ""),
 xlab = expression(theta), ylab = "posterior density")
 lines(theta, dbeta(theta, ab[1], ab[2]), lty = 2, col = "grey")
 abline(v = y / n, col = "red")
}
par(mfrow = c(1, 1))

Conjugacy 2 · Gamma–Poisson

For $y_1, \dots, y_n \mid \lambda \sim \text{Pois}(\lambda)$ and $\lambda \sim \text{Gamma}(\alpha, \beta)$ (shape-rate),

Again, the prior plays the role of pseudo-observations: $\alpha$ events accumulated over $\beta$ units of effort. After observing $n$ new units and $\sum y_i$ events, you just add them in. Note that the bookkeeping differs from the Beta–Binomial case: here the prior's weight is carried by $\beta$ alone — the pseudo-effort that your $n$ new units are added to — not by $\alpha+\beta$. Each conjugate family has its own currency for prior strength; read it off the update rule rather than assuming one formula covers them all.

Worked example: a parasitologist counts $y = (3, 1, 4, 2, 0, 5, 2, 3, 1, 4)$ ticks on 10 birds (total 25 in 10 birds). Prior $\lambda \sim \text{Gamma}(2, 1)$ (weakly informative around mean 2 per bird).

y <- c(3, 1, 4, 2, 0, 5, 2, 3, 1, 4); n <- length(y)
alpha0 <- 2; beta0 <- 1
alpha1 <- alpha0 + sum(y)
beta1 <- beta0 + n
c(post_mean = alpha1 / beta1,
 post_sd = sqrt(alpha1 / beta1^2),
 MLE = mean(y)) # compare to MLE

lam <- seq(0, 6, length.out = 400)
plot(lam, dgamma(lam, alpha1, beta1), type = "l", lwd = 2,
 xlab = expression(lambda), ylab = "posterior")
lines(lam, dgamma(lam, alpha0, beta0), lty = 2, col = "grey")
abline(v = mean(y), col = "red")

Why conjugacy is not enough — and what MCMC buys you

Conjugacy gives elegant closed forms for one-parameter problems. Real models have:

• Nonlinear functions of parameters (logistic growth: $\mu = r - r/K \cdot x$ depends on the ratio $r/K$).
• Many parameters interacting (random effects, hierarchical means and variances).
• Awkward likelihoods (mixtures, occupancy, capture-recapture, state-space).

For these models the marginal $p(y) = \int p(y,\theta)\,d\theta$ is intractable. MCMC sidesteps the integral by generating dependent samples from the unnormalized posterior; the long-run frequency of those samples is the posterior itself. The two algorithms you will meet today and tomorrow:

The Metropolis algorithm in 30 lines of R

Same Beta–Binomial spotted-owl example, but pretending we don't know the closed form. The Metropolis recipe at each iteration:

1. Propose $\theta^\star \sim \text{Normal}(\theta_{\text{current}}, s)$ (symmetric proposal).
2. Compute the acceptance ratio $\alpha = \min\!\left(1, \dfrac{p(\theta^\star)\, L(\theta^\star \mid y)}{p(\theta_{\text{current}})\, L(\theta_{\text{current}} \mid y)}\right)$.
3. Draw $u \sim \text{Uniform}(0,1)$. If $u < \alpha$, accept; otherwise stay.

y <- 9; n <- 12; alpha0 <- 2; beta0 <- 2

log_post <- function(p) {
 if (p <= 0 || p >= 1) return(-Inf)
 dbeta(p, alpha0, beta0, log = TRUE) +
 dbinom(y, size = n, prob = p, log = TRUE)
}

niter <- 20000; chain <- numeric(niter); chain[1] <- 0.5
accept <- 0; s <- 0.35 # proposal SD
set.seed(2)
for (i in 2:niter) {
 prop <- rnorm(1, mean = chain[i - 1], sd = s)
 log_alpha <- log_post(prop) - log_post(chain[i - 1])
 if (log(runif(1)) < log_alpha) { chain[i] <- prop; accept <- accept + 1 }
 else chain[i] <- chain[i - 1]
}

burn <- 2000
cat("acceptance rate =", round(accept / niter, 3), "\n")
hist(chain[-(1:burn)], breaks = 60, freq = FALSE, main = "Posterior of p",
 xlab = "p", col = "grey90", border = "white")
# Overlay the closed-form Beta(alpha0+y, beta0+n-y) posterior
curve(dbeta(x, alpha0 + y, beta0 + n - y), col = "red", lwd = 2, add = TRUE)

The main event · logistic growth with Nimble

We now fit the ecological example the morning's lm() couldn't honestly answer: posterior distributions for the logistic-growth parameters $r$, $K$, and $\sigma$, plus derived quantities $K/2$ (population size at maximum growth) and $\text{MSY} = rK/4$ (theoretical maximum sustainable yield).

The model is the per-capita form of $dN/dt = rN(1 - N/K)$:

The complete Nimble workflow

library(nimble); library(MCMCvis); library(HDInterval); library(coda)
Logistic <- read.csv("data/Logistic.csv", header = TRUE)
Logistic <- Logistic[order(Logistic$PopulationSize), ]
N_grid <- seq(0, 1500, by = 10); N_grid[1] <- 1

logistic_code <- nimbleCode({
 # priors
 K ~ dunif(0, 4000)
 r ~ dunif(0, 2)
 sigma ~ dunif(0, 2)
 tau <- 1 / sigma^2 # BUGS-style precision

 # likelihood
 for (i in 1:n) {
 mu[i] <- r - (r / K) * x[i]
 y[i] ~ dnorm(mu[i], tau)
 }

 # derived quantities
 N_at_maxdNdt <- K / 2
 MSY <- r * K / 4

 for (j in 1:nN) {
 dNdt[j] <- r * N[j] * (1 - N[j] / K)
 }
})

# Nimble splits the JAGS "data" list into:
# constants : sizes, indices, fixed covariates, prediction grids
# data : observed stochastic nodes (y[i] here)
constants <- list(
 n = nrow(Logistic),
 x = as.double(Logistic$PopulationSize),
 N = N_grid,
 nN = length(N_grid)
)
data <- list(y = as.double(Logistic$GrowthRate))

inits_list <- list(
 list(K = 1500, r = 0.20, sigma = 1.0),
 list(K = 1000, r = 0.15, sigma = 0.1),
 list(K = 900, r = 0.30, sigma = 0.01)
)

set.seed(10)
samples <- nimbleMCMC(
 code = logistic_code,
 constants = constants,
 data = data,
 inits = inits_list,
 monitors = c("K", "r", "sigma", "mu", "N_at_maxdNdt", "MSY", "dNdt"),
 nchains = 3,
 nburnin = 5000,
 niter = 25000,
 thin = 1,
 samplesAsCodaMCMC = TRUE
)

MCMCtrace(samples, params = c("K", "r", "sigma"), pdf = FALSE, Rhat = TRUE, n.eff = TRUE)
MCMCsummary(samples, params = c("K", "r", "sigma", "N_at_maxdNdt", "MSY"), round = 3)

# P(r > 0.22 | data)
r_chain <- MCMCchains(samples, params = "r")
mean(r_chain > 0.22)

Two kinds of credible interval: BCI vs HPDI

A 95% Bayesian credible interval (BCI, equal-tailed) cuts 2.5% off each tail of the posterior. A 95% highest posterior density interval (HPDI) is the shortest interval containing 95% of the posterior — by construction, every point inside has higher posterior density than every point outside.

• For symmetric posteriors (a typical sample-size-large $r$ posterior here), BCI and HPDI are nearly identical.
• For skewed posteriors (often $\sigma$, $K$ when the prior bound bites, ratios of parameters), HPDI is shorter and shifted toward the mode.
• HPDI is the right interval when you want to say "the 95% most-credible values of $\theta$."
• BCI is the right interval when you want a quantile-symmetric summary (e.g. for forecasts where the median matters).

BCI <- MCMCpstr(samples, params = "mu",
 func = function(x) quantile(x, c(.025, .50, .975)))
HPDI <- MCMCpstr(samples, params = "mu",
 func = function(x) hdi(x, .95))

plot(constants$x, data$y, pch = 19,
 xlab = "Population size N", ylab = "Per-capita growth (1/yr)")
lines(constants$x, BCI$mu[, 2], lwd = 2)
lines(constants$x, BCI$mu[, 1], lty = 2, col = "red")
lines(constants$x, BCI$mu[, 3], lty = 2, col = "red")
lines(constants$x, HPDI$mu[, 1], lty = 2, col = "blue")
lines(constants$x, HPDI$mu[, 2], lty = 2, col = "blue")

Bonus walk-through 8 · Bayesian binomial estimation of survival

This example shows beta-binomial conjugacy from end to end, with the prior justification spelled out and a comparison between (i) the closed-form posterior, (ii) a Nimble fit, and (iii) Metropolis. All three give the same answer — they have to, because the model is the same.

Setting. You banded 25 juvenile Western Bluebirds in spring; 18 were resighted alive the following spring. The biological literature suggests annual juvenile survival in this species is "around 0.45 with a fair amount of variability." We encode this prior as Beta(9, 11), which has mean $9/20 = 0.45$ and an effective sample size of 20 trials.

y <- 18; n <- 25; alpha0 <- 9; beta0 <- 11
alpha1 <- alpha0 + y # 27
beta1 <- beta0 + n - y # 18
c(mean = alpha1 / (alpha1 + beta1),
 median = qbeta(0.5, alpha1, beta1),
 BCI_lo = qbeta(0.025, alpha1, beta1),
 BCI_hi = qbeta(0.975, alpha1, beta1))

surv_code <- nimbleCode({
 phi ~ dbeta(9, 11) # prior
 y ~ dbinom(phi, n) # likelihood
})

set.seed(7)
surv_samps <- nimbleMCMC(
 code = surv_code,
 constants = list(n = 25),
 data = list(y = 18),
 inits = list(list(phi = 0.5), list(phi = 0.3), list(phi = 0.7)),
 monitors = "phi",
 nchains = 3, nburnin = 2000, niter = 10000,
 samplesAsCodaMCMC = TRUE)
MCMCsummary(surv_samps, round = 3)

# P(phi > 0.5 | data)?
mean(MCMCchains(surv_samps, "phi") > 0.5)
# closed-form check:
1 - pbeta(0.5, alpha1, beta1)

Bonus walk-through 9 · Bayesian Normal regression, with posterior predictive checks

Body mass as a function of body length in Pacific giant salamanders. We will (i) think carefully about priors on the scale of the data, (ii) fit the regression in Nimble, and (iii) do a posterior predictive check by simulating replicate datasets from the posterior and comparing them to the real one. This is the workflow you should use on every Bayesian model you fit.

set.seed(21)
n <- 120
length_mm <- runif(n, 80, 250)
true_int <- -35
true_slope <- 0.32
true_sigma <- 6.0
mass_g <- true_int + true_slope * length_mm + rnorm(n, 0, true_sigma)
plot(length_mm, mass_g, pch = 19, xlab = "length (mm)", ylab = "mass (g)")

# Prior reasoning:
# * intercept: at length=0 the mean mass should be plausibly small; N(0, 100^2) is mild.
# * slope: a mm of length should add at most a few grams; N(0, 1) is very wide on this scale.
# * sigma: residual SD must be positive and on the same scale as mass (grams).
# half-N(0, 20) keeps it positive with most mass below 50.

reg_code <- nimbleCode({
 # priors
 beta0 ~ dnorm(0, sd = 100)
 beta1 ~ dnorm(0, sd = 1)
 sigma ~ T(dnorm(0, sd = 20), 0, ) # half-normal via truncation

 # likelihood + posterior predictive draws
 for (i in 1:n) {
 mu[i] <- beta0 + beta1 * length_mm[i]
 mass[i] ~ dnorm(mu[i], sd = sigma)
 mass_rep[i] ~ dnorm(mu[i], sd = sigma) # replicate dataset
 }

 # No discrepancy statistics in here. sd(mass) is a constant — putting it in the
 # model re-computes it every iteration and monitors a number that never moves.
 # We compute discrepancies in R below from mu and mass_rep, so we can try
 # several of them without refitting.
})

constants <- list(n = n, length_mm = length_mm)
data <- list(mass = mass_g)
inits <- list(
 list(beta0 = 0, beta1 = 0, sigma = 10),
 list(beta0 = -50, beta1 = 0.5, sigma = 5),
 list(beta0 = 20, beta1 = -0.1, sigma = 15))

reg_samps <- nimbleMCMC(
 code = reg_code,
 constants = constants,
 data = data,
 inits = inits,
 monitors = c("beta0", "beta1", "sigma", "mu", "mass_rep"),
 nchains = 3, nburnin = 3000, niter = 15000,
 samplesAsCodaMCMC = TRUE)
MCMCsummary(reg_samps, params = c("beta0", "beta1", "sigma"), round = 3)

mu_post  <- MCMCchains(reg_samps, "mu")        # fitted mass for each salamander
rep_post <- MCMCchains(reg_samps, "mass_rep")  # a fake dataset, one per draw

# Two ways to summarise "how far are the points from the line?":
#   SSR          = the total squared distance
#   biggest miss = the single worst salamander
# Compute each on the real data and on the fake data, at every posterior draw.
S <- nrow(mu_post)
obs_ssr <- numeric(S)
rep_ssr <- numeric(S)
obs_miss <- numeric(S)
rep_miss <- numeric(S)

for (s in 1:S) {
 e_obs <- mass_g - mu_post[s, ] # residuals, real data
 e_rep <- rep_post[s, ] - mu_post[s, ] # residuals, fake data

 obs_ssr[s] <- sum(e_obs^2)
 rep_ssr[s] <- sum(e_rep^2)
 obs_miss[s] <- max(abs(e_obs))
 rep_miss[s] <- max(abs(e_rep))
}

# Bayesian p-value: how often is the fake data as bad as the real data?
mean(rep_ssr >= obs_ssr) # sum of squared residuals
mean(rep_miss >= obs_miss) # biggest miss

# Where does the real dataset sit among the fake ones?
hist(rep_miss, xlab = "biggest miss in a fake dataset (g)", main = "")
abline(v = mean(obs_miss), col = "red", lwd = 2) # real data: inside the pile, or off in the tail?

# Visual PPC: overlay 100 replicate datasets on the data
mass_rep_chain <- MCMCchains(reg_samps, "mass_rep")
plot(length_mm, mass_g, pch = 19, xlab = "length (mm)", ylab = "mass (g)")
for (i in sample(nrow(mass_rep_chain), 100)) {
 points(length_mm, mass_rep_chain[i, ], col = rgb(0, 0, 0, 0.04), pch = 19, cex = 0.5)
}
points(length_mm, mass_g, pch = 19, col = "red")

Exercises

Eight numbered exercises spanning the day. Work them with a partner; check the answer key only after you've drafted your own solution. Difficulty pills indicate rough time commitment.

library(stats4)
nll <- function(p) -(30 * log(p) + 20 * log(1 - p))
fit <- mle(nll, start = list(p = 0.5),
 method = "L-BFGS-B", lower = 0.001, upper = 0.999)
summary(fit); confint(fit)

library(stats4)
TheData <- read.csv("data/vonbert.csv")
Age <- TheData$age; Len <- TheData$length

profile_kappa <- function(k) {
 nll <- function(Linf, sigma) {
 pred <- Linf * (1 - exp(-k * Age))
 -sum(dnorm(Len, pred, sigma, log = TRUE))
 }
 fit <- mle(nll, start = list(Linf = 100, sigma = 10),
 method = "L-BFGS-B", lower = c(10, 0.01))
 -logLik(fit)[1]
}
grid <- seq(0.05, 0.45, by = 0.005)
prof <- sapply(grid, profile_kappa)
best <- min(prof)
plot(grid, prof - best, type = "l", lwd = 2, xlab = expression(kappa),
 ylab = expression(-Delta * log~L))
abline(h = 1.92, lty = 2, col = "red")

library(HDInterval)
post <- rbeta(1e5, 25, 25)
c(mean = mean(post), median = median(post),
 quantile(post, c(.025, .975)),
 hpdi = hdi(post, .95))

y <- c(3,1,4,2,0,5,2,3,1,4)
log_post <- function(log_lam) {
 lam <- exp(log_lam)
 dgamma(lam, 2, 1, log = TRUE) + sum(dpois(y, lam, log = TRUE)) + log_lam # Jacobian
}
ll <- numeric(3e4); ll[1] <- log(2); acc <- 0; s <- 0.30
for (i in 2:length(ll)) {
 prop <- rnorm(1, ll[i - 1], s)
 if (log(runif(1)) < log_post(prop) - log_post(ll[i - 1])) { ll[i] <- prop; acc <- acc + 1 } else ll[i] <- ll[i - 1]
}
lam <- exp(ll[-(1:3000)])
hist(lam, breaks = 60, freq = FALSE)
curve(dgamma(x, 2 + sum(y), 1 + length(y)), add = TRUE, col = "red", lwd = 2)
acc / length(ll)