Bayesian Inference

Agricultural Statistics & Data Science · Statistical theory core · also Bayesian Statistics, Bayesian Methods (PhD)

11units
11chapters
5,343words
100practice questions

Inference through prior and posterior: conjugate analysis, prior choice, credible intervals, Bayes factors and decision theory, hierarchical models, and computation by Markov-chain Monte Carlo, with agricultural uses in variety trials, genomic prediction and small-area estimation.

In the exams: Not an explicit unit in AICE-JRF Ag Stat; PhD course in most programmes; GATE ST Bayes' theorem (light); UPSC ISS decision theory (present)

Syllabus

  1. Bayes' theorem for parameters: prior, likelihood and posterior
  2. Conjugate families
  3. Non-informative and weakly informative priors
  4. Posterior summaries and credible intervals
  5. Bayesian hypothesis testing and Bayes factors
  6. Decision theory: loss functions and Bayes risk
  7. Hierarchical (multilevel) models
  8. Markov chain Monte Carlo: Gibbs sampling and Metropolis-Hastings
  9. Convergence diagnostics
  10. Posterior predictive checks and model comparison (DIC, WAIC)
  11. Empirical Bayes

Study manual

Bayesian inference treats unknown parameters as random quantities with probability distributions that express what is known about them. A prior distribution states knowledge before the data; the likelihood says how probable the data are for each parameter value; Bayes' theorem combines them into a posterior distribution, from which every inference follows: estimates, intervals, probabilities of hypotheses, predictions and decisions. The approach answers directly the questions practitioners ask, such as 'what is the probability that the new variety out-yields the check by at least 5%?', which frequentist p-values do not.

Agricultural research has become one of its main fields of application. Breeders use Bayesian and empirical-Bayes methods (BLUP and its genomic versions) to predict breeding values; variety-testing programmes borrow strength across trials and seasons through hierarchical models; crop-yield estimation for small areas combines survey and satellite data with Fay-Herriot-type shrinkage; plant-protection services update pest-risk forecasts as monitoring data arrive; and agricultural economists build decision models in which uncertainty about yields and prices enters through posterior distributions.

This manual covers Bayes' theorem for parameters; conjugate families; non-informative and weakly informative priors; posterior summaries and credible intervals; Bayes factors and testing; decision theory; hierarchical models; Markov chain Monte Carlo (Gibbs sampling and Metropolis-Hastings); convergence diagnostics; posterior predictive checks and information criteria (DIC, WAIC); and empirical Bayes. It builds on Probability Theory (conditional probability and Bayes' theorem for events), Probability Distributions and Sampling Distributions (the beta, gamma, normal and exponential-family forms on which conjugacy rests) and Statistical Inference I and II (likelihood, sufficiency, Fisher information and the frequentist tests with which Bayesian answers are compared).

Bayesian inference is not an explicit unit of the ICAR AICE-JRF Agricultural Statistics paper, but it is a PhD course in most agricultural-statistics programmes and appears in GATE Statistics (Bayes' theorem and posterior computations), in the decision-theory part of the Indian Statistical Service syllabus, and increasingly in research-institute interviews and in the methods sections of breeding and survey papers.

1. Prior, likelihood and posterior

Bayes' theorem for parameters

Let θ be an unknown parameter with prior density π(θ), and let the data x have likelihood f(x | θ). Bayes' theorem gives the posterior density π(θ | x) = f(x | θ)π(θ)/m(x), where m(x) = ∫ f(x | θ)π(θ) dθ is the marginal (prior predictive) density of the data. Since m(x) does not depend on θ, the posterior is proportional to likelihood times prior: π(θ | x) ∝ f(x | θ)π(θ). For discrete θ the integral becomes a sum, exactly as in Bayes' theorem for events (see Probability Theory).

The posterior depends on the data only through the likelihood, and hence only through a sufficient statistic. It obeys the likelihood principle: two experiments giving proportional likelihoods give the same posterior, whatever their sampling plans or stopping rules. That is a sharp contrast with frequentist p-values, which depend on outcomes that could have occurred but did not.

Bayesian updating is sequential. With independent data sets x₁ and x₂, the posterior after x₁ serves as the prior for x₂, and the final posterior is the same as from analysing both together. A plant-protection service can therefore update the posterior probability of a pest outbreak week by week as trap catches arrive.

A discrete example and the role of the prior

A seed lot is either good (germination 0.90) or poor (0.70); from past records, 80% of lots are good. Ten seeds are tested and 7 germinate. The likelihoods are C(10, 7)(0.9)⁷(0.1)³ = 0.0574 for a good lot and C(10, 7)(0.7)⁷(0.3)³ = 0.2668 for a poor one. The posterior probability that the lot is good is 0.8 × 0.0574/(0.8 × 0.0574 + 0.2 × 0.2668) = 0.0459/0.1093 = 0.462. The data have moved the probability from 0.80 to 0.46: the evidence favoured 'poor' by a likelihood ratio of about 4.6, which more than offset prior odds of 4 to 1 for 'good'.

The prior matters most when data are few and matters less as data accumulate: under regularity conditions the posterior concentrates around the true value and is approximately normal with mean near the MLE and variance near the inverse Fisher information (the Bernstein-von Mises theorem), whatever smooth prior with positive density at the true value was used. In practice the prior should be stated, justified from previous trials, expert knowledge or deliberately vague choices, and checked by sensitivity analysis: does the conclusion change under other reasonable priors?

The posterior predictive distribution for a future observation y is p(y | x) = ∫ f(y | θ)π(θ | x) dθ. It accounts for both the randomness of y and the uncertainty about θ, and so is wider than the plug-in distribution f(y | θ̂). Predictions for a farmer's next season, or for the yield of a new trial site, should use it.

For the beta-binomial example in Chapter 2, the probability that the next seed germinates is simply the posterior mean of θ; for the normal model it is a normal distribution whose variance is σ² plus the posterior variance of μ. These predictive distributions are also the basis of the model checks in Chapter 10: a model whose predictions do not resemble the data it was fitted to cannot be trusted to predict new data.

Remember

2. Conjugate families

Beta-binomial and gamma-Poisson

A family of priors is conjugate to a likelihood if the posterior belongs to the same family. Conjugacy gives closed-form posteriors whose parameters update by simple arithmetic, and it exists for every exponential-family likelihood (see the distributions manual).

Binomial data with a beta prior: if θ ~ beta(a, b) and x successes are seen in n trials, the posterior is beta(a + x, b + n − x). The prior acts like a + b earlier trials with a successes. The posterior mean is (a + x)/(a + b + n), a weighted average of the prior mean a/(a + b) and the sample proportion x/n with weights (a + b) and n. With a uniform beta(1, 1) prior and 7 germinations in 10 seeds, the posterior is beta(8, 4) with mean 8/12 = 0.667. The posterior predictive distribution of successes in m further trials is beta-binomial.

Poisson counts with a gamma prior: if λ ~ gamma(α, β) (shape α, rate β) and counts x₁, ..., xₙ are observed, the posterior is gamma(α + Σxᵢ, β + n), with mean (α + Σxᵢ)/(β + n). The prior acts like β earlier observation periods containing α events. The posterior predictive distribution of a future count is negative binomial, which is why over-dispersed counts arise naturally when Poisson rates vary.

Other standard pairs: the exponential or gamma likelihood with a gamma prior on the rate; the multinomial with a Dirichlet prior, whose parameters add to the category counts; the geometric and negative binomial with a beta prior; and the uniform (0, θ) with a Pareto prior.

Normal models

Normal data with known variance σ² and a normal prior on the mean: if μ ~ N(μ₀, τ₀²) and x₁, ..., xₙ ~ N(μ, σ²), then μ | x ~ N(μₙ, τₙ²) with precision 1/τₙ² = 1/τ₀² + n/σ² and mean μₙ = τₙ²(μ₀/τ₀² + nx̄/σ²). Precisions add, and the posterior mean is a precision-weighted average of the prior mean and the sample mean. Example: prior belief about a variety's mean yield is N(4.0, 0.5²) t/ha; nine plots with σ = 0.6 give x̄ = 4.6. Prior precision 4, data precision 9/0.36 = 25, posterior precision 29, posterior variance 0.0345 (SD 0.186), posterior mean (4 × 4.0 + 25 × 4.6)/29 = 4.517 t/ha.

With unknown variance, the conjugate prior for (μ, σ²) is normal-inverse-gamma (equivalently normal-gamma on the precision): σ² has an inverse-gamma marginal and μ given σ² is normal. The marginal posterior of μ is then a Student-t distribution, and with the standard non-informative prior π(μ, σ²) ∝ 1/σ² it reproduces the frequentist t-interval exactly, a useful bridge between the two schools. For the normal linear model the same structure gives normal-inverse-gamma posteriors for the coefficients and the error variance.

Conjugate priors are convenient but not obligatory. With modern computation (Chapter 8) any prior can be used, and the choice should reflect knowledge, not algebra. Mixtures of conjugate priors keep tractability while allowing more flexible prior shapes, such as a prior that puts some weight on 'no effect' and some on a range of effects.

Remember

3. Non-informative and weakly informative priors

Flat, Jeffreys and reference priors

When little is known, or when a 'default' analysis is wanted, priors are chosen to let the data dominate. A flat (uniform) prior is the simplest choice, but 'flat' depends on the parameterisation: a uniform prior on a probability p is not uniform on the log-odds log[p/(1 − p)]. Many flat priors on unbounded parameters are improper (they do not integrate to 1). An improper prior is acceptable only if the resulting posterior is proper, which must be checked; for example, a flat prior on the log of a variance component in a hierarchical model can give an improper posterior.

Jeffreys' prior, π(θ) ∝ √I(θ), with I the Fisher information, is invariant under reparameterisation: the same rule applied on any scale gives equivalent posteriors. It gives beta(1/2, 1/2) for a binomial proportion, π(λ) ∝ λ^(−1/2) for a Poisson mean, a flat prior for a normal mean, and π(σ) ∝ 1/σ (flat on log σ) for a normal standard deviation. For multiparameter problems Jeffreys' rule can behave poorly, and reference priors (Bernardo, Berger) refine it by treating parameters of interest and nuisance parameters differently.

Default priors often give credible intervals with good frequentist coverage (matching priors), which reassures users of both schools. They are not 'objective' in an absolute sense: every prior is a modelling choice.

Weakly informative priors

A weakly informative prior deliberately contains less information than is actually available but enough to rule out absurd values and stabilise computation. Examples are a normal prior with a wide but finite scale on a regression coefficient measured in standardised units, and a half-Cauchy or half-normal prior on a standard deviation in a hierarchical model, which keeps the posterior away from implausibly large values while allowing values near zero; Gelman (2006) argued for half-Cauchy and related priors on group-level standard deviations instead of the inverse-gamma(ε, ε) priors that were once common and can be surprisingly informative when the group variance is small.

In agriculture, prior information is often genuine and substantial: heritabilities of a trait from many earlier studies, typical yield responses to nitrogen, coefficients of variation of trials at a station. Using it through informative priors, especially in small trials, improves estimates; the price is the obligation to justify the prior and to show how much the conclusion depends on it. A prior effective sample size, such as a + b for a beta prior, helps readers judge how strong a prior is relative to the data.

Prior elicitation from experts asks for quantities they can judge, such as a median and a 90% range for a yield response, and fits a distribution to them; experience shows that experts tend to be over-confident, so elicited priors are often widened.

Remember

4. Posterior summaries and credible intervals

Point summaries

The posterior is the full answer; summaries are for communication. The posterior mean, median and mode are the Bayes estimates under squared-error, absolute-error and 0-1 loss respectively (Chapter 6). For a beta(8, 4) posterior the mean is 8/12 = 0.667, the mode (a − 1)/(a + b − 2) = 7/10 = 0.70, which equals the MLE under a uniform prior, and the median lies between them. For skewed posteriors, such as those for variances or for rare-event probabilities, the three can differ substantially, and the median is often the most stable summary. Posterior standard deviations play the role of standard errors.

Posterior probabilities of statements are often the most useful summary: P(θ > 0 | x) for a treatment effect, P(new variety yields at least 5% more than the check | x), or P(infestation above the economic threshold | x). These are direct answers to practical questions and can feed straight into decisions.

Credible intervals

A 100(1 − α)% credible interval is any interval C with P(θ ∈ C | x) = 1 − α. Two constructions are standard. The equal-tailed interval runs from the α/2 to the 1 − α/2 posterior quantile; it is invariant under monotone transformations of θ (the interval for log θ is the log of the interval for θ). The highest posterior density (HPD) interval contains the values with the highest posterior density, so every point inside is more probable than every point outside; it is the shortest interval with the given probability, but it is not invariant to transformations, and for multimodal posteriors it can be a union of intervals. For symmetric unimodal posteriors the two coincide.

The interpretation differs from a confidence interval: a 95% credible interval contains θ with posterior probability 0.95 given the data and the model, which is the statement users usually want. In large samples with diffuse priors, credible and confidence intervals nearly coincide numerically (Bernstein-von Mises). In small samples, or with informative priors, they can differ, and the difference is the effect of the prior information.

Example: for the normal-model yield example of Chapter 2, the posterior is N(4.517, 0.186²), so the 95% credible interval is 4.517 ± 1.96 × 0.186 = 4.15 to 4.88 t/ha, and P(μ > 4.2 | x) = P(Z > (4.2 − 4.517)/0.186) = P(Z > −1.70) ≈ 0.96.

Credible intervals computed from MCMC output (Chapter 8) are read from the sorted sample: the 2.5% and 97.5% sample quantiles for the equal-tailed interval, or the shortest interval containing 95% of the draws for the HPD interval, with Monte Carlo error that depends on the effective sample size.

Remember

5. Bayes factors and hypothesis testing

Posterior odds and the Bayes factor

To compare hypotheses H₀ and H₁ (or models M₀ and M₁), Bayesian inference assigns prior probabilities and computes posterior probabilities. The posterior odds equal the prior odds times the Bayes factor: P(H₁ | x)/P(H₀ | x) = [P(H₁)/P(H₀)] × B₁₀, with B₁₀ = m₁(x)/m₀(x), where mᵢ(x) = ∫ f(x | θ, Hᵢ)πᵢ(θ) dθ is the marginal likelihood under Hᵢ. For simple hypotheses the Bayes factor is just the likelihood ratio; for composite ones it averages the likelihood over the prior within each hypothesis, rather than maximising it as the likelihood-ratio test does.

Scales for interpreting Bayes factors exist; Kass and Raftery (1995) suggested that B₁₀ between 1 and 3 is barely worth mentioning, 3 to 20 positive, 20 to 150 strong and above 150 very strong evidence for H₁, while Jeffreys' original scale uses different cut-points. Such labels are conventions; the Bayes factor itself, and the posterior probability it implies under stated prior odds, are what should be reported.

In the seed-lot example of Chapter 1 the Bayes factor for 'poor' against 'good' is 0.2668/0.0574 = 4.65; with prior odds of 1 to 4 for 'poor' the posterior odds are 1.16, i.e. P(poor | x) = 0.54.

Properties and pitfalls

Bayes factors automatically penalise complexity: a model with diffuse priors on extra parameters spreads its prior predictive probability thinly and receives a lower marginal likelihood unless the data demand the extra parameters (an 'Occam's razor' effect). The Schwarz criterion gives a large-sample approximation: 2 log B₁₀ ≈ BIC₀ − BIC₁, so a BIC difference of 6 corresponds to a Bayes factor of about e³ ≈ 20.

This sensitivity to the prior is also the main pitfall. For testing a point null θ = θ₀ against a composite alternative, the Bayes factor depends heavily on how spread out the prior under H₁ is; making it ever more diffuse drives B₁₀ towards zero, so H₀ is favoured whatever the data. Improper priors cannot be used under H₁ for this reason (the arbitrary constant does not cancel). Lindley's paradox is the extreme form: with a large sample, a result significant at p = 0.01 in a frequentist test can yield a Bayes factor favouring H₀, because the evidence is judged against a diffuse alternative. Intrinsic and fractional Bayes factors, and reporting Bayes factors for a range of priors, are the remedies.

For interval hypotheses, such as 'the new variety is at least 5% better', the posterior probability of the interval is the direct and robust answer, and is often more useful than a point-null test. For estimation of an effect that is certainly not exactly zero, credible intervals usually say more than any test.

In variety testing, posterior probabilities of superiority are now reported in some programmes alongside or instead of significance tests: the probability that a candidate out-yields the check across the target region, computed from a hierarchical model of multi-location data, is a statement that a release committee can weigh directly against the costs of releasing a variety that turns out worse.

Remember

6. Decision theory

Loss, risk and Bayes rules

Statistical decision theory formalises choosing an action a when the state of nature θ is uncertain. A loss function L(θ, a) states the cost of action a when θ is true. For a decision rule δ(x), the frequentist risk is R(θ, δ) = E_θ[L(θ, δ(X))], an average over data for fixed θ. The Bayes risk is r(π, δ) = ∫ R(θ, δ)π(θ) dθ, the risk averaged over the prior. A Bayes rule minimises the Bayes risk, and it can be found by minimising, for each observed x, the posterior expected loss ∫ L(θ, a)π(θ | x) dθ.

For estimation: squared-error loss (θ − a)² gives the posterior mean; absolute-error loss |θ − a| the posterior median; 0-1 loss (in the limit) the posterior mode; and asymmetric losses give posterior quantiles. Under the loss k₁(θ − a) for under-estimation and k₂(a − θ) for over-estimation, the Bayes estimate is the k₁/(k₁ + k₂) posterior quantile. An agronomist who loses three times as much from under-applying fertiliser as from over-applying would choose the 75% quantile of the posterior for the optimal rate, not its mean.

For testing: with loss c₀ for wrongly rejecting H₀ (Type I) and c₁ for wrongly accepting it (Type II), and zero loss for correct decisions, the Bayes rule rejects H₀ when P(H₀ | x)c₀ < P(H₁ | x)c₁, i.e. when the posterior probability of H₀ is below c₁/(c₀ + c₁). With equal losses, reject when P(H₀ | x) < 1/2.

Admissibility, minimax and practice

A rule is inadmissible if another rule has risk no larger for every θ and smaller for some θ. Bayes rules with proper priors are admissible under mild conditions (when the Bayes rule is unique or the risk functions are continuous). A minimax rule minimises the maximum risk over θ; it is a cautious criterion, and a Bayes rule whose risk is constant in θ is minimax. Stein's discovery that the sample mean of three or more normal means is inadmissible under total squared-error loss (Chapter 11) showed that intuitive estimators can be beaten by shrinkage, a result with deep links to empirical Bayes.

In agricultural practice, decision analysis is used for spray-or-not decisions given a posterior distribution of pest density and an economic threshold, for choosing fertiliser rates under yield-response uncertainty, for deciding whether to release a variety given posterior probabilities of superiority across environments, and for valuing additional information (the expected value of sample information tells whether another season of trials is worth its cost). Losses are economic quantities; they vary with prices and farm circumstances, and an analysis should state the losses it assumed.

Remember

7. Hierarchical models

Exchangeability and partial pooling

Many agricultural data sets have a grouped structure: plots within trials, trials within locations, animals within herds, farms within villages, varieties within a testing programme. A hierarchical (multilevel) model gives each group its own parameter and treats the group parameters as draws from a common population distribution whose parameters (hyperparameters) are also estimated. For variety means across trials, yᵢ ~ N(θᵢ, σᵢ²) for trial i, θᵢ ~ N(μ, τ²), with priors on μ and τ.

The justification is exchangeability: if nothing distinguishes the groups before seeing the data, their parameters are exchangeable, and de Finetti's theorem says an exchangeable sequence can be represented as conditionally independent given some common parameter. Hierarchical models are the Bayesian expression of the random-effects models in the mixed-models subject.

The posterior mean of each group parameter, given the hyperparameters, is a compromise between the group's own estimate and the overall mean: E(θᵢ | y, μ, τ) = (1 − Bᵢ)yᵢ + Bᵢμ, with shrinkage factor Bᵢ = σᵢ²/(σᵢ² + τ²). Groups with noisy estimates (large σᵢ², few plots) are pulled strongly towards μ; precise groups keep nearly their own estimates. This is partial pooling: between no pooling (each group alone) and complete pooling (one common value).

A worked example and practical points

Suppose a trait's true trial means vary with τ = 0.4 t/ha around μ = 5.0 t/ha. Trial A has mean 6.2 with standard error 0.6 (few plots); trial B has mean 6.2 with standard error 0.2. For A, B_A = 0.36/(0.36 + 0.16) = 0.692 and the shrunken estimate is 0.308 × 6.2 + 0.692 × 5.0 = 5.37; for B, B_B = 0.04/0.20 = 0.20 and the estimate is 0.8 × 6.2 + 0.2 × 5.0 = 5.96. The same observed mean earns different credibility according to its precision, which is exactly how breeders treat single high-yielding results from small trials.

Practical points. The hyperparameters, especially τ, need priors that are not unintentionally informative (half-Cauchy or half-normal priors on τ are common). With few groups, τ is poorly estimated, and full Bayesian analysis propagates that uncertainty where plug-in methods do not. Hierarchical models extend naturally to several levels (plots, trials, years), to genotype-by-environment structures, to spatial correlation among neighbouring fields, and to non-normal responses (binomial germination, Poisson counts) through generalised linear mixed models. The connection to BLUP is direct: best linear unbiased predictors of random effects are the posterior means under normal hierarchical models with known variance components.

Remember

8. Markov chain Monte Carlo

Why simulation, and the Gibbs sampler

For most realistic models the posterior has no closed form and its normalising constant m(x) is an intractable integral. Markov chain Monte Carlo (MCMC) constructs a Markov chain whose stationary distribution is the posterior; after an initial burn-in, the chain's states are draws (dependent, but valid for averaging) from the posterior. Any posterior summary is then estimated by the corresponding sample summary: means by averages, quantiles by sample quantiles, probabilities by proportions. MCMC needs the posterior only up to a constant, which is why it unlocked Bayesian practice from about 1990.

The Gibbs sampler updates one parameter (or block) at a time by drawing from its full conditional distribution given the current values of all others: θ₁ from π(θ₁ | θ₂, ..., x), then θ₂ from π(θ₂ | θ₁, θ₃, ..., x), and so on. With conditionally conjugate priors these full conditionals are standard distributions. For the normal model with unknown mean and variance and semi-conjugate priors, the sampler alternates a normal draw for μ given σ² and an inverse-gamma draw for σ² given μ. Data augmentation adds latent variables (missing plots, latent classes, latent normal variables behind probit responses) to make the full conditionals simple. Gibbs sampling mixes slowly when parameters are strongly correlated in the posterior, because each coordinate moves only a little given the others.

Metropolis-Hastings and modern samplers

The Metropolis-Hastings (MH) algorithm needs no full conditionals. From the current value θ, propose θ* from a proposal density q(θ* | θ) and accept it with probability α = min{1, [π(θ* | x)q(θ | θ*)]/[π(θ | x)q(θ* | θ)]}; otherwise stay at θ. The normalising constant cancels in the ratio. With a symmetric proposal, such as a random walk θ* = θ + ε with ε ~ N(0, s²), the acceptance probability reduces to min{1, π(θ* | x)/π(θ | x)} (the original Metropolis algorithm). The chain satisfies detailed balance with respect to the posterior, which ensures that the posterior is its stationary distribution.

Tuning matters. Proposals that are too small are almost always accepted but move slowly; proposals that are too large are rarely accepted. Theory for random-walk proposals in high dimensions (Roberts, Gelman and Gilks, 1997) suggests acceptance rates near 0.234 are efficient, with higher rates (around 0.44) for one-dimensional updates. The Gibbs sampler is the special case of MH in which the proposal is the full conditional and every proposal is accepted.

Modern software uses Hamiltonian Monte Carlo (HMC) and its adaptive version, the no-U-turn sampler, which use gradients of the log posterior to make long, efficient moves and cope well with correlated parameters; tools built on them are now standard for hierarchical models in agricultural research. Whatever the sampler, its output must be checked before use (Chapter 9).

Remember

9. Convergence diagnostics

Has the chain converged, and how much information does it carry?

MCMC output is only useful if the chain has reached its stationary distribution and explored it adequately. No diagnostic can prove convergence; diagnostics can only reveal failures. Standard practice is to run several chains (four is common) from dispersed starting values, discard an initial burn-in, and examine the following.

Trace plots of each parameter against iteration should look like stationary 'hairy caterpillars', with all chains overlapping; trends, long excursions or chains stuck in different regions signal problems. The Gelman-Rubin potential scale reduction factor R̂ compares between-chain and within-chain variance; values near 1 indicate that the chains agree. The traditional rule was R̂ < 1.1; more recent recommendations (Vehtari and colleagues, 2021) use rank-normalised R̂ with a threshold of 1.01. Geweke's diagnostic compares means from the early and late parts of a single chain with a z-score.

Because successive MCMC draws are autocorrelated, N draws carry less information than N independent ones. The effective sample size is ESS = N/(1 + 2Σρₖ), where ρₖ is the lag-k autocorrelation. With N = 10,000 draws and a sum of autocorrelations of 4.5, ESS = 10,000/10 = 1,000. The Monte Carlo standard error of a posterior mean estimate is the posterior SD divided by √ESS: with posterior SD 0.5 and ESS 1,000, MCSE = 0.016, small enough for most purposes. An ESS of a few hundred per parameter is commonly considered adequate for means and central intervals, more for tail quantiles.

Remedies for poor mixing

When diagnostics fail, the options are: run longer; reparameterise (centring or non-centring hierarchical models, which changes posterior correlations dramatically when group-level variances are small); block correlated parameters so they are updated together; improve the proposal; switch to HMC; or simplify an over-parameterised model whose parameters are not identified by the data. Thinning (keeping every k-th draw) reduces storage but does not increase information; it is not a cure for autocorrelation. Divergent transitions reported by HMC samplers indicate regions of the posterior the sampler cannot explore reliably and should be investigated, not ignored. Reporting R̂, ESS and the number of chains and iterations in a paper's methods section lets readers judge the computation.

A practical example: a hierarchical model for 40 varieties tested in 12 trials is fitted with four chains of 2,000 iterations each after 1,000 warm-up iterations. All R̂ values are below 1.01 except that for the between-trial standard deviation, which is 1.08, and its ESS is only 90. The trace plot shows one chain wandering near zero. The remedy is a non-centred parameterisation of the trial effects (writing each effect as the standard deviation times a standard normal variable), after which R̂ falls to 1.00 and the ESS rises above 1,000. Reporting the posterior of the variety effects from the first run would have understated uncertainty about the trial-to-trial variation, which drives how strongly variety means are shrunk.

Remember

10. Posterior predictive checks and model comparison

Posterior predictive checks

A model can fit the parameters well yet be wrong. Posterior predictive checking simulates replicated data y_rep from the posterior predictive distribution and compares them with the observed data through test quantities T(y), such as the proportion of zeros in pest counts, the maximum yield, the variance-to-mean ratio, or the lag-1 autocorrelation of residuals along a field. If the observed T(y) is extreme relative to the replicated T(y_rep), the model misses that feature. The posterior predictive p-value is P(T(y_rep) ≥ T(y) | y); values near 0 or 1 signal misfit. For example, if a Poisson model for aphid counts produces replicated data with far fewer zeros than observed, a zero-inflated or negative binomial model is needed.

Posterior predictive p-values are not uniformly distributed under a correct model (they tend to cluster near 0.5, because the data are used twice), so they are diagnostics, not formal tests. Graphical comparisons of observed and replicated data are often more informative than any single number.

Information criteria and cross-validation

For comparing models by predictive accuracy, Bayesian information criteria estimate out-of-sample deviance. The deviance information criterion (Spiegelhalter and colleagues, 2002) is DIC = D̄ + p_D, where D̄ is the posterior mean of the deviance −2 log f(y | θ) and the effective number of parameters is p_D = D̄ − D(θ̄), the posterior mean deviance minus the deviance at the posterior mean. Equivalently DIC = D(θ̄) + 2p_D. Smaller DIC indicates better expected predictive performance. DIC is easy to compute from MCMC output but can be unreliable for strongly non-normal posteriors, mixtures and hierarchical models where p_D can even be negative.

The Watanabe-Akaike (widely applicable) information criterion uses the log pointwise predictive density, lppd = Σᵢ log[(1/S)Σₛ f(yᵢ | θˢ)], penalised by an effective number of parameters p_WAIC (the sum over observations of the posterior variance of log f(yᵢ | θ)); on the deviance scale WAIC = −2(lppd − p_WAIC). It is fully Bayesian, uses the whole posterior and is asymptotically equivalent to leave-one-out cross-validation; Pareto-smoothed importance-sampling LOO (PSIS-LOO) estimates LOO directly from MCMC output and flags observations for which the estimate is unreliable.

These criteria answer 'which model predicts new data better?', whereas Bayes factors answer 'which model is more probable given the priors?'. The two can disagree, and neither replaces checking the model against the data. In agricultural applications, cross-validation that respects the structure of prediction (leaving out whole trials, seasons or locations, not single plots) gives the honest estimate of how a model will predict at a new site.

Remember

11. Empirical Bayes

Estimating the prior from the data

Empirical Bayes (EB) uses the data to estimate the prior, or its hyperparameters, and then applies Bayes' theorem with the estimated prior. In parametric EB, the marginal distribution of the data, m(x | η) = ∫ f(x | θ)π(θ | η) dθ, is used to estimate the hyperparameters η by maximum likelihood or moments; the posterior for each unit is then computed with η̂ plugged in. In the normal hierarchical model, EB estimates μ and τ² from the spread of the group means and gives shrunken estimates (1 − B̂ᵢ)yᵢ + B̂ᵢμ̂: exactly the hierarchical estimates of Chapter 7 with the hyperparameters fixed at their estimates.

Robbins' nonparametric EB for Poisson counts shows the idea without any assumed prior: the Bayes estimate E(θ | x) = (x + 1)m(x + 1)/m(x), and replacing the marginal probabilities by the observed frequencies gives an estimate computed from the data alone. If, among many fields, 40 had one pest colony and 20 had two, the EB estimate of the rate for a field showing one colony is 2 × 20/40 = 1.0.

The James-Stein estimator is the classic link between EB and frequentist theory. For p ≥ 3 independent normal means with unit variance, the estimator [1 − (p − 2)/‖x‖²]x has smaller total mean squared error than x itself for every value of the means: the usual estimator is inadmissible. Efron and Morris showed that James-Stein is an empirical Bayes estimator under a normal prior centred at zero; shrinking towards the grand mean instead uses p − 3 in place of p − 2.

Uses and limitations

Empirical Bayes is the working method behind much of applied agricultural statistics: small-area estimation of crop yields (Fay-Herriot models combine direct survey estimates for each district with a regression on auxiliary data such as satellite indices, shrinking noisy district estimates towards the regression prediction); genetic evaluation (BLUP with estimated variance components is EB); variety-trial analysis in which variety means are shrunk according to trial precision; and false discovery rate methods for marker screens, which can be read as EB estimates of the proportion of null hypotheses.

Its main limitation is that plugging in η̂ treats estimated hyperparameters as known, so posterior variances and intervals are too narrow, especially when the number of groups is small. Corrections (Morris's adjustments, parametric bootstrap, or a full hierarchical Bayesian analysis with a prior on η) address this. With many groups the difference between EB and full Bayes is small; with few groups full Bayes is preferable. EB also uses the data twice, once to choose the prior and once to update it, which is why it is best viewed as an approximation to a hierarchical Bayesian analysis rather than as a separate philosophy.

Remember

What the exams ask

How it is used, across a career

As a student

A student meets Bayesian inference as a PhD course and in GATE ST and ISS decision-theory items. The concrete skills are: deriving posteriors for beta-binomial, gamma-Poisson and normal models; computing posterior means, modes and credible intervals; writing Bayes estimators under different losses; computing a Bayes factor for simple hypotheses; explaining shrinkage in a hierarchical model; writing a Gibbs sampler for a simple model; and reading R̂, ESS and DIC or WAIC from software output. A thesis using hierarchical models or genomic prediction needs all of these.

Five years in

A statistician five years into work at a breeding programme or crop-estimation agency uses Bayesian tools daily: fitting hierarchical models to multi-trial variety data, shrinking district yield estimates with Fay-Herriot models, choosing and justifying priors for small trials, running and diagnosing MCMC or HMC fits, checking models by posterior predictive simulation, comparing models by LOO or WAIC, and reporting posterior probabilities of superiority that decision makers can use.

Twenty-five years in

A senior statistician with twenty-five years of experience decides where Bayesian methods enter institutional practice: official small-area crop statistics, genetic-evaluation systems, variety-release rules based on probabilities of superiority, and decision analyses for pest management and input recommendations. Such a statistician sets standards for prior justification, sensitivity analysis and computational reporting, and explains to policy users what posterior probabilities mean and what they depend on.

PhD, twenty-five years in

A PhD with long experience develops Bayesian methodology for agricultural problems: hierarchical spatial and spatio-temporal models for field trials and surveys, Bayesian genomic prediction with sparse or structured priors, efficient samplers for large mixed models, prior elicitation from accumulated trial databases, and theoretical properties such as posterior consistency and frequentist coverage of credible intervals. Such a researcher supervises doctoral work and referees the growing Bayesian literature in agricultural science.

Written from

Proofs are sketched rather than given. Not covered in depth: Bayesian nonparametrics (Dirichlet processes), variational inference and integrated nested Laplace approximation (INLA), sequential Monte Carlo, Bayesian model averaging, reversible-jump MCMC, spatial Bayesian models in detail and Bayesian genomic prediction models (BayesA, BayesB and relatives), which belong to specialised subjects. Interpretation scales for Bayes factors and thresholds for R̂ and effective sample size are conventions that differ between authors and have changed over time; the manual states which convention it quotes. No software is prescribed.

Original practice questions written to the syllabus and pattern of the exams named in each question's examTags. They are not past papers and do not reproduce any board's questions. Every numerical answer was computed with a script; the values needed are stated in the question.

Practise it, and prove it

These notes are free to read here, with no account.

What the Laisha app adds:

Free for ten days, then a subscription. We would rather write that here than have you find it later.

Open Laisha →