Maximum Likelihood Estimation Based on Adaptive Quadrature

Quadrature methods, like the Laplace approximation, approximate integrals. If you choose METHOD=QUAD for a generalized linear mixed model, the GLIMMIX procedure approximates the marginal log likelihood with an adaptive Gauss-Hermite quadrature rule. Gaussian quadrature is particularly well suited to numerically evaluate integrals against probability measures (Lange 1999, Ch. 16). And Gauss-Hermite quadrature is appropriate when the density has kernel and integration extends over the real line, as is the case for the normal distribution. Suppose that is a probability density function and the function is to be integrated against it. Then the quadrature rule is

where denotes the number of quadrature points, the are the quadrature weights, and the are the abscissas. The Gaussian quadrature chooses abscissas in areas of high density, and if is continuous, the quadrature rule is exact if is a polynomial of up to degree . In the generalized linear mixed model the roles of and are played by the conditional distribution of the data given the random effects, and the random-effects distribution, respectively. Quadrature abscissas and weights are those of the standard Gauss-Hermite quadrature (Golub and Welsch 1969; see also Table 25.10 of Abramowitz and Stegun 1972; Evans 1993).

A numerical integration rule is called adaptive when it uses a variable step size to control the error of the approximation. For example, an adaptive trapezoidal rule uses serial splitting of intervals at midpoints until a desired tolerance is achieved. The quadrature rule in the GLIMMIX procedure is adaptive in the following sense: if you do not specify the number of quadrature points (nodes) with the QPOINTS= suboption of the METHOD=QUAD option, then the number of quadrature points is determined by evaluating the log likelihood at the starting values at a successively larger number of nodes until a tolerance is met (for more details see the text under the heading "Starting Values" in the next section). Furthermore, the GLIMMIX procedure centers and scales the quadrature points by using the empirical Bayes estimates (EBEs) of the random effects and the Hessian (second derivative) matrix from the EBE suboptimization. This centering and scaling improves the likelihood approximation by placing the abscissas according to the density function of the random effects. It is not, however, adaptiveness in the previously stated sense.

Objective Function

Let denote the vector of fixed-effects parameters and the vector of covariance parameters. For quadrature estimation in the GLIMMIX procedure, includes the G-side parameters and a possible scale parameter , provided that the conditional distribution of the data contains such a scale parameter. is the vector of the G-side parameters. The marginal distribution of the data for subject in a mixed model can be expressed as

Suppose denotes the number of quadrature points in each dimension (for each random effect) and denotes the number of random effects. For each subject, obtain the empirical Bayes estimates of as the vector that minimizes

If are the standard abscissas for Gauss-Hermite quadrature, and is a point on the -dimensional quadrature grid, then the centered and scaled abscissas are

As for the Laplace approximation, is the second derivative matrix with respect to the random effects,

These centered and scaled abscissas, along with the Gauss-Hermite quadrature weights , are used to construct the -dimensional integral by a sequence of one-dimensional rules

The right-hand side of this expression, properly accumulated across subjects, is the objective function for adaptive quadrature estimation in the GLIMMIX procedure.

   proc glimmix method=quad(qpoints=1);

   proc glimmix method=quad(qpoints=5);
      class A id;
      model y = / dist=negbin;
      random A / subject=id;
   run;

   proc glimmix method=quad(qpoints=5);
      class A B id;
      model y = / dist=negbin;
      random A A*B / subject=id;
   run;