Latent-Variable Augmentation for Binary GLMMs

Why Introduce a Latent Response?

A binary observation contributes a Bernoulli likelihood to a generalized linear mixed model. For a probit model, that likelihood involves the normal distribution function, so it is not a Gaussian likelihood in the regression coefficients and random effects. Latent-variable augmentation gives each binary observation an unobserved continuous response whose threshold determines the recorded outcome. Conditional on those latent responses, the regression part of the model is Gaussian, making several posterior updates easier to derive and sample.

The construction below is specifically for a probit binary GLMM. It is not a claim that a logistic likelihood is Gaussian after the same construction. The logistic model has different augmentation methods, including the Pólya-Gamma method considered separately.

Definition: For observation $i$ in group $j$, let $y_{ij}\in\{0,1\}$ be the observed outcome and let $w_{ij}$ be a latent continuous response. In a probit GLMM, $$ w_{ij}=x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_j+\epsilon_{ij},\qquad \epsilon_{ij}\sim N(0,1),\qquad y_{ij}=\mathbb{1}(w_{ij}>0). $$ Here $x_{ij}$ and $z_{ij}$ are known fixed-effect and random-effect design vectors, $\beta$ is the fixed-effect coefficient vector, and $b_j$ is the group-specific random-effect vector. A common assumption is $b_j\sim N(0,D)$, independently across groups, where $D$ is a covariance matrix.

Conditional on $\beta$ and $b_j$, define $\eta_{ij}=x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_j$. Then $w_{ij}\sim N(\eta_{ij},1)$, and thresholding gives $$ P(y_{ij}=1\mid \beta,b_j)=P(w_{ij}>0\mid \beta,b_j)=\Phi(\eta_{ij}). $$ Thus the observed-data likelihood contribution is $\Phi(\eta_{ij})$ when $y_{ij}=1$, and $1-\Phi(\eta_{ij})=\Phi(-\eta_{ij})$ when $y_{ij}=0$. The augmentation retains this likelihood after the latent responses are integrated out; it does not change the model for the observed data.

Assumptions and the Conditional Updates

The basic Gibbs construction assumes observations are conditionally independent given the regression parameters and group effects, the latent residuals are independent standard normal variables, and the outcome is determined by the stated zero threshold. It also needs specified priors. For the standard Gaussian updates, take $\beta\sim N(m_0,V_0)$, with positive-definite $V_0$, independent of the random effects and their covariance parameters. These assumptions must match the scientific model; augmentation is a computational representation, not a repair for a misspecified link or dependence structure.

Given the observed $y_{ij}$ and current parameters, each latent response is a normal draw restricted to the half-line consistent with the observation: $$ w_{ij}\mid y_{ij},\beta,b_j\sim \begin{cases} N(\eta_{ij},1)\text{ truncated to }(0,\infty),&y_{ij}=1,\\ N(\eta_{ij},1)\text{ truncated to }(-\infty,0],&y_{ij}=0. \end{cases} $$ The truncation is essential: drawing an unrestricted normal value and then retaining the observed binary outcome would not produce the conditional distribution of $w_{ij}$ given $y_{ij}$.

Once the latent responses are included, the regression equation is a normal linear model. Stack the observations into vectors $w$ and $\eta$, and let $X$ and $Z$ denote the corresponding design matrices. Conditional on the group effects, $w-Zb=X\beta+\epsilon$, where $\epsilon\sim N(0,I)$. Combining this likelihood with the normal prior gives $$ V_\beta=(V_0^{-1}+X^{\mathsf T}X)^{-1},\qquad m_\beta=V_\beta\{V_0^{-1}m_0+X^{\mathsf T}(w-Zb)\}, $$ so $\beta\mid w,b\sim N(m_\beta,V_\beta)$.

For a group with design matrix $Z_j$, the same argument gives, conditional on $\beta,D$, $$ V_j=(D^{-1}+Z_j^{\mathsf T}Z_j)^{-1},\qquad m_j=V_j Z_j^{\mathsf T}(w_j-X_j\beta), $$ and $b_j\mid w_j,\beta,D\sim N(m_j,V_j)$. Here $w_j$ is the vector of latent responses in group $j$, and $X_j$ is its fixed-effect design matrix. If a variance component has a conjugate prior, its update may also remain conjugate conditional on the current random effects; otherwise its update needs a suitable nonconjugate sampling step.

Why the Gaussian Updates Follow

The key calculation is completing the square in the product of a Gaussian likelihood and a Gaussian prior. Write $r=w-Zb$. Up to a factor that does not depend on $\beta$, the conditional density is $$ \exp\left\{-\frac{1}{2}(r-X\beta)^{\mathsf T}(r-X\beta) -\frac{1}{2}(\beta-m_0)^{\mathsf T}V_0^{-1}(\beta-m_0)\right\}. $$ The terms quadratic in $\beta$ have precision matrix $X^{\mathsf T}X+V_0^{-1}$. The linear terms are $\beta^{\mathsf T}(X^{\mathsf T}r+V_0^{-1}m_0)$. Completing the square therefore gives precision $V_\beta^{-1}=X^{\mathsf T}X+V_0^{-1}$ and mean $V_\beta(X^{\mathsf T}r+V_0^{-1}m_0)$, exactly the conditional distribution above. Applying the same calculation to each $b_j$, with prior precision $D^{-1}$, gives the group-effect update.

Worked Example: Drawing a Latent Response

Suppose one hypothetical observation has $y=1$ and current linear predictor $\eta=0.5$. Under the model, $w\sim N(0.5,1)$ before conditioning on the observed outcome. Because $y=1$ means $w>0$, the conditional draw must come from this normal distribution truncated above zero.

The probability of observing $y=1$ at this predictor is $\Phi(0.5)\approx0.6915$. To illustrate an inverse-CDF draw, take a hypothetical uniform draw $u=0.5$. The conditional probability scale for a lower-truncated normal is $$ \Phi(-0.5)+u\{1-\Phi(-0.5)\} \approx 0.3085+0.5(0.6915)=0.6543. $$ The corresponding standard normal quantile is approximately $\Phi^{-1}(0.6543)=0.397$, so the latent draw is $w=0.5+0.397=0.897$. It satisfies $w>0$, as required by the observed $y=1$.

This is one valid conditional draw, not an estimate of the latent response. In computation, tail-stable truncated-normal algorithms may be needed when the truncation point lies far into a distribution's tail.

Worked Example: Updating a Fixed Intercept

Consider a hypothetical model with an intercept $\beta$, no random effects, and two latent responses currently equal to $1.2$ and $-0.4$. The associated observed outcomes are $1$ and $0$, respectively, so both latent values agree with their thresholds. Let the prior be $\beta\sim N(0,1)$. The design matrix is a two-row column of ones, and the latent-response model is $w_i=\beta+\epsilon_i$.

The posterior precision is $1+2=3$, since the prior precision is $1$ and there are two unit-variance observations. The posterior variance is $1/3$, and the posterior mean is $$ \frac{1(0)+1.2+(-0.4)}{3}=\frac{0.8}{3}=0.2667. $$ Therefore $\beta\mid w\sim N(0.2667,\,0.3333)$. The second argument is the variance, not the standard deviation. This is a conditional update given the current latent responses; it is not the posterior distribution based only on the two binary outcomes.

Worked Example: Updating a Group Intercept

Suppose a hypothetical group has two observations, both with fixed-effect predictor $\beta=0.3$. The current latent responses are $w_1=0.8$ and $w_2=1.1$. Let its random intercept have prior $b_j\sim N(0,1)$, so $D=1$. Both latent responses are positive and would correspond to binary outcomes of one.

After subtracting the fixed predictor, the residual values are $0.5$ and $0.8$. The random-intercept posterior precision is $1+2=3$, giving variance $1/3$, and its mean is $$ \frac{0.5+0.8}{3}=\frac{1.3}{3}=0.4333. $$ Thus $b_j\mid w,\beta,D\sim N(0.4333,\,0.3333)$. The update combines the two residuals with the zero-centered prior; it does not treat the group intercept as an unconstrained fixed parameter.

A Gibbs Sweep, Diagnostics, and Limitations

A basic Gibbs sweep cycles through the conditional distributions. With the model and priors specified, one possible order is to draw all latent responses given the binary outcomes and current parameters, draw $\beta$ given the latent responses and group effects, draw each $b_j$ given the latent responses and current covariance parameters, and then update covariance parameters from their full conditionals. Repeating these steps produces dependent draws whose stationary distribution is the joint posterior, provided the updates are valid and the chain is run appropriately. Conjugacy for the Gaussian regression blocks does not guarantee that every block, or the entire chain, will mix quickly.

For example, in a random-intercept model with $b_j\sim N(0,\tau^2)$ and independent inverse-gamma prior $\tau^2\sim\operatorname{IG}(a_0,c_0)$, using density proportional to $v^{-(a_0+1)}\exp(-c_0/v)$, the conditional update is $$ \tau^2\mid b\sim \operatorname{IG}\left(a_0+\frac{J}{2},\ c_0+\frac{1}{2}\sum_{j=1}^{J}b_j^2\right). $$ This illustrates a conjugate variance update conditional on the current group effects. Other prior choices need not yield this inverse-gamma conditional, so one should derive the update from the chosen prior rather than assume conjugacy.

Several checks help identify problems:

Threshold consistency: every sampled $w_{ij}$ must be positive when $y_{ij}=1$, and nonpositive when $y_{ij}=0$. Violations indicate an incorrect truncation step.
Tail computation: when $\eta_{ij}$ strongly conflicts with the observed outcome, the required truncated normal may be an extreme-tail distribution. Naive CDF subtraction or inverse-CDF code can lose numerical precision.
Mixing: latent responses, fixed effects, and random effects can be strongly dependent in their joint posterior. Autocorrelated draws or slow movement in variance components can persist even though every conditional update is exact. Trace plots and other chain diagnostics are useful; the conditional formulas alone do not establish adequate mixing.
Model scale: the residual variance is fixed at one to identify the probit scale. Estimating it freely without an identifying restriction changes the model and leaves the latent scale confounded with the regression coefficients.
Link and interpretation: the coefficients belong to a probit linear predictor. They are not logistic regression coefficients, and changing the link changes the likelihood and interpretation.

The main computational alternative is to sample parameters using the observed binary likelihood directly, for example with a Metropolis-Hastings step. That avoids latent responses but can require careful proposals for correlated parameter blocks. For a logistic random-effects model, Pólya-Gamma augmentation is a distinct method designed for the logistic likelihood; the probit truncated-normal construction above should not be substituted for it. Latent probit augmentation is attractive when Gaussian conditional updates are valuable, but it introduces extra variables and may mix poorly in some data and parameter regimes. It also does not remove the need to choose defensible priors or assess the fitted model.

Key takeaway: In a probit GLMM, thresholding a unit-variance latent normal response reproduces the binary likelihood. Given those latent responses, normal priors yield Gaussian updates for fixed and random effects; the latent-response step itself is a truncated-normal draw.

Check Your Understanding

Answer each question, then check the solutions below.

In a probit GLMM, what restriction does $y_{ij}=0$ impose on the latent response $w_{ij}$?
For a fixed-intercept model with prior $\beta\sim N(0,1)$ and one current latent response $w=0.6$, what are the conditional posterior mean and variance of $\beta$?
Why must the latent response be drawn from a truncated normal rather than an unrestricted normal?
What assumption fixes the latent-response residual variance at one, and why does it matter?
Does the same thresholded-normal augmentation automatically provide Gaussian conditional updates for a logistic GLMM? Explain.

Solutions. (1) It requires $w_{ij}\leq0$. (2) The posterior precision is $1+1=2$, so the variance is $1/2$; the mean is $(0+0.6)/2=0.3000$. (3) The conditional distribution must respect the observed binary outcome, which records on which side of zero the latent response lies. (4) The probit scale is identified by fixing the residual variance to one; otherwise its scale is confounded with the coefficients and latent-response scale. (5) No. The construction derives from a normal latent response and gives the probit likelihood. Logistic models require a method matched to the logistic likelihood, such as Pólya-Gamma augmentation.

Latent-Variable Augmentation for Binary GLMMs

What You'll Learn