Main Content

nlmefit

Fit nonlinear mixed-effects estimation

    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

    example

    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.

    example

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

    example

    Examples

    collapse all

    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

    ϕ1X1exp(ϕ2X2V)+ϕ3X3.

    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])

    Figure contains 6 axes objects. Axes object 1 with title beta indexOf 1 baseline contains an object of type line. This object represents Rep 1. Axes object 2 with title beta indexOf 2 baseline contains an object of type line. This object represents Rep 1. Axes object 3 with title beta indexOf 3 baseline contains an object of type line. This object represents Rep 1. Axes object 4 with title Psi indexOf 11 baseline contains an object of type line. This object represents Rep 1. Axes object 5 with title Psi indexOf 22 baseline contains an object of type line. This object represents Rep 1. Axes object 6 with title Psi indexOf 33 baseline contains an object of type line. This object represents Rep 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.

    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 ϕ1X1exp(ϕ2X2V)+ϕ3X3.

    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)

    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]
    
    

    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.

    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=ϕ1exp(-ϕ2t)+ϕ3exp(-ϕ4t),

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

    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

    Figure contains an axes object. The axes object with title Indomethacin Elimination, xlabel Time (hours), ylabel Concentration (mcg/ml) contains 12 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 1, fitted curve 2, fitted curve 3, fitted curve 4, fitted curve 5, fitted curve 6.

    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

    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

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

    Data Types: single | double

    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

    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

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

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

    yfitFitted 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

    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

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

    expand all

    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

    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

    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

    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

    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

    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

    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

    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

    expand all

    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

    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

    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

    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 NameFormulastats.errorparam
    "constant"

    y=f+aε

    a
    "proportional"

    y=f+bfε

    b
    "combined"

    y=f+(a+bf)ε

    [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

    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, specified as a statistics options structure returned by the statset function. The structure contains the following fields.

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

    StreamsStreams used by nlmefit when generating random numbers, specified as a RandStream object. The default is the current global stream.
    TolXTermination tolerance for the fixed- and random-effects estimates. The default is 1e-4.

    Example: Options=statset("nlmefit")

    Data Types: struct

    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

    φ^=Aβ+Bbφ=g(φ^)

    where φ and φ^ 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.

    ValueFunction
    0

    φ^=φ

    1

    φ^=log(φ)

    2

    φ^=probit(φ)

    3

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

    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

    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.

    ValueDescription
    "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

    collapse all

    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

    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

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

    FieldDescription
    loglMaximized loglikelihood for the fitted model. This field is empty if the LogLikMethod name-value argument is "none" during fitting.
    rmseSquare root of the estimated error variance. nlmefit calculates rmse on the log scale when the ErrorModel name-value argument is "exponential" during fitting.
    errorparamEstimated parameters of the error model specified by the ErrorModel name-value argument during fitting.
    aicAkaike 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.

    sebetaStandard 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 yypop, 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 yyind, 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

    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