nlmefit - Nonlinear mixed-effects estimation

Syntax

beta = nlmefit(X,y,group,V,fun,beta0) [beta,PSI] = nlmefit(...) [beta,PSI,stats] = nlmefit(...) [beta,PSI,stats,B] = nlmefit(...) [...] = nlmefit(...,param1,val1,param2,val2,...)

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.

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 function 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(...) 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(...) also returns stats, a structure with fields:

logl — The maximized log-likelihood for the fitted model
mse — The estimated error variance for the fitted model
aic — The Akaike information criterion for the fitted model
bic — The Bayesian information criterion for the fitted model
sebeta — The standard errors for beta
dfe — The error degrees of freedom for the model

[beta,PSI,stats,B] = nlmefit(...) 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.

[...] = nlmefit(...,param1,val1,param2,val2,...) specifies additional parameter name/value pairs.

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 between 1 and 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 between 1 and 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'`.
`'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.
`'Options'`	A structure of the form returned by `statset`. `nlmefit` uses the following `statset` parameters: `'TolX'` — Termination tolerance on the estimated fixed and random effects. The default is `1e-4`. `'TolFun'` — Termination tolerance on the log-likelihood function. The default is `1e-4`. `'MaxIter'` — Maximum number of iterations allowed. The default is `200`. `'Display'` — Level of iterative display during estimation. Choices are: `'off'` — Displays no information. This is the default. `'iter'` — Displays iterative output to the command window. `'final'` — Displays the final output. `'FunValCheck'` — Check for invalid values, such as `NaN` or `Inf`, from `modelfun`. Choices are `'on'` and `'off'`. The default is `'on'`. `'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.
`'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 only specify `'fminunc'` 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.7609
  346.2518
PSI1 =
  962.1533         0         0
         0    0.0000         0
         0         0  297.9931
stats1 = 
      logl: -131.5457
       mse: 59.7881
       aic: 277.0913
       bic: 287.9788
    sebeta: NaN
       dfe: 28

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.3190
  723.7610
  346.2527
PSI2 =
  962.0491         0
         0  298.1869
stats2 = 
      logl: -131.5457
       mse: 59.7881
       aic: 275.0913
       bic: 284.4234
    sebeta: NaN
       dfe: 29
b2 =
  -28.5254   31.6061  -36.5071   39.0738   -5.6475
   10.0034   -0.7633    6.0080   -9.4630   -5.7853

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.

nlmefit - Nonlinear mixed-effects estimation

Syntax

Description

Examples

References

See Also

Recommended Products