nlmefitsa

Fit nonlinear mixed-effects model with stochastic EM algorithm

collapse all in page

Syntax

beta = nlmefitsa(X,y,group,V,modelfun,beta0)

beta = nlmefitsa(X,y,group,V,modelfun,beta0,Name=Value)

[beta,psi] = nlmefitsa(___)

[beta,psi,stats] = nlmefitsa(___)

[beta,psi,stats,b] = nlmefitsa(___)

Description

beta = nlmefitsa(X,y,group,V,modelfun,beta0) uses the stochastic approximation expectation-maximization (SAEM) algorithm to fit the nonlinear mixed-effects regression model modelfun to the data in X, y, group, and V and returns the fixed-effects estimates in beta. nlmefitsa uses the initial fixed-effects values in beta0 to fit the model.

For more information about the SAEM algorithm, see Algorithms.

example

beta = nlmefitsa(X,y,group,V,modelfun,beta0,Name=Value) specifies additional options using one or more name-value arguments. For example, you can specify design matrices for the fixed and random effects, and specify the method for approximating the loglikelihood.

[beta,psi] = nlmefitsa(___) additionally returns a matrix of covariance estimates for the random effects, using any of the input argument combinations in the previous syntaxes.

[beta,psi,stats] = nlmefitsa(___) additionally returns a structure that contains model statistics, including the maximized loglikelihood, root mean squared error, Akaike information criterion (AIC), and Bayesian information criterion (BIC) for the fitted model.

example

[beta,psi,stats,b] = nlmefitsa(___) additionally returns the random-effects estimates.

example

Examples

collapse all

Specify Group Predictors

Open Live Script

Load nonlinear sample data and display the predictor, response, and grouping variables.

load nonlineardata.mat
[X,y,group]

ans = 30×5

    8.1472    0.7060   75.1267  573.4851    1.0000
    9.0579    0.0318   25.5095  188.3748    1.0000
    1.2699    0.2769   50.5957  356.7075    1.0000
    9.1338    0.0462   69.9077  499.6050    1.0000
    6.3236    0.0971   89.0903  631.6939    1.0000
    0.9754    0.8235   95.9291  679.1466    1.0000
    2.7850    0.6948   54.7216  398.8715    1.0000
    5.4688    0.3171   13.8624  109.1202    1.0000
    9.5751    0.9502   14.9294  207.5047    1.0000
    9.6489    0.0344   25.7508  190.7724    1.0000
    1.5761    0.4387   84.0717  593.2222    1.0000
    9.7059    0.3816   25.4282  203.1922    1.0000
    9.5717    0.7655   81.4285  634.8833    1.0000
    4.8538    0.7952   24.3525  205.9043    1.0000
    8.0028    0.1869   92.9264  663.2529    1.0000
      ⋮

Display the group predictor variables.

X is a matrix of predictor data, and y is a vector containing the response variable. group contains data for the grouping variable, and v contains the group predictor variables.

Create an anonymous nonlinear function that accepts a vector of coefficients, a matrix of predictor data, and a vector of group predictor variables.

model = @(phi,xfun,vfun)(phi(1).*xfun(:,1).*exp(phi(2).*xfun(:,2)./vfun)+phi(3).*xfun(:,3))

model = function_handle with value:
    @(phi,xfun,vfun)(phi(1).*xfun(:,1).*exp(phi(2).*xfun(:,2)./vfun)+phi(3).*xfun(:,3))

model is a handle for a function given by the formula

$ϕ_{1} X_{1} \exp (\frac{ϕ_{2} X_{2}}{V}) + ϕ_{3} X_{3}$ .

Fit model to the data in X, y, group, and v using the nlmefitsa function. Specify a vector of ones as the initial estimate for the fixed-effects coefficients.

beta = nlmefitsa(X,y,group,V,model,[1 1 1])

beta = 3×1

    1.0008
    4.9980
    6.9999

The output shows the progress of estimating the fixed effects and the elements of the random-effects covariance matrix, as well as the final fixed-effects estimates beta. In each plot, the horizontal axis shows the iteration step, and the vertical axis shows the value of the estimation.

Specify Output Function

Open Live Script

Load some nonlinear sample data.

load nonlineardata.mat;

X is a matrix of predictor data and y is a vector containing the response variable. group and V, respectively contain data for a grouping variable and group predictor.

Create an anonymous nonlinear function that accepts a vector of coefficients, a matrix of predictor data, and a vector of group predictor data.

model = @(phi,xfun,vfun)(phi(1).*xfun(:,1).*exp(phi(2).*xfun(:,2)./vfun)+phi(3).*xfun(:,3));

model is a handle for a function given by the formula $ϕ_{1} X_{1} \exp (\frac{ϕ_{2} X_{2}}{V}) + ϕ_{3} X_{3}$ .

Define an output function for nlmefitsa. For more information about the form of the output function, see the OutputFcn field description for the Options name-value argument.

function stop = outputFunction(beta,status,state)

    stop = 0;

    hold on
    plot3(status.iteration,beta(2),beta(1),"mo")

    state = string(state);

    if state=="done"
        stop=1;
    end

end

outputFunction is a function that plots the iteration number for the fitting algorithm together with the first and second fixed effects. outputFunction returns 1 when the fitting algorithm has completed its final iteration.

Use the statset function to create an options structure for nlmefitsa that uses outputFunction as its output function.

default_opts=statset("nlmefitsa");
opts = statset(default_opts,OutputFcn=@outputFunction);

opts is a statistics options structure that contains options for the stochastic expectation maximization fitting algorithm.

Create a figure and define axes for outputFunction to plot into. Fit model to the predictor data in X and the response data in y using the options in opts.

figure
ax = axes(view=[12,10]);
xlabel("Iteration")
ylabel("beta(2)")
zlabel("beta(1)")
box on

[beta,psi,stats] = nlmefitsa(X,y,group,V,model,[1 1 1],Options=opts)

Figure contains an axes object. The axes object with xlabel Iteration, ylabel beta(2) contains 168 objects of type line.

beta = 3×1

    1.0008
    4.9980
    6.9999

psi = 3×3
10^-4 ×

    0.0415         0         0
         0    0.2912         0
         0         0    0.0004

stats = struct with fields:
          logl: []
           aic: []
           bic: []
        sebeta: []
           dfe: 23
          covb: []
    errorparam: 0.0139
          rmse: 0.0012
          ires: [30×1 double]
          pres: [30×1 double]
         iwres: [30×1 double]
         pwres: [30×1 double]
         cwres: [30×1 double]

nlmefitsa calls outputFunction after every iteration of the fitting algorithm. The figure shows that the fixed effects beta(1) and beta(2) are near 1 when the iteration number is near 0. This is consistent with the initial values for beta. As the iteration number increases, beta(1) makes a large jump up to around 3.3 before converging to 1. As beta(1) converges to 1, beta(2) converges to 5. The output argument beta contains the final values for the fixed effects. psi and stats, respectively, contain the covariance matrix for the random effects and additional statistics about the fit.

Transform Parameters and Plot the Fitted Model

Open Live Script

Load the indomethacin data set.

load indomethacin

The variables time, concentration, and subject contain time series data for the blood concentration of the drug indomethacin in six patients.

Create an anonymous nonlinear function that accepts a vector of coefficients and a vector of predictor variables.

model = @(phi,t)(phi(1).*exp(-phi(2).*t)+phi(3).*exp(-phi(4).*t));

model is a handle for a function given by the formula

$c = ϕ_{1} \exp (- ϕ_{2} t) + ϕ_{3} \exp (- ϕ_{4} t)$ ,

where $c$ is the concentration of indomethacin, $ϕ_{i}$ for $i = 1, 2, 3, 4$ are coefficients, and $t$ is time. The function does not contain group-specific predictor variables because the formula does not include them.

Fit the model to the data using time as the predictor variable, subject as the grouping variable, and concentration as the response. Specify a log transformation function for the second and fourth coefficients.

phi0 = [1 1 1 1];
xform = [0 1 0 1];
[beta,psi,stats,b] = nlmefitsa(time,concentration, ...
   subject,[],model,phi0,ParamTransform=xform)

psi = 4×4

    0.0585         0         0         0
         0    0.0248         0         0
         0         0    0.5068         0
         0         0         0    0.0139

stats = struct with fields:
          logl: []
           aic: []
           bic: []
        sebeta: []
           dfe: 57
          covb: []
    errorparam: 0.0811
          rmse: 0.0772
          ires: [66×1 double]
          pres: [66×1 double]
         iwres: [66×1 double]
         pwres: [66×1 double]
         cwres: [66×1 double]

b = 4×6

   -0.2302   -0.0033    0.1625    0.1774   -0.3334    0.1129
    0.0363   -0.1502    0.0071    0.0471    0.0068   -0.0481
   -0.7631   -0.0553    0.8780   -0.8120    0.5429    0.1695
   -0.0030   -0.0223    0.0192   -0.0830    0.0505   -0.0066

The output argument beta contains the fixed effects for the model, and b contains the random effects. The plots show the progress of the Monte Carlo simulation used to fit the coefficients. The maximum likelihood estimates for beta and the random-effects covariance matrix psi converge after about 300 iterations.

Plot the sample data together with the model, using only the fixed effects in beta for the model coefficients. Use the gscatter function to color code the data according to the subject. To reverse the log transformation on the second and fourth coefficients, take their exponentials using the exp function.

figure
hold on

gscatter(time,concentration,subject);

phi = [beta(1),exp(beta(2)),beta(3),exp(beta(4))];
tt = linspace(0,8);
cc = model(phi,tt);
plot(tt,cc,LineWidth=2,Color="k")
 
legend("Subject 1","Subject 2","Subject 3",...
    "Subject 4","Subject 5","Subject 6","Fitted curve")
xlabel("Time (hours)")
ylabel("Concentration (mcg/ml)")
title("Indomethacin Elimination")
hold off

Figure contains an axes object. The axes object with title Indomethacin Elimination, xlabel Time (hours), ylabel Concentration (mcg/ml) contains 7 objects of type line. One or more of the lines displays its values using only markers These objects represent Subject 1, Subject 2, Subject 3, Subject 4, Subject 5, Subject 6, Fitted curve.

The plot shows that the blood concentration of indomethacin decreases over eight hours, and the fitted model passes through the bulk of the data.

Plot the data together with the model again, using both the fixed effects and the random effects in b for the model coefficients. For each subject, plot the data and the fitted curve in the same color.

figure
hold on

h = gscatter(time,concentration,subject);

for j=1:6
    phir = [beta(1)+b(1,j),exp(beta(2)+b(2,j)), ...
            beta(3)+b(3,j),exp(beta(4)+b(4,j))];
    ccr = model(phir,tt);
    col = h(j).Color;
    plot(tt,ccr,Color=col)
end

legend("Subject 1","Subject 2","Subject 3",...
    "Subject 4","Subject 5","Subject 6",...
    "Fitted curve 1","Fitted curve 2","Fitted curve 3",...
    "Fitted curve 4","Fitted curve 5","Fitted curve 6")
xlabel("Time (hours)")
ylabel("Concentration (mcg/ml)")
title("Indomethacin Elimination")

hold off

The plot shows that, for each subject, the fitted curve follows the bulk of the data more closely than the curve for the fixed-effects model in the previous figure. This result suggests that the random effects improve the fit of the model.

Input Arguments

collapse all

`X` — Predictor variables
numeric matrix

Predictor variables, specified as an n-by-h numeric matrix, where n is the number of observations and h is the number of predictor variables.

Data Types: single | double

`y` — Response variable
numeric vector

Response variable, specified as a numeric vector with n elements, where n is the number of observations.

Data Types: single | double

`group` — Grouping variable
numeric vector | categorical vector | string vector | cell array of character vectors

Grouping variable, specified as a numeric, categorical, or string vector, or a cell array of character vectors.

Example: ["subject1" "subject4" "subject3"…"subject3" "subject2"]

Example: [ones(50,1);2*ones(50,1);3*ones(50,1)]

Data Types: double | single | string | cell | categorical

`V` — Group predictor variables
numeric matrix | cell array

Group predictor variables, specified as a numeric matrix or cell array. Group predictor variables are variables that have the same value for all observations in a group.

If the group-specific predictor variables are the same size across groups, specify V as an m-by-g matrix, where m is the number of groups and g is the number of group predictor variables.
If the group-specific predictor variables vary in size across groups, specify V as an m-by-g cell array.
If you do not have group predictor variables, specify V as [].

Example: [1 2; 2 4; 5 2; 7 1; 4 5; 3 6]

Data Types: single | double | cell

`modelfun` — Model function
function handle

Model function, specified as a function handle. If V is empty, the function can have the form yfit = modelfun(phi,xfun). Otherwise, it must have the form yfit = modelfun(phi,xfun,vfun).

Argument	Description
`phi`	An array of model coefficients. The size of `phi` depends on the value of the `Vectorization` name-value argument. `"SinglePhi"` — `phi` is a vector of coefficients for the observations in `xfun`. This value is the default when you do not specify `Vectorization`. `"SingleGroup"` or `"Full"` — `phi` is a vector or matrix of coefficients. If `phi` is a vector, the coefficients apply to all observations in `xfun`. If `phi` is a matrix, it has the same number of rows as `xfun`. Each row of `phi` corresponds to an observation in `xfun`, and each column corresponds to a coefficient.
`xfun`	An array of predictor variables. The size of `xfun` depends on the value of the `Vectorization` name-value argument. `"SinglePhi"` — `xfun` is a vector containing a single observation from `X`, or a matrix containing observations for a single group. This value is the default when you do not specify `Vectorization`. `"SingleGroup"` — `xfun` contains observations from a single group. `"Full"` — `xfun` represents `X` or a subset of observations from `X`.
`vfun`	An array of group predictor variables. The size of `vfun` depends on the value of the `Vectorization` name-value argument. `"SinglePhi"` — `vfun` is a vector or matrix. When `vfun` is a vector, it represents the row of `V` corresponding to the group for `xfun`. When `vfun` is a matrix, it has the same number of rows as `xfun`. Each row of `vfun` corresponds to a row of `xfun`. This value is the default when you do not specify `Vectorization`. `"SingleGroup"` — `vfun` is a vector representing the row in `V` that corresponds to the group for `xfun`. `"Full"` — `vfun` is a vector or matrix. When `vfun` is a vector, the predictors apply to all observations in `xfun`. When `vfun` is a matrix, it has the same number of rows as `xfun`. Each row of `vfun` corresponds to a row of `xfun`. If `V` is an empty array, `nlmefitsa` calls the model function with only two inputs.
`yfit`	Fitted response values. `yfit` contains an element for each row in `xfun`.

For improved performance, allow modelfun to accept predictor variables for multiple observations and specify the Vectorization name-value argument as "SingleGroup" or "Full".

Example: model = @(phi,xfun,vfun)(phi(:,1).*exp(vfun(:,1).*xfun(:,1))+log(phi(:,2).*xfun(:,2)./vfun(:,1))); beta = nlmefitsa(X,y,group,V,model,[1 1 1]);

Data Types: function_handle

`beta0` — Initial fixed-effects values
numeric vector | numeric matrix

Initial fixed-effects values, specified as a numeric vector or matrix.

If beta0 is a vector, its size must be equal to the number of fixed effects.
If beta0 is a matrix, it must have as many rows as the number of model coefficients. For each column in beta0, nlmefitsa repeats the estimation of the fixed effects using the column entries as the initial fixed-effects values.

Example: [1 1 1 1]

Data Types: single | double

Name-Value Arguments

expand all

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Example: nlmefitsa(X,y,group,V,model,[1 1 1],Replicates=2,Vectorization="SingleGroup") specifies to estimate the model function parameters twice using the same starting point, and to pass the sample data to the model function in batches corresponding to each group.

Fixed and Random Effects

expand all

`FEParamsSelect` — Fixed-effects coefficients
numeric vector | logical vector

Fixed-effects coefficients, specified as a numeric or logical vector. This argument indicates which coefficients in the model function input argument phi contain a fixed effect. For more information about the model function, see modelfun. By default, all the elements of phi contain fixed effects.

When you specify FEParamsSelect, you cannot specify FEConstDesign or FEGroupDesign.

Example: FEParamsSelect=[1 0 1 1]

Data Types: single | double | logical

`FEConstDesign` — Fixed-effects design matrix
design matrix

Fixed-effects design matrix, specified as a p-by-f design matrix. p is the number of coefficients in the model function input argument phi, and f is the number of fixed effects. For more information, see modelfun and Mixed-Effects Model Hierarchy.

When you specify FEConstDesign, you cannot specify FEParamsSelect or FEGroupDesign.

Example: FEConstDesign=[0.1 0.5; 3 0.7; 0.8 2; 4 8] specifies a design matrix for four model parameters and two fixed effects.

Data Types: single | double

`FEGroupDesign` — Fixed-effects group design matrices
p-by-f-by-m array of design matrices

Fixed-effects group design matrices, specified as a p-by-f-by-m array of design matrices. p is the number of coefficients in the model function input argument PHI, f is the number of fixed effects, and m is the number of groups specified by group. FEGroupDesign specifies a different p-by-f design matrix for each of the m groups.

When you specify FEGroupDesign, you cannot specify FEConstDesign or FEParamsSelect.

Example: A(:,:,1) = [0.1 0.5; 3 0.7; 0.8 2]; A(:,:,2) = ones(3,2); FEGroupDesign=A specifies a design matrix for three model parameters, two fixed effects, and two groups.

Data Types: single | double

`REParamsSelect` — Random-effects coefficients
numeric vector | logical vector

Random-effects coefficients, specified as a numeric or logical vector. This argument indicates which coefficients in the model function input argument PHI contain a random effect. For more information about the model function, see modelfun. By default, all the elements of PHI contain random effects.

When you specify REParamsSelect, you cannot specify REConstDesign.

Example: REParamsSelect=[0 0 1 1]

Data Types: single | double | logical

`REConstDesign` — Random-effects design matrix
numeric vector | logical vector

Random-effects design matrix, specified as a p-by-r design matrix. p is the number of coefficients in the model function input argument phi, and r is the number of random effects. For more information, see modelfun and Mixed-Effects Model Hierarchy.

When you specify REConstDesign, you cannot specify REParamsSelect.

Example: REConstDesign=[1 0 2 1; 3 1 0 0; 0 0 1 0] specifies a design matrix for three parameters and four random effects.

Data Types: single | double

Iterative Algorithm

expand all

`Vectorization` — Size of model function inputs
`"SinglePhi"` (default) | string scalar | character vector

Size of the model function inputs, specified as a string scalar or character vector with one of the following values. For more information about the model function inputs, see modelfun.

Value Description

Value	Description
`"SinglePhi"`	`phi` is a vector, `xfun` is a row vector or matrix, and `vfun` is a vector or matrix. If `xfun` is a matrix, the observations must come from a single group. The model function calculates `yfit` for a single set of model parameters at a time.
`"SingleGroup"`	`phi` is a vector or matrix, `xfun` is a matrix, and `vfun` is a vector. The observations in `xfun` must come from a single group. The model function calculates `yfit` for a single group at a time.
`"Full"`	`phi` is a vector or matrix, `xfun` is a matrix, and `vfun` is a vector or matrix. The observations in `xfun` can come from different groups. The model function calculates `yfit` for multiple observations that can come from different groups.

"SinglePhi"

phi is a vector, xfun is a row vector or matrix, and vfun is a vector or matrix. If xfun is a matrix, the observations must come from a single group.

The model function calculates yfit for a single set of model parameters at a time.

"SingleGroup"

phi is a vector or matrix, xfun is a matrix, and vfun is a vector. The observations in xfun must come from a single group.

The model function calculates yfit for a single group at a time.

"Full"

phi is a vector or matrix, xfun is a matrix, and vfun is a vector or matrix. The observations in xfun can come from different groups.

The model function calculates yfit for multiple observations that can come from different groups.

If V is an empty array, nlmefitsa calls the model function with two inputs.

Example: Vectorization="SingleGroup"

Data Types: string | char

`CovPattern` — Random-effects covariance matrix pattern
r-by-r identity matrix (default) | r-by-r numeric matrix | r-by-r logical matrix | 1-by-r numeric vector

Random-effects covariance matrix pattern, specified as an r-by-r identity, numeric, or logical matrix, or a 1-by-r vector of integers, where r is the number of random effects. nlmefitsa calculates variance and covariance estimates for the random effects specified by CovPattern. For more information, see psi.

When CovPattern is a matrix, each off-diagonal element corresponds to a pair of random effects, and each element along the diagonal corresponds to the variance for a single random effect. nlmefitsa calculates variance and covariance estimates for the random effects corresponding to the nonzero elements of CovPattern. Zero elements in CovPattern indicate that the variances and covariances are constrained to zero.
If CovPattern is not a row-column permutation of a block diagonal matrix, nlmefitsa automatically adds nonzero elements to produce such a pattern.
When CovPattern is a vector, elements with equal values define groups of random effects. In addition to variances, nlmefitsa calculates covariance estimates for pairs of random effects within groups, and constrains covariances across groups to zero.

The default value for CovPattern is an r-by-r identity matrix, which corresponds to uncorrelated random effects.

Example: CovPattern=ones(3)

Example: CovPattern=[1 1 2]

Data Types: single | double | logical

`Cov0` — Initial covariance values
r-by-r positive definite matrix

Initial covariance values, specified as an r-by-r positive definite matrix, where r is the number of random effects. A positive definite matrix is symmetric with all nonnegative eigenvalues. If you do not specify Cov0, nlmefitsa calculates Cov0 using beta0.

Example: Cov0=[4 1 2; 1 2 0; 2 0 3]

Data Types: single | double

`ComputeStdErrors` — Flag to calculate standard errors
`false` or `0` (default) | `true` or `1`

Flag to calculate standard errors, specified as a numeric or logical 0 (false) or 1 (true). nlmefitsa returns the standard errors in the stats output argument.

Example: ComputeStdErrors=true

Data Types: single | double | logical

`ErrorModel` — Model for error term
`"constant"` (default) | string scalar | character vector

Model for the error term, specified as a string scalar or character vector containing the model name. nlmefitsa stores the parameter values for the error model in the errorparam field of the stats output argument.

Model Name	Formula	`stats.errorparam`
`"constant"`	$y = f + a ε$	a
`"proportional"`	$y = f + b f ε$	b
`"combined"`	$y = f + (a + b f) ε$	[a,b]
`"exponential"`	$y = \exp (a ε)$	a

In the above table, $y$ is the response variable y, $f$ is the fitted model, and $a$ and $b$ are model parameters.

Example: ErrorModel="exponential"

Data Types: string | char

`ErrorParameters` — Initial parameters for error model
numeric scalar (default) | 1-by-2 numeric vector | 2-by-1 numeric vector

Initial parameters for the error model, specified as a numeric scalar, 1-by-2 numeric vector, or 2-by-1 numeric vector.

If ErrorModel is "combined", ErrorParameters must be a two-element vector. Otherwise, ErrorParameters must be a scalar.

Example: ErrorParameters=0.75

Example: ErrorParameters=[0.3,0.5]

Data Types: single | double

`LogLikMethod` — Loglikelihood approximation method
`"none"` (default) | string scalar | character vector

Loglikelihood approximation method, specified as a string scalar or character vector containing one of the following values.

Value	Description
`"none"`	Omit the loglikelihood approximation.
`"is"`	Importance sampling
`"gq"`	Gaussian quadrature
`"lin"`	Linearization

Example: LogLikMethod="lin"

Data Types: string | char

`NBurnIn` — Number of burn-in iterations
`5` (default) | nonnegative integer

Number of burn-in iterations for the Markov chain Monte Carlo (MCMC) algorithm, specified as a nonnegative integer. nlmefitsa uses the MCMC algorithm during fitting to update the expected value of the loglikelihood function. During the burn-in iterations, nlmefitsa does not recalculate the parameter estimates for the loglikelihood function.

For more information about the nlmefitsa fitting algorithm, see Algorithms.

Example: NBurnIn=10

Data Types: single | double

`NChains` — Number of chains
positive integer

Number of chains nlmefitsa simulates for each iteration of the fitting algorithm, specified as a positive integer. During an iteration, nlmefitsa calculates NChains coefficient vectors for each group. If you do not specify NChains, the function calculates a value that results in approximately 100 groups across all chains.

For more information about the nlmefitsa fitting algorithm, see Algorithms.

Example: NChains=5

Data Types: single | double

`NIterations` — Number of iterations
`[150 150 100]` (default) | positive integer scalar | 1-by-3 vector of positive integers | 3-by-1 vector of positive integers

Number of iterations for each phase of the fitting algorithm, specified as a positive integer scalar, a 1-by-3 vector of positive integers, or a 3-by-1 vector of positive integers. NIterations determines the number of iterations for the following:

Simulated annealing in the first step
Full step size calculation in the second step
Reduced step size calculation in the third step

When NIterations is a vector, each element contains the number of iterations for the step corresponding to its index. When NIterations is a scalar, nlmefitsa distributes the iterations across the three phases of the fitting algorithm in the same proportions as the default.

For more information about the nlmefitsa fitting algorithm, see Algorithms.

Example: NIterations=500

Data Types: single | double

`NMCMCIterations` — Number of MCMC iterations
`[2 2 2]` (default) | positive integer scalar | 1-by-3 vector of positive integers | 3-by-1 vector of positive integers

Number of Markov chain Monte Carlo (MCMC) iterations in each simulation step of the fitting algorithm, specified as a positive integer scalar, a 1-by-3 vector of positive integers, or a 3-by-1 vector of positive integers. NMCMCIterations determines the number of iterations for the following:

Full multivariate update in the first step
Single coordinate update in the second step
Multiple coordinate update in the third step

When NMCMCIterations is a vector, each element contains the number of iterations for the step corresponding to its index. When NMCMCIterations is a scalar, nlmefitsa expands the scalar to a three-element vector with elements equal to the scalar.

For more information about the nlmefitsa fitting algorithm, see Algorithms.

Example: NMCMCIterations=[5,2,4]

Data Types: single | double

`OptimFun` — Optimization function
`"fminsearch"` (default) | `"fminunc"` (requires Optimization Toolbox™)

Optimization function for maximizing the likelihood function, specified as "fminsearch" or "fminunc".

"fminsearch" — Uses a direct search method involving function evaluations only. For more information, see fminsearch.
"fminunc" — Uses gradient methods. This option, which requires Optimization Toolbox, is generally more efficient. For more information, see fminunc (Optimization Toolbox).

For more information about the nlmefitsa fitting algorithm, see Algorithms.

Example: OptimFun="fminunc"

Data Types: string | char

`Options` — Options
statistics options structure

Options, specified as a statistics options structure returned by the statset function. The structure contains the following fields.

Field	Description
`DerivStep`	Relative difference used in the finite difference gradient calculation, specified as a positive scalar, or as a vector of the same length as the `PHI` input argument for the model function. The default value is `eps^(1/3)`. For more information about the model function, see `modelfun`.
`Display`	Option to display algorithm information, specified as one of the following: `"off"` (default) — Display no information. `"final"` — Display information about the final iteration of the algorithm that maximizes the likelihood function. `"iter"` — Display information about each iteration of the algorithm that maximizes the likelihood function.
`FunValCheck`	Flag to check for invalid values returned by the model function, specified as one of the following. `"on"` (default) — Check for invalid values, including `NaN` or `Inf`. `"off"` — Do not check for invalid values.
`OutputFcn`	Output function, specified as a function handle or a cell array of function handles. `nlmefitsa` calls all output functions after each iteration of the fitting algorithm described in Algorithms, and stops iterating when `OutputFcn` returns a logical `1`. `OutputFcn` must have the form `stop = outputfcn(beta,status,state)`: `beta` — Current fixed effects, specified as a numeric vector. `status` — Fitting status, specified as a structure that includes the following fields: `iteration` — Current iteration for the fitting algorithm, specified as a nonnegative integer. `Psi` — Current random-effects covariance matrix, specified as a nonnegative numeric matrix. `status` includes additional fields not used by `OutputFcn`. `state` — Type of iteration, specified as one of the following: `'init'` — Initial iteration. `'iter'` — Intermediate iteration. `'done'` — `OutputFcn` finished iterating. `stop` — Flag to stop the fitting algorithm, specified as a numeric logical scalar. When `stop` is `1` for at least one output function, `nlmefitsa` stops the fitting algorithm. The default output function plots the progress of the fitting algorithm for estimating the fixed effects, and plots the random-effects variances.
`Streams`	Streams used by `nlmefitsa` when generating random numbers, specified as a `RandStream` object. The default is the current global stream.
`TolX`	Termination tolerance for the fixed- and random-effects estimates. The default is `1e-4`.

Example: Options=statset("nlmefitsa")

Data Types: struct

`ParamTransform` — Coefficient transformation function
vector of zeros (default) | integers between 0 and 3

Coefficient transformation function for the model function, specified as a vector of integers between 0 and 3, inclusive. ParamTransform must have the same number of elements as the model function input argument PHI. For more information about the model function, see modelfun.

The coefficients for the model function are given by

$\begin{array}{l} \hat{φ} = A β + B b \\ φ = g (\hat{φ}) \end{array}$

where $φ$ and $\hat{φ}$ are the original and transformed model coefficients, respectively. $A$ and $B$ are the design matrices for the fixed and random effects $β$ and $b$ , respectively. The value of ParamTransform indicates the form of the transformation function $g$ according to the following scheme.

Value	Function
`0`	$\hat{φ} = φ$
`1`	$\hat{φ} = \log (φ)$
`2`	$\hat{φ} = probit (φ)$
`3`	$\hat{φ} = logit (φ)$

For more information about the fixed-effects design matrix, see FEConstDesign and FEGroupDesign. For more information about the random-effects design matrix, see REConstDesign.

Example: ParamTransform=[0 1 0]

Data Types: single | double

`Replicates` — Number of estimates
positive integer

Number of estimates for the fixed and random effects, specified as a positive integer.

When beta0 is a vector, nlmefitsa estimates the fixed and random effects Replicates number of times using the initial values in beta0. In this case, the default value for Replicates is 1.
When beta0 is a matrix, nlmefitsa estimates the fixed- and random-effects for each column in beta0. In this case, Replicates must be equal to the number of columns in beta0.

Example: Replicates=3

Data Types: single | double

Output Arguments

collapse all

`beta` — Fixed-effects coefficient estimates
numeric vector | numeric matrix

Fixed-effects coefficient estimates, returned as a numeric vector or an f-by-k matrix. f is the number of fixed effects and k is the number of estimates specified by Replicates.

Data Types: single | double

`psi` — Random-effects covariance matrices
numeric matrix | array of numeric matrices

Random-effects covariance matrices, returned as an r-by-r numeric matrix or an r-by-r-by-k array of numeric matrices. r is the number of random effects and k is the number of estimates specified by Replicates.

Data Types: single | double

`stats` — Fitting information
structure

Fitting information, returned as a structure with the following fields.

Field	Description
`logl`	Maximized loglikelihood for the fitted model. This field is empty if the `LogLikMethod` name-value argument is `"none"` during fitting.
`rmse`	Square root of the estimated error variance. `nlmefitsa` calculates `rmse` on the log scale when the `ErrorModel` name-value argument is `"exponential"` during fitting.
`errorparam`	Estimated parameters of the error model specified by the `ErrorModel` name-value argument during fitting.
`aic`	Akaike information criterion (AIC). This field is empty if the `LogLikMethod` is name-value argument is `"none"` during fitting. `nlmefitsa` calculates the AIC as `aic = -2logl+2numParam`, where `logl` is the maximized loglikelihood for the fitted model, and `numParam` is the number of fitting parameters. The fitting parameters include the degrees of freedom for the random-effects covariance matrix, the number of fixed effects, and the number of parameters in the error model.
`bic`	Bayesian information criterion (BIC). `nlmefitsa` calculates the BIC as `bic = -2logl+log(m)numParam`, where `logl` is the maximized loglikelihood for the fitted model, `m` is the number of groups specified by the `group` input argument, and `numParam` is the number of fitting parameters. The fitting parameters include the degrees of freedom for the random-effects covariance matrix, the number of fixed effects, and the number of parameters in the error model. This field is empty if the `LogLikMethod` name-value argument is `"none"` during fitting.
`sebeta`	Standard errors for the fixed-effects estimates in `beta`. This field is empty if the `ComputeStdErrors` name-value argument is `false` during fitting.
`covb`	Estimated covariance matrix for the fixed-effects estimates in `beta`. This field is empty if the `ComputeStdErrors` name-value argument is `false` during fitting.
`dfe`	Error degrees of freedom
`pres`	Population residuals, given by `y` – `y_pop`, where `y_pop` are the predicted response values for the population. `nlmefitsa` calculates `y_pop` by passing `X` to the fitted model with the random effects set to zero.
`ires`	Individual residuals, given by `y` – `y_ind`, where `y_ind` are the predicted response values for the individual observations. `nlmefitsa` calculates `y_ind` by passing `X` to the fitted model.
`pwres`	Population weighted residuals
`cwres`	Conditional weighted residuals
`iwres`	Individual weighted residuals

Data Types: struct

`b` — Random-effects estimates
numeric vector | numeric matrix | array of numeric matrices

Random-effects estimates, returned as an r-by-m numeric matrix or an r-by-m-by-k array of numeric matrices. r is the number of random effects, m is the number of groups specified by group, and k is the number of replicates specified by Replicates.

Data Types: single | double

Algorithms

By default, nlmefitsa fits a model where each model coefficient is the sum of a fixed effect and a random effect. To fit a model with a different number of fixed and random effects, use the name-value arguments in the Fixed and Random Effects section.

To estimate the parameters of a nonlinear mixed-effects model, nlmefitsa uses an iterative algorithm that includes a random Monte Carlo simulation. To specify parameters for the iterative algorithm, use the name-value arguments in the Iterative Algorithm section. Examine the plot of simulation results to ensure that the simulation has converged. You can also compare the results of calling nlmefitsa multiple times with different starting values, or use the Replicates name-value argument to perform multiple simulations.

During fitting, nlmefitsa calculates maximum likelihood estimates, which are parameter values that maximize a likelihood function. nlmefitsa uses the following likelihood function:

$p (y | β, σ^{2}, \sum) = \int p (y | β, b, σ^{2}) p (b | \sum) d b,$

where

$y$ is the response data.
$β$ is the vector of population coefficients.
$σ^{2}$ is the residual variance.
$\sum$ is the covariance matrix for the random effects.
$b$ is the set of unobserved random effects.

Each pdf $p$ on the right side of the formula is a normal (Gaussian) likelihood function that might depend on covariates.

nlmefitsa uses the stochastic approximation expectation-maximization (SAEM) algorithm to calculate the maximum likelihood estimates [1]. The SAEM algorithm takes these steps:

Simulation — Use a Markov chain Monte Carlo algorithm together with the posterior density function $p (b | y, \sum)$ and the current parameter estimates to simulate values for $b$ . Calculate the loglikelihood for each simulated $b$ , and then take the average of the loglikelihoods.
Stochastic approximation — Update the expected value of the loglikelihood function by moving its current value toward the average calculated in the previous step.
Maximization — Choose new estimates for $β$ , $σ^{2}$ , and $\sum$ that maximize the loglikelihood function given the simulated values of the random effects.

References

[1] Delyon, B., M. Lavielle, and E. Moulines. "Convergence of a Stochastic Approximation Version of the EM Algorithm." Annals of Statistics, 27, 94–128, 1999.

[2] Mentré, F., and M. Lavielle. "Stochastic EM Algorithms in Population PKPD Analyses." American Conference on Pharmacometrics, 2008.

Version History

Introduced in R2010a

nlmefitsa

Syntax

Description

Examples

Specify Group Predictors

Specify Output Function

Transform Parameters and Plot the Fitted Model

Input Arguments

X — Predictor variables numeric matrix

y — Response variable numeric vector

group — Grouping variable numeric vector | categorical vector | string vector | cell array of character vectors

V — Group predictor variables numeric matrix | cell array

modelfun — Model function function handle

beta0 — Initial fixed-effects values numeric vector | numeric matrix

Name-Value Arguments

Fixed and Random Effects

FEParamsSelect — Fixed-effects coefficients numeric vector | logical vector

FEConstDesign — Fixed-effects design matrix design matrix

FEGroupDesign — Fixed-effects group design matrices p-by-f-by-m array of design matrices

REParamsSelect — Random-effects coefficients numeric vector | logical vector

REConstDesign — Random-effects design matrix numeric vector | logical vector

Iterative Algorithm

Vectorization — Size of model function inputs "SinglePhi" (default) | string scalar | character vector

CovPattern — Random-effects covariance matrix pattern r-by-r identity matrix (default) | r-by-r numeric matrix | r-by-r logical matrix | 1-by-r numeric vector

Cov0 — Initial covariance values r-by-r positive definite matrix

ComputeStdErrors — Flag to calculate standard errors false or 0 (default) | true or 1

ErrorModel — Model for error term "constant" (default) | string scalar | character vector

ErrorParameters — Initial parameters for error model numeric scalar (default) | 1-by-2 numeric vector | 2-by-1 numeric vector

LogLikMethod — Loglikelihood approximation method "none" (default) | string scalar | character vector

NBurnIn — Number of burn-in iterations 5 (default) | nonnegative integer

NChains — Number of chains positive integer

NIterations — Number of iterations [150 150 100] (default) | positive integer scalar | 1-by-3 vector of positive integers | 3-by-1 vector of positive integers

NMCMCIterations — Number of MCMC iterations [2 2 2] (default) | positive integer scalar | 1-by-3 vector of positive integers | 3-by-1 vector of positive integers

OptimFun — Optimization function "fminsearch" (default) | "fminunc" (requires Optimization Toolbox™)

Options — Options statistics options structure

ParamTransform — Coefficient transformation function vector of zeros (default) | integers between 0 and 3

Replicates — Number of estimates positive integer

Output Arguments

beta — Fixed-effects coefficient estimates numeric vector | numeric matrix

psi — Random-effects covariance matrices numeric matrix | array of numeric matrices

stats — Fitting information structure

b — Random-effects estimates numeric vector | numeric matrix | array of numeric matrices

Algorithms

References

Version History

See Also

Topics

`X` — Predictor variables
numeric matrix

`y` — Response variable
numeric vector

`group` — Grouping variable
numeric vector | categorical vector | string vector | cell array of character vectors

`V` — Group predictor variables
numeric matrix | cell array

`modelfun` — Model function
function handle

`beta0` — Initial fixed-effects values
numeric vector | numeric matrix

`FEParamsSelect` — Fixed-effects coefficients
numeric vector | logical vector

`FEConstDesign` — Fixed-effects design matrix
design matrix

`FEGroupDesign` — Fixed-effects group design matrices
p-by-f-by-m array of design matrices

`REParamsSelect` — Random-effects coefficients
numeric vector | logical vector

`REConstDesign` — Random-effects design matrix
numeric vector | logical vector

`Vectorization` — Size of model function inputs
`"SinglePhi"` (default) | string scalar | character vector

`CovPattern` — Random-effects covariance matrix pattern
r-by-r identity matrix (default) | r-by-r numeric matrix | r-by-r logical matrix | 1-by-r numeric vector

`Cov0` — Initial covariance values
r-by-r positive definite matrix

`ComputeStdErrors` — Flag to calculate standard errors
`false` or `0` (default) | `true` or `1`

`ErrorModel` — Model for error term
`"constant"` (default) | string scalar | character vector

`ErrorParameters` — Initial parameters for error model
numeric scalar (default) | 1-by-2 numeric vector | 2-by-1 numeric vector

`LogLikMethod` — Loglikelihood approximation method
`"none"` (default) | string scalar | character vector

`NBurnIn` — Number of burn-in iterations
`5` (default) | nonnegative integer

`NChains` — Number of chains
positive integer

`NIterations` — Number of iterations
`[150 150 100]` (default) | positive integer scalar | 1-by-3 vector of positive integers | 3-by-1 vector of positive integers

`NMCMCIterations` — Number of MCMC iterations
`[2 2 2]` (default) | positive integer scalar | 1-by-3 vector of positive integers | 3-by-1 vector of positive integers

`OptimFun` — Optimization function
`"fminsearch"` (default) | `"fminunc"` (requires Optimization Toolbox™)

`Options` — Options
statistics options structure

`ParamTransform` — Coefficient transformation function
vector of zeros (default) | integers between 0 and 3

`Replicates` — Number of estimates
positive integer

`beta` — Fixed-effects coefficient estimates
numeric vector | numeric matrix

`psi` — Random-effects covariance matrices
numeric matrix | array of numeric matrices

`stats` — Fitting information
structure

`b` — Random-effects estimates
numeric vector | numeric matrix | array of numeric matrices