## Integral Approximations

An important part of the marginal maximum likelihood method described
previously is the computation of the integral over the random effects.
The default method in PROC NLMIXED for computing this integral
is adaptive Gaussian quadrature as described in Pinheiro and
Bates (1995). Another approximation method is the first-order
method of Beal and Sheiner (1982, 1988).
A description of these two methods follows.
*Adaptive Gaussian Quadrature*

A quadrature method approximates a given integral by a weighted sum over
predefined abscissas for the random effects. A good approximation
can usually be obtained with an adequate number of quadrature points
as well as appropriate centering and scaling of the abscissas.
Adaptive Gaussian quadrature for the integral over *u*_{i} centers the
integral at the empirical Bayes estimate of *u*_{i},
defined as the vector that minimizes

with and set equal to their current estimates. The final
Hessian matrix from this optimization can be used to scale
the quadrature abscissas.
Suppose (*z*_{j}, *w*_{j}; *j* = 1, ... ,*p*) denote the standard Gauss-Hermite
abscissas and weights (Golub and Welsch 1969, or Table 25.10 of
Abramowitz and Stegun 1972). The adaptive Gaussian quadrature
integral approximation is as follows.

where *r* is the dimension of *u*_{i}, is the
Hessian matrix from the empirical Bayes minimization,
*z*_{j1, ... ,jr} is a vector with elements
(*z*_{j1}, ... ,*z*_{jr}), and

PROC NLMIXED selects the number of quadrature points adaptively by
evaluating the log likelihood function at the starting values of the
parameters until two successive evaluations have a relative difference
less than the value of the QTOL= option. The specific search sequence
is described under the QFAC= option. Using the QPOINTS= option, you
can adjust the number of quadrature points *p* to obtain different
levels of accuracy. Setting *p*=1 results in the Laplacian
approximation as described in Beal and Sheiner (1992), Wolfinger
(1993), Vonesh (1992, 1996), Vonesh and Chinchilli (1997), and
Wolfinger and Lin (1997).

The NOAD option in the PROC NLMIXED statement requests nonadaptive
Gaussian quadrature. Here all are set equal to zero, and
the Cholesky root of the estimated variance matrix of the random
effects is substituted for in the
preceding expression for *a*_{j1, ... ,jr}. The NOADSCALE option
requests the same scaling substitution but with the empirical Bayes
.

PROC NLMIXED computes the derivatives of the adaptive Gaussian
quadrature approximation when carrying out the default dual
quasi-Newton optimization.

*First-Order Method*

Another integral approximation available in PROC NLMIXED is the
first-order method of Beal and Sheiner (1982, 1988) and Sheiner and
Beal (1985). This approximation is used only in the case where
is normal, that is,

where *n*_{i} is the dimension of *y*_{i}, *R*_{i} is a diagonal
variance matrix, and *m*_{i} is the conditional mean vector of *y*_{i}.
The first-order approximation is obtained by expanding
with a one-term Taylor series expansion about
*u*_{i} = 0, resulting in the approximation

where is the Jacobian matrix
evaluated at *u*_{i} = 0.
Assuming that is normal with mean 0 and variance
matrix , the first-order integral approximation is computable
in closed form after completing the square:

where . The resulting approximation for is then
minimized over to obtain the first-order
estimates. PROC NLMIXED uses finite-difference derivatives of the
first-order integral approximation when carrying out the default
dual quasi-Newton optimization.

Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.