Bayesian computation in R
Most ecologists meet Bayesian inference through a modelling language: you write the model, press go, and read the summary. That works until it does not, and when it stops working the error message is about the sampler rather than about the ecology. A model that is perfectly well specified can sample badly. A chain that looks calm can be exploring a fraction of the posterior. An approximation can be excellent for the parameters and useless for the quantity you actually wanted.
These tutorials build the machinery by hand in base R so that those symptoms become readable. Nothing here is a substitute for the modelling languages, and the tutorials say so. The claim is narrower: once you have written a Metropolis step, watched a leapfrog trajectory get repaired by its accept step, and seen a funnel defeat a sampler that was fine a moment earlier, the diagnostics stop being ritual and start being information. Every model here is simulated, so the posterior the sampler should have found is known while it runs.
The sampler, written out
Start where the machinery is simplest and the answer is checkable against something exact.
- Metropolis-Hastings from scratch is the whole idea in a handful of lines: propose, compare, accept or stay, and then the interval the chain gives you set against the one the asymptotics would have handed over.
- Gibbs sampling with conjugate updates derives the full conditionals for a mean and a variance and checks the sampled posterior against the closed form, which is the only honest way to know a sampler is right rather than plausible.
- MCMC convergence diagnostics from scratch builds the split potential scale reduction factor and effective sample size from their definitions, then runs them against a chain that passes the eye test and should not.
- Bayesian hierarchical models with MCMC samples a partially pooled model by hand and follows the posterior predictive check through, so that shrinkage stops being a slogan and becomes something you can watch happen.
- Bayesian model comparison: WAIC, LOO and DIC computes the out-of-sample scores from the log-likelihood matrix, and is clear about which of them are estimating the same thing and where they part company.
- Choosing priors: flat is not uninformative pushes a flat prior through the transformation that matters and shows how strongly it can speak on the scale you actually care about.
- Bayesian logistic regression under separation takes the case where the maximum likelihood estimate runs off to infinity and shows what a weakly informative prior is really doing when it rescues the fit.
- Bayesian occupancy with latent-state MCMC introduces the move the rest of this page depends on: stop integrating the hidden state out, sample it instead, and the awkward likelihood becomes a sequence of easy updates.
Gradients, approximations and the diagnostics that matter
This is the layer that explains what a modern sampler is doing, and what its warnings are trying to tell you.
- Hamiltonian Monte Carlo from scratch writes the leapfrog integrator, shows the accept step repairing the integration error rather than hiding it, and prices the whole thing in the currency that counts, which is effective samples per gradient evaluation.
- The Laplace approximation in R puts a Gaussian at the mode, collects the marginal likelihood that falls out for free, and then works through the named ways this cheap and widely used approximation misleads.
- Reparameterisation and Neal’s funnel writes one model two ways and finds the crossover where the standard advice reverses, which is the most useful thing to know about a geometry that has stalled.
- Checking a Bayesian computation is the audit: an exact benchmark, simulation-based calibration, what a divergence is evidence of, and the Monte Carlo standard error that tells you when to stop running.
The latent-variable move on classic ecological models
Each of these belongs to its own ecological topic as much as to this one. They are gathered here because they are the same computational idea applied four times, and reading them together makes the idea portable.
- Bayesian Cormack-Jolly-Seber survival samples the latent alive state directly, which turns the usual likelihood gymnastics into bookkeeping and makes the survival and recapture parameters separable by eye.
- Bayesian N-mixture abundance puts latent site abundance on a discrete grid and updates it conjugately, so the notorious identifiability problem of the model becomes visible in the posterior rather than hidden in an optimiser.
- Bayesian state-space movement recovers a whole latent path with forward filtering and backward sampling, which is the cleanest demonstration on the site of why sampling the hidden quantity beats integrating it away.
- Bayesian closed capture-recapture handles an unknown population size by parameter-expanded data augmentation, turning a parameter whose dimension is unknown into one that is merely large.
When there is no likelihood to write
Some ecological models are easy to simulate forward and hopeless to write a density for. That is not a reason to abandon inference, but it changes what the diagnostics can promise.
- Approximate Bayesian computation from scratch fits such a model by comparing simulated data with the data in hand, and is explicit about what the tolerance is buying and what it is quietly costing.
- Choosing summary statistics for ABC makes the central point of the method: in ABC the summaries are the data, and choosing them badly costs more than any tolerance setting.
- ABC-MCMC and sequential ABC samplers measures what a better sampler can and cannot win, and finds the ceiling set by the ratio of prior volume to posterior volume.
- Checking an ABC analysis draws the line the method makes easy to blur: rank calibration audits the sampler, not the model, and it will happily pass an analysis that has learned nothing at all.
What these tutorials assume
Base R and one plotting package throughout. Every sampler, every diagnostic and every approximation is written out rather than called, and each is checked against something whose answer is known independently, either a conjugate posterior, a closed form, or a simulation whose truth was chosen before the fitting started. Familiarity with likelihood and with the models being fitted is assumed; familiarity with any probabilistic programming language is not. Where a real modelling language would be the sensible production choice, the tutorials name it and move on.
Where this connects
The computation exists to serve models that live elsewhere on this site, and reads better alongside them.
- Occupancy modelling in R is where the latent presence state comes from, and where the likelihood version of the same model is fitted and diagnosed.
- Capture-recapture in R supplies the survival and abundance models that two of the latent-state tutorials rewrite in Bayesian form.
- Movement ecology in R is the home of the state-space path model, including the tag error and behavioural state questions the sampler is ultimately serving.
- Population models in R provides the demographic machinery that integrated and hierarchical population analyses put a posterior on.
- Reproducible workflows in R matters more here than almost anywhere else on this site: a sampler is stochastic, slow and version-sensitive, so the seed, the lockfile and the pipeline are what let anyone reproduce a posterior rather than merely a similar one.