sbiofitmixed

Fit nonlinear mixed-effects model (requires Statistics and Machine Learning Toolbox software)

Syntax

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo)

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo,dosing)

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo,dosing,functionName)

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo,dosing,functionName,opt)

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo,dosing,functionName,opt,variants)

fitResults = sbiofitmixed(_,'UseParallel',tf_parallel)

fitResults = sbiofitmixed(_,'ProgressPlot',tf_progress)

[fitResults,simDataI,simDataP] = sbiofitmixed(_)

Description

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo) performs nonlinear mixed-effects estimation using the SimBiology^® model sm and returns a NLMEResults object fitResults.

grpData is a groupedData object specifying the data to fit. responseMap defines the mapping between the model components and response data in grpData. covEstiminfo is a CovariateModel object or an array of estimatedInfo objects that defines the parameters to be estimated.

If the model contains active doses and variants, they are applied during the simulation.

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo,dosing) uses the dosing information specified by a matrix of SimBiology dose objects dosing instead of using the active doses of the model sm if there are any.

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo,dosing,functionName) uses the estimation function specified by functionName that must be either 'nlmefit' or 'nlmefitsa'.

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo,dosing,functionName,opt) uses the additional options specified by opt for the estimation function functionName.

fitResults = sbiofitmixed(sm,grpData,responseMap,covEstiminfo,dosing,functionName,opt,variants) applies variant objects specified as variants instead of using any active variants of the model.

fitResults = sbiofitmixed(_,'UseParallel',tf_parallel) specifies whether to estimate parameters in parallel if Parallel Computing Toolbox™ is available.

fitResults = sbiofitmixed(_,'ProgressPlot',tf_progress) specifies whether to show the progress of parameter estimation.

[fitResults,simDataI,simDataP] = sbiofitmixed(_) returns a vector of results objects fitResults, vector of simulation results simDataI using individual-specific parameter estimates, and vector of simulation results simDataP using population parameter estimates.

Note

• sbiofitmixed unifies sbionlmefit and sbionlmefitsa estimation functions. Use sbiofitmixed to perform nonlinear mixed-effects modeling and estimation.

• sbiofitmixed simulates the model using a SimFunction object, which automatically accelerates simulations by default. Hence it is not necessary to run sbioaccelerate before you call sbiofitmixed.

Examples

Fit a One-Compartment PK Model to the Phenobarbital Data

This example uses data collected on 59 preterm infants given phenobarbital during the first 16 days after birth [1]. Each infant received an initial dose followed by one or more sustaining doses by intravenous bolus administration. A total of between 1 and 6 concentration measurements were obtained from each infant at times other than dose times, for a total of 155 measurements. Infant weights and APGAR scores (a measure of newborn health) were also recorded.

Load the data.

load pheno.mat ds

Convert the dataset to a groupedData object, a container for holding tabular data that is divided into groups. It can automatically identify commonly used variable names as the grouping variable or independent (time) variable. Display the properties of the data and confirm that GroupVariableName and IndependentVariableName are correctly identified as 'ID' and 'TIME', respectively.

data = groupedData(ds);
data.Properties

ans = struct with fields:
                Description: ''
                   UserData: []
             DimensionNames: {'Observations'  'Variables'}
              VariableNames: {'ID'  'TIME'  'DOSE'  'WEIGHT'  'APGAR'  'CONC'}
       VariableDescriptions: {}
              VariableUnits: {}
         VariableContinuity: []
                   RowNames: {}
           CustomProperties: [1x1 matlab.tabular.CustomProperties]
          GroupVariableName: 'ID'
    IndependentVariableName: 'TIME'

Create a simple one-compartment PK model with bolus dosing and linear clearance to fit such data. Use the PKModelDesign object to construct the model. Each compartment is defined by a name, dosing type, a clearance type, and whether or not the dosing requires a lag parameter. After constructing the model, you can also get a PKModelMap object map that lists the names of species and parameters in the model that are most relevant for fitting.

pkmd = PKModelDesign;
addCompartment(pkmd,'Central','DosingType','Bolus',...
                    'EliminationType','linear-clearance',...
                    'HasResponseVariable',true,'HasLag',false);
[onecomp, map] = pkmd.construct;

Describe the experimentally measured response by mapping the appropriate model component to the response variable. In other words, indicate which species in the model corresponds to which response variable in the data. The PKModelMap property Observed indicates that the relevant species in the model is Drug_Central, which represents the drug concentration in the system. The relevant data variable is CONC, which you visualized previously.

map.Observed

ans = 1x1 cell array
    {'Drug_Central'}

Map the Drug_Central species to the CONC variable.

responseMap = 'Drug_Central = CONC';

The parameters to estimate in this model are the volume of the central compartment Central and the clearance rate Cl_Central. The PKModelMap property Estimated lists these relevant parameters. The underlying algorithm of sbiofit assumes parameters are normally distributed, but this assumption may not be true for biological parameters that are constrained to be positive, such as volume and clearance. Specify a log transform for the estimated parameters so that the transformed parameters follow a normal distribution. Use an estimatedInfo object to define such transforms and initial values (optional).

map.Estimated

ans = 2x1 cell
    {'Central'   }
    {'Cl_Central'}

Define such estimated parameters, appropriate transformations, and initial values.

estimatedParams = estimatedInfo({'log(Central)','log(Cl_Central)'},'InitialValue',[1 1]);

Each infant received a different schedule of dosing. The amount of drug is listed in the data variable DOSE. To specify these dosing during fitting, create dose objects from the data. These objects use the property TargetName to specify which species in the model receives the dose. In this example, the target species is Drug_Central, as listed by the PKModelMap property Dosed.

map.Dosed

ans = 1x1 cell array
    {'Drug_Central'}

Create a sample dose with this target name and then use the createDoses method of groupedData object data to generate doses for each infant based on the dosing data DOSE.

sampleDose = sbiodose('sample','TargetName','Drug_Central');
doses = createDoses(data,'DOSE','',sampleDose);

Fit the model.

[nlmeResults,simI,simP] = sbiofitmixed(onecomp,data,responseMap,estimatedParams,doses,'nlmefit');

Visualize the fitted results using individual-specific parameter estimates.

fig1 = plot(nlmeResults,'ParameterType','individual');

% Resize the figure.
fig1.hFig.Position(:) = [100 100 1200 800];

Visualize the fitted results using population parameter estimates.

fig2 = plot(nlmeResults,'ParameterType','population'); 

% Resize the figure.
fig2.hFig.Position(:) = [100 100 1200 800];

Input Arguments

sm — SimBiology model

SimBiology model object

SimBiology model, specified as a SimBiology model object. The active configset object of the model contains solver settings for simulation. Any active doses and variants are applied to the model during simulation unless specified otherwise using the dosing and variants input arguments, respectively.

grpData — Data to fit

groupedData object

Data to fit, specified as a groupedData object.

The name of the time variable must be defined in the IndependentVariableName property of grpData. For instance, if the time variable name is 'TIME', then specify it as follows.

grpData.Properties.IndependentVariableName = 'TIME';

grpData must have at least two groups, and the name of grouping variable name must be defined in the GroupVariableName property of grpData. For example, if the grouping variable name is 'GROUP', then specify it as follows.

grpData.Properties.GroupVariableName = 'GROUP';

A group usually refers to a set of measurements that represent a single time course, often corresponding to a particular individual or experimental condition.

Note

sbiofitmixed uses the categorical function to identify groups. If any group values are converted to the same value by categorical, then those observations are treated as belonging to the same group. For instance, if some observations have no group information (that is, empty character vector), then categorical converts empty character vectors to <undefined>, and these observations are treated as one group.

responseMap — Mapping information of model components to response data

character vector | string | string vector | cell array of character vectors

Mapping information of model components to grpData, specified as a character vector, string, string vector, or cell array of character vectors.

Each character vector or string is an equation-like expression, similar to assignment rules in SimBiology. It contains the name (or qualified name) of a quantity (species, compartment, or parameter) or an observable object in the model sm, followed by the character '=' and the name of a variable in grpData. For clarity, white spaces are allowed between names and '='.

For example, if you have the concentration data 'CONC' in grpData for a species 'Drug_Central', you can specify it as follows.

responseMap = 'Drug_Central = CONC';

To name a species unambiguously, use the qualified name, which includes the name of the compartment. To name a reaction-scoped parameter, use the reaction name to qualify the parameter.

If the model component name or grpData variable name is not a valid MATLAB^® variable name, surround it by square brackets, such as:

responseMap = '[Central 1].Drug = [Central 1 Conc]';

If a variable name itself contains square brackets, you cannot use it in the expression to define the mapping information.

An error is issued if any (qualified) name matches two components of the same type. However, you can use a (qualified) name that matches two components of different types, and the function first finds the species with the given name, followed by compartments and then parameters.

covEstiminfo — Estimated parameters

vector of estimatedInfo objects | CovariateModel object

Estimated parameters, specified as a vector of estimatedInfo objects or a CovariateModel object that defines the estimated parameters in the model sm, their initial estimates (optional), and their relationship to group-specific covariates included in grpData (optional). If this is a vector of estimatedInfo objects, then no covariates are used, and all parameters are estimated with group-specific random effects.

You can also specify parameter transformations if necessary. Supported transforms are log, logit, and probit. For details, see EstimatedInfo object and CovariateModel object.

If covEstiminfo is a vector of estimatedInfo objects, the CategoryVariableName property of each of these objects is ignored.

dosing — Dosing information

[] | 2-D matrix of dose objects

Dosing information, specified as an empty array or 2-D matrix of dose objects (ScheduleDose object or RepeatDose object).

• If empty, no doses are applied during simulation, even if the model has active doses.

• If not empty, the matrix must have a single row or one row per group in the input data. If it has a single row, same doses are applied to all groups during simulation. If it has multiple rows, each row is applied to a separate group, in the same order as the groups appear in the input data.

• Multiple columns are allowed so that you can apply multiple dose objects to each group. All doses in a column must reference the same components in the model. For example, doses in the same column must have the same dose target (TargetName). If you parameterize any property of a dose, then all doses within the column must have this property parameterized to the same parameter. If some groups require more doses than others, then fill in the matrix with default (dummy) doses.

• A default dose has default values for all properties, except for the Name property. Create a default dose as follows.

  d1 = sbiodose('d1');

• In addition to manually constructing dose objects using sbiodose, if grpData has dosing information, you can use the createDoses method to construct doses.

functionName — Estimation function name

character vector | string

Estimation function name, specified as a character vector or string. Choices are 'nlmefit' or 'nlmefitsa'. For the summary supported methods and fitting options, see Supported Methods for Parameter Estimation in SimBiology.

opt — Options specific to estimation function

struct

Options specific to the estimation function, specified as a structure. The structure contains fields and default values that are the name-value pair arguments accepted by nlmefit (Statistics and Machine Learning Toolbox) and nlmefitsa (Statistics and Machine Learning Toolbox), except the following that are not supported.

• 'FEConstDesign'

• 'FEGroupDesign’

• 'FEObsDesign'

• 'FEParamsSelect'

• 'ParamTransform'

• 'REConstDesign'

• 'REGroupDesign'

• 'REObsDesign'

• 'Vectorization'

'REParamsSelect' is only supported when covEstiminfo is a vector of estimatedInfo objects.

Use the statset (Statistics and Machine Learning Toolbox) function only to set the 'Options' field of the opt structure as follows.

opt.Options = statset('Display','iter','TolX',1e-3,'TolFun',1e-3);

For other supported name-value pair arguments (see nlmefit (Statistics and Machine Learning Toolbox) and nlmefitsa (Statistics and Machine Learning Toolbox)), set them as follows.

opt.ErrorModel = 'proportional';
opt.ApproximationType = 'LME';

variants — Variants

[] | vector of variant objects

Variants, specified as an empty array or vector of variant objects. If empty, no variants are used even if the model has active variants.

tf_parallel — 'UseParallel' option

false (default) | true

'UseParallel' option, specified as true or false. If true, and Parallel Computing Toolbox is available, the function performs parameter estimation in parallel.

tf_progress — 'ProgressPlot' option

false (default) | true

'ProgressPlot' option, specified as true or false. If true, a new figure opens containg plots that show the progress of parameter estimation. Specifically, the plots show the values of fixed effects parameters (theta), the estimates of the variance parameters, that is, the diagonal elements of the covariance matrix of the random effects (Ψ), and the log-likelihood. For details, see Progress Plot.

Output Arguments

fitResults — Estimation results

NLMEResults object

Estimation results, returned as an NLMEResults object.

simDataI — Simulation results

vector of SimData objects

Simulation results, returned as a vector of SimData objects representing simulation results for each group (or individual) using fixed-effect and random-effect estimates (individual-specific parameter estimates).

The states reported in simDataI are the states that were included in the responseMap input argument as well as any other states listed in the StatesToLog property of the runtime options (RuntimeOptions) of the SimBiology model sm.

simDataP — Simulation results

vector of SimData objects

Simulation results, returned as a vector of SimData objects representing simulation results for each group (or individual) using only fixed-effect estimates (population parameter estimates).

The states reported in simDataP are the states that were included in the responseMap input argument as well as any other states listed in the StatesToLog property of the runtime options (RuntimeOptions) of the SimBiology model sm.

References

[1] Grasela Jr, T.H., Donn, S.M. (1985) Neonatal population pharmacokinetics of phenobarbital derived from routine clinical data. Dev Pharmacol Ther. 8(6), 374–83.

Extended Capabilities

Automatic Parallel Support

Accelerate code by automatically running computation in parallel using Parallel Computing Toolbox™.

To run in parallel, set 'UseParallel' to true.

For more information, see the 'UseParallel' name-value pair argument.