Drawing Correlated Fixed Effects for Simulation

Why Draw Fixed Effects Jointly?

An NLME fit gives estimates of population fixed-effect parameters, along with an estimated covariance matrix that describes uncertainty in those estimates. When simulation is intended to reflect that parameter uncertainty, drawing each estimate separately is generally not enough: estimates from the same fit may move together or in opposite directions. A simulation that omits this dependence can misrepresent uncertainty in predicted quantities, even if each parameter's marginal standard error is correct.

The task here is to draw a vector of fixed-effect parameters jointly, using its estimated covariance matrix. These are draws representing uncertainty about the fitted fixed effects under an approximation. They are not draws of subject-level random effects, and the fixed-effect covariance matrix is not the random-effect covariance matrix.

Definition: Let $\hat{\boldsymbol{\theta}}$ be a vector of estimated NLME fixed effects, in a specified parameter order, and let $\widehat{\mathbf V}_{\theta}$ be its estimated covariance matrix. A correlated fixed-effect draw is generated from the approximate distribution $$ \boldsymbol{\theta}^{*}\sim N\!\left(\hat{\boldsymbol{\theta}},\widehat{\mathbf V}_{\theta}\right). $$ The draw has expected value $\hat{\boldsymbol{\theta}}$ and covariance $\widehat{\mathbf V}_{\theta}$ under this simulation approximation.

Assumptions and the Cholesky Construction

The multivariate normal approximation is most defensible when the fixed-effect estimator is approximately normal near the fitted value. In practice, this is more plausible when the model is identifiable, the estimate is not near a parameter boundary, and the likelihood or estimating procedure has behaved regularly enough for its covariance estimate to be informative. The approximation is conditional on the fitted model and the estimated covariance matrix: it does not account for uncertainty in that covariance matrix itself.

For $p$ fixed effects, $\hat{\boldsymbol{\theta}}$ has $p$ entries and $\widehat{\mathbf V}_{\theta}$ is a $p\times p$ matrix. It should be symmetric and positive definite for the ordinary Cholesky construction. Write its Cholesky factor as a lower-triangular matrix $\mathbf L$, satisfying $\mathbf L\mathbf L^{T}=\widehat{\mathbf V}_{\theta}$. Generate a vector $\mathbf z$ of $p$ independent standard normal values, then transform it as follows:

\boldsymbol{\theta}^{*} = \hat{\boldsymbol{\theta}}+\mathbf L\mathbf z, \qquad \mathbf z\sim N(\mathbf 0,\mathbf I).

This construction has the desired mean because $E(\mathbf z)=\mathbf 0$, so $E(\boldsymbol{\theta}^{*})=\hat{\boldsymbol{\theta}}$. Its covariance follows from the covariance rule for a linear transformation:

\operatorname{Cov}(\boldsymbol{\theta}^{*}) = \operatorname{Cov}(\mathbf L\mathbf z) = \mathbf L\operatorname{Cov}(\mathbf z)\mathbf L^{T} = \mathbf L\mathbf L^{T} = \widehat{\mathbf V}_{\theta}.

Thus the off-diagonal entries of the fitted covariance matrix are reproduced by the joint draws. In particular, a positive covariance tends to produce draws in which the corresponding estimates are above or below their fitted values together; a negative covariance tends to produce opposite-direction deviations. The degree of linear dependence is summarized by the correlation, obtained by dividing the covariance by the product of the two standard errors.

Worked Example: A Two-Parameter Draw With Positive Correlation

Suppose a hypothetical NLME fit estimates clearance $CL=5.0$ L/h and volume $V=40.0$ L. In the order $(CL,V)$, suppose the estimated covariance matrix is

\widehat{\mathbf V}_{\theta} = \begin{pmatrix} 0.16 & 0.96\\ 0.96 & 16.00 \end{pmatrix}.

The standard errors are $\sqrt{0.16}=0.40$ L/h for clearance and $\sqrt{16.00}=4.00$ L for volume. Their estimated correlation is $0.96/(0.40\times4.00)=0.60$. A Cholesky factor is

\mathbf L= \begin{pmatrix} 0.40 & 0\\ 2.40 & 3.20 \end{pmatrix}, \qquad \mathbf L\mathbf L^{T} = \begin{pmatrix} 0.16 & 0.96\\ 0.96 & 16.00 \end{pmatrix}.

For an illustrative standard-normal vector $\mathbf z=(1.00,-0.50)^{T}$, the simulated deviation is

\mathbf L\mathbf z = \begin{pmatrix} 0.40(1.00)\\ 2.40(1.00)+3.20(-0.50) \end{pmatrix} = \begin{pmatrix} 0.40\\ 0.80 \end{pmatrix}.

Adding this to the estimate gives $CL^{*}=5.40$ L/h and $V^{*}=40.80$ L. This is one hypothetical draw, not a revised estimate. The two positive deviations are consistent with the positive fitted correlation; other draws can have deviations in opposite directions.

Why Dependence Matters for Predictions

A predicted quantity can depend on several fixed effects at once. If its local dependence on those effects is approximately linear, the covariance matrix determines how much uncertainty is associated with their combined contribution. For a linear combination $a^{T}\boldsymbol{\theta}$, its variance under the joint approximation is $a^{T}\widehat{\mathbf V}_{\theta}a$. Dropping off-diagonal covariances changes this variance unless the relevant covariances are zero or happen to cancel in that particular combination.

Worked Example: The Effect of Covariance on a Linear Prediction

Suppose two fixed effects, $\beta_0$ and $\beta_1$, both have units of a response measure, and a model prediction at covariate value $x=2$ is $\beta_0+2\beta_1$. For this hypothetical fit, let the standard errors be $0.40$ and $0.50$, respectively, and let their covariance be $0.12$. The variance of the prediction from joint draws is

\begin{aligned} \operatorname{Var}(\beta_0+2\beta_1) &= \operatorname{Var}(\beta_0)+4\operatorname{Var}(\beta_1) +4\operatorname{Cov}(\beta_0,\beta_1)\\ &= 0.40^2+4(0.50^2)+4(0.12)\\ &= 0.16+1.00+0.48=1.64. \end{aligned}

The prediction standard deviation is $\sqrt{1.64}=1.281$ response units. If the two effects were incorrectly drawn independently, the covariance term would be omitted: the variance would be $0.16+1.00=1.16$, and the standard deviation would be $\sqrt{1.16}=1.077$ response units. In this example, independent draws understate the uncertainty in this prediction because the covariance is positive and the prediction uses positive weights for both effects. A different covariance sign or prediction contrast could instead make independent draws overstate uncertainty.

Generating and Checking a Set of Draws

For many simulation replicates, repeat the same operation with a fresh independent standard-normal vector each time, while keeping the parameter order and covariance matrix fixed. Each simulated vector must use the same ordering for $\hat{\boldsymbol{\theta}}$, the rows and columns of $\widehat{\mathbf V}_{\theta}$, and the entries of $\mathbf z$. A mismatch can produce numerically plausible values with the wrong dependence structure.

Specify the parameter order.
Record the fixed-effect names and units in one order, and verify that the covariance matrix uses that same order in both dimensions.

Check and factor the covariance matrix.
Confirm symmetry and positive definiteness, then calculate a Cholesky factor $\mathbf L$ such that $\mathbf L\mathbf L^{T}=\widehat{\mathbf V}_{\theta}$.

Generate joint draws.
For each replicate, generate $p$ independent standard normal values and compute $\hat{\boldsymbol{\theta}}+\mathbf L\mathbf z$.

Check the simulated sample.
Across many draws, compare empirical means, variances, and covariances with the fitted values. Investigate discrepancies beyond ordinary Monte Carlo variation.

The empirical mean and covariance from a finite set of draws will not equal their targets exactly. They should approach $\hat{\boldsymbol{\theta}}$ and $\widehat{\mathbf V}_{\theta}$ as the number of draws increases, but random differences remain. Also check whether draws violate model constraints, such as positivity. Such violations can signal that the normal approximation is unsuitable on the original parameter scale; they should not be silently discarded or treated as valid parameter values.

Worked Example: A Three-Parameter Draw and Matrix Check

Consider three hypothetical fixed effects with estimate $\hat{\boldsymbol{\theta}}=(1.20,0.80,2.00)^{T}$, where each component is in its model-defined parameter units. Suppose the factor is

\mathbf L= \begin{pmatrix} 0.20&0&0\\ 0.10&0.30&0\\ -0.05&0.12&0.25 \end{pmatrix}.

Multiplying the factor by its transpose gives the implied covariance matrix:

\widehat{\mathbf V}_{\theta} = \begin{pmatrix} 0.0400&0.0200&-0.0100\\ 0.0200&0.1000&0.0310\\ -0.0100&0.0310&0.0794 \end{pmatrix}.

For $\mathbf z=(0.50,-1.00,2.00)^{T}$, the deviation is \[ \mathbf L\mathbf z = \begin{pmatrix} 0.20(0.50)\\ 0.10(0.50)+0.30(-1.00)\\ -0.05(0.50)+0.12(-1.00)+0.25(2.00) \end{pmatrix} = \begin{pmatrix} 0.10\\ -0.25\\ 0.355 \end{pmatrix}. \] The resulting draw is $(1.30,0.55,2.355)^{T}$. The positive covariance between the second and third effects is $0.0310$, while the covariance between the first and third is $-0.0100$. Those signs are part of the joint uncertainty information being carried into the simulation.

Alternatives, Diagnostics, and Limits

The most direct alternative is to draw each fixed effect independently from a normal distribution using its own estimate and standard error. This is appropriate only when the covariance matrix is diagonal, or when ignoring its off-diagonal entries is a deliberate approximation justified for the purpose at hand. Otherwise it imposes zero correlation and can distort uncertainty in model predictions. A multivariate-normal routine that accepts a mean vector and covariance matrix is another valid implementation; it should produce the same target distribution as the Cholesky construction, up to implementation and Monte Carlo differences.

Cholesky factorization requires a positive-definite covariance matrix. If it fails, first check for transcription errors, inconsistent parameter ordering, asymmetry from numerical rounding, or a poorly estimated covariance matrix. A singular or nearly singular matrix may arise when parameters are weakly identified or redundant. Do not automatically modify the matrix until it factors: any adjustment changes the target distribution and needs a statistical and computational justification.

Useful diagnostics include checking that the factor reconstructs the covariance matrix, that simulated marginal standard deviations are close to the fitted standard errors, and that simulated correlations are close to the fitted correlations. Also inspect draw ranges for implausible or impermissible values. These checks can reveal coding errors or an unsuitable approximation, but passing them does not establish that the fitted NLME model is adequate.

The method treats the estimated covariance matrix as known and uses a normal approximation for the estimator. It does not capture uncertainty in the covariance estimate, non-normal sampling distributions, or all consequences of parameter boundaries. It is also distinct from drawing random effects for new subjects: the latter represents between-subject variability, whereas these draws represent uncertainty in the estimated fixed effects. The appropriate simulation depends on which source of uncertainty is intended.

Common pitfall: “The standard errors are correct, so independent normal draws are sufficient.” This is wrong when fixed-effect estimates are correlated: matching marginal standard errors does not reproduce the joint distribution. Use the full covariance matrix, and check the simulated cross-covariances as well as the marginal variances.

Check Your Understanding

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

A two-effect covariance matrix has diagonal entries $0.09$ and $0.25$, and off-diagonal entries $-0.06$. What are the standard errors and correlation?
If $\mathbf L\mathbf L^{T}=\widehat{\mathbf V}_{\theta}$, why does adding $\mathbf L\mathbf z$ to the estimate preserve the target covariance?
For a prediction $\beta_0+2\beta_1$, the two variances are $0.16$ and $0.25$, and the covariance is $-0.05$. Calculate the joint variance and compare it with the variance under independent draws.
A Cholesky factorization fails for the supplied covariance matrix. Name two checks to make before changing the matrix.

Solutions. (1) The standard errors are $\sqrt{0.09}=0.30$ and $\sqrt{0.25}=0.50$. The correlation is $-0.06/(0.30\times0.50)=-0.40$. (2) Since $\operatorname{Cov}(\mathbf z)=\mathbf I$, the transformed covariance is $\mathbf L\mathbf I\mathbf L^{T}=\widehat{\mathbf V}_{\theta}$; adding a constant estimate does not change covariance. (3) The joint variance is $0.16+4(0.25)+4(-0.05)=0.96$. Independent draws give $0.16+4(0.25)=1.16$, which is larger here because the negative covariance reduces the variance of this positively weighted combination. (4) Verify that the parameter order matches across the estimate and both covariance-matrix dimensions, and check symmetry, transcription, and positive definiteness. A failed factorization should prompt diagnosis rather than an automatic matrix alteration.