nlmefit

Fit nonlinear mixed-effects estimation

Syntax

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

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

[beta,psi] = nlmefit(___)

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

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

Description

beta = nlmefit(X,y,group,V,modelfun,beta0) fits the nonlinear mixed-effects regression model modelfun to the data in X and Y, and returns the fixed-effects estimates in beta

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

[beta,psi] = nlmefit(___) additionally returns a matrix of covariance estimates for the random effects.

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

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

Examples

Specify Group Predictors

Load some nonlinear sample data and display the 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
      ⋮

V

V = 2×1

     2
     3

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 = 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

$${\varphi }_1{X}_1\exp \left(\frac{{\varphi }_2{X}_2}{V}\right)+{\varphi }_3{X}_3$$.

Use the stochastic EM algorithm to fit model to the data in X, y, group, and V. Use 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-effect 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

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 $${\varphi }_1{X}_1\exp \left(\frac{{\varphi }_2{X}_2}{V}\right)+{\varphi }_3{X}_3$$.

Define an output function for nlmefit. 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("nlmefit");
opts = statset(default_opts,OutputFcn=@outputFunction);

opts is a statistics options structure that contains options for the stocastic 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)

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]

nlmefit 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 Fitted Model

Load the indomethacin data set.

load indomethacin

The variables time, concentration, and subject contain timeseries 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 data.

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={\varphi }_1\exp (-{\varphi }_{2}t)+{\varphi }_3\exp (-{\varphi }_{4}t)$$,

where $$c$$ is the concentration of indomethacin, $${\varphi }_i$$ for $$i=1,2,3,4$$ are coefficients, and $$t$$ is time.The function does not contain group-specific predictor variables because they are not present in the formula.

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] = nlmefit(time,concentration, ...
   subject,[],model,phi0,ParamTransform=xform)

beta = 4×1

    0.4606
   -1.3459
    2.8277
    0.7729

psi = 4×4

    0.0124         0         0         0
         0    0.0000         0         0
         0         0    0.3264         0
         0         0         0    0.0250

stats = struct with fields:
           dfe: 57
          logl: 54.5882
           mse: 0.0066
          rmse: 0.0787
    errorparam: 0.0815
           aic: -91.1765
           bic: -93.0506
          covb: [4×4 double]
        sebeta: [0.1092 0.2244 0.2558 0.1066]
          ires: [66×1 double]
          pres: [66×1 double]
         iwres: [66×1 double]
         pwres: [66×1 double]
         cwres: [66×1 double]

b = 4×6

   -0.1111    0.0871    0.0670    0.0115   -0.1315    0.0769
         0         0         0         0         0         0
   -0.7396   -0.0704    0.8004   -0.5654    0.4127    0.1624
    0.0279    0.0287    0.0040   -0.2315    0.1984   -0.0276

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

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

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, nlmefit 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 = nlmefit(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, nlmefit 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

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: nlmefit(X,y,group,v,model,[1 1 1],Replicates=2,Vectorization="SingleGroup") specifies for nlmefit 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

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, FEGroupDesign, or FEObsDesign.

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, FEGroupDesign, or FEObsDesign.

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, FEParamsSelect, or FEObsDesign.

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

FEObsDesign — Fixed-effects observation design matrices

p-by-f-by-n array of design matrices

Fixed-effects observation design matrices, specified as p-by-f-by-n 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 n is the number of observations specified by X and y. FEObsDesign specifies a different p-by-f design matrix for each of the n observations.

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

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

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, REGroupDesign, or REObsDesign.

Example: FEParamsSelect=[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, REGroupDesign, or REObsDesign.

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

REGroupDesign — Random-effects group design matrices

p-by-r-by-m array of design matrices

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

When you specify REGroupDesign, you cannot specify REConstDesign, REParamsSelect, or REObsDesign.

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

Data Types: single | double

REObsDesign — Random-effects observation design matrices

p-by-r-by-n array of design matrices

Random-effects observation design matrices, specified as p-by-r-by-n array of design matrices. p is the number of coefficients in the model function input argument PHI, r is the number of random effects, and n is the number of observations specified by X and y. REObsDesign specifies a different p-by-r design matrix for each of the n observations.

When you specify REObsDesign, you cannot specify REGroupDesign, REConstDesign, or REParamsSelect.

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

Data Types: single | double

Iterative Algorithm

ApproximationType — Method for approximating likelihood

"LME" (default) | "RELME" | "FO" | "FOCE"

Method for approximating the likelihood of the model, specified as one of the following:

• "LME" — Likelihood for the linear mixed-effects model at the current conditional estimates of beta and b

• "RELME" — Restricted likelihood for the linear mixed-effects model at the current conditional estimates of beta and b

• "FO" — First-order Laplacian approximation without random effects

• "FOCE" — First-order Laplacian approximation at the conditional estimates of b

The conditional estimates for beta and b are calculated using the current estimates of the parameters in the covariance matrix for b.

Example: ApproximationType="FO"

Data Types: string | char

CovParameterization — Parametrization for scaled covariance matrix

"logm" (default) | "chol"

Parametrization for the scaled covariance matrix, specified as one of the following:

• "logm" — Matrix logarithm

• "chol" — Cholesky factorization

Example: CovParameterization="chol"

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. nlmefit 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. nlmefit 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, nlmefit 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, nlmefit 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

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. nlmefit 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\epsilon$$	a
"proportional"	$$y=f+bf\epsilon$$	b
"combined"	$$y=f+(a+bf)\epsilon$$	[a,b]
"exponential"	$$y=\exp (a\epsilon )$$	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

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 nlmefit 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. nlmefit 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, nlmefit 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 nlmefit 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("nlmefit")

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}\widehat{\phi }=A\beta +Bb\\ \phi =g(\widehat{\phi })\end{array}$$

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

Value	Function
0	$$\widehat{\phi }=\phi$$
1	$$\widehat{\phi }=\log (\phi )$$
2	$$\widehat{\phi }=\text{probit}(\phi )$$
3	$$\widehat{\phi }=\text{logit}(\phi )$$

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

RefineBeta0 — Flag to refine initial fixed-effects values

"on" (default) | "off"

Flag to refine initial fixed-effects values, specified as "on" or "off". When RefineBeta0 is "on", nlmefit first fits the model without random effects, then replaces beta0 with the fitted fixed-effects.

Example: RefineBeta0="on"

Data Types: string | char

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
"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, nlmefit calls the model function with two inputs.

Example: Vectorization="SingleGroup"

Data Types: string | char

Output Arguments

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. nlmefit 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. nlmefit calculates the AIC as aic = -2*logl+2*numParam, 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). nlmefit calculates the BIC as bic = -2*logl+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 – ypop, where ypop are the predicted response values for the population. nlmefit calculates ypop by passing X to the fitted model with the random effects set to zero.
ires	Individual residuals, given by y – yind, where yind are the predicted response values for the individual observations. nlmefit calculates yind 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

nlmefit fits the model by maximizing an estimate of the marginal likelihood with the random effects integrated out. During fitting, the function makes the following assumptions:

• The random effects come from a multivariate normal distribution and are independent between groups.

• The observation errors are independent, come from the same normal distribution, and are independent of the random effects. To change this default assumption, specify the ErrorModel name-value argument.

References

[1] Lindstrom, M. J., and D. M. Bates. “Nonlinear mixed-effects models for repeated measures data.” Biometrics. Vol. 46, 1990, pp. 673–687.

[2] Davidian, M., and D. M. Giltinan. Nonlinear Models for Repeated Measurements Data. New York: Chapman & Hall, 1995.

[3] Pinheiro, J. C., and D. M. Bates. “Approximations to the log-likelihood function in the nonlinear mixed-effects model.” Journal of Computational and Graphical Statistics. Vol. 4, 1995, pp. 12–35.

[4] Demidenko, E. Mixed Models: Theory and Applications. Hoboken, NJ: John Wiley & Sons, Inc., 2004.

Version History

Introduced in R2008b