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 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 $${\epsilon}_{ij}\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}={\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 {\rm N}(0,{\sigma}_{\eta}^{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}={\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({x}_{ij},{p}_{i})+{\epsilon}_{ij}\\ {p}_{i}={A}_{i}\theta +{B}_{i}{\eta}_{i}\text{}\\ {\epsilon}_{ij}\text{\hspace{0.17em}}\sim N(0,{\sigma}^{2})\text{}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\\ {\eta}_{i}\sim N(0,\Psi )\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 |
σ^{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.
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.
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 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 a complete workflow example, see Modeling the Population Pharmacokinetics of Phenobarbital in Neonates.
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(0,\Psi )$$. 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)$$.
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;
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.
During the Nonlinear Mixed-Effects Modeling Workflow, you can optionally specify an error model using a structure.
options.ErrorModel = 'proportional';
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.
SimBiology estimates the parameters of a nonlinear mixed-effects model by maximizing a likelihood function. The function can be described as:
$$p(y|\theta ,{\sigma}^{2},\Psi )={\displaystyle \int p(}y|\theta ,\eta ,{\sigma}^{2})p(\eta |\Psi )\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(y|\theta ,{\sigma}^{2},\Psi )$$ is the marginal density of y, $$p(y|\theta ,\eta ,{\sigma}^{2})$$ is the conditional density of y given the random effects η, and the prior distribution of η is $$p(\eta |\Psi )$$.
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.
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.