## 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}_{ij}=\text{​}\text{\hspace{0.17em}}\frac{{D}_{i}}{{V}_{}}{e}^{-{k}_{i}{t}_{ij}}+a{\epsilon }_{ij}$,

where yij is the jth response of the ith individual, Di is the dose administered to the ith individual, V is the population mean volume of distribution, a is an error parameter, and ${\epsilon }_{ij}\sim N\left(0,1\right)$, representing some measurement error. The elimination rate parameter (ki) depends on the clearance and volume of the central compartment ${k}_{i}=\frac{C{l}_{i}}{{V}_{}}$. Both ki and Cli 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}={\theta }_{1}+{\eta }_{i}$, where θ1 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 ${\eta }_{i}\sim Ν\left(0,{\sigma }_{\eta }^{2}\right)$.

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 wi is the weight of the ith individual, then the model becomes $C{l}_{i}={\theta }_{1}+{\theta }_{2}*{w}_{i}+{\eta }_{i}$, where θ2 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}_{ij}=f\left({x}_{ij},{p}_{i}\right)+{\epsilon }_{ij}\\ {p}_{i}={A}_{i}\theta +{B}_{i}{\eta }_{i}\text{​}\\ {\epsilon }_{ij}\text{\hspace{0.17em}}\sim N\left(0,{\sigma }^{2}\right)\text{​}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\\ {\eta }_{i}\sim N\left(0,\Psi \right)\end{array}$

 yij Data vector of individual-specific response values f General, real-valued function of pi and xij xij Data matrix of individual-specific predictor values pi Vector of individual-specific model parameters θ Vector of fixed effects, modeling population parameters ηi Vector of multivariate normally distributed individual-specific random effects Ai Individual-specific design matrix for combining fixed effects Bi 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 σ2 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.

1. Convert the data to the `groupedData` format.

2. Define dosing data. For details, see Doses in SimBiology Models.

3. Create a structural model (one-, two-, or multicompartment model). For details, see Create Pharmacokinetic Models.

4. Create a covariate model to define parameter-covariate relationships if any. For details, see Specify a Covariate Model.

5. 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).

6. 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).

7. (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.

8. Estimate parameters using `sbiofitmixed`, which performs Maximum Likelihood Estimation.

9. (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 a complete workflow example, see Modeling the Population Pharmacokinetics of Phenobarbital in Neonates.

### Specify a Covariate Model

When specifying a nonlinear mixed-effects model, you define parameter-covariate relationship using a covariate model (`CovariateModel object`). 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}={\theta }_{1}+{\theta }_{2}*{w}_{i}+{\eta }_{1i}$,

${V}_{i}={\theta }_{3}+{\eta }_{2i}$,

${k}_{i}={\theta }_{4}+{\eta }_{3i}$,

where θs represent fixed effects and ηs represent random effects, which are normally distributed $\left(\begin{array}{l}{\eta }_{1i}\\ {\eta }_{2i}\\ {\eta }_{3i}\end{array}\right)\sim MVN\left(0,\Psi \right)$. By default, the random effects are uncorrelated. So $\Psi =\left(\begin{array}{ccc}{\sigma }_{1}^{2}& 0& 0\\ 0& {\sigma }_{2}^{2}& 0\\ 0& 0& {\sigma }_{3}^{2}\end{array}\right)$.

1. Construct an empty `CovariateModel` object.

`covModel = CovariateModel;`

2. 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'```
3. 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;```

4. 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 η1, η2, and η3, and there is a covariance σ12 between η1 and η2. 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 $\left(\begin{array}{l}1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}0\\ 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}0\\ 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}1\end{array}\right)$ represents $\left(\begin{array}{ccc}{\sigma }_{1}^{2}& {\sigma }_{12}& 0\\ {\sigma }_{21}& {\sigma }_{2}^{2}& 0\\ 0& 0& {\sigma }_{3}^{2}\end{array}\right)$.

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\left(y|\theta ,{\sigma }^{2},\Psi \right)=\int p\left(y|\theta ,\eta ,{\sigma }^{2}\right)p\left(\eta |\Psi \right)\text{\hspace{0.17em}}d\eta$,

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

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.

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