Chapter Contents
Chapter Contents
The NLMIXED Procedure

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 ui centers the integral at the empirical Bayes estimate of ui, defined as the vector \hat{u_i} that minimizes
-\log [ p(y_i | X_i,\phi, u_i) q(u_i | \xi) ]
with \phi and \xi set equal to their current estimates. The final Hessian matrix from this optimization can be used to scale the quadrature abscissas.

Suppose (zj, wj; 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.

eqn{ \int p(y_i | X_i, \phi, u_i) q(u_i | \xi) d u_i \approx } \ & & 2^{r/2} | \...
 .... ,j_r}) q(a_{j_1, ... ,j_r} | \xi) 
 \prod_{k=1}^r w_{j_k} \exp{z^2_{j_k}} ] \

where r is the dimension of ui, \Gamma(X_i,\theta) is the Hessian matrix from the empirical Bayes minimization, zj1, ... ,jr is a vector with elements (zj1, ... ,zjr), and

a_{j_1, ... ,j_r} = \hat{u}_i + 
 2^{1/2} \Gamma(X_i,\theta)^{-1/2} z_{j_1, ... ,j_r}

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 \hat{u_i} are set equal to zero, and the Cholesky root of the estimated variance matrix of the random effects is substituted for \Gamma(X_i,\theta)^{-1/2} in the preceding expression for aj1, ... ,jr. The NOADSCALE option requests the same scaling substitution but with the empirical Bayes \hat{u_i}.

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 p(y_i | X_i, \phi, u_i) is normal, that is,
eqn{ p(y_i | X_i, \phi, u_i) = 
 (2 \pi)^{-n_i/2} | R_i(X_i,\phi) | ^{-1/2} } \ ...
 ..._i - m_i(X_i,\phi,u_i) ]^T R_i(X_i,\phi)^{-1} 
 [ y_i - m_i(X_i,\phi,u_i) ] \}
where ni is the dimension of yi, Ri is a diagonal variance matrix, and mi is the conditional mean vector of yi.

The first-order approximation is obtained by expanding m(X_i, \phi, u_i) with a one-term Taylor series expansion about ui = 0, resulting in the approximation

p(y_i | X_i, \phi, u_i) & \approx & 
 (2 \pi)^{-n_i/2} | R_i(X_i,\phi) | ^{-1/2}...
 ...i(X_i,\phi)^{-1} . \ & &  . [ y_i - m_i(X_i,\phi,0) - Z_i(X_i,\phi) 
 u_i ] )
where Z_i(X_i,\phi) is the Jacobian matrix \partial m_i(X_i,\phi,u_i) / \partial u_i evaluated at ui = 0.

Assuming that q(u_i | \xi) is normal with mean 0 and variance matrix G(\xi), the first-order integral approximation is computable in closed form after completing the square:

eqn{ \int p(y_i | X_i, \phi, u_i) q(u_i | \xi) d u_i \approx
 (2 \pi)^{-n_i/2} |...
 ...[ y_i - m_i(X_i,\phi,0) ]^T V_i(X_i,\theta)^{-1} 
 [ y_i - m_i(X_i,\phi,0) ] )

where V_i(X_i,\theta) = Z_i(X_i,\phi) G(\xi) Z_i(X_i,\phi)^T +
R_i(X_i,\phi). The resulting approximation for f(\theta) is then minimized over \theta = (\phi, \xi) 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.

Chapter Contents
Chapter Contents

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