Random Slopes and Multidimensional Likelihood Integrals

Why a Random Slope Adds an Integration Dimension

A random intercept lets each cluster have its own baseline level. A random slope also lets the effect of a predictor vary across clusters. Both effects are latent: the observed outcomes provide information about them, but they are not observed directly. To obtain a likelihood for the population parameters alone, the cluster-specific effects must be integrated out.

With one Gaussian random intercept, that marginalization is a one-dimensional integral. Adding one Gaussian random slope makes the latent-effects vector two-dimensional, so each cluster’s likelihood involves integration over both coordinates. The observations within a cluster share the same intercept and slope; the slope is not redrawn for each observation.

Definition: For cluster $i$, let $b_i=(b_{0i},b_{1i})^T$ be a Gaussian random-effects vector, where $b_{0i}$ is the random intercept and $b_{1i}$ is the random slope for a scalar predictor $z_{ij}$. The linear predictor for observation $j$ is $$ \eta_{ij}=x_{ij}^T\beta+b_{0i}+b_{1i}z_{ij}. $$ Here $x_{ij}$ is the fixed-effects design vector and $\beta$ its coefficient vector. The two-dimensional cluster marginal likelihood is the integral of the conditional response likelihood, multiplied by the joint density of $b_i$, over $(b_{0i},b_{1i})\in\mathbb{R}^2$.

Model and Assumptions

Suppose there are $m$ independent clusters, with $n_i$ observations in cluster $i$. Conditional on $b_i$, the responses within cluster $i$ are independent and have specified densities or probability mass functions $f(y_{ij}\mid b_i;\beta,\phi)$. The distribution of $b_i$ is Gaussian with mean zero and a positive-definite covariance matrix $D$. Random-effects vectors from different clusters are independent. These assumptions imply that observations in a cluster may be dependent marginally, because they share $b_i$, even though they are independent conditional on it.

Write $\phi_2(b_i;0,D)$ for the bivariate Gaussian density. For a positive-definite $D$, its standard form is $$ \phi_2(b_i;0,D)=\frac{1}{2\pi |D|^{1/2}}\exp\left(-\frac{1}{2}b_i^TD^{-1}b_i\right). $$ The response model, link function, and any exposure or trial-count terms determine the conditional likelihood. They do not change the fact that the random-effects integration is over two coordinates.

Deriving the Two-Dimensional Integrand

The conditional likelihood for the observations in cluster $i$, at a specified value of its random effects, is $$ L_i(y_i\mid b_i;\beta,\phi)=\prod_{j=1}^{n_i} f(y_{ij}\mid b_i;\beta,\phi). $$ The joint density of the responses and random effects factors into this conditional likelihood times the random-effects density. Integrating the joint density over all possible values of the unobserved $b_i$ gives the marginal likelihood:

L_i(\beta,\phi,D) =\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \left\{\prod_{j=1}^{n_i} f(y_{ij}\mid b_{0i},b_{1i};\beta,\phi)\right\} \phi_2\!\left((b_{0i},b_{1i})^T;0,D\right) \,db_{1i}\,db_{0i}.

Thus the integrand is the product of two pieces: the likelihood of the cluster’s observed data at the proposed random-effects values, and the Gaussian density assigned to those values. The order of integration is immaterial when the integral exists; the displayed order simply makes the two coordinates explicit. Since the clusters are independent after integrating their separate random-effects vectors, the full likelihood is

L(\beta,\phi,D)=\prod_{i=1}^{m} L_i(\beta,\phi,D).

This is a marginal likelihood, not the conditional likelihood evaluated at $b_i=0$. Setting $b_i=0$ uses only one point in the plane; marginalization averages the conditional likelihood over the full Gaussian distribution of both effects.

Key takeaway: A random intercept and one random slope produce a two-dimensional integral because the shared latent vector has two components. The number of observations in a cluster changes the product inside the integral, not the dimension of the integral.

Worked Integrands

Worked Example: Binary Outcomes With a Logistic Random Slope

Consider one hypothetical cluster with three binary outcomes $y=(0,1,1)$ observed at $z=(-1,0,1)$. Suppose $$ \operatorname{logit}(p_{ij})=-0.2+0.8z_{ij}+b_{0i}+b_{1i}z_{ij}, $$ where the random-effects vector has covariance $D=\operatorname{diag}(0.25,0.16)$. Conditional on the random effects, assume the three outcomes are independent Bernoulli variables. These assumptions give a valid cluster likelihood and a bivariate Gaussian density.

At the point $b_{0i}=0.1,\ b_{1i}=-0.3$, the three linear predictors are $$ \eta_1=-0.2-0.8+0.1+0.3=-0.6,\qquad \eta_2=-0.2+0.1=-0.1,\qquad \eta_3=-0.2+0.8+0.1-0.3=0.4. $$

Applying $p=\exp(\eta)/\{1+\exp(\eta)\}$ gives $p_1\approx0.3543$, $p_2\approx0.4750$, and $p_3\approx0.5987$. The conditional likelihood at this point is $$ (1-p_1)p_2p_3 =(1-0.3543)(0.4750)(0.5987) \approx0.1836. $$ The random-effects density at this point is $$ \phi_2(b;0,D) =\frac{1}{2\pi(0.5)(0.4)} \exp\left[-\frac12\left\{\left(\frac{0.1}{0.5}\right)^2+ \left(\frac{-0.3}{0.4}\right)^2\right\}\right] \approx0.5888. $$ Consequently, the full integrand at this point is approximately $0.1836(0.5888)=0.1081$. This is the integrand’s value at one pair of random-effects values, not the marginal likelihood; the latter requires integrating the integrand over the plane.

Worked Example: Poisson Counts With Exposure

For a hypothetical cluster, suppose counts $y=(0,2,1)$ are observed at $z=(0,1,2)$, with exposures $e=(1,2,1.5)$. Use the conditional Poisson model $$ Y_{ij}\mid b_i\sim\operatorname{Poisson}(\lambda_{ij}),\qquad \log(\lambda_{ij})=\log(e_{ij})-0.2+0.3z_{ij}+b_{0i}+b_{1i}z_{ij}. $$ Assume conditional independence within the cluster and Gaussian random effects with a positive-definite covariance matrix. The exposure enters as an offset, and both latent effects are shared across the three observations.

At $b_{0i}=0.1,\ b_{1i}=0.2$, the three log rates apart from exposure are $-0.1,\ 0.4,\ 0.9$. Thus $$ \lambda_1=1\exp(-0.1)\approx0.9048,\qquad \lambda_2=2\exp(0.4)\approx2.9836,\qquad \lambda_3=1.5\exp(0.9)\approx3.6894. $$

The conditional likelihood is $$ \frac{e^{-0.9048}(0.9048)^0}{0!} \frac{e^{-2.9836}(2.9836)^2}{2!} \frac{e^{-3.6894}(3.6894)^1}{1!} \approx0.008402. $$ The numerator and factorial terms are the Poisson probabilities for $0,2,1$, respectively. The cluster integrand at this random-effects point is therefore $0.008402\,\phi_2((0.1,0.2)^T;0,D)$. The density factor must still be included before integrating; the conditional likelihood alone is not the marginal likelihood.

Worked Example: Grouped Binomial Responses

Suppose a hypothetical cluster contains two groups of Bernoulli trials. At $z=0$, there are $n_1=3$ trials and $y_1=2$ successes; at $z=1$, there are $n_2=2$ trials and $y_2=1$ success. Use a logistic model with fixed coefficients $\beta_0=-0.2,\ \beta_1=0.3$, and random intercept and slope. Conditional on the random effects, assume the trial counts are binomial and the two group outcomes are independent.

At $b_0=0.1,\ b_1=0.3$, the logits are $-0.2+0.1=-0.1$ at $z=0$ and $-0.2+0.3+0.1+0.3=0.5$ at $z=1$. Therefore $p_1=\operatorname{logit}^{-1}(-0.1)\approx0.4750$ and $p_2=\operatorname{logit}^{-1}(0.5)\approx0.6225$. The grouped conditional likelihood, including the binomial coefficients, is $$ {3\choose2}(0.4750)^2(0.5250) {2\choose1}(0.6225)(0.3775) =3(0.4750)^2(0.5250)\times2(0.6225)(0.3775) \approx0.1670. $$ The corresponding integrand is this conditional likelihood multiplied by $\phi_2((0.1,0.3)^T;0,D)$. Integrating that product over both random effects gives the group’s marginal likelihood.

Why the Dimension and the Product Matter

The derivation follows directly from the shared-effects structure. Conditional independence gives a product over observations, while the common $b_i$ remains the same in every factor. The Gaussian density is also joint: it weights a pair $(b_{0i},b_{1i})$, allowing the two effects to be associated when $D$ has nonzero off-diagonal entries. Whether the covariance matrix has zero or nonzero covariance does not change the two-dimensional domain.

A useful comparison is a random-intercept-only model. Its latent-effects vector has one component, so the analogous marginal likelihood is a single integral. With an intercept and one slope, it is a double integral. More generally, a model with $q$ jointly Gaussian random effects has a $q$-dimensional integral. This dimension is determined by the latent vector, not by the cluster size.

Another distinction is between conditional and marginal analysis. A conditional GLM likelihood treats the random effects as given; the GLMM marginal likelihood accounts for their distribution by integration. GEE instead specifies a marginal mean and uses a working correlation for inference, rather than this random-effects likelihood. These approaches target different model structures and should not be treated as interchangeable calculations.

Diagnostics, Limitations, and Common Errors

Before evaluating or approximating the integral, check that the response likelihood matches the outcome type, the linear predictor uses the intended slope covariate, and the random-effects density is a proper joint Gaussian density. The covariance matrix must be positive definite. Also verify that the same $b_{0i}$ and $b_{1i}$ appear in every row of cluster $i$, and that the response factors correspond to the observed outcomes, exposures, or trial counts.

Data structure can make the model difficult to estimate even when the integral is correctly defined. If $z_{ij}$ barely varies within clusters, the observations may provide little information for separating a cluster-specific intercept from a cluster-specific slope. Sparse clusters can also yield weak information about variation in slopes. These are identification and information concerns; they do not turn the two-dimensional integral into a one-dimensional one. Numerical evaluation may be challenging when likelihood values are very small or concentrated in a narrow region, so implementations need stable integration methods and suitable diagnostics.

Wrong: Multiply the separate marginal probabilities for each observation, integrating a different random slope for every row. Correction: The cluster shares one random-effects vector; multiply the conditional response factors first, then integrate that product once over the cluster’s vector.
Wrong: Call the conditional likelihood at $b_i=0$ the marginal likelihood. Correction: The marginal likelihood integrates over all $(b_{0i},b_{1i})$, with the Gaussian density as a weight.
Wrong: Omit the Gaussian density and integrate only the response likelihood over the plane. Correction: The density is essential: it defines the distribution over which the random effects are marginalized.
Wrong: Assume that more observations mean a higher-dimensional integral. Correction: More observations add factors inside the integrand; the number of random effects determines its dimension.

Laplace approximation and Gaussian quadrature are important alternatives to direct evaluation when exact integration is unavailable, but they approximate this same marginal likelihood. They do not change the integrand’s definition or remove the need to specify the shared effects and their joint density correctly.

Check Your Understanding

Answer each question, then compare with the solutions below.

A model has a random intercept and random slope for one predictor. What is the latent-effects vector, and over what domain is each cluster likelihood integrated?
Why is the conditional response likelihood a product across observations, but not a product of separately marginalized observation likelihoods?
For one cluster with two Poisson counts, write the general form of its integrand if the log mean includes exposure, a fixed predictor, a random intercept, and a random slope.
A random-intercept-only model is changed to include one random slope. What changes in the marginal likelihood, and what does not change when the cluster gains additional observations?

Solutions. (1) The vector is $b_i=(b_{0i},b_{1i})^T$, and the integration domain is $\mathbb{R}^2$, or equivalently the full $(b_{0i},b_{1i})$ plane. (2) Conditional independence gives a product when the same cluster effects are fixed. Those effects are shared across observations, so the product must be formed before integrating the joint density once; separate integrations would incorrectly assign separate effects to the observations. (3) If $\lambda_{ij}=e_{ij}\exp(x_{ij}^T\beta+b_{0i}+b_{1i}z_{ij})$, the integrand is $\prod_{j=1}^{2}\{e^{-\lambda_{ij}}\lambda_{ij}^{y_{ij}}/y_{ij}!\}\,\phi_2((b_{0i},b_{1i})^T;0,D)$. (4) The integral changes from one dimension to two and gains a random-slope contribution in each observation’s linear predictor. Additional observations add conditional-likelihood factors; they do not increase the integration dimension.