Nonlinear Mixed-Effects Modeling

What Is a Nonlinear Mixed-Effects Model?

A mixed-effects model is a statistical model that incorporates both fixed effects and random effects. Fixed effects are population parameters assumed to be the same each time data is collected, and random effects are random variables associated with each sample (individual) from a population. Mixed-effects models work with small sample sizes and sparse data sets, and are often used to make inferences on features underlying profiles of repeated measurements from a group of individuals from a population of interest.

As with all regression models, their purpose is to describe a response variable as a function of the predictor (independent) variables. Mixed-effects models, however, recognize correlations within sample subgroups, providing a reasonable compromise between ignoring data groups entirely, thereby losing valuable information, and fitting each group separately, which requires significantly more data points.

For instance, consider population pharmacokinetic data that involve the administration of a drug to several individuals and the subsequent observation of drug concentration for each individual, and the objective is to make a broader inference on population-wide parameters while considering individual variations. The nonlinear function often used for such data is an exponential function since many drugs once distributed in a patient are eliminated in an exponential fashion. Thus the measured drug concentration of an individual can be described as:

$y_{i j} = \frac{D_{i}}{V_{}} e^{- k_{i} t_{i j}} + a ε_{i j}$ ,

where y_ij is the jth response of the ith individual, D_i is the dose administered to the ith individual, V is the population mean volume of distribution, a is an error parameter, and $ε_{i j} \sim N (0, 1)$ , representing some measurement error. The elimination rate parameter (k_i) depends on the clearance and volume of the central compartment $k_{i} = \frac{C l_{i}}{V_{}}$ . Both k_i and Cl_i are for the ith patient, meaning they are patient-specific parameters.

To account for variations between individuals, assume that the clearance is a random variable depending on individuals, varying around the population mean. For the ith individual, $C l_{i} = θ_{1} + η_{i}$ , where θ₁ is the fixed effect (population parameter for the clearance) and η_i is the random effect, that is, the deviation of the ith individual from the mean clearance of the population $η_{i} \sim Ν (0, σ_{η}^{2})$ .

If you have any individual-specific covariates such as weight w that linearly relate to the clearance, you can try explaining some of the between-individual differences. For example, if w_i is the weight of the ith individual, then the model becomes $C l_{i} = θ_{1} + θ_{2} * w_{i} + η_{i}$ , where θ₂ is the fixed effect of weight on clearance.

A general nonlinear mixed-effects (NLME) model with constant variance is as follows:

$\begin{array}{l} y_{i j} = f (x_{i j}, p_{i}) + ε_{i j} \\ p_{i} = A_{i} θ + B_{i} η_{i} \\ ε_{i j} \sim N (0, σ^{2}) \\ η_{i} \sim N (0, Ψ) \end{array}$

y_ij	Data vector of individual-specific response values
f	General, real-valued function of p_i and x_ij
x_ij	Data matrix of individual-specific predictor values
p_i	Vector of individual-specific model parameters
θ	Vector of fixed effects, modeling population parameters
η_i	Vector of multivariate normally distributed individual-specific random effects
A_i	Individual-specific design matrix for combining fixed effects
B_i	Individual-specific design matrix for combining random effects
ε_ij	Vector of group-specific errors, assumed to be independent, identically, normally distributed, and independent of η_i
Ψ	Covariance matrix for the random effects
σ²	Error variance, assumed to be constant across observations

In addition to the constant error model, there are other error models such as proportional, exponential, and combined error models. For details, see Error Models.

Nonlinear Mixed-Effects Modeling Workflow

SimBiology^® lets you estimate fixed effects θs and random effects ηs as well as the covariance matrix of random effects Ψ. However, you cannot alter A and B design matrices since they are automatically determined from the covariate model you specify. Use the sbiofitmixed function to estimate nonlinear mixed-effects parameters. These steps show one of the workflows you can use at the command line.

Import data.
Convert the data to the groupedData format.
Define dosing data. For details, see Doses in SimBiology Models.
Create a structural model (one-, two-, or multicompartment model). For details, see Create Pharmacokinetic Models.
Create a covariate model to define parameter-covariate relationships if any. For details, see Specify a Covariate Model.
Map the response variable from data to the model component. For example, if you have the measured drug concentration data for the central compartment, then map it to the drug species in the central compartment (typically the Drug_Central species).
Specify parameters to estimate using the EstimatedInfo object. It lets you optionally specify parameter transformations, initial values, and parameter bounds. Supported transforms are log, probit, logit, and none (no transform).
(Optional) You can also specify an error model. The default model is the constant error model. For instance, you can change it to the proportional error model if you assume the measurement error is proportional to the response data. See Specify an Error Model.
Estimate parameters using fitproblem or sbiofitmixed, which performs Maximum Likelihood Estimation.
(Optional) If you have a large, complex model, the estimation might take longer. SimBiology lets you check the status of fitting as it progresses. See Obtain the Fitting Status.

For workflow examples, see:

Specify a Covariate Model

When specifying a nonlinear mixed-effects model, you define parameter-covariate relationship using a covariate model (CovariateModel). For example, suppose you have PK profile data for multiple individuals and are estimating three parameters (clearance Cl, compartment volume V, and elimination rate k) that have both fixed and random effects. Assume the clearance Cl has a correlation with a covariate variable weight (w) of each individual. Each parameter can be described as a linear combination of fixed and random effects.

$C l_{i} = θ_{1} + θ_{2} * w_{i} + η_{1 i}$ ,

$V_{i} = θ_{3} + η_{2 i}$ ,

$k_{i} = θ_{4} + η_{3 i}$ ,

where θs represent fixed effects and ηs represent random effects, which are normally distributed $(\begin{array}{l} η_{1 i} \\ η_{2 i} \\ η_{3 i} \end{array}) \sim M V N (0, Ψ)$ . By default, the random effects are uncorrelated. So $Ψ = (\begin{matrix} σ_{1}^{2} & 0 & 0 \\ 0 & σ_{2}^{2} & 0 \\ 0 & 0 & σ_{3}^{2} \end{matrix})$ .

Construct an empty CovariateModel object.
```
covModel = CovariateModel;
```
Set the Expression property to define the relationships between parameters (Cl, V, and k) and covariate (w). You must use theta as a prefix for all fixed effects and eta for random effects.
```
covModel.Expression = {'Cl = theta1 + theta2*w + eta1','V = theta3 + eta2','k = theta4 + eta3'};
```
The FixedEffectNames property displays the fixed effects defined in the model.
```
covModel.FixedEffectNames
```
```
ans = 

    'theta1'
    'theta3'
    'theta4'
    'theta2'
```
The FixedEffectDescription property displays which fixed effects correspond to which parameter. For instance, theta1 is the fixed effect for the Cl parameter, and theta2 is the fixed effect for the weight covariate that has a correlation with Cl parameter, denoted as Cl/w.
```
covModel.FixedEffectDescription
```
```
ans = 

    'Cl'
    'V'
    'k'
    'Cl/w'
```

Specify initial estimates for the fixed effects. Create a structure containing initial estimates using the constructDefaultFixedEffectValues function.

initialEstimates = constructDefaultFixedEffectValues(covModel)

initialEstimates = 

    theta1: 0
    theta2: 0
    theta3: 0
    theta4: 0

initialEstimates.theta1 = 1.20;
initialEstimates.theta2 = 0.30;
initialEstimates.theta3 = 0.90;
initialEstimates.theta4 = 0.10;

Set the initial estimates to the FixedEffectValues property.
```
covModel.FixedEffectValues = initialEstimates;
```

Specify a Covariance Pattern Among Random Effects

By default, sbiofitmixed assumes no covariance among random effects, that is, a diagonal covariance matrix is used. Suppose you have η₁, η₂, and η₃, and there is a covariance σ₁₂ between η₁ and η₂. You can indicate this using a pattern matrix where 1 indicates a variance or covariance parameter which is estimated by sbiofitmixed. For instance, a pattern matrix $(\begin{array}{l} 1 1 0 \\ 1 1 0 \\ 0 0 1 \end{array})$ represents $(\begin{matrix} σ_{1}^{2} & σ_{12} & 0 \\ σ_{21} & σ_{2}^{2} & 0 \\ 0 & 0 & σ_{3}^{2} \end{matrix})$ .

Define such a pattern using an options struct.

options.CovPattern = [1 1 0;1 1 0;0 0 1];

Then you can use options as an input argument for sbiofitmixed. For a complete workflow, see Nonlinear Mixed-Effects Modeling Workflow.

Specify an Error Model

During the Nonlinear Mixed-Effects Modeling Workflow, you can optionally specify an error model using a structure.

options.ErrorModel = 'proportional';

Then you can use options as one of the input arguments when you run sbiofitmixed.

Supported error models are constant (default), proportional, combined, and exponential models. For details, see Error Models.

Maximum Likelihood Estimation

SimBiology estimates the parameters of a nonlinear mixed-effects model by maximizing a likelihood function. The function can be described as:

$p (y | θ, σ^{2}, Ψ) = \int p (y | θ, η, σ^{2}) p (η | Ψ) d η$ ,

where y is the response data, θ is the vector of fixed effects, σ² is the error variance, Ψ is the covariance matrix for random effects, and η is the vector of unobserved random effects. $p (y | θ, σ^{2}, Ψ)$ is the marginal density of y, $p (y | θ, η, σ^{2})$ is the conditional density of y given the random effects η, and the prior distribution of η is $p (η | Ψ)$ .

This integral contains a nonlinear function of the fixed effects and variance parameters that you want to maximize. Typically for nonlinear models, the integral does not have a closed form, and needs to be solved numerically, which involves simulating the function at each time step of an optimization algorithm. Therefore, the estimation can take a long time for complex models, and initial values of parameters might play an important role for successful convergence. SimBiology provides these iterative algorithms to solve the integral and maximize the likelihood if you have Statistics and Machine Learning Toolbox™.

LME — Use the likelihood for the linear mixed-effects model at the current conditional estimates of θ and η. This is the default.
RELME — Use the restricted likelihood for the linear mixed-effects model at the current conditional estimates of θ and η.
FO — First-order (Laplacian) approximation without random effects.
FOCE — First-order (Laplacian) approximation at the conditional estimates of θ.
stochastic EM — Use the Expectation-Maximization (EM) algorithm in which the E step is replaced by a stochastic procedure.

For a complete workflow, see Nonlinear Mixed-Effects Modeling Workflow.

Obtain the Fitting Status

During the estimation of mixed-effects parameters of a large and complex model that may take a longer time, you may want to obtain the status of fitting as it progresses. Set 'ProgressPlot' to true when you run sbiofitmixed to display the progress of the fitting. For details, see Progress Plot.