Constructing the Integrand for a Binary Random-Intercept GLMM

From conditional outcomes to a cluster integrand

In a binary random-intercept GLMM, observations in a cluster share a latent intercept. Conditional on that intercept, their outcomes are modeled with Bernoulli probabilities; marginally, the shared intercept induces dependence among them. To form the likelihood for the model parameters, we combine the cluster’s conditional likelihood with the random-intercept density and integrate over the possible intercept values.

The focus here is the integrand being integrated—not an approximation to the integral or an algorithm for finding a value that makes the integrand large. Getting the integrand right requires keeping three pieces distinct: the conditional Bernoulli probabilities, the fact that one intercept is shared within a cluster, and the density of that intercept.

Definition: For clusters $i=1,\ldots,M$, let $Y_{ij}\in\{0,1\}$ be the outcome for observation $j=1,\ldots,n_i$. Let $\mathbf{x}_{ij}$ be a fixed covariate vector and $\boldsymbol{\beta}$ its coefficient vector. A binary random-intercept logit GLMM specifies $$ Y_{ij}\mid b_i,\mathbf{x}_{ij}\sim \operatorname{Bernoulli}(p_{ij}),\qquad \operatorname{logit}(p_{ij})=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_i,\qquad b_i\sim N(0,\sigma_b^2). $$ The random intercepts are independent across clusters. Conditional on the random intercepts and covariates, all observations are independent. The variance satisfies $\sigma_b^2>0$.

Write $\eta_{ij}(b)=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b$. Then $p_{ij}(b)=\exp\{\eta_{ij}(b)\}/[1+\exp\{\eta_{ij}(b)\}]$. This is a conditional probability: it describes the outcome probability at a specified value of the cluster intercept. The marginal likelihood contribution for a cluster requires averaging its entire conditional likelihood over the distribution of its one shared intercept.

Deriving the integrand and its log form

For a Bernoulli outcome $y$, the probability mass at $y\in\{0,1\}$ is $p^y(1-p)^{1-y}$. Under the logit link, substituting $p=\exp(\eta)/[1+\exp(\eta)]$ gives

p^y(1-p)^{1-y} = \left(\frac{\exp(\eta)}{1+\exp(\eta)}\right)^y \left(\frac{1}{1+\exp(\eta)}\right)^{1-y} = \frac{\exp(y\eta)}{1+\exp(\eta)}.

Thus the conditional likelihood for cluster $i$, at a specified value $b$, is the product of these Bernoulli probabilities over the observations in that cluster. The density of $b_i$ is $f(b)=\exp[-b^2/(2\sigma_b^2)]/(\sqrt{2\pi}\sigma_b)$. Multiplying the conditional likelihood by this density gives the integrand; integrating it over the real line gives the marginal cluster likelihood.

Definition: The one-dimensional marginal-likelihood integrand for cluster $i$ is $$ g_i(b) = \left\{\prod_{j=1}^{n_i} \frac{\exp\!\left[Y_{ij}\eta_{ij}(b)\right]} {1+\exp\!\left[\eta_{ij}(b)\right]}\right\} \frac{1}{\sqrt{2\pi}\sigma_b} \exp\!\left(-\frac{b^2}{2\sigma_b^2}\right), \qquad -\infty<b<\infty. $$ The cluster’s marginal likelihood contribution is $L_i=\int_{-\infty}^{\infty}g_i(b)\,db$. Under the stated independence assumptions, the full marginal likelihood is $\prod_{i=1}^{M}L_i$.

The same construction is often easier to inspect on the log scale. Taking the logarithm of the Bernoulli probability contributes $Y_{ij}\eta_{ij}(b)-\log[1+\exp\{\eta_{ij}(b)\}]$. Taking the logarithm of the Gaussian density contributes a quadratic term and a normalization term. Therefore, defining $h_i(b)=\log g_i(b)$,

h_i(b) = \sum_{j=1}^{n_i} \left\{ Y_{ij}\eta_{ij}(b) -\log\!\left[1+\exp\{\eta_{ij}(b)\}\right] \right\} -\frac{b^2}{2\sigma_b^2} -\log(\sqrt{2\pi}\sigma_b).

This is the log of the integrand, not the log of the marginal likelihood. In particular, the marginal log-likelihood contribution is $\log\{\int \exp[h_i(b)]\,db\}$, which is generally not the same as $\int h_i(b)\,db$. The integration is over the common random intercept $b$, and the normalizing constant in its density is part of the integrand.

Worked Example: Three Outcomes in an Intercept-Only Cluster

Consider one hypothetical cluster with outcomes $Y=(1,0,1)$, an intercept-only fixed effect $\beta_0=\log 2$, and $\sigma_b=1$. At $b=0$, each linear predictor is $\eta=\log 2$, so each conditional success probability is $2/3$. The conditional likelihood is $$ (2/3)(1-2/3)(2/3)=\frac{4}{27}=0.1481. $$

The Gaussian density at zero is $\exp(0)/\sqrt{2\pi}=0.3989$. The integrand at this point is the product: $$ g(0)=\frac{4}{27}\times 0.3989=0.05910. $$

On the log scale, the conditional log-likelihood is $2\log(2/3)+\log(1/3)=\log(4/27)=-1.9095$. The log density is $-\log(\sqrt{2\pi})=-0.9189$, so $h(0)=-1.9095-0.9189=-2.8285$, and $\exp\{h(0)\}=0.05910$, up to rounding. This value is the integrand at $b=0$, not the marginal likelihood: the latter requires integrating over all $b$.

Reading and evaluating the log-scale expression

The sum in $h_i(b)$ accounts for every observed outcome in the cluster. If $Y_{ij}=1$, its contribution is $\eta_{ij}(b)-\log[1+\exp\{\eta_{ij}(b)\}]$, the log of the conditional success probability. If $Y_{ij}=0$, its contribution is $-\log[1+\exp\{\eta_{ij}(b)\}]$, the log of the conditional failure probability. The Gaussian terms then weight that conditional likelihood according to how plausible $b$ is under the assumed random-effects distribution.

When evaluating the expression numerically, $\log[1+\exp(\eta)]$ is sometimes called the softplus function. Directly computing $\exp(\eta)$ can overflow for very large positive $\eta$, even when the desired log probability is well defined. Stable software implementations evaluate this function without forming an unnecessarily large exponential. The mathematical expression remains the same.

Worked Example: Different Covariates and a Nonzero Intercept Value

Suppose a hypothetical cluster has two observations with $Y_1=0$, $Y_2=1$, covariates $x_1=0$, $x_2=1$, and fixed-effects model $\eta_j(b)=\beta_0+\beta_1x_j+b$, where $\beta_0=-0.4$, $\beta_1=0.8$, and $\sigma_b=0.7$. Evaluate the integrand at $b=0.2$.

The linear predictors are $\eta_1=-0.4+0+0.2=-0.2$ and $\eta_2=-0.4+0.8+0.2=0.6$. The observed conditional probabilities are $P(Y_1=0\mid b)=1/[1+\exp(-0.2)]=0.5498$ and $P(Y_2=1\mid b)=\exp(0.6)/[1+\exp(0.6)]=0.6457$. Their product is approximately $0.3550$.

The random-intercept density at $0.2$ is $$ \frac{1}{0.7\sqrt{2\pi}}\exp\!\left(-\frac{0.2^2}{2(0.7^2)}\right) = \frac{1}{1.7546}\exp(-0.04082) \approx 0.5471. $$ Thus $g(0.2)\approx 0.3550\times0.5471=0.1942$.

The log-scale calculation gives the same result. The Bernoulli terms sum to $-\log[1+\exp(-0.2)]+0.6-\log[1+\exp(0.6)]\approx-1.0353$. The log density is $-\log(0.7\sqrt{2\pi})-0.2^2/[2(0.7^2)]\approx-0.6031$. Hence $h(0.2)\approx-1.6384$, and $\exp\{h(0.2)\}\approx0.1942$.

One shared intercept means one integral per cluster

The product over $j$ belongs inside the integral because a single value $b_i$ affects all observations in cluster $i$. The cluster contribution is therefore an average of the joint conditional likelihood across possible shared intercepts. Integrating each observation separately and multiplying those integrals would instead give each observation its own independent random intercept. That is a different model and generally a different likelihood.

Under the assumptions above, different clusters have independent random intercepts and conditionally independent outcomes. Integrating each cluster’s intercept therefore produces a marginal likelihood $L_i$ for that cluster, and the full likelihood is the product of those contributions. Within a cluster, however, the observations are generally dependent after $b_i$ is integrated out. The integrand construction preserves that dependence.

Worked Example: Identical Covariates and the Grouped-Count Constant

Suppose a hypothetical cluster contains four individually recorded outcomes $1,1,0,0$, all with the same intercept-only linear predictor $\eta(b)=b$. Let $\sigma_b=1$, and evaluate at $b=\log 2$. Then $p=2/3$, and the likelihood of this particular ordered set of four outcomes is $$ (2/3)^2(1/3)^2=\frac{4}{81}=0.04938. $$

At $b=\log 2$, the normal density is $$ \frac{1}{\sqrt{2\pi}}\exp\!\left(-\frac{(\log 2)^2}{2}\right) =0.3989\exp(-0.2402)\approx0.3137. $$ Therefore $g(\log 2)\approx(4/81)(0.3137)=0.01549$. On the log scale, the Bernoulli contribution is $2\log(2/3)+2\log(1/3)=\log(4/81)\approx-3.0082$; adding the log density, approximately $-1.1592$, gives $h(\log 2)\approx-4.1673$.

If the same data are represented only as a grouped count of two successes out of four, the binomial probability for that count includes the coefficient $\binom{4}{2}=6$. That coefficient is absent from the likelihood of the particular ordered outcome sequence above. Either representation can be used, but its likelihood must match the way the data are represented; omitting or adding the coefficient inconsistently changes the integrand by a factor of six.

Checks, alternatives, and limitations

Before an integrand is used in a likelihood calculation, check that it matches the data structure and model assumptions. In particular, verify that outcomes are coded as zero or one, that the fixed part of each linear predictor uses the intended covariates, and that every observation in a cluster shares the same $b$. The Gaussian density must use the stated variance and standard deviation consistently: its exponent contains $2\sigma_b^2$, while its normalizing factor contains $\sigma_b$.

Sign check. A success contributes $\eta-\log(1+\exp\eta)$; a failure contributes $-\log(1+\exp\eta)$. Reversing these terms swaps success and failure probabilities.
Density check. The normal density is a density over $b$, so the factor $1/(\sqrt{2\pi}\sigma_b)$ is required. Leaving it out changes the likelihood as a function of $\sigma_b$.
Dependence check. There is one integration variable per cluster, not one per observation, for a random-intercept model.
Scale check. $h_i(b)$ is a log integrand. It must be exponentiated before it is integrated on the original $b$ scale.

A tempting alternative is to set $b_i=0$ and use the resulting conditional likelihood. This evaluates the conditional model at the mean of the random-effect distribution; it does not average the likelihood over that distribution and is not generally the marginal likelihood. A GLM without random effects makes a different modeling assumption, while GEE targets population-averaged inference through a different estimating approach. Neither should be substituted for the random-effects integral without considering the inferential target and assumptions. Laplace approximation and numerical quadrature are methods for evaluating integrals such as $L_i$; they do not change the integrand that the model specifies.

The construction also has limits. It assumes a Gaussian random-intercept distribution, a logit conditional mean, and conditional Bernoulli independence. If the random-effect distribution, link, or dependence structure differs, the density or conditional likelihood—and therefore the integrand—must change. A correct expression is necessary for estimation, but it does not establish that the model is substantively adequate.

Key takeaway: For each cluster, multiply all conditional Bernoulli probabilities evaluated at the same intercept $b$ by the Gaussian density of $b$. The resulting one-dimensional function $g_i(b)$ is integrated over the real line to obtain that cluster’s marginal likelihood contribution.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

In a cluster with three observations, how many random-intercept integration variables appear in its marginal likelihood contribution, and why?
For a failure $Y_{ij}=0$, write its contribution to the log integrand from the Bernoulli likelihood.
Explain why $\exp\{h_i(0)\}$ is generally not the marginal likelihood contribution $L_i$.
If a model uses $b_i\sim N(0,\sigma_b^2)$, identify the correct normalizing factor and quadratic exponent in its density.

Solutions. (1) One variable, $b_i$, because the random intercept is shared by all observations in the cluster; their conditional likelihoods are multiplied before that common intercept is integrated out. (2) The contribution is $-\log[1+\exp\{\eta_{ij}(b)\}]$, the log of the conditional failure probability. (3) $\exp\{h_i(0)\}=g_i(0)$ is the integrand at one point. The marginal contribution is the area $\int_{-\infty}^{\infty}g_i(b)\,db$, not that single point value. (4) The density is $\{\sqrt{2\pi}\sigma_b\}^{-1}\exp[-b_i^2/(2\sigma_b^2)]$; its normalizing factor uses the standard deviation $\sigma_b$, and its quadratic exponent uses the variance $\sigma_b^2$.