Quadrature with Random-Intercept and Random-Slope Covariance

Why covariance changes the quadrature calculation

A random intercept shifts a cluster’s linear predictor, while a random slope changes how that predictor responds to a covariate. When these effects are correlated, their joint distribution is tilted: clusters with larger intercepts may tend to have smaller slopes, or vice versa. Quadrature must account for that joint distribution, not integrate the two effects as if they were independent.

The key operation is to express the two correlated effects as a linear transformation of independent standard normal coordinates. Quadrature is then applied in those standard normal coordinates, with the transformation carrying the covariance into each evaluation point. In adaptive quadrature, a second transformation scales and orients the nodes to match the local curvature of the cluster’s integrand.

Definition: For cluster $j$, let $b_j=(b_{0j},b_{1j})^\mathsf{T}$ denote a random intercept and random slope, with $$ b_j\sim N_2(0,\Sigma),\qquad \Sigma= \begin{pmatrix} \sigma_0^2&\sigma_{01}\\ \sigma_{01}&\sigma_1^2 \end{pmatrix}. $$ For observation $i$ in that cluster, write the linear predictor as $\eta_{ij}=x_{ij}^\mathsf{T}\beta+b_{0j}+z_{ij}b_{1j}$, where $x_{ij}$ is the fixed-effect design vector and $z_{ij}$ is the covariate with a random slope. The conditional response distribution and link function determine the conditional likelihood $L_j(b_j)$. The cluster’s marginal likelihood contribution is $I_j=\int L_j(b_j)\phi_2(b_j;0,\Sigma)\,db_j$.

Here $\phi_2(\,\cdot\,;0,\Sigma)$ is the bivariate normal density. The standard setup assumes clusters are independent given the model parameters, responses within a cluster are conditionally independent given its random effects, and the random effects have the stated mean-zero Gaussian distribution. The covariance matrix must be positive definite for the ordinary two-dimensional Cholesky transformation below.

Covariance transformation and nonadaptive nodes

Let $\Sigma=LL^\mathsf{T}$, where $L$ is its lower-triangular Cholesky factor, and let $u_j\sim N_2(0,I)$. Then $b_j=Lu_j$ has covariance $\Sigma$. Writing $L$ explicitly shows how covariance enters every quadrature point:

L= \begin{pmatrix} \sigma_0&0\\ \sigma_{01}/\sigma_0& \sqrt{\sigma_1^2-\sigma_{01}^2/\sigma_0^2} \end{pmatrix}, \qquad \begin{pmatrix}b_{0j}\\b_{1j}\end{pmatrix} = L\begin{pmatrix}u_{0j}\\u_{1j}\end{pmatrix}.

The diagonal entries are real and positive when $\Sigma$ is positive definite. In particular, a nonzero covariance makes $b_{1j}$ depend on both standard normal coordinates, even though $u_{0j}$ and $u_{1j}$ are independent. This is the random-intercept-and-slope application of Cholesky whitening.

For a one-dimensional Gauss–Hermite rule with nodes $t_k$ and weights $w_k$ satisfying $\int e^{-t^2}f(t)\,dt\approx\sum_k w_k f(t_k)$, a standard-normal expectation uses nodes $u_k=\sqrt{2}\,t_k$ and normalized weights $w_k/\sqrt{\pi}$. The two-dimensional tensor rule uses every pair of nodes, with product weights. Equivalently, transform each pair to $b=L u$, evaluate the conditional likelihood there, and average using those normalized product weights.

Worked Example: Turning Correlated Nodes into Random Effects

Suppose a hypothetical model has $\sigma_0=0.60$, $\sigma_1=0.25$, and correlation $\rho=-0.40$. Thus $\sigma_{01}=\rho\sigma_0\sigma_1=-0.0600$, and $$ \Sigma=\begin{pmatrix}0.3600&-0.0600\\-0.0600&0.0625\end{pmatrix}. $$ The Cholesky factor is $$ L=\begin{pmatrix}0.6000&0\\-0.1000&\sqrt{0.0525}\end{pmatrix} \approx\begin{pmatrix}0.6000&0\\-0.1000&0.2291\end{pmatrix}. $$ Indeed, the lower-right entry is $\sqrt{0.0625-(-0.0600)^2/0.3600}=\sqrt{0.0525}$, and multiplying $LL^\mathsf{T}$ gives the stated covariance matrix.

Use the three-point standard-normal rule, whose nodes are $0,\pm\sqrt{3}$ and whose normalized weights are $2/3,1/6,1/6$, respectively. At the node $u=(\sqrt{3},0)^\mathsf{T}$, the corresponding effects are $$ b_0=0.6000\sqrt{3}=1.0392,\qquad b_1=-0.1000\sqrt{3}=-0.1732. $$ For an observation with random-slope covariate $z=2$, their contribution to the linear predictor is $b_0+2b_1=1.0392-0.3464=0.6928$.

If, hypothetically, the fixed part of this observation’s linear predictor is $-1.000$, then $\eta=-1.000+0.6928=-0.3072$. Under a Bernoulli-logit model, the conditional probability of response 1 at this node is $1/(1+\exp(0.3072))\approx0.4238$. The node’s weight in the two-dimensional rule is $(1/6)(2/3)=1/9$, so its contribution to the marginal probability for a response of 1 is $(1/9)(0.4238)\approx0.0471$. This is one term in the quadrature sum, not the full marginal probability. The negative covariance is reflected in the negative slope effect at this node.

Adaptive scaling around the cluster mode

Nonadaptive quadrature places nodes around zero in the whitened coordinates. Adaptive quadrature instead places and scales the nodes using the shape of the cluster’s integrand. Define $h_j(u)=\log L_j(Lu)$ and $$ g_j(u)=h_j(u)-\tfrac12u^\mathsf{T}u. $$ Apart from the constant normal-density factor, the cluster integral is the integral of $\exp(g_j(u))$. Let $\hat u_j$ be a local maximum of $g_j$, and define $H_j=-\nabla^2g_j(\hat u_j)$. The adaptive construction requires a finite mode and positive-definite $H_j$; it uses a matrix $A_j$ satisfying $A_jA_j^\mathsf{T}=H_j^{-1}$. A Cholesky factor of $H_j^{-1}$ is a convenient choice.

The following change of variables derives both the transformed nodes and the determinant factor. Write $u=\hat u_j+\sqrt{2}A_jt$. The local quadratic approximation is $g_j(u)\approx g_j(\hat u_j)-\tfrac12(u-\hat u_j)^\mathsf{T}H_j(u-\hat u_j)=g_j(\hat u_j)-t^\mathsf{T}t$. The Jacobian determinant is $2^{d/2}|A_j|$ in dimension $d$. Factoring out the Gaussian weight $e^{-t^\mathsf{T}t}$ therefore gives:

I_j=(2\pi)^{-d/2}\int e^{g_j(u)}\,du \ \approx\ e^{g_j(\hat u_j)}|A_j|\pi^{-d/2} \sum_{k_1,\ldots,k_d} \left(\prod_{r=1}^d w_{k_r}\right) \exp\left\{ g_j(u_{k_1,\ldots,k_d})-g_j(\hat u_j) +\sum_{r=1}^d t_{k_r}^2 \right\}, \qquad u_{k_1,\ldots,k_d}=\hat u_j+\sqrt{2}A_jt.

Here $t=(t_{k_1},\ldots,t_{k_d})^\mathsf{T}$, and the $w_k$ are the weights for the $e^{-t^2}$ convention. The correction term $\sum_r t_{k_r}^2$ cancels the Gaussian decay already represented by the quadrature weight. Omitting it, changing the determinant, or mixing a standard-normal rule’s weights with this $e^{-t^2}$ formula changes the approximation.

Worked Example: Checking an Adaptive Mode and Scale

Consider a hypothetical random-intercept-and-slope logistic cluster with $\Sigma=I$, fixed linear predictor zero, and four observations: two at $z=-1$ and two at $z=1$, with one response 0 and one response 1 at each covariate value. In whitened coordinates, the random part of the predictor is $u_0+zu_1$. At $u=0$, all four response probabilities are $1/2$.

For a Bernoulli-logit observation, the log-likelihood score with respect to $u$ is $(y-p)(1,z)^\mathsf{T}$. At each value of $z$, the paired responses 0 and 1 have opposite score contributions, so the total score is zero. The standard-normal prior also has score zero at $u=0$. Thus $\hat u=0$ satisfies the stationarity equation.

The negative Hessian of the log posterior kernel is $I+\sum_i p_i(1-p_i)a_i a_i^\mathsf{T}$, where $a_i=(1,z_i)^\mathsf{T}$. Since $p_i(1-p_i)=1/4$, and the two observations at each of $z=-1$ and $z=1$ give $\sum_i a_i a_i^\mathsf{T}=4I$, it follows that $H=I+\tfrac14(4I)=2I$. Therefore $A=H^{-1/2}=I/\sqrt{2}$ is consistent with the mode, and $AA^\mathsf{T}=I/2=H^{-1}$; also $|A|=1/2$.

For a three-point rule, the transformed coordinate is $u=\sqrt{2}At=t$. At the tensor-rule center $t=(0,0)^\mathsf{T}$, $g(\hat u)=-4\log 2$, and the exponential correction is 1. The center term in the adaptive approximation is consequently $$ e^{-4\log 2}\,(1/2)\pi^{-1}\,(\sqrt{\pi}\,2/3)^2 =\frac{1}{16}\cdot\frac12\cdot\frac{4}{9} =\frac{1}{72}\approx0.01389. $$

As a check away from the center, take $t=(\sqrt{3/2},0)^\mathsf{T}$. Then $u=t$, and $p=\operatorname{logit}^{-1}(\sqrt{3/2})\approx0.7729$ for the $z=1$ observations. Each covariate level contributes a response pair with likelihood $p(1-p)$, so $L(u)=[p(1-p)]^2$. The correction exponent is $$ g(u)-g(0)+t^\mathsf{T}t =2\log\left(\frac{p(1-p)}{1/4}\right)-\frac{3}{4}+\frac{3}{2} \approx0.04266. $$ Thus the exponential correction is about $1.0436$. The mode, curvature, and evaluated likelihood are all computed from the same specified cluster.

Practical checks, alternatives, and limitations

For a correlated random intercept and slope, there are two transformations to keep distinct. First, $b=Lu$ maps independent normal coordinates to effects with the model covariance. In adaptive quadrature, $u=\hat u+\sqrt{2}At$ then maps quadrature nodes to coordinates centered and scaled for the cluster integrand. The actual effects at an adaptive node are $b=L(\hat u+\sqrt{2}At)$. Applying $L$ twice, or using $A$ in place of $L$, incorrectly changes the covariance or local scale.

Useful numerical checks include verifying $LL^\mathsf{T}=\Sigma$, checking the score $\nabla g(\hat u)$ is near zero, confirming $H$ is positive definite, and checking $AA^\mathsf{T}$ against $H^{-1}$. Also verify that increasing quadrature order stabilizes the cluster integral and that the computation uses the same nodes, weights, and scaling convention throughout. A mode that is not stationary or a scale matrix inconsistent with the Hessian makes the adaptive rule internally incoherent, even if the resulting sum is finite.

Tensor-product quadrature is direct and accurate when the dimension and node count are manageable, but with $K$ nodes per coordinate it requires $K^2$ evaluations for this pair of effects, and much more as the random-effects dimension grows. Sparse grids reduce the number of evaluations in some moderate-dimensional settings; Laplace approximation is cheaper but may be less accurate when the integrand is skewed or poorly represented by a local Gaussian approximation. These are alternatives, not interchangeable formulas: each has a different approximation error and computational cost.

Adaptive quadrature remains a finite-node approximation. It can be sensitive to low quadrature order, weak identification of the slope variance, near-singular covariance, and integrands with substantial skewness or multiple important modes. A positive-definite local Hessian supports the local scaling, but does not prove that a single mode captures the whole integral. If $\Sigma$ is singular, the ordinary full-rank Cholesky transformation is unavailable; the effective lower-dimensional support must be handled rather than forcing an invalid inverse or square root.

Common mistake: “The intercept and slope are separate effects, so I can use independent quadrature nodes for them and ignore their covariance.” This is wrong when $\sigma_{01}\ne0$: independent standard-normal coordinates are appropriate only before transforming by $L$. The corrected calculation evaluates the likelihood at $b=L u$, preserving the covariance at every node. In adaptive quadrature, the covariance transform $L$ and local scale transform $A$ are both needed, and the latter must satisfy $AA^\mathsf{T}=H^{-1}$.

Check Your Understanding

Answer each question, then compare with the solutions.

If $\sigma_{01}=0$, what simplifies in the Cholesky map from $u$ to $b$?
Why does the adaptive quadrature formula include $+\sum_r t_{k_r}^2$ in its exponential correction?
For $H=4I$ in two dimensions, give a valid symmetric scale matrix $A$ and its determinant.
What must be checked before using a proposed adaptive mode and scale matrix?

Solutions. (1) The Cholesky factor is diagonal, $L=\operatorname{diag}(\sigma_0,\sigma_1)$, so each random effect is scaled from its own independent standard-normal coordinate. (2) The quadrature weights already contain $e^{-t^\mathsf{T}t}$; adding $t^\mathsf{T}t$ in the correction accounts for that Gaussian weight when approximating the remaining integrand. (3) $A=I/2$, since $AA^\mathsf{T}=I/4=H^{-1}$, and $|A|=1/4$. (4) Check that the mode has near-zero score, that the negative Hessian there is positive definite, and that $AA^\mathsf{T}=H^{-1}$; also keep the covariance and quadrature transformations distinct.

Quadrature with Random-Intercept and Random-Slope Covariance

What You'll Learn