Simulating Parameters With Lognormal Random Effects

Why Lognormal Random Effects Need Careful Interpretation

In an NLME model, an individual parameter such as clearance or volume may need to be positive. A normal random effect can enforce positivity by entering the parameter through an exponential transformation. The resulting parameter is lognormally distributed: its distribution is asymmetric, and its typical value is not generally its arithmetic mean. These distinctions matter when generating individual parameters from a fitted model and summarising the resulting simulated population.

This tutorial focuses on the distribution induced by that transformation. It treats the population parameters as fixed at their fitted values and simulates variability between individuals. It does not include uncertainty in the fitted population parameters or residual error in observations.

The model assumes that \(\theta_{\mathrm{pop}}\) is positive, that the random effect is normally distributed with mean zero and variance \(\omega^2\), and that the exponential relationship is the specified model for individual variability. With \(\omega>0\), every simulated \(\theta_i\) is positive, and the distribution is generally right-skewed. If \(\omega=0\), there is no between-individual variability and every individual has \(\theta_i=\theta_{\mathrm{pop}}\).

Median, Mean, and Variability

Because \(\eta_i\) has median zero, the median individual parameter is \(\theta_{\mathrm{pop}}\). In this model, the population typical value is therefore the median of the individual-parameter distribution, not its arithmetic mean. The mean is larger whenever \(\omega>0\), a consequence of the curvature of the exponential function.

The coefficient of variation (CV) is dimensionless: it is the standard deviation divided by the arithmetic mean. If a CV is reported as a percentage, multiply the displayed ratio by 100. Notice that the CV depends on \(\omega\), not on the size or units of \(\theta_{\mathrm{pop}}\). The mean and standard deviation, in contrast, have the same units as the individual parameter.

Why These Formulas Hold

For a normal random variable \(Z\sim N(0,\omega^2)\), its moment-generating identity is \(E\{\exp(tZ)\}=\exp(t^2\omega^2/2)\). Applying this identity to \(\theta_i=\theta_{\mathrm{pop}}\exp(\eta_i)\) gives the mean directly. Applying it to \(\theta_i^2\) gives the second moment; subtracting the squared mean gives the variance.

Taking the square root of the variance and dividing by the mean cancels \(\theta_{\mathrm{pop}}\) and leaves \(\sqrt{\exp(\omega^2)-1}\). For the median, the exponential is strictly increasing, so the median of \(\exp(\eta_i)\) is \(\exp(0)=1\); multiplying by \(\theta_{\mathrm{pop}}\) gives the stated result. More generally, the \(q\)-quantile, for \(0<q<1\), is

where \(\Phi^{-1}(q)\) is the \(q\)-quantile of the standard normal distribution. Thus a central 95% interval for individual parameters under the model is \(\theta_{\mathrm{pop}}\exp(-1.96\omega)\) to \(\theta_{\mathrm{pop}}\exp(1.96\omega)\). It describes the central 95% of the modelled individual-parameter distribution, not uncertainty about the fitted population value.

Simulating and Checking Parameters

For a single parameter, simulation at fixed fitted population values has a simple sequence: draw \(\eta_i\) from \(N(0,\omega^2)\), then calculate \(\theta_i=\theta_{\mathrm{pop}}\exp(\eta_i)\). For multiple parameters, the random effects may be drawn jointly from a multivariate normal distribution with the fitted covariance structure, then transformed component by component. The transformation preserves positivity, but the correlations on the parameter scale need not equal the correlations on the random-effect scale.

For example, if two log-scale random effects have standard deviations \(\omega_1,\omega_2\) and correlation \(\rho\), then their transformed parameters have correlation

This follows by calculating the covariance of the exponentials using the joint normal moment-generating function and dividing by their standard deviations. It is a reminder to distinguish the random-effect correlation from the correlation between simulated parameter values.

Useful checks on a simulation include verifying that every transformed parameter is positive, comparing the simulated median with \(\theta_{\mathrm{pop}}\), and comparing the simulated mean and CV with \(\theta_{\mathrm{pop}}\exp(\omega^2/2)\) and \(\sqrt{\exp(\omega^2)-1}\). A finite simulation will not match these theoretical values exactly. Larger samples generally reduce Monte Carlo fluctuation, but a discrepancy that is large relative to sampling variation may indicate an incorrect transformation, an incorrect value of \(\omega\), or a coding error.

Choosing and Interpreting the Distribution

The lognormal model is useful when individual parameters must be positive and variability is naturally described by multiplicative departures from a typical value. A random effect of \(0.1\) multiplies the typical value by \(\exp(0.1)\), while a random effect of \(-0.1\) multiplies it by \(\exp(-0.1)\). These are reciprocal factors, not equal additive changes. By contrast, an additive normal model \(\theta_i=\theta_{\mathrm{pop}}+\eta_i\) describes additive departures and can assign positive probability to negative parameter values. That may be unsuitable for a parameter that must remain positive.

The lognormal model is not automatically appropriate for every positive parameter. Its right-skewness can produce a long upper tail, and the model can imply very large values when \(\omega\) is large. The fitted distribution and the scientific plausibility of its implied range should be examined rather than relying only on the fact that exponentiation guarantees positivity.