Day 4 · Plain-language summary

Likelihood, priors, posteriors, and MCMC, in normal English

Read this before tackling the Day 4 lab. It explains every idea (maximum likelihood, Bayes' rule, MCMC, Nimble) with as little math as possible.

What we did today, in one paragraph

We learned two ways of asking the question "what parameter values do my data suggest?" In the morning we did maximum likelihood, picking the parameter values that make the data we actually saw as plausible as possible. In the afternoon we did Bayes, combining those same plausibilities (the likelihood) with what we knew before (the prior) to get a full probability distribution over the parameters (the posterior). We saw how to compute posteriors in two cases by hand (Beta–Binomial, Gamma–Poisson) and how to use Markov chain Monte Carlo, implemented in Nimble, to compute posteriors for any model.

What is a likelihood?

A likelihood is a number, well, a function, that tells you how compatible your data are with a particular guess about the parameters. It uses the same formula as the probability of seeing the data, but it is read differently:

Probability: Hold the parameter fixed, ask "how likely is each possible dataset?" Sum/integrate over datasets.
Likelihood: Hold the data fixed (at what you actually observed), ask "how compatible is each possible parameter value with this dataset?" Do not sum/integrate over parameters, because it isn't a probability over parameters.

A likelihood is not a probability about parameters. That sentence is the most important conceptual hurdle in classical statistics. If you want a probability about parameters, you need Bayes (afternoon).

What is a maximum-likelihood estimate?

The MLE is the parameter value that makes the data as plausible as possible. If you watched 50 cars and 30 ran the red light, the MLE for the probability of running the light is 30/50 = 0.6. If you weighed 100 fish and want $L_\infty$ for a Von Bertalanffy growth curve, the MLE is the value of $L_\infty$ that minimizes the squared distance between the curve and the points (when the noise is Gaussian, exactly the case where MLE and least-squares give the same answer).

For complicated models the MLE doesn't have a closed form. We use a computer to numerically optimize the log-likelihood. R's mle() and optim() functions do this for us; they also give back standard errors based on the curvature of the log-likelihood at the maximum (the second derivative, the "Hessian").

Why ecologists love MLE. It gives a single best estimate, a standard error, and a CI. It is computationally cheap for moderate-size problems. It plays well with information criteria (AIC) and likelihood-ratio tests. Most published ecological models are MLE fits. Anything that says lm, glm, or nlme is using MLE under the hood.

Two things you can do once you have a likelihood

The morning lab spends a lot of time on two tools that fall straight out of the likelihood, and both are easier than they sound. One compares models. The other builds an honest confidence interval for a single parameter.

The likelihood-ratio test: is the bigger model worth it?

Suppose you fit two models to the same data. The simple one says one growth curve describes every fish in the file. The complicated one gives each sex its own curve, two extra parameters. The complicated model will always fit at least as well, because it contains the simple one as a special case. It can always just set the two curves equal. So "it fit better" proves nothing on its own. The real question is whether it fit enough better to justify the two extra parameters.

The likelihood-ratio test (LRT) answers that. Take the log-likelihood of each fit, the number that says how well each model explains the data, and compute

$\Lambda = 2\,\big(\ell_{\text{complicated}} - \ell_{\text{simple}}\big).$

Big $\Lambda$ means the extra parameters bought you a lot. The useful fact is that we know how big $\Lambda$ gets by luck alone when the simple model is actually true. It follows a $\chi^2$ distribution with degrees of freedom equal to the number of extra parameters. So you compare your $\Lambda$ to that reference and read off a p-value. In the Day 4 growth example $\Lambda = 75.4$ on 2 df, which gives $p \approx 4\times10^{-17}$, so luck is not a plausible explanation and the sexes really do grow differently.

The one rule you must not break: the models have to be nested. Nested means the simple model is a special case of the complicated one, reachable by fixing some of its parameters. One-curve-for-all is nested inside a-curve-per-sex (set the two curves equal). A Poisson model is not nested inside a lognormal one. If the models are not nested the $\chi^2$ reference is simply wrong and the p-value is meaningless. That is precisely the gap AIC fills: AIC ranks any set of models fit to the same data, nested or not, but it gives you a ranking, not a test, and never a p-value. Use the LRT when you want a yes/no on a specific added effect; use AIC when you want to sort a collection of candidates.

Profile likelihood: a confidence interval that does not lie

The quick way to get a 95% CI is "estimate $\pm$ 1.96 standard errors." That recipe quietly assumes the likelihood is a symmetric hill around its peak. For a mean of a large sample, fine. For nonlinear parameters like a growth coefficient $\kappa$, a rate, or a variance, the hill is usually lopsided, steep on one side, with a long shallow shoulder on the other, and the symmetric interval then puts its bounds in the wrong places. It can even hand you a negative lower bound for a quantity that cannot be negative.

The profile likelihood fixes that by simply looking at the actual shape of the hill. Pick a candidate value for the parameter you care about, say $\kappa = 0.22$. Freeze it there, and re-fit everything else, letting all the other parameters move to their best values given that $\kappa$. Record how much worse the fit got. Do that across a grid of $\kappa$ values and you have traced a curve. The 95% CI is every value whose fit is within 1.92 log-likelihood units of the best one. (1.92 is not arbitrary. It is half of $\chi^2_{1,\,0.95} = 3.84$, the same LRT arithmetic, run one parameter at a time.) In R, confint() on an mle() fit does exactly this for you.

The words "re-fit everything else" are the whole idea, and the commonest way to get this wrong is to skip them. If you hold the other parameters frozen at their best-fit values while you slide $\kappa$, you are not profiling, you are taking a slice, and the interval comes out too narrow. Day 4's lab shows this happening on purpose. A grid over $\kappa$ with $L_\infty$ nailed to 100 reports $\kappa \approx 0.15$, while letting $L_\infty$ move returns $\kappa = 0.259$ on the same fish.

The Bayesian echo. Both tools exist because maximum likelihood gives you a point and then has to work to describe the uncertainty around it. In the afternoon you get a whole posterior distribution, so the "interval" is just a question you ask of the draws you already have. That is the same trade-off you will meet again with derived quantities such as MSY, awkward for likelihood and free for Bayes.

From prior to posterior, in pictures

Bayes does something the MLE cannot. It gives you a distribution over the parameter, not just a point. Three ingredients:

Prior $p(\theta)$: what you believed about $\theta$ before seeing the data. Often "I don't know much", a wide, weakly informative distribution.
Likelihood $p(y \mid \theta)$: exactly the same likelihood as the morning.
Posterior $p(\theta \mid y)$: what you believe now. By Bayes' rule, $\text{posterior} \propto \text{prior} \times \text{likelihood}$.

If the prior is uniform and the likelihood is peaked, the posterior looks like the likelihood, and the posterior mode is the MLE. If the prior is strong, it pulls (or "shrinks") the posterior toward what you previously believed. As you get more data, the likelihood gets sharper and dominates, and the prior becomes less and less important.

Priors are a feature, not a bug. They let you encode real prior knowledge (e.g. "survival in this species is biologically plausible only between 0.2 and 0.8") and they regularize unstable estimates. They are also the thing reviewers question hardest, so always justify them.

The big Bayesian flip

Frequentist statistics computes things like $P(\text{data} \mid H_0)$, the famous p-value, and then asks you to translate that into a verdict on a hypothesis. Bayes computes the thing you actually wanted in the first place:

Frequentist: "If the null hypothesis is true, the data I observed (or anything more extreme) would happen in 3% of repetitions."

Bayesian: "Given the data I observed, there is a 97% probability that the effect is positive."

Notice that the Bayesian sentence puts the probability where you wanted it all along, on the hypothesis. That is the payoff for buying into the prior + posterior framework.

What does MCMC actually do?

For real models (anything with more than one or two parameters, or a non-conjugate likelihood) we cannot write the posterior down in closed form. Markov chain Monte Carlo is a clever algorithm that generates a long sequence of random samples whose long-run frequency is the posterior. You don't compute the posterior; you draw from it.

The simplest version (Metropolis):

Start somewhere.
Propose a small random move.
If the move improves the posterior, accept. If it makes the posterior worse, accept with some probability (the ratio of new to old posterior).
Otherwise stay put.
Repeat thousands of times. Throw out the first chunk ("burn-in"). The rest is a sample from the posterior.

Once you have those samples, every summary you want is a one-liner:

posterior mean of $r$? mean(r_samples)
95% credible interval? quantile(r_samples, c(.025, .975))
probability $r > 0.22$? mean(r_samples > 0.22)
posterior of $K/2$? Compute $K/2$ for every sample of $K$. That is the posterior of $K/2$.

Where does Nimble fit in?

You don't want to write Metropolis by hand for every new model. Nimble (and JAGS, Stan, etc.) take a BUGS-style description of your model, priors and likelihood, and automatically:

build the computational graph,
pick a sensible sampler for each parameter,
run several independent chains from different starting values (nimbleMCMC() runs them one after another by default; true parallel execution takes extra setup),
return posterior samples for whatever quantities you ask to monitor.

Why Nimble specifically (instead of JAGS)? Nimble's BUGS dialect is nearly identical to JAGS (most JAGS code runs unchanged inside a nimbleCode({...}) block), but it is much more extensible. You can write custom samplers, custom distributions, and compile to fast C++ all from R. Whole subfields (ecological state-space, capture-recapture, spatial models) have migrated to Nimble for this reason.

JAGS → Nimble in three changes:

cat("model {…}", file="x.txt") → nimbleCode({ … })
One data list → two lists: constants (sample sizes, indices, fixed covariates) and data (the observed random variables).
jags.model() + update() + coda.samples() → one call to nimbleMCMC(code=, constants=, data=, inits=, monitors=, nchains=, nburnin=, niter=).

The model body (dnorm, dunif, dbeta, dbinom, dpois, loops, derived quantities) translates verbatim.

Common confusions to avoid

"The MLE is unbiased." Not always! MLEs are consistent (right answer with infinite data) but can be biased in finite samples. Variance components in mixed models, for example, have MLEs that are systematically too small, which is why REML exists.

"A 95% confidence interval has a 95% chance of containing the true parameter." No. A 95% CI is a procedure that, in repeated sampling, would cover the true parameter 95% of the time. Once you have a specific CI from a specific dataset, it either contains the parameter or it doesn't, and no probability statement applies. The Bayesian credible interval does support that natural reading; that's part of why people switch.

"My MCMC has run 50,000 iterations, so it must have converged." Iteration count alone tells you nothing. Look at trace plots (multiple chains should overlap), $\hat R$ values (should be very close to 1), and effective sample sizes (should be at least a few hundred for the parameters you care about). MCMC failure is silent, so diagnose every run.

"The prior is subjective, the likelihood is objective." Both are modeling choices. Picking Gaussian errors instead of t-distributed errors is a modeling choice; picking a uniform prior over $K$ instead of a log-uniform prior is a modeling choice. Justify both.

Looking ahead to Day 5

Tomorrow we take the Bayesian machinery you just learned and apply it to harder, more realistic models, hierarchical (random-effects) models, occupancy models, capture-recapture, and state-space models. Same Nimble syntax, deeper biology. The morning's maximum likelihood and the afternoon's conjugacy and logistic-growth fit you did today are the building blocks; everything tomorrow plugs together pieces you already understand.

Things to keep in your back pocket going into Day 5:

How to write a nimbleCode block.
The split between constants and data.
What MCMCsummary and MCMCtrace tell you.
How to compute any tail probability from posterior samples.
The difference between a BCI and an HPDI.