Integrating Correlated Random Effects

Why Change Coordinates Before Integrating?

In a shared random-effects joint model, the longitudinal measurements and event-time outcome for a subject are conditionally linked through latent subject-specific effects. The subject’s marginal likelihood is obtained by integrating its conditional contribution over the distribution of those effects. When a random intercept and random slope are correlated, their joint density contains that covariance structure. A change of variables can instead express the effects as linear combinations of independent standard-normal variables. This gives an equivalent integral with a simpler, fixed reference density.

The change is useful for writing and checking the likelihood integral, and for later numerical calculations. It does not make the conditional likelihood simpler by itself: wherever the original model uses a random intercept or slope, the transformed values must be substituted into that model.

Model and Transformation

For subject $i$, let $b_i=(b_{0i},b_{1i})^\mathsf{T}$ denote the random intercept and random slope. Assume the effects are independent across subjects and have a bivariate normal distribution with mean zero and positive-definite covariance matrix $D$. Write that matrix in terms of standard deviations $\sigma_0>0,\sigma_1>0$ and correlation $-1<\rho<1$:

b_i\sim N_2(0,D),\qquad D= \begin{pmatrix} \sigma_0^2 & \rho\sigma_0\sigma_1\\ \rho\sigma_0\sigma_1 & \sigma_1^2 \end{pmatrix}.

Here $\sigma_0$ is measured in the units of the longitudinal outcome, while $\sigma_1$ is measured in outcome units per unit of time if the slope is per unit of time. The correlation $\rho$ is dimensionless. Positive definiteness requires $\sigma_0,\sigma_1>0$ and $|\rho|<1$.

Definition: Let $z_i=(z_{0i},z_{1i})^\mathsf{T}$, where $z_{0i}$ and $z_{1i}$ are independent $N(0,1)$ variables. The standardized representation of the random effects is $$ b_{0i}=\sigma_0z_{0i},\qquad b_{1i}=\sigma_1\left\{\rho z_{0i}+\sqrt{1-\rho^2}\,z_{1i}\right\}. $$ Equivalently, $b_i=Lz_i$, where $$ L= \begin{pmatrix} \sigma_0 & 0\\ \rho\sigma_1 & \sigma_1\sqrt{1-\rho^2} \end{pmatrix}. $$

The same $z_{0i}$ contributes to both effects, inducing their correlation; the independent $z_{1i}$ supplies the remaining slope variation. In particular, $LL^\mathsf{T}=D$, so the transformed vector has the required covariance. The construction assumes a bivariate normal random-effects distribution; it is not merely a way to reproduce a correlation while retaining an arbitrary distribution.

How the Marginal-Likelihood Integral Changes

Let $y_i$ denote subject $i$’s longitudinal data, $T_i$ the observed event or censoring time, and $\delta_i$ the event indicator. Write $L_i(b_i;\theta)$ for the conditional subject contribution, with parameters $\theta$. It includes the longitudinal and event-time terms appropriate to the model and observation scheme. The marginal contribution integrates this conditional contribution against the random-effects density:

L_i^{\mathrm{marg}}(\theta) = \int_{\mathbb R^2} L_i(b_i;\theta)\,\phi_2(b_i;0,D)\,db_i.

Here $\phi_2(\,\cdot\,;0,D)$ is the bivariate normal density. Under the change $b_i=Lz_i$, the absolute Jacobian determinant is $|L|=\sigma_0\sigma_1\sqrt{1-\rho^2}$. Since $|D|=|L|^2$, the density and the Jacobian cancel their covariance-dependent normalizing factors. The result is an integral under the standard bivariate normal density:

\begin{aligned} L_i^{\mathrm{marg}}(\theta) &=\int_{\mathbb R^2}L_i(Lz_i;\theta)\, \phi_2(Lz_i;0,D)\,|L|\,dz_i\\ &=\int_{\mathbb R^2}L_i(Lz_i;\theta)\, \frac{1}{2\pi}\exp\left\{-\frac{z_{0i}^2+z_{1i}^2}{2}\right\}\,dz_i. \end{aligned}

The cancellation follows from $D=LL^\mathsf{T}$: the quadratic form $b_i^\mathsf{T}D^{-1}b_i$ becomes $z_i^\mathsf{T}z_i$, and the normal density’s factor $1/\sqrt{|D|}$ is $1/|L|$. Multiplication by the change-of-variables factor $|L|$ removes it. The conditional contribution remains $L_i(Lz_i;\theta)$, so the covariance parameters still affect the integral through the transformed random effects.

Worked Example: Transforming a Pair of Standard-Normal Values

Suppose the random intercept has standard deviation $\sigma_0=2$ outcome units, the random slope has standard deviation $\sigma_1=0.5$ outcome units per year, and their correlation is $\rho=0.6$. For $z_{0i}=1$ and $z_{1i}=-1$, calculate the corresponding random effects.

The transformation gives $b_{0i}=2(1)=2$ outcome units. For the slope,

b_{1i}=0.5\{0.6(1)+\sqrt{1-0.6^2}(-1)\} =0.5(0.6-0.8)=-0.1

outcome units per year. Thus this subject has an intercept deviation of 2 units and a slope deviation of $-0.1$ units per year. At time $t=4$ years, their combined contribution to the linear predictor is $b_{0i}+4b_{1i}=2+4(-0.1)=1.6$ outcome units. This calculation checks the substitution into a trajectory; it does not by itself determine the subject’s likelihood.

Worked Example: Checking the Covariance and Density Transformation

Use the same $\sigma_0=2,\sigma_1=0.5,\rho=0.6$. First calculate the covariance matrix. Then, for $b_i=(1,0.2)^\mathsf{T}$, calculate its standardized coordinates and check the density-times-Jacobian identity.

The covariance matrix is

D= \begin{pmatrix} 4 & 0.6\\ 0.6 & 0.25 \end{pmatrix}, \qquad |D|=4(0.25)-0.6^2=0.64.

The off-diagonal entry is $\rho\sigma_0\sigma_1=0.6$, and the determinant is positive, as required. The inverse transformation gives

z_{0i}=\frac{b_{0i}}{\sigma_0}=\frac{1}{2}=0.5,\qquad z_{1i}=\frac{b_{1i}/\sigma_1-\rho z_{0i}}{\sqrt{1-\rho^2}} =\frac{0.4-0.3}{0.8}=0.125.

The Jacobian magnitude is $|L|=2(0.5)(0.8)=0.8$. At this point, the correlated density is

\phi_2(b_i;0,D) =\frac{1}{2\pi(0.8)} \exp\left\{-\frac{0.5^2+0.125^2}{2}\right\} \approx 0.1742.

Multiplying by the Jacobian gives $0.1742(0.8)\approx0.1394$. Directly, the standard bivariate normal density at $(0.5,0.125)$ is $\frac{1}{2\pi}\exp\{-(0.5^2+0.125^2)/2\}\approx0.1394$. The agreement verifies the density transformation, to rounding.

Worked Example: An Exact Integral Check

Consider a simple hypothetical conditional contribution $L_i(b_i)=\exp(-b_{0i}^2/2)$, with $\sigma_0=2$. This is a bounded function of the intercept alone, so it can be used to check the change of variables without specifying a full joint-model likelihood. The slope distribution and its correlation with the intercept do not affect this particular integral because the contribution does not depend on the slope.

In the original coordinates, $b_{0i}\sim N(0,4)$. Integrating the normal density multiplied by this contribution gives

E\{\exp(-b_{0i}^2/2)\} =\int_{-\infty}^{\infty}\exp(-b^2/2) \frac{1}{\sqrt{2\pi(4)}}\exp\left(-\frac{b^2}{2(4)}\right)\,db =\frac{1}{\sqrt{5}} \approx0.4472.

In standardized coordinates $b_{0i}=2z_{0i}$, the same expectation is

\int_{\mathbb R^2}\exp\{-(2z_{0i})^2/2\} \frac{1}{2\pi}\exp\left\{-\frac{z_{0i}^2+z_{1i}^2}{2}\right\} \,dz_{0i}\,dz_{1i} =\frac{1}{\sqrt{5}}\cdot1 \approx0.4472.

The first factor is the one-dimensional normal integral for $z_{0i}$; the second is the integral of the standard-normal density for $z_{1i}$, equal to 1. Both coordinate systems therefore give the same marginal contribution.

Assumptions, Checks, and Failure Modes

The transformation is valid when the random intercept and slope follow the stated bivariate normal distribution and its covariance matrix is positive definite. The subject contribution must be integrated over that distribution, and the substitution must be applied consistently wherever the random effects enter the longitudinal model or event hazard. Conditional on the random effects and included covariates, the subject contribution must also be the one specified by the joint model; a coordinate change cannot repair an incorrectly specified conditional model.

Useful implementation checks include verifying $LL^\mathsf{T}=D$, confirming the inverse transformation recovers the original $b_i$, and checking that the transformed density integrates to one. Simulated standard-normal pairs should produce sample variances near 1 and sample correlation near 0; after transformation, the corresponding random effects should have variances near $\sigma_0^2,\sigma_1^2$ and correlation near $\rho$, allowing for simulation error. In an actual likelihood calculation, compare the original and transformed integrands at matching values using the Jacobian identity.

Reversing the correlation term: using $b_{1i}=\sigma_1\{\rho z_{1i}+\sqrt{1-\rho^2}z_{0i}\}$ with $b_{0i}=\sigma_0z_{0i}$ generally produces the wrong covariance. Check the cross-product: the specified construction gives $\operatorname{Cov}(b_{0i},b_{1i})=\sigma_0\sigma_1\rho$.
Omitting the Jacobian in a density calculation: the density in $b$-coordinates is not simply the standard-normal density evaluated at $z$. The Jacobian is essential for equality of the two integrals, even though it cancels the covariance factor in the transformed measure.
Using $\rho=\pm1$ as an ordinary case: then $|D|=0$, $L$ is singular, and this two-dimensional density transformation is not valid. The model has a degenerate random-effects distribution and needs a lower-dimensional treatment.
Dropping covariance dependence: after transformation, the reference density no longer contains $\sigma_0,\sigma_1,\rho$, but $L_i(Lz_i;\theta)$ does. Treating the entire integral as independent of the covariance parameters gives an incorrect likelihood.

What the Transformation Does—and Does Not—Do

One can integrate directly in correlated $b$-coordinates using $\phi_2(b;0,D)$, or integrate in $z$-coordinates using the independent standard-normal density and substituting $b=Lz$ into the conditional contribution. These are mathematically equivalent representations. The transformed representation separates the reference distribution from its covariance parameters, which is convenient when setting up a numerical integral. The direct representation may be easier to interpret in terms of subject-specific intercepts and slopes.

Neither representation removes the need to evaluate the marginal integral, and neither guarantees accurate numerical results. Strong curvature in the conditional contribution, weak information about a variance component, or a correlation close to a boundary can still make computation difficult. The transformation also does not justify assuming independent random intercepts and slopes: independence applies to $z_{0i}$ and $z_{1i}$, not generally to $b_{0i}$ and $b_{1i}$. The integral representation is the necessary setup; the choice and accuracy of a numerical integration method are separate issues.

Key takeaway: For a positive-definite bivariate normal intercept–slope distribution, write $b_i=Lz_i$ with independent standard-normal coordinates. The change of variables converts the marginal contribution to an integral against a standard bivariate normal density, while preserving all covariance effects through the transformed conditional contribution $L_i(Lz_i;\theta)$.

Check Your Understanding

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

If $\rho=0$, what does the transformation reduce to, and what does that imply about the random intercept and slope?
For $\sigma_0=3,\sigma_1=2,\rho=-0.5$, calculate $b_{0i}$ and $b_{1i}$ when $z_{0i}=1,z_{1i}=1$.
Why does the transformed marginal-likelihood integral still depend on $\sigma_0,\sigma_1,\rho$?
What fails in the stated two-dimensional change of variables if $\rho=1$?
State one numerical check that can detect an incorrect implementation of the transformation.

Solutions. (1) It becomes $b_{0i}=\sigma_0z_{0i}$ and $b_{1i}=\sigma_1z_{1i}$; the effects are independent under the assumed bivariate normal model because their covariance is zero. (2) $b_{0i}=3$, and $b_{1i}=2\{-0.5+\sqrt{0.75}\}\approx2(-0.5+0.8660)=0.7321$, in the respective intercept and slope units. (3) The parameters enter through $b_i=Lz_i$ inside the conditional contribution, even though the density of $z_i$ is fixed. (4) The covariance matrix is singular, $|D|=0$, and the transformation matrix has zero determinant, so there is no ordinary two-dimensional density or invertible change of variables. (5) Verify numerically that $LL^\mathsf{T}=D$, or transform simulated standard-normal pairs and check that their empirical covariance approaches $D$.