Fit nonlinear mixedeffects model with stochastic EM algorithm
[
BETA
,PSI
,STATS
,B
]
= nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0)
[BETA
,PSI
,STATS
,B
]
= nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0,'Name
',Value
)
[
fits a nonlinear
mixedeffects regression model and returns estimates of the fixed
effects in BETA
,PSI
,STATS
,B
]
= nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0)BETA
. By default, nlmefitsa
fits
a model where each model parameter is the sum of a corresponding fixed
and random effect, and the covariance matrix of the random effects
is diagonal, i.e., uncorrelated random effects.
The BETA
, PSI
, and other
values this function returns are the result of a random (Monte Carlo)
simulation designed to converge to the maximum likelihood estimates
of the parameters. Because the results are random, it is advisable
to examine the plot of simulation to results to be sure that the simulation
has converged. It may also be helpful to run the function multiple
times, using multiple starting values, or use the 'Replicates'
parameter
to perform multiple simulations.
[
accepts
one or more commaseparated parameter name/value pairs. Specify BETA
,PSI
,STATS
,B
]
= nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0,'Name
',Value
)Name
inside
single quotes.
Definitions:
In the following list of arguments, the following variable definitions apply:
n — number of observations
h — number of predictor variables
m — number of groups
g — number of groupspecific predictor variables
p — number of parameters
f — number of fixed effects

An nbyh matrix of n observations on h predictor variables. 

An nby1 vector of responses. 

A grouping variable indicating to which of m groups
each observation belongs. 

An mbyg matrix of g groupspecific
predictor variables for each of the m groups in
the data. These are predictor values that take on the same value for
all observations in a group. Rows of 

A handle to a function that accepts predictor values and model
parameters, and returns fitted values.


An fby1 vector with initial estimates for
the f fixed effects. By default, f is
equal to the number of model parameters p. 
By default, nlmefitsa
fits a model where
each model parameter is the sum of a corresponding fixed and random
effect. Use the following parameter name/value pairs to fit a model
with a different number of or dependence on fixed or random effects.
Use at most one parameter name with an 'FE'
prefix
and one parameter name with an 'RE'
prefix. Note
that some choices change the way nlmefitsa
calls MODELFUN
,
as described further below.

A vector specifying which elements of the model parameter vector 

A pbyf design matrix 

A pbyfbym array specifying a different pbyf fixed effects design matrix for each of the m groups. 

A vector specifying which elements of the model parameter vector 

A pbyr design matrix 
The default model is equivalent to setting both FEConstDesign
and REConstDesign
to eye(p)
,
or to setting both FEParamsSelect
and REParamsSelect
to
1:p.
Additional optional parameter name/value pairs control the iterative algorithm used to maximize the likelihood:

Specifies an rbyr logical
or numeric matrix Alternatively, specify 

Initial value for the covariance matrix 



A character vector specifying the form of the error term. Default
is
If this parameter is given, the output


A scalar or twoelement vector specifying starting values for
parameters of the error model. This specifies the a, b,
or [a b] values depending on
the 

Specifies the method for approximating the loglikelihood. Choices are:


Number of initial burnin iterations during which the parameter estimates are not recomputed. Default is 5. 

Number c of "chains" simulated. Default is 1. Setting c>1 causes c simulated coefficient vectors to be computed for each group during each iteration. Default depends on the data, and is chosen to provide about 100 groups across all chains. 

Number of iterations. This can be a scalar or a threeelement vector. Controls how many iterations are performed for each of three phases of the algorithm:
Default is 

Number of Markov Chain Monte Carlo (MCMC) iterations. This can be a scalar or a threeelement vector. Controls how many of three different types of MCMC updates are performed during each phase of the main iteration:
Default is 

Either 

A structure created by a call to


A vector of pvalues specifying a transformation
function XB = ADESIGN*BETA + BDESIGN*B PHI = f(XB) PHI :


Number 

Determines the possible sizes of the
The default for 

Estimates of the fixed effects 

An rbyr estimated covariance matrix for the random effects. By default, r is equal to the number of model parameters p. 

A structure with the following fields:

In order to estimate the parameters of a nonlinear mixed effects model, we would like to choose the parameter values that maximize a likelihood function. These values are called the maximum likelihood estimates. The likelihood function can be written in the form
$$p\left(y\beta ,{\sigma}^{2},\Sigma \right)={\displaystyle \int p\left(y\beta ,b,{\sigma}^{2}\right)p\left(b\Sigma \right)db}$$
y is the response data
β is the vector of population coefficients
σ^{2} is the residual variance
∑ is the covariance matrix for the random effects
b is the set of unobserved random effects
Each p() function on the righthandside is a normal (Gaussian) likelihood function that may depend on covariates.
Since the integral does not have a closed form, it is difficult to find parameters that maximize it. Delyon, Lavielle, and Moulines [1] proposed to find the maximum likelihood estimates using an ExpectationMaximization (EM) algorithm in which the E step is replaced by a stochastic procedure. They called their algorithm SAEM, for Stochastic Approximation EM. They demonstrated that this algorithm has desirable theoretical properties, including convergence under practical conditions and convergence to a local maximum of the likelihood function. Their proposal involves three steps:
Simulation: Generate simulated values of the random effects b from the posterior density p(bΣ) given the current parameter estimates.
Stochastic approximation: Update the expected value of the loglikelihood function by taking its value from the previous step, and moving part way toward the average value of the loglikelihood calculated from the simulated random effects.
Maximization step: Choose new parameter estimates to maximize the loglikelihood function given the simulated values of the random effects.
[1] Delyon, B., M. Lavielle, and E. Moulines, Convergence of a stochastic approximation version of the EM algorithm, Annals of Statistics, 27, 94128, 1999.
[2] Mentré, France, and Marc Lavielle, Stochastic EM algorithms in population PKPD analyses, 2008.