Past exam of the mathematics course of the University of Cambridge / 2010 / iii / Paper 36 (Applied Bayesian Statistics)

If ${\pi }_0=\mathbb{P}(H_0)$ and ${\pi }_1=\mathbb{P}(H_1)=1-{\pi }_0$, the prior odds on the null hypothesis areThese odds concern which model is true before seeing the observations. They are separate from the prior distribution of the success probability within the alternative model.

Use the orientation that favours the null hypothesis:For a composite model the denominator is its Bayesian model evidence, obtained by integrating the likelihood against the model's prior distribution. This Bayes factor exceeds one when the observed data have greater marginal probability under $H_0$. By Bayes' theorem,Thus the posterior odds equal the prior odds multiplied by the Bayes factor. The reciprocal orientation instead gives $B_{10}=1/B_{01}$.

Conditional on the success probability, the count has binomial distribution. Integrating its likelihood under the uniform prior distribution, or equivalently a beta distribution with parameters $(1,1)$, givesThe prior predictive distribution is therefore uniform over the possible counts. Excluding the single point $\theta =1/2$ does not alter this integral because a continuous prior distribution gives that point zero probability. This calculation is the marginalization underlying Beta-binomial conjugacy.

Under the null hypothesis, $Y\sim \mathrm{Bin}(n,1/2)$, with expected value $n/2$ and variance $n/4$. Approximating its probability at an integer by the normal distribution density over a unit-width cell givesDividing by the prior predictive distribution from the preceding part gives the uniform-alternative binomial Bayes factorThe second approximation replaces $n+1$ by $n$. Both are large-sample approximations for counts in the central range, rather than exact identities at extreme counts.

The normal distribution standardized statistic for the binomial distribution count isBecause the alternative is two-sided, the approximate p-value is $2\{1-\Phi (3)\}\approx 0.0027$. A continuity correction replaces $3$ by $(149.5)/50=2.99$ and gives approximately $0.0028$, with the same conclusion. Thus the result is statistically significant at both the 5% and 1% levels. The p-value is a null tail probability, not a posterior probability that the null hypothesis is true.

Substituting the standardized deviation into the Bayes factor approximation givesThis is only a slight preference for the alternative, despite the small p-value. The two procedures measure different things: the p-value sums a tail under the null hypothesis, while the Bayes factor compares probabilities of the observed count under both models. The uniform alternative prior distribution spreads its probability over many success probabilities that fit the observations poorly. Only a narrow range around $0.515$ has large likelihood, and integrating over the whole alternative penalizes this diffuse prediction.

The prefactor makes the distinction especially clear. At a fixed standardized deviation $z$, the p-value stays approximately $2\{1-\Phi (|z|)\}$, whereas $B_{01}\approx \sqrt{2n/\pi }{e}^{-z^2/2}$ eventually grows with $n$. A fixed diffuse alternative can consequently favour the point null even while a conventional test rejects it. This is the same prior-width mechanism illustrated by the Jeffreys-Lindley paradox for a Gaussian point null.

An alternative that allows small departures from chance, while assigning little probability to near-perfect prediction or near-perfect failure, is more plausible than a uniform prior distribution. For example, a symmetric beta distribution $\mathrm{Beta}(a,a)$ with $a>1$ is centered at $1/2$, withChoosing $a=50$ gives a prior standard deviation of about $0.05$ and concentrates on modest deviations. Its concentration should represent substantive beliefs specified before examining this experiment, rather than being tuned to the observed count. A continuous prior distribution centered at $1/2$ still places no point mass at $H_0$; the separate model prior probability supplies that mass.

No. A Bayes factor favouring the alternative means $B_{01}<1$, but the posterior odds are $B_{01}{O}_{01}$. They are less than one only if $B_{01}<1/O_{01}$. For example, prior odds of $1000:1$ on the null hypothesis, combined with a Bayes factor of $100:1$ in favour of the alternative, still give posterior odds of $10:1$ on the null hypothesis. An unusual claim can require substantial evidence to overcome low initial plausibility.

Take the exposure $E_i$ as known and positive. The Poisson distribution log-likelihood, up to a parameter-independent constant, is ${\ell }_i({\lambda }_i)=y_i\log {\lambda }_i-E_i{\lambda }_i$. Its first two derivatives areSince $\mathbb{E}(Y_i)=E_i{\lambda }_i$, the Fisher information isThe Jeffreys prior is proportional to its square root. Removing the constant $\sqrt{E_i}$ givesThis is an improper prior: its integral diverges at infinity.

Multiplying the Poisson distribution likelihood by the Jeffreys prior gives the posterior density kernelThus, in shape-rate convention,It is a proper posterior distribution even when $y_i=0$, provided $E_i>0$. The posterior mean and posterior variance are $(y_i+1/2)/E_i$ and $(y_i+1/2)/E_i^2$. This is Poisson-gamma conjugacy with unequal exposures, with the zero-rate prior interpreted as an improper kernel.

A proper gamma approximation to a Poisson Jeffreys prior uses shape $1/2$ and a small positive rate $\epsilon$:Its density is proportional to ${\lambda }_i^{-1/2}{e}^{-\epsilon {\lambda }_i}$, so where $\epsilon {\lambda }_i$ is small its kernel approximates the Jeffreys prior. After observing the data, Poisson-gamma conjugacy gives $\mathrm{Gamma}(y_i+1/2,E_i+\epsilon )$, which converges to the preceding posterior distribution as $\epsilon \downarrow 0$. This is an approximation of an improper kernel on the relevant parameter range, not a weak limit to a normalized Jeffreys probability distribution. The shape must be $1/2$; a small-shape gamma prior would approximate a different power of ${\lambda }_i$.

Under independent area counts and independent area prior distributions, the joint posterior distribution factors into the gamma posteriors. For each simulation $s=1,\dots ,S$, draw one joint vector by independently generatingRank the risks within that joint draw, taking rank one to mean highest risk:The empirical posterior rank distribution and highest-risk probability areIndependent gamma draws make this direct Monte Carlo integration, so Markov chain Monte Carlo is unnecessary here. Continuous posteriors give ties probability zero, and the highest-risk probabilities sum to one. The same ranking calculation can be applied to dependent joint draws from a richer hierarchical Bayesian model; ranking separate posterior means would discard the uncertainty that the program is intended to measure.

For a shape-rate gamma distribution $\mathrm{Gamma}(a,b)$, the mean and variance conditions requireThe first equation gives $a=b$, and the second then gives $a=b=16$. HenceThis proper prior distribution is centered on the national reference risk, while allowing genuine between-area variation.

The new prior density contributes the factor ${\lambda }_i^{15}{e}^{-16{\lambda }_i}$. Multiplication by the Poisson distribution likelihood givesIts posterior mean displays the resulting shrinkage:The observed risk ratio $y_i/E_i$ is pulled towards one, most strongly for small exposures. The posterior variance is $(y_i+16)/(E_i+16{)}^2$.

The informative gamma distribution prior distribution discourages extreme risks caused by small denominators. Its influence is strongest where the exposure is small, while large-exposure areas are predominantly determined by their likelihood. Consequently an area previously appearing near the top because of a noisy high raw risk ratio may have its highest-risk probability reduced, and the posterior rank distributions of poorly observed areas tend to overlap more rather than suggesting a precise league table.

There is no universal claim that every rank distribution becomes wider or narrower: unequal exposures produce unequal shrinkage and can change the ordering, while the new prior also changes the uncertainty in the risks. In the limiting case of negligible information for every area, all area posteriors approach the same independent continuous prior distribution. Symmetry then gives each area probability $1/I$ of occupying each rank. This illustrates why modest differences between weakly observed areas should not yield confident rankings.

The posterior probability $q_i=\mathbb{P}({\lambda }_i>1\mid y_i)$ answers a natural directional question: how likely is this area's risk to exceed the reference risk? It accounts for uncertainty rather than declaring an increase from the point estimate alone. However, it does not distinguish a negligible increase from a clinically important one, and it combines prior beliefs with evidence from the data. With an informative prior distribution, a large probability may already have been present before observation, or a genuinely surprising count may be strongly moderated by the prior.

To measure the evidence supplied by the observations, compare posterior and prior odds. If $q_0=\mathbb{P}({\lambda }_i>1)$ under the chosen proper prior distribution, the Bayes factor for increased versus nonincreased risk isThis follows from the posterior odds identity when the two models use the original prior conditioned respectively on ${\lambda }_i>1$ and ${\lambda }_i\le 1$. For the specified prior, $q_0$ is the upper tail at one of $\mathrm{Gamma}(16,16)$; centering its mean at one does not make this tail probability exactly $1/2$. To address seriousness itself, choose a substantive excess $c>0$ and report the posterior probability of ${\lambda }_i>1+c$, together with an interval for the size of the excess. An evidence ratio and a meaningful risk threshold answer different, complementary questions.

The logarithm is strictly increasing, so maximizing the likelihood function is equivalent to maximizing its logarithm. Multiplication by $-2$ reverses that ordering. Thus, whenever the maximum is attained,Here the Bayesian deviance includes the likelihood constants. Parameter-independent constants do not alter the minimizer, but a common convention must be retained when comparing criteria between models.

Under the usual regular, identifiable, fixed-dimensional setting, the maximum-likelihood estimator is an interior point with vanishing log-likelihood gradient, and the observed information is positive definite. DefineThe second-order Taylor expansion of the Bayesian deviance isWith a prior density approximately constant over the likelihood's concentration region, the posterior density is therefore proportional toRecognizing the multivariate normal density givesTypically $J$ grows at rate $n$ and the posterior concentrates in a region of radius $n^{-1/2}$, making higher-order terms negligible under the regularity assumptions. The local flatness assumption concerns this region, not necessarily the entire parameter space. Boundary maxima, unidentified directions, or persistent multiple modes can invalidate this approximation.

Conditional on the observed data, the preceding multivariate normal distribution approximation implies thatThe deviance expansion consequently givesThe sum of squares of $p$ independent standard normal distribution variables has chi-squared distribution with $p$ degrees of freedom. Hence, as a statement about the posterior distribution of the Bayesian deviance,The data are fixed here and the parameter varies according to its posterior; this is not a claim about resampling the data.

Let $\overline{\theta }=\mathbb{E}(\theta \mid y)$ and $\overline{D}=\mathbb{E}\{D(\theta )\mid y\}$. The usual effective parameter count in DIC isWith the locally flat prior distribution, $\overline{\theta }\approx \hat{\theta }$ and the mean of the chi-squared distribution excess is $p$. Thus $\overline{D}\approx D(\hat{\theta })+p$, while $D(\overline{\theta })\approx D(\hat{\theta })$, giving $p_D\approx p$.

More generally, the quadratic posterior deviance moments give the effective count. Within the same approximation, let the posterior covariance matrix be $\Sigma$. The expectation of a quadratic form giveswhereas evaluation at $\overline{\theta }$ omits the trace term. Therefore $p_D\approx \mathrm{tr}(J\Sigma )$. For a flat prior, $\Sigma \approx J^{-1}$ and the trace is $p$; informative prior shrinkage can reduce this effective count.

Using the same Bayesian deviance convention, the deviance information criterion and Akaike information criterion areBoth combine a lack-of-fit term with a complexity penalty; smaller values are preferred among comparable models. Under the locally flat prior distribution and the regular multivariate normal distribution approximation, $D(\overline{\theta })\approx D(\hat{\theta })$ and $p_D\approx p$, soAll fitted parameters counted in the likelihood, including unknown dispersion parameters, belong in $p$. The equivalence need not hold for informative priors, nonregular models, or different choices of the likelihood's latent-variable representation.

For retained Markov chain Monte Carlo draws ${\theta }^{(1)},\dots ,{\theta }^{(S)}$, calculate both the parameter average $\hat{m}=S^{-1}{\sum }_s{\theta }^{(s)}$ and the average deviance ${\hat{D}}_{\mathrm{av}}=S^{-1}{\sum }_{s}D({\theta }^{(s)})$. Two estimates areThe first uses the approximate equality of the posterior mean and maximum-likelihood estimator. The second subtracts the mean $p$ of the chi-squared distribution deviance excess. These are different operations: evaluation after averaging versus averaging after evaluation. Both rely on a well-mixed chain and the regular locally flat-prior approximation. The smallest sampled deviance is an upper bound on the true minimum, but merely taking that minimum can be inefficient in high dimension because few draws approach the mode closely.

The alternative effective count isUnder the locally flat-prior approximation, the constant $D(\hat{\theta })$ contributes no variance, and the remaining chi-squared distribution term has variance $2p$. Hence $p_V\approx p$, the same calibration as the usual effective parameter count in DIC. It can be estimated by half the sample variance of deviance draws, without evaluating deviance at a parameter average.

This is a calibrated alternative, not a general identity with $p_D$. For example, with a quadratic deviance and a multivariate normal distribution posterior of mean $m$ and covariance $\Sigma$, putting $b=m-\hat{\theta }$ givesThey agree when $b=0$ and $\Sigma =J^{-1}$, but informative priors can change both the covariance and the displacement of the posterior mean. Their equality should therefore not be assumed outside the stated regime.

The common normal distribution for the child-specific random intercepts produces partial pooling. Children with little information borrow strength from the common population, while well-observed children retain more of their individual estimate. For fixed slopes and hyperparameters, let ${\overline{r}}_i$ average the three responses after subtracting their slope terms. Multiplying the normal likelihood by the normal intercept prior distribution giveswith conditional variance $(3/{\sigma }^2+1/{\tau }^2{)}^{-1}$. Thus small between-child variation produces stronger shrinkage towards $\delta$. The common mean and variation are themselves learned from all children in this hierarchical Bayesian model.

Center the log-time predictor and the baseline log measurement near their sample means, so the intercept describes a typical predictor configuration. This reduces intercept-slope posterior correlation and can improve Markov chain Monte Carlo mixing. If exactly the original model is required, transform the intercept prior consistently when changing coordinates rather than silently assigning a different independent prior.

For weakly informed child effects, a non-centered Gaussian random-effect parameterization can also help:This retains the same conditional prior distribution for ${\alpha }_i$ but removes its prior-scale dependence from the sampled $z_i$. It can reduce strong dependence between the effects and $\tau$. A centered parameterization can be better when each child's effect is very accurately observed, so improvement is not automatic. Blocking strongly correlated coefficients is another option. Check the resulting chains with Markov chain Monte Carlo convergence diagnostics and effective sample size of a Markov chain, rather than judging convergence from the number of iterations alone.

The bounded uniform prior distributions on the slopes and population intercept are proper and broad on their stated scale. They allow either sign but are not invariant to changes of units. The residual precision parameter $q=1/{\sigma }^2$ has a proper shape-rate $\mathrm{Gamma}(0.001,0.001)$ prior. Equivalently, ${\sigma }^2$ has inverse-gamma distribution with the same shape and scale. This is a commonly used diffuse choice, but it is not literally noninformative and its tail behaviour merits sensitivity analysis.

The prior for $\tau$ is uniform on $(0,100)$, so the induced density of $v={\tau }^2$ isIt is integrable at zero, unlike the usual scale-invariant variance kernel $\pi (v)\propto 1/v$. The latter kernel, often called the standard variance Jeffreys prior, would cause an improper posterior from a log-uniform random-effect scale prior here. After integrating the child intercepts, each child's response vector has a normal likelihood with covariance matrix ${\sigma }^2{I}_3+v11^T$. For fixed positive ${\sigma }^2$, this covariance remains nonsingular at $v=0$, and the marginal likelihood has a positive finite limit there. On a compact set of the other parameters its positivity givesConsequently the scale-invariant kernel cannot be normalized to a posterior distribution. A proper prior on the standard deviation avoids this problem. This familiar scale prior is not necessarily the actual Jeffreys prior for an additive variance component derived from the observed-data likelihood; the latter can be finite at zero.

The baseline assay is a noisy measurement just like the follow-up assays. On the log scale, a normal distribution for its error is a reasonable approximation if variability is approximately multiplicative on the original titre scale. The latent ${\mu }_{0i}$ represents the true baseline log-titre, withUsing a common variance is a simplifying assumption justified if baseline and follow-up measurements have comparable assay variability; it should not be treated as a fact guaranteed by the study design. Treating $y_{0i}$ as exact ignores measurement error in a predictor and can lead to attenuation bias from classical measurement error and understated uncertainty.

Model 2 is an errors-in-variables model: the follow-up regression depends on a latent true baseline rather than substituting its noisy observation. A common normal distribution for the true baselines provides partial pooling and allows uncertainty about each baseline to propagate to the regression coefficients and predictions. Conditional on the shared parameters, the child-specific factorization isThe Directed acyclic graph below includes the shared parameters, the two local latent variables, both observed measurements, and the deterministic regression mean. The outer plate repeats over children and the inner plate over follow-up visits. Filled nodes are observed; rectangular mean nodes are deterministic.

One complete adapted BUGS specification uses the non-centered form of the random intercept, retaining the original priors and adding proper broad priors for $\theta$ and $\psi$:

for (i in 1:N) {
  z[i] ~ dnorm(0, 1)
  alpha[i] <- delta + tau*z[i]
  truebase[i] ~ dnorm(theta, baseprec)
  baseline[i] ~ dnorm(truebase[i], noiseprec)
  for (j in 1:3) {
    followmean[i,j] <- alpha[i] + beta*log(time[i,j]) + gamma*truebase[i]
    followup[i,j] ~ dnorm(followmean[i,j], noiseprec)
  }
}
noiseprec ~ dgamma(0.001, 0.001)
beta ~ dunif(-100, 100)
gamma ~ dunif(-100, 100)
delta ~ dunif(-100, 100)
tau ~ dunif(0, 100)
theta ~ dunif(-100, 100)
psi ~ dunif(0, 100)
baseprec <- 1/(psi*psi)

Here $N$ is the number of children, baseline and followup contain the observed data, and dnorm uses precision rather than variance. The priors on the newly introduced population baseline mean and standard deviation are explicit choices, not consequences of the likelihood; sensible alternatives can be checked for sensitivity. The code uses the noisy baseline only in its own measurement likelihood, and uses truebase in every follow-up mean.

The PDF output places $-1$ inside the 95% credible interval for the time slope, whose endpoints are approximately $-1.329$ and $-0.7955$, and places $1$ inside the baseline slope interval, approximately $0.7915$ to $1.231$. Relative to their posterior standard deviations, the proposed fixed values are onlyfrom the posterior means. Thus the simpler choice $\beta =-1,\gamma =1$ is compatible with the fitted coefficient summaries and has a natural proportional-decay interpretation. It is a reasonable constrained model to compare, not proof of the exact equalities: separate marginal intervals do not establish a joint constraint, and a continuous posterior distribution assigns a specified point probability zero.

The reported Monte Carlo errors measure numerical uncertainty in the estimated posterior means, not uncertainty about the coefficients themselves. In particular, the baseline slope's Monte Carlo error is about 9% of its posterior standard deviation, appreciably larger than the corresponding proportion for the time slope. Longer or better-mixed chains may be needed for precise summaries.

The reported deviance information criterion differs byModel 3 has a mean Bayesian deviance worse by only $0.2$, but its effective parameter count in DIC is lower by $2.1$. The small fit cost is more than offset by the reduced complexity penalty, so the criterion gives a mild preference to Model 3. A difference this small should not be presented as decisive evidence, especially without the Monte Carlo error of the criterion difference.

The counts $143.6$ and $141.5$ are effective rather than literal parameter numbers. The likelihood includes child-specific intercepts and uncertain true baselines as well as global coefficients. With 106 of each local quantity and seven global parameters, Model 2 has 219 stochastic unknowns in a natural full representation, but partial pooling and other prior information reduce their effective contribution. Fixing two slopes lowers the effective count by about two, with a further small change from the fitted posterior. There is no reason for $p_D$ to be an integer, or to equal just the number of population-level regression coefficients.

For illustration, the deviances evaluated at posterior means implied by the table are $\overline{D}-p_D=984.5$ for Model 2 and $986.8$ for Model 3. These differ from the average deviances because averaging a nonlinear deviance is not the same as evaluating it at an averaged parameter. This interpretation presumes the same observed-data likelihood convention in both fits; integrating local variables out before forming the deviance would define a different comparison. The DIC difference is not a Bayes factor.

Let $T_{0i}=e^{{\mu }_{0i}}$ denote the child's latent true baseline titre. With the constrained slopes, the conditional mean of log follow-up titre isExponentiating this log-scale center givesThus the conditional median, or geometric-mean titre, expressed as a fraction of true baseline titre, decays inversely with elapsed time; the child-specific intercept supplies the proportionality constant. It is not an exponential decay in time, and the formula is a follow-up model for $t>0$, not an extrapolation to the immunisation instant $t=0$.

For the arithmetic conditional mean, the lognormal distribution correction givesThe correction changes the constant but not the inverse-time dependence when the residual variance is constant. Averaging over random intercepts similarly preserves the time factor whenever the required expectation is finite.