nlmefit - Nonlinear mixed-effects estimation

Syntax

Description

beta = nlmefit(X,y,group,V,fun,beta0) fits a nonlinear mixed-effects regression model and returns estimates of the fixed effects in beta. By default, nlmefit fits a model in which each parameter is the sum of a fixed and a random effect, and the random effects are uncorrelated (their covariance matrix is diagonal).

X is an n-by-h matrix of n observations on h predictors.

y is an n-by-1 vector of responses.

group is a grouping variable indicating m groups in the observations. group is a categorical variable, a numeric vector, a character matrix with rows for group names, or a cell array of strings. For more information on grouping variables, see Grouped Data.

V is an m-by-g matrix or cell array of g group-specific predictors. These are predictors that take the same value for all observations in a group. The rows of V are assigned to groups using grp2idx, according to the order specified by grp2idx(group). Use a cell array for V if group predictors vary in size across groups. Use [] for V if there are no group-specific predictors.

fun is a handle to a function that accepts predictor values and model parameters and returns fitted values. fun has the form

yfit = modelfun(PHI,XFUN,VFUN)

The arguments are:

• PHI — A 1-by-p vector of model parameters.

• XFUN — A k-by-h array of predictors, where:

  • k = 1 if XFUN is a single row of X.

  • k = n_i if XFUN contains the rows of X for a single group of size n_i.

  • k = n if XFUN contains all rows of X.

• VFUN — Group-specific predictors given by one of:

  • A 1-by-g vector corresponding to a single group and a single row of V.

  • An n-by-g array, where the jth row is V(I,:) if the jth observation is in group I.

  If V is empty, nlmefit calls modelfun with only two inputs.

• yfit — A k-by-1 vector of fitted values

When either PHI or VFUN contains a single row, it corresponds to all rows in the other two input arguments.

Note   If modelfun can compute yfit for more than one vector of model parameters per call, use the 'Vectorization' parameter (described later) for improved performance.

beta0 is a q-by-1 vector with initial estimates for q fixed effects. By default, q is the number of model parameters p.

nlmefit fits the model by maximizing an approximation to the marginal likelihood with random effects integrated out, assuming that:

• Random effects are multivariate normally distributed and independent between groups.

• Observation errors are independent, identically normally distributed, and independent of the random effects.

[beta,PSI] = nlmefit(X,y,group,V,fun,beta0) also returns PSI, an r-by-r estimated covariance matrix for the random effects. By default, r is equal to the number of model parameters p.

[beta,PSI,stats] = nlmefit(X,y,group,V,fun,beta0) also returns stats, a structure with fields:

• dfe — The error degrees of freedom for the model

• logl — The maximized log-likelihood for the fitted model

• rmse — The square root of the estimated error variance (computed on the log scale for the exponential error model)

• errorparam — The estimated parameters of the error variance model

• aic — The Akaike information criterion, calculated as aic = -2 * logl + 2 * numParam, where numParam is the number of fitting parameters, including the degree of freedom for covariance matrix of the random effects, the number of fixed effects and the number of parameters of the error model, and logl is a field in the stats structure

• bic — The Bayesian information criterion, calculated as bic = –2*logl + log(M) * numParam

  • M is the number of groups.

  • numParam and logl are defined as in aic.

  Note that some literature suggests that the computation of bic should be , bic = –2*logl + log(N) * numParam, where N is the number of observations.

• covb — The estimated covariance matrix of the parameter estimates

• sebeta — The standard errors for beta

• ires — The population residuals (y-y_population), where y_population is the individual predicted values

• pres — The population residuals (y-y_population), where y_population is the population predicted values

• iwres — The individual weighted residuals

• pwres — The population weighted residuals

• cwres — The conditional weighted residuals

[beta,PSI,stats,B] = nlmefit(X,y,group,V,fun,beta0) also returns B, an r-by-m matrix of estimated random effects for the m groups. By default, r is equal to the number of model parameters p.

[beta,PSI,stats,B] = nlmefit(X,y,group,V,fun,beta0,'Name',value) specifies one or more optional parameter name/value pairs. Specify Name inside single quotes.

Use the following parameters to fit a model different from the default. (The default model is obtained by setting both FEConstDesign and REConstDesign to eye(p), or by setting both FEParamsSelect and REParamsSelect to 1:p.) Use at most one parameter with an 'FE' prefix and one parameter with an 'RE' prefix. The nlmefit function requires you to specify at least one fixed effect and one random effect.

Parameter	Value
FEParamsSelect	A vector specifying which elements of the parameter vector PHI include a fixed effect, given as a numeric vector of indices from 1 to p or as a 1-by-p logical vector. If q is the specified number of elements, then the model includes q fixed effects.
FEConstDesign	A p-by-q design matrix ADESIGN, where ADESIGN*beta are the fixed components of the p elements of PHI.
FEGroupDesign	A p-by-q-by-m array specifying a different p-by-q fixed-effects design matrix for each of the m groups.
FEObsDesign	A p-by-q-by-n array specifying a different p-by-q fixed-effects design matrix for each of the n observations.
REParamsSelect	A vector specifying which elements of the parameter vector PHI include a random effect, given as a numeric vector of indices from 1 to p or as a 1-by-p logical vector. The model includes r random effects, where r is the specified number of elements.
REConstDesign	A p-by-r design matrix BDESIGN, where BDESIGN*B are the random components of the p elements of PHI.
REGroupDesign	A p-by-r-by-m array specifying a different p-by-r random-effects design matrix for each of m groups.
REObsDesign	A p-by-r-by-n array specifying a different p-by-r random-effects design matrix for each of n observations.

Use the following parameters to control the iterative algorithm for maximizing the likelihood:

Parameter	Value
RefineBeta0	Determines whether nlmefit makes an initial refinement of beta0 by first fitting modelfun without random effects and replacing beta0 with beta. Choices are 'on' and 'off'. The default value is 'on'.
ErrorModel	A string specifying the form of the error term. Default is 'constant'. Each model defines the error using a standard normal (Gaussian) variable e, the function value f, and one or two parameters a and b. Choices are: 'constant': y = f + a*e 'proportional': y = f + b*f*e 'combined': y = f + (a+b*f)*e 'exponential': y = f*exp(a*e), or equivalently log(y) = log(f) + a*e If this parameter is given, the output stats.errorparam field has the value a for 'constant' and 'exponential' b for 'proportional' [a b] for 'combined'
ApproximationType	The method used to approximate the likelihood of the model. Choices are: 'LME' — Use the likelihood for the linear mixed-effects model at the current conditional estimates of beta and B. This is the default. 'RELME' — Use the 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.
Vectorization	Indicates acceptable sizes for the PHI, XFUN, and VFUN input arguments to modelfun. Choices are: 'SinglePhi' — modelfun can only accept a single set of model parameters at a time, so PHI must be a single row vector in each call. nlmefit calls modelfun in a loop, if necessary, with a single PHI vector and with XFUN containing rows for a single observation or group at a time. VFUN may be a single row that applies to all rows of XFUN, or a matrix with rows corresponding to rows in XFUN. This is the default. 'SingleGroup' — modelfun can only accept inputs corresponding to a single group in the data, so XFUN must contain rows of X from a single group in each call. Depending on the model, PHI is a single row that applies to the entire group or a matrix with one row for each observation. VFUN is a single row. 'Full' — modelfun can accept inputs for multiple parameter vectors and multiple groups in the data. Either PHI or VFUN may be a single row that applies to all rows of XFUN or a matrix with rows corresponding to rows in XFUN. This option can improve performance by reducing the number of calls to modelfun, but may require modelfun to perform singleton expansion on PHI or V.
CovParameterization	Specifies the parameterization used internally for the scaled covariance matrix. Choices are 'chol' for the Cholesky factorization or 'logm' the matrix logarithm. The default is 'logm'.
CovPattern	Specifies an r-by-r logical or numeric matrix P that defines the pattern of the random-effects covariance matrix PSI. nlmefit estimates the variances along the diagonal of PSI and the covariances specified by nonzeroes in the off-diagonal elements of P. Covariances corresponding to zero off-diagonal elements in P are constrained to be zero. If P does not specify a row-column permutation of a block diagonal matrix, nlmefit adds nonzero elements to P as needed. The default value of P is eye(r), corresponding to uncorrelated random effects. Alternatively, P may be a 1-by-r vector containing values in 1:r, with equal values specifying groups of random effects. In this case, nlmefit estimates covariances only within groups, and constrains covariances across groups to be zero.
ParamTransform	A vector of P values specifying a transformation function f() for each of the P parameters: XB = ADESIGN*BETA + BDESIGN*B PHI = f(XB)Each element of the vector must be one of the following integer codes specifying the transformation for the corresponding value of PHI: 0: PHI = XB (default for all parameters) 1: log(PHI) = XB 2: probit(PHI) = XB 3: logit(PHI) = XB
Options	A structure of the form returned by statset. nlmefit uses the following statset parameters: 'DerivStep' — Relative difference used in finite difference gradient calculation. May be a scalar, or a vector whose length is the number of model parameters p. The default is eps^(1/3). 'Display' — Level of iterative display during estimation. Choices are: 'off' (default) — Displays no information 'final' — Displays information after the final iteration 'iter' — Displays information at each iteration 'FunValCheck' — Check for invalid values, such as NaN or Inf, from modelfun. Choices are 'on' and 'off'. The default is 'on'. 'MaxIter' — Maximum number of iterations allowed. The default is 200. 'OutputFcn' — Function handle specified using @, a cell array with function handles or an empty array (default). The solver calls all output functions after each iteration. 'TolFun' — Termination tolerance on the log-likelihood function. The default is 1e-4. 'TolX' — Termination tolerance on the estimated fixed and random effects. The default is 1e-4.
OptimFun	Specifies the optimization function used in maximizing the likelihood. Choices are 'fminsearch' to use fminsearch or 'fminunc' to use fminunc. The default is 'fminsearch'. You can specify 'fminunc' only if Optimization Toolbox software is installed.

Examples

Display data on the growth of five orange trees:

CIRC = [30 58 87 115 120 142 145;
        33 69 111 156 172 203 203;
        30 51 75 108 115 139 140;
        32 62 112 167 179 209 214;
        30 49 81 125 142 174 177];
time = [118 484 664 1004 1231 1372 1582];

h = plot(time,CIRC','o','LineWidth',2);
xlabel('Time (days)')
ylabel('Circumference (mm)')
title('{\bf Orange Tree Growth}')
legend([repmat('Tree ',5,1),num2str((1:5)')],...
       'Location','NW')
grid on
hold on

Use an anonymous function to specify a logistic growth model:

model=@(PHI,t)(PHI(:,1))./(1+exp(-(t-PHI(:,2))./PHI(:,3)));

Fit the model using nlmefit with default settings (that is, assuming each parameter is the sum of a fixed and a random effect, with no correlation among the random effects):

TIME = repmat(time,5,1);
NUMS = repmat((1:5)',size(time));

beta0 = [100 100 100];
[beta1,PSI1,stats1] = nlmefit(TIME(:),CIRC(:),NUMS(:),...
                              [],model,beta0)
beta1 =
  191.3189
  723.7608
  346.2517

PSI1 =
  962.1534         0         0
         0    0.0000         0
         0         0  297.9881

stats1 = 
           dfe: 28
          logl: -131.5457
           mse: 59.7882
          rmse: 7.9016
    errorparam: 7.7323
           aic: 277.0913
           bic: 274.3574
          covb: [3x3 double]
        sebeta: [15.2249 33.1579 26.8235]
          ires: [35x1 double]
          pres: [35x1 double]
         iwres: [35x1 double]
         pwres: [35x1 double]
         cwres: [35x1 double]

The negligible variance of the second random effect, PSI1(2,2), suggests that it can be removed to simplify the model:

[beta2,PSI2,stats2,b2] = nlmefit(TIME(:),CIRC(:),...
    NUMS(:),[],model,beta0,'REParamsSelect',[1 3])

beta2 =
  191.3194
  723.7628
  346.2548

PSI2 =
  962.2114         0
         0  298.3989

stats2 = 
           dfe: 29
          logl: -131.5456
           mse: 59.7851
          rmse: 7.7640
    errorparam: 7.7321
           aic: 275.0913
           bic: 272.7479
          covb: [3x3 double]
        sebeta: [15.2252 33.1572 26.8246]
          ires: [35x1 double]
          pres: [35x1 double]
         iwres: [35x1 double]
         pwres: [35x1 double]
         cwres: [35x1 double]

b2 =
  -28.5250   31.6063  -36.5070   39.0735   -5.6479
   10.0097   -0.7638    6.0117   -9.4685   -5.7892

The log-likelihood logl is unaffected, and both the Akaike and Bayesian information criteria (aic and bic) are reduced, supporting the decision to drop the second random effect from the model.

Use the estimated fixed effects in beta2 and the estimated random effects for each tree in b2 to plot the model through the data:

PHI = repmat(beta2,1,5) + ...          % Fixed effects
      [b2(1,:);zeros(1,5);b2(2,:)];    % Random effects

colors = get(h,'Color');
tplot = 0:0.1:1600;
for I = 1:5
  fitted_model=@(t)(PHI(1,I))./(1+exp(-(t-PHI(2,I))./ ... 
PHI(3,I)));
  plot(tplot,fitted_model(tplot),'Color',colors{I}, ...
			'LineWidth',2)
end

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.

See Also

nlinfit | nlmefitsa | nlpredci

How To

• Mixed-Effects Models
• Grouped Data

© 1984-2012- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS