Worked Bayesian Random-Intercept Logistic GLMM with Prior Sensitivity

From a hierarchical model to a posterior fit

A Bayesian GLMM fit combines a sampling model for the observed outcomes with probability distributions for the unknown coefficients and random effects. In a random-intercept logistic model, the intercept for each group shifts that group’s log odds relative to the population-level intercept. The key inferential question in this tutorial is how much variation among groups is supported by the data, and how much the answer changes when the prior scale for that variation changes.

The examples use hypothetical clinic data. The outcome is whether a patient experiences a specified event, recorded as $Y_{ij}=1$ for an event and $Y_{ij}=0$ otherwise, where $i$ indexes patients and $j$ indexes clinics. No covariates are included, so the example isolates the random-intercept variance component.

Model and assumptions: For patient $i$ in clinic $j$, assume $Y_{ij}\mid b_j,\alpha$ are conditionally independent Bernoulli outcomes with probability $p_{ij}$, where $$ \operatorname{logit}(p_{ij})=\log\left(\frac{p_{ij}}{1-p_{ij}}\right)=\alpha+b_j. $$ The clinic effects are independent and identically distributed as $b_j\mid\sigma\sim N(0,\sigma^2)$, with $\sigma\geq 0$. The model assumes a common population intercept, normally distributed clinic effects, and conditional independence of outcomes given the effects. It also assumes that the recorded groups and outcomes are correctly specified and that the data contain no unmodeled covariates needed for the intended comparison.

A simple prior specification is $\alpha\sim N(0,2.5^2)$ and $\sigma\sim\operatorname{HalfNormal}(A)$, independently, where $A$ is the prior scale and a half-normal distribution is the distribution of the absolute value of a $N(0,A^2)$ random variable. These priors are proper: they integrate to one over their respective parameter spaces. The scale $A$ is stated in log-odds units because $\sigma$ is the standard deviation of clinic effects on the log-odds scale.

For computation, one may use the centered effects $b_j$, or the noncentered representation $b_j=\sigma z_j$, with independent $z_j\sim N(0,1)$. These are two parameterizations of the same model, not different scientific assumptions. The earlier material on centered and noncentered parameterizations explains why their sampling behavior can differ.

Likelihood and posterior target

If clinic $j$ has $n_j$ patients and $k_j$ events, its conditional event probability is $p_j=\operatorname{logit}^{-1}(\alpha+b_j)$. The grouped-data likelihood is a product of binomial likelihoods. The binomial coefficient does not depend on the parameters, but retaining it makes the sampling model explicit.

L(\alpha,\mathbf b)= \prod_{j=1}^{J}{n_j\choose k_j} \left[\operatorname{logit}^{-1}(\alpha+b_j)\right]^{k_j} \left[1-\operatorname{logit}^{-1}(\alpha+b_j)\right]^{n_j-k_j}.

With the stated independent priors, the posterior density is proportional to the likelihood times the prior densities. Under the noncentered representation, the target density for $(\alpha,\sigma,\mathbf z)$ is

\pi(\alpha,\sigma,\mathbf z\mid \mathbf y) \propto L(\alpha,\sigma\mathbf z)\, \phi(\alpha;0,2.5^2)\, \pi_A(\sigma)\, \prod_{j=1}^{J}\phi(z_j;0,1), \qquad \sigma\geq 0,

where $\phi(\,\cdot\,;\mu,v)$ denotes a normal density with mean $\mu$ and variance $v$, and $\pi_A$ is the half-normal density with scale $A$. MCMC produces dependent draws whose stationary distribution is this posterior. Posterior summaries are calculated from those draws after warmup and adaptation; they are not obtained by treating the draws as independent observations.

Derivation of the noncentered target: In the centered model, the conditional density of the $J$ effects contains the factor $\sigma^{-J}\prod_j\phi(b_j/\sigma;0,1)$. Substituting $b_j=\sigma z_j$ changes the volume element by $d\mathbf b=\sigma^Jd\mathbf z$. The factors $\sigma^{-J}$ and $\sigma^J$ cancel, leaving $\prod_j\phi(z_j;0,1)$ in the noncentered target. Thus the transformation preserves the model’s probability distribution while changing the coordinates used by the sampler.

Worked Example: Writing the clinic likelihood

Suppose eight hypothetical clinics each contribute ten patients. Their event counts are $1,2,3,4,6,7,8,9$, respectively. The total is $40$ events among $80$ patients, so the observed overall event proportion is $40/80=0.500$. This overall proportion is descriptive; it does not by itself estimate the between-clinic standard deviation.

For the first clinic, $n_1=10$ and $k_1=1$, so its likelihood contribution conditional on its random intercept is $$ {10\choose 1}p_1(1-p_1)^9, \qquad p_1=\operatorname{logit}^{-1}(\alpha+b_1). $$ For the last clinic, the contribution is ${10\choose 9}p_8^9(1-p_8)$. Multiplying the eight such contributions gives the grouped likelihood. For example, if $\alpha=-0.10$ and the first clinic’s effect were $b_1=0.70$, its log odds would be $-0.10+0.70=0.60$, and $p_1=1/(1+\exp(-0.60))\approx0.6457$. The corresponding expected event count among ten patients would be $10(0.6457)=6.457$. These values illustrate the model calculation; they are not fitted estimates for this data set.

The interpretation of $b_1=0.70$ is also specific: holding the population intercept fixed, the clinic’s conditional odds are multiplied by $\exp(0.70)\approx2.014$ relative to a clinic with effect zero. It is not a 70-percentage-point increase in event probability.

Summarizing variance-component draws

The standard deviation $\sigma$ and variance $\sigma^2$ describe the same random-effect distribution on different scales. For each retained posterior draw $\sigma^{(m)}$, calculate $(\sigma^{(m)})^2$ if a variance-scale summary is needed. Because squaring is increasing for nonnegative values, quantiles transform directly: the posterior median of $\sigma^2$ is the square of the posterior median of $\sigma$, and an interval with endpoints $L$ and $U$ for $\sigma$ becomes an interval with endpoints $L^2$ and $U^2$ for $\sigma^2$. Posterior means do not generally transform this way.

Worked Example: Converting a posterior interval to the variance scale

Suppose a hypothetical MCMC fit gives a posterior median of $\sigma=1.20$ and a central 90% credible interval of $0.55$ to $2.30$, in log-odds units. The median on the variance scale is $$ (1.20)^2=1.44. $$ The corresponding 90% interval for $\sigma^2$ is $$ (0.55^2,\;2.30^2)=(0.3025,\;5.290). $$ Thus, under the specified model and prior, the posterior places 90% probability on variance values between $0.3025$ and $5.290$ log-odds-squared units. This is a posterior probability statement conditional on the model and prior; it is not a claim that the fixed interval has a 90% probability of containing the parameter.

The interval is wide, which signals substantial uncertainty about the amount of clinic heterogeneity. A single estimate such as the median should therefore be read together with the interval and with the prior-sensitivity analysis.

Posterior probabilities can be estimated by counting draws in a specified region. If $M$ retained draws are available, the estimate of $P(\sigma>c\mid\mathbf y)$ is the fraction of those draws above $c$. This gives a direct summary for a substantively chosen threshold, but the threshold must be expressed on the correct scale and justified by the application.

Worked Example: Estimating posterior probability from draws

Suppose four chains provide $M=16{,}000$ retained draws in total, and $14{,}240$ of them have $\sigma>0.50$. The empirical posterior probability is $$ \widehat{P}(\sigma>0.50\mid\mathbf y)=\frac{14{,}240}{16{,}000}=0.890. $$ For the hypothetical clinics, this says that the fitted posterior assigns approximately 89.0% probability to a between-clinic standard deviation greater than $0.50$ log-odds units. The estimate is meaningful only if the chains adequately explore the posterior; a large raw draw count cannot compensate for poor mixing or nonconvergence.

Checking chain behavior before interpreting the fit

A posterior summary is trustworthy only to the extent that the MCMC has explored the target distribution. For this example, suppose four chains were run from dispersed, plausible initial values, with warmup discarded. The hypothetical retained output has overlapping trace plots for $\sigma$, no persistent chain separation or long-term drift, $\widehat R=1.002$, and an effective sample size of 6,200 for $\sigma$. These checks are mutually supportive: the traces show movement across similar ranges, the between-chain and within-chain comparison is close to one, and the effective sample size indicates that the 16,000 correlated draws contain information comparable to about 6,200 independent draws for this summary.

If the posterior standard deviation of the retained $\sigma$ draws is about $0.75$, an approximate Monte Carlo standard error for their mean is $0.75/\sqrt{6200}\approx0.0095$. This measures simulation error in estimating the posterior mean, not the posterior uncertainty about $\sigma$; the latter is represented by the spread of the posterior itself. These diagnostics do not establish that the model fits the data. They assess computation, and substantive model checking remains a separate task.

Warning signs include chains trapped in different regions, a trace with a sustained trend, poor effective sample size for $\sigma$, or a noticeably elevated $\widehat R$. If these occur, posterior intervals and tail probabilities should not be reported as settled results. Reconsider parameterization, run length, initialization, or sampler settings, then repeat the checks. The earlier material on mixing, effective sample size, and Monte Carlo error gives fuller treatments of those diagnostics.

Prior-scale sensitivity

A half-normal prior puts more probability near zero than at large values, while its scale $A$ determines how quickly that probability decreases. Changing $A$ changes the prior information about plausible between-clinic variation. Sensitivity analysis fits the same likelihood under alternative defensible values of $A$, then compares relevant posterior summaries. It does not select a prior by seeking the most convenient conclusion.

For the hypothetical eight-clinic data, consider three half-normal priors for $\sigma$, with scales $A=0.50$, $1.50$, and $3.00$ log-odds units. The following table gives illustrative posterior output from well-mixed chains; these numbers are hypothetical and are provided to demonstrate interpretation.

Half-normal scale $A$	Posterior median of $\sigma$	Central 90% credible interval	$P(\sigma>0.50\mid\mathbf y)$
0.50	0.62	0.27 to 1.16	0.63
1.50	1.19	0.55 to 2.27	0.89
3.00	1.38	0.61 to 2.76	0.92

Worked Example: Comparing prior scales

Under the $A=0.50$ prior, the posterior median is $0.62$, compared with $1.19$ under $A=1.50$. The difference is $1.19-0.62=0.57$ log-odds units. The posterior probability above the threshold also changes, from $0.63$ to $0.89$, a difference of $0.26$. The broad-prior fit gives a median of $1.38$, only $1.38-1.19=0.19$ higher than the $A=1.50$ fit, with a similar posterior probability above $0.50$ ($0.92$ versus $0.89$).

For these illustrative results, the narrow prior materially pulls the variance component toward smaller values, while the two broader priors give more similar summaries. The appropriate conclusion is that the posterior assessment of heterogeneity is sensitive to the prior scale, particularly when the prior is relatively concentrated. One should report the specified prior and the sensitivity, rather than presenting one summary as if the data alone determined it.

The likelihood-based GLMM approach is a main alternative: it estimates model parameters by likelihood methods rather than combining likelihood with explicit parameter priors. A Bayesian fit is useful here because it yields a full posterior distribution for $\sigma$ and permits direct probability statements about thresholds. Those summaries depend on the prior, especially when the number of groups is small or information about heterogeneity is weak. Conversely, a single likelihood estimate does not express prior sensitivity in the same way. Neither approach removes the need to state the model assumptions and assess whether the data support the intended analysis.

Common pitfalls and limitations

Equating good chain diagnostics with a good model. Converged chains indicate that the sampler explored the specified posterior; they do not show that the normal random-intercept distribution or conditional-independence assumption is appropriate.
Calling prior sensitivity a convergence problem. If chains mix well under each prior but posterior summaries differ across prior scales, the difference reflects the model-plus-prior analyses, not necessarily faulty computation.
Reporting a variance without its scale or uncertainty. State whether the summary is for $\sigma$ or $\sigma^2$, include units, and report a posterior interval or other uncertainty summary.
Interpreting $\sigma$ as a probability difference. It is a standard deviation on the log-odds scale; the probability difference depends on the intercept and starting probability.

With only eight clinics, the data provide limited information about the shape and spread of the random-effect distribution. The posterior can therefore be broad and prior-sensitive even when each clinic contributes several patients. More patient observations within the same small set of clinics do not automatically provide the same information about between-clinic variation as observing more independent clinics. The results also depend on the logistic link, the assumed normal distribution of clinic effects, and the absence of omitted predictors or other dependence structures.

Key takeaway: A Bayesian random-intercept logistic GLMM fit is the combination of a conditional Bernoulli likelihood, a distribution for group effects, and explicit priors. Interpret variance-component summaries only after checking chain behavior, and compare plausible prior scales to show how much the posterior conclusion depends on prior specification.

Check Your Understanding

Answer each question, then compare with the brief solutions.

If a clinic effect is $b_j=0.40$, what is its conditional odds multiplier relative to a clinic with effect zero?
A posterior interval for $\sigma$ is $0.40$ to $1.60$. Give the corresponding interval for $\sigma^2$.
Out of 8,000 retained draws, 6,000 have $\sigma>0.30$. Estimate the posterior probability above that threshold.
Chains have $\widehat R$ close to one, but their traces show long periods in separate ranges. Is the posterior summary ready to report? Explain.
Two prior scales give materially different posterior medians, while both fits have satisfactory chain diagnostics. What does this indicate?

Solutions. (1) The multiplier is $\exp(0.40)\approx1.492$, so the conditional odds are about 1.49 times those for a zero-effect clinic. (2) Square the endpoints: $(0.40^2,1.60^2)=(0.16,2.56)$. (3) The estimate is $6{,}000/8{,}000=0.750$, or 75.0%. (4) No. Chain separation suggests inadequate exploration despite the reported $\widehat R$; investigate mixing and rerun or adjust the computation before interpreting summaries. (5) The posterior conclusion is prior-sensitive over the scales examined; report that sensitivity rather than attributing the difference to chain failure or treating either summary as prior-free.