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

The binomial likelihood is $p(x\mid \theta )=\left(\genfrac{}{}{0ex}{}{m}{x}\right){\theta }^x(1-\theta {)}^{m-x}$. Multiplication by the beta distribution prior density and Bayes' theorem give the posterior density kernel ${\theta }^{a+x-1}(1-\theta {)}^{b+m-x-1}$. Normalizing by the beta function yields Beta-binomial conjugacy:This is proper for $a,b>0$ and $0\le x\le m$, including counts at either boundary.

For the $m$ independent Bernoulli trials, the log likelihood function, apart from its constant, is $\ell (\theta )=x\log \theta +(m-x)\log (1-\theta )$. Its negative expected second derivative, the Fisher information, isFor $m>0$, the Jeffreys prior is proportional to its square root. Since $\mathrm{B}(1/2,1/2)=\pi$, its normalized form isThe Jeffreys prior is invariant under smooth one-to-one reparameterization as a measure. If $\varphi =h(\theta )$, then $I_{\varphi }(\varphi )=I_{\theta }(\theta )(d\theta /d\varphi {)}^2$, so $\sqrt{I_{\varphi }}d\varphi =\sqrt{I_{\theta }}|d\theta |$. Thus recomputing it in the new coordinate gives precisely the transformed prior, rather than a new prior. In several dimensions the corresponding density is $\sqrt{\det I(\theta )}$, with the same Jacobian transformation property.

Put $\alpha =a+x$, $\beta =b+m-x$, and $S=\alpha +\beta =m+a+b$. Integrate the future binomial distribution against the current beta distribution posterior density:This is the beta-binomial distribution. For positive integer $a,b$, the gamma function values reduce to factorials, givingUsing $(S-1)!=(S-1)(S-2)!$ and $(S+n-1)!=(S+n-1)(S+n-2)!$ rewrites this asSubstitution of $\alpha ,\beta ,S$ gives exactly the two required factors. The probability is zero outside $0\le y\le n$.

With a uniform beta distribution prior, the posterior predictive distribution for another batch of the same size isBut $\mathrm{B}(x+1,m-x+1)=1/[(m+1)(\genfrac{}{}{0ex}{}{m}{x})]$. HenceThe displayed expression is symmetric in $x,y$. This is the uniform-count predictive property: the joint batch counts are exchangeable and both marginal counts are uniform. Conditional reversal is not a general property of arbitrary priors; it holds here because those marginal probabilities coincide.

Before collecting any observations, a uniform success-chance prior distribution givesThus the uniform-count predictive property assigns equal probability to every possible number of successes in a fixed-size batch under the prior predictive distribution. It is uniform over counts, not over all ordered binary sequences. After observing $x$ successes in $m$ trials, its familiar one-step posterior predictive probability is $(x+1)/(m+2)$; the count-uniform property refers to prediction before those data are observed.

For the uniform prior, write $s=x+y$ for the pooled success count in $m+n$ trials. Conditional on $s$, exchangeability makes all placements of the successes equally likely. The successes in the second batch therefore have a hypergeometric distribution, and the factor in question isThis is hypergeometric allocation of exchangeable Bernoulli counts: sample $n$ positions without replacement from a population of $m+n$ positions containing $s$ successes. The uniform-count predictive property gives $\mathrm{Pr}(X+Y=s)=1/(m+n+1)$ and $\mathrm{Pr}(X=x)=1/(m+1)$, soAs $y$ varies with $x$ fixed, $s=x+y$ also varies. Therefore $B$ alone is not a normalized hypergeometric distribution over $y$; its normalization factor is essential.

Under the discrete uniform distribution, the observation has probability $1/N$ when its value is an available label, and zero otherwise. Therefore the likelihood function isThis is the discrete uniform endpoint posterior model's likelihood; all endpoints smaller than the observed label are excluded.

On its admissible support the likelihood function $1/N$ strictly decreases. Its maximum-likelihood estimate is consequently the smallest admissible endpoint:

A constant improper prior would give the formal posterior density kernel $N^{-1}{1}_{\{N\ge y\}}$. However,There is no normalized posterior distribution under this prior. The failure is not simply that the prior is improper: some improper priors give proper posteriors. Here posterior propriety itself fails after this single observation, so probabilities, quantiles and a posterior mean cannot be assigned to that kernel.

Write $H_k={\sum }_{j=1}^{k}1/j$, with $H_0=0$. The truncated uniform prior gives the posterior distributionIts posterior mean isEach term in the numerator of the expectation is $N/N=1$, which explains the exact count $M-y+1$.

Keep $y$ fixed. Integral bounds for the decreasing function $1/t$ give $H_M=\log M+O(1)$, so $H_M-H_{y-1}=\log M+O_y(1)$. ConsequentlyThe posterior mean thus grows like $M/\log M$ instead of approaching a data-determined limit. Even the posterior distribution mass below a fixed endpoint $n$ is $(H_n-H_{y-1})/(H_M-H_{y-1})\to 0$. Hence increasing the supposedly innocuous prior cutoff moves the inference to larger and larger endpoints. In an integral approximation its median is of order $\sqrt{yM}$, also depending strongly on the cutoff. This explains the sensitivity of the entire inference, not just of its posterior mean.

Multiplying the reciprocal improper prior by the likelihood function givesThus this discrete uniform endpoint posterior is proper. Its exact cumulative probability is $S_y^{-1}{\sum }_{k=y}^n{k}^{-2}$. Integral approximation gives $S_y\approx {\int }_y^{\infty }{t}^{-2}dt=1/y$ and ${\sum }_{k=y}^n{k}^{-2}\approx {\int }_y^n{t}^{-2}dt=1/y-1/n$, soThese approximations are useful when $y$ is moderately large and are not exact discrete identities. For example the exact untruncated median at $y=100$ is $199$. Although its quantiles are finite, its posterior mean is infinite, since its expectation numerator is the divergent sum ${\sum }_{N=y}^{\infty }1/N$.

The WinBUGS model uses a finite truncation of the reciprocal prior. The numbered lines have the following roles:

Combining that categorical distribution observation likelihood with the truncated reciprocal prior producesThe Markov chain Monte Carlo output approximates this posterior. The actual PDF separates the observation line from the preceding comment; it must remain an active stochastic statement.

The posterior distribution has a long right tail. Its posterior mean is pulled upward by relatively rare large endpoints and is sensitive to the artificial upper cutoff; its untruncated posterior mean does not exist. The median is a central quantile that stays finite as the cutoff is removed and is much less sensitive to extreme draws. It is also the Bayes estimator under absolute-error loss, whereas a mean corresponds to squared-error loss when the relevant expectation exists.

The median better represents the typical posterior endpoint in this strongly right-skewed distribution. It need not be preferable for every decision problem, since the loss function matters. For comparison, exact summation of the truncated kernel at cutoff $5000$ gives median $195$ and mean about $397.67$; the displayed simulation summaries are Monte Carlo estimates, not exact posterior calculations.

At several doses the three replicate plates differ much more than the square-root-of-mean scale expected from a Poisson distribution. For example the group at dose $100$ has sample mean about $42.67$ and sample variance about $274.33$, rather than variance near $42.67$. Its counts span a range of $33$, compared with a Poisson standard deviation about $6.53$. The dose-$33$ and dose-$1000$ groups likewise have variances above their means. Thus extra plate-to-plate variation is plausible. With only three replicates per dose this is evidence suggesting overdispersion, rather than a precise determination of it; not every group has excess sample variance.

Let ${\eta }_i=\alpha +\beta \log (x_i+10)+\gamma x_i$. The Poisson-lognormal random-effect model assigns each plate a random intensity ${\mu }_{ij}=e^{{\eta }_i+{\lambda }_{ij}}$. Conditional on that intensity its Poisson distribution has mean and variance both ${\mu }_{ij}$, but marginalizing over the normal distribution of the random effect givesThe law of total variance yieldsThis exceeds the mean whenever $\tau >0$. The independent plate effects supply additional heterogeneity without inducing dependence between plates conditional on the global parameters. Also $e^{{\eta }_i}$ is the median random intensity, not its marginal mean; the factor $e^{{\tau }^2/2}$ matters for population predictions.

Centre and scale both predictors, particularly the raw dose whose scale is much larger than that of $\log (x_i+10)$. For example use $z_{1i}=[\log (x_i+10)-c_1]/s_1$, $z_{2i}=(x_i-c_2)/s_2$ and the linear predictor ${\eta }_i=a_0+b_1{z}_{1i}+b_2{z}_{2i}$. This reduces the intercept-slope correlations and numerical scale disparities; further orthogonalizing correlated columns or updating coefficients jointly can improve Markov chain Monte Carlo mixing. The exact old prior can be retained through the corresponding coefficient transformation; assigning new independent priors to transformed coefficients would instead change the prior.

For weakly informed plate effects or small $\tau$, use the non-centred Gaussian random-effect parameterisation ${\lambda }_{ij}=\tau z_{ij}$ with $z_{ij}\sim N(0,1)$. This often reduces dependence between the scale and local effects; it is not guaranteed to outperform the centred version when plate effects are strongly informed. Sensibly calibrated proper priors also avoid the enormous initial log intensities permitted by the stated coefficient ranges. Use dispersed chains and Markov chain Monte Carlo convergence diagnostics to assess stationarity, agreement between chains and effective sample size of a Markov chain. Reparameterization can improve mixing; simply retaining more iterations does not remove poor mixing.

The regression coefficients have independent proper uniform distributions on very wide finite intervals. These are flat in the coefficients, not flat in fitted rates, and can be highly informative on the rate scale. The scale $\tau$ is uniform on $(0,100)$, and WinBUGS receives the corresponding precision $1/{\tau }^2$ for each normal distribution random effect. Writing $v={\tau }^2$, the induced prior density isIt is proper and integrable at zero.

The usual scale-invariant variance prior $p(v)\propto 1/v$ would instead cause an improper posterior from a log-uniform random-effect scale prior. After integrating out the plate effects, the likelihood function tends as $v\downarrow 0$ to the ordinary positive Poisson likelihood. On any compact interior set of coefficient values this limit is positive and bounded away from zero. Hence the posterior normalizing integral contains a divergent factor ${\int }_0^{\epsilon }dv/v$. There is no posterior distribution under that prior, even if formal full conditional distributions appear usable. Thus a proper prior with integrable mass near zero is needed here. The reciprocal-scale rule is the standard pure-scale Jeffreys prior; it should not be confused with the model-specific information prior for an additive variance component.

Delete the plate random effects from the linear predictor and remove their stochastic nodes and the unused scale/precision nodes. The observation layer then remains a Poisson distribution withFor example the essential WinBUGS observation block can be written asRetain the coefficient priors. This is the $\tau =0$ sampling model implemented directly; attempting to keep 1/(tau*tau) while fixing tau to zero would divide by zero.

Under the common observation-level Bayesian deviance convention, DIC is $\overline{D}+p_D$, with the effective parameter count in DIC defined by $p_D=\overline{D}-D(\overline{\psi })$. The no-random-effect model's effective count is near three, corresponding to its three fitted regression coefficients. With plate random effects, the effective count rises to $13.6$: the local effects add flexibility, but their partial pooling means they do not count as eighteen fully independent unrestricted fitted parameters. The effective count need not be an integer or equal the raw parameter count; the scale parameter controls that shrinkage.

Adding the plate effects reduces the average Bayesian deviance by $28.6$ at an effective-complexity cost of $10.7$. Therefore its DIC is lower byThe reported criterion favours the overdispersed model. This is a fit-complexity comparison, not a posterior probability or Bayes factor. In a hierarchical Bayesian model, the interpretation depends on whether deviance is conditional on local effects or has marginalized them out; these figures use the supplied observation-level convention.

For each plate use its conditional fitted intensity, including the fitted random effect. At a posterior draw $s$ form the Pearson residualA simpler plot uses ${\hat{\mu }}_{ij}$ in this expression, but posterior draws also display uncertainty in the fitted values. Look for residuals centred around zero, unexpectedly large tails or isolated outliers, and patterns with dose, fitted intensity or plate grouping. For known conditional parameters the Poisson distribution gives mean zero and variance one for this residual; approximate normality is reasonable only at adequate count sizes. Fitting the same data and estimating latent effects changes its calibration, so compare the diagnostics with residuals from conditional posterior predictive checks, rather than treating them as independent exact standard normals.

This checks the plate-level Poisson sampling layer given the effects and can reveal remaining unmodelled variation. It does not by itself establish the population dose-response curve: fitted plate effects can absorb departures from that curve.

Use conditional versus marginal posterior predictive checks to replicate new plates, rather than recycling the fitted plate effects. For each joint posterior draw of $(\alpha ,\beta ,\gamma ,\tau )$, draw fresh independent ${\lambda }_{ij}^{\mathrm{rep}}\sim N(0,{\tau }^2)$ and thenUse the same doses and three-plate design, giving the posterior predictive distribution for repeated experiments. Compare observed dose means and their trend, the high-dose downturn, within-dose spreads, and outliers with the replicated curves and intervals. For example compare a discrepancywith $T(y^{\mathrm{rep}};\psi )$ draw by draw. The resulting posterior predictive tail fraction or graphical rank assesses whether the observed trend is unusual under this model. Fresh plate effects are needed to assess the population dose-response assumption. Reusing the fitted effects would mainly check the conditional sampling layer and could hide systematic departures from the proposed curve.

Compute the Bayesian model evidence for each model by integrating its sampling density against its parameter prior:The posterior odds are the prior odds multiplied by the Bayes factor:The parameter priors used for evidence must be normalized; arbitrary normalization constants in model-specific improper priors would make this comparison undefined.

Within the positive component, normal-normal conjugacy givesCompleting the square yieldsHere positive denotes membership in the variable-effect component, not the condition ${\theta }_i>0$; its normal prior permits effects of either sign.

Integrating the positive component's normal distribution prior gives its prior predictive distribution:The negative component gives $Y_i\mid -\sim N(0,1)$. Let $f_1(y)={\varphi }_{1+V}(y)$ and $f_0(y)={\varphi }_1(y)$ denote those centred normal densities. The point-null mixture prior assigns component odds $q/(1-q)$, so the posterior odds of the positive component areIts posterior probability is $O_i/(1+O_i)$. Evidence depends on $y_i^2$, so large effects of either sign favour the positive component. This is the marginal distribution before observing $y_i$, not the posterior predictive law for a second observation conditional on it.

Under the independent-gene version of the model, integrating both component labels and effects givesThis product uses independence of the gene effects and observation errors conditional on the shared parameters. Maximize its logarithm over $0\le q\le 1$, for example with a one-dimensional bounded optimizer or a grid followed by refinement. Its derivatives areThus the log likelihood function is concave. An interior maximum-likelihood estimate solves the score equation; otherwise it is at zero or one, determined by the endpoint score signs. An alternative expectation-maximization algorithm iteratesFor $V=0$, the components coincide and $q$ is not identifiable; for $V>0$ the full observation values carry information about the mixing fraction.

Use prior elicitation to make the qualitative judgments quantitative: ask whether $10\%$ represents a mean, median or mode, and what probability “unlikely above $15\%$” should mean, for example a $5\%$ upper-tail probability. A beta distribution on $(0,1)$ is a convenient family. If $0.1$ is taken as the mean, write its parameters as $a_q=0.1\kappa$, $b_q=0.9\kappa$, then choose $\kappa$ so that the desired tail probability above $0.15$ is obtained. Check its quantiles and prior predictive distribution for plausibility, and assess sensitivity to reasonable alternative elicited values. The information suggests an informative prior near $0.1$, but does not uniquely determine its parameters.

Introduce independent component indicators $z_i\mid q\sim \mathrm{Bernoulli}(q)$, and auxiliary slab effects $u_i\sim N(0,V)$, independent of the indicators. Set ${\theta }_i=z_i{u}_i$. Then inactive genes have exactly zero effect, while active genes have the stipulated normal distribution. With $V>0$ supplied as known data and elicited positive beta shapes, one suitable WinBUGS model isThe normal arguments are precisions, not variances. Supply all observed y values, monitor q, theta and optionally z, and assess Markov chain Monte Carlo convergence diagnostics. The proper auxiliary prior remains defined for inactive genes, where $u_i$ is not informed by the likelihood. This represents the exact point-null mixture prior, rather than replacing its point mass by a narrow continuous spike. Conditional on the indicators, $q\mid z$ has beta distribution $\mathrm{Beta}(a_q+{\sum }_i{z}_i,b_q+N-{\sum }_i{z}_i)$, providing a useful check on the update.

Let $p_0=\overline{\Phi }(2)$ and $p_1=\overline{\Phi }(2/\sqrt{1+V})$, where $\overline{\Phi }$ is the standard normal distribution upper-tail probability. The mixture model givesFor $V>0$, threshold inference for a mixture proportion therefore givesThis is a method-of-moments estimate when it lies in $[0,1]$. More formally, if $K$ of the $N$ independent genes exceed the threshold, then $K\sim \mathrm{Binomial}(N,(1-q)p_0+q{p}_1)$, and its constrained maximum-likelihood estimate clips the displayed value to $[0,1]$. No numerical estimate is determined without the given value of $V$.

Here $p_0\approx 0.02275$. Matching $15\%$ exactly requires $p_1\ge 0.15$, equivalently $V\ge [2/{\Phi }^{-1}(0.85){]}^2-1\approx 2.724$. For smaller positive $V$, the maximum is at $q=1$, and a large-sample exceedance fraction this high suggests a mismatch with the prescribed component variances. Sampling fluctuation is still possible; an observed fraction need not equal its expectation. If $V=0$, the denominator vanishes and this summary gives no information about $q$. Because $p_1<1/2$ for finite $V$, matching $15\%$ in expectation also forces $q>(0.15-p_0)/(0.5-p_0)\approx 0.267$. Thus the summary suggests tension with the earlier prior concentrated near $0.1$, under this mean-zero alternative model. The event is the one-sided exceedance $Y_i>2$, not $|Y_i|>2$.