Nonlinear mixed-effects models in Julia

Linear mixed-effects models

A linear mixed-effects model is characterized by the distribution of two vector-valued random variables: the $n$-dimensional response, $Y$, and the $q$-dimensional random effects vector, $B$. The unconditional distribution of $B$ is multivariate normal $$ B\sim N(0,\Sigma_\theta) $$ as is the conditional distribution of $Y$ given $B=b$ $$ (Y|B=b)\sim N(X\beta+Zb, \sigma^2 I) $$

In the MixedModels package for Julia we represent the covariance matrix in the unconditional distribution of $B$ as $$ \Sigma_\theta=\sigma^2\Lambda_\theta\Lambda_\theta' $$ where $\Lambda_\theta$ is the $q\times q$ relative covariance factor.

For given values of $\theta$ and $\beta$ we solve a penalized linear least squares problem of the form $$ r^2_{\beta,\theta}=\min_u \|(y-X\beta)-Z\Lambda_\theta u\|^2 + \|u\|^2 $$ for which we compute the sparse Cholesky factorization $$ L_\theta L_\theta' = \Lambda_\theta'Z'Z\Lambda_\theta + I $$ where $L_\theta$ is a lower-triangular sparse matrix. Because $L_\theta$ is triangular, we can easily evaluate its determinant as the product of its diagonal elements. By construction these diagonal elements are positive and the log-determinant, $\log(|\Lambda_\theta'Z'Z\Lambda_\theta + I|)=2\log(|L_\theta|)$ is easily evaluated.

The log-likelihood for a linear mixed-effects model, $\ell(\beta,\theta|y)$, can be expressed as $$ -2\ell(\beta,\theta|y) = \log(|\Lambda_\theta'Z'Z\Lambda_\theta + I|)+n\left[1+\log\left(\frac{2\pi r^2_{\beta,\theta}}{n}\right)\right] $$

Nonlinear mixed-effects models

The formulation of a nonlinear mixed-effects model (NLMM) is nearly the same as that of an LMM except that the mean of the conditional distribution $Y|B=b$ is a nonlinear function of the parameters $\beta$ and the random effects, $b$ or, equivalently, the spherical random effects $u$ where $b=\Lambda_\theta u$.

The nonlinear model function, $f$, is a function of a $k$-dimensional model parameter, $\phi$ and covariates. In our formulation we write the vec of the $n\times k$ model parameter matrix, $\Phi$ as a linear predictor $$ vec(\Phi)=X\beta+Z\Lambda_\theta u $$ and evaluate $$ \mu_{Y|U=u}=f(T,\Phi) $$ where $T$ is a matrix of covariate values, time in the case of a population pharmacokinetic model.

The penalized linear least squares (PLS) problem solved for each evaluation of the objective in the case of a linear mixed-effects model is replaced by a penalized nonlinear least squares (PNLS) problem $$ r^2_{\beta,\theta}=\min_u \|y-\mu_{Y|U=u}\|^2+\|u\|^2 $$ in an NLMM.
The optimizer, $\tilde{u}=\arg\min_u r^2_{\beta,\theta}$, is the mode of the distribution of $U$ given $Y=y$ and the current values of the parameters, $\beta$ and $\theta$.

In general there will not be a closed-form expression for the log-likelihood. However, the Laplace approximation to the log-likelihood is exactly the same expression as the log-likelihood for an LMM. The Laplace approximation is equivalent to a one-point adaptive Gauss-Hermite approximation to the integral that defines the log-likelihood. This approximation is adaptive in the sense that the Gauss-Hermite approximation is taken at the mode of the conditional distribution, not at the mean of the unconditional distribution.

A simple example - first-order kinetics, 1 compartment, single bolus unit dose

One of the simplest pharmacokinetics models is that for a one compartment model with a single bolus dose. The predicted concentration at time $t$ is given by $$ c(t)=V\exp(-Kt) $$ where V is the volume of distribution and K is the elimination rate constant. We represent an instance of such a model with the user-defined type logsd1.

Fitting the model

Simulated data from a one compartment model with a single bolus dose are provided in a compressed csv file in the data directory of the NLreg package. We read these data as a DataFrame object. The head function shows us the first few rows.

We create a nonlinear regression model from these data as

The SimpleNLMM type represents a nonlinear mixed-effects model with a single grouping factor for the random effects (ID in this case), and random-effects for each of the nonlinear model parameters. The third and fourth arguments to the external constructor are the initial values for the template $\Lambda$ matrix and for the fixed-effects parameters $\beta$.

The prss! function evaluates the penalized residual sum of squares, $r^2_{\beta,\theta}(u)$ at the initial parameter estimates. The fixed-effects parameter, $\beta$, is the parameter vector from the nonlinear least squares fit, the template for the relative covariance factor, $\Lambda$, is initialized to the identity and the random-effects, $u$, are initialized to zero.

The individual values of the nonlinear model parameters, $\phi_i,i=1,\dots,200$ are stored as the phi element of the model structure.