Recall from the section Notation for the Generalized Linear Mixed Model that
where and . Following Wolfinger and O’Connell (1993), a firstorder Taylor series of about and yields
where
is a diagonal matrix of derivatives of the conditional mean evaluated at the expansion locus. Rearranging terms yields the expression
The left side is the expected value, conditional on , of
and
You can thus consider the model
which is a linear mixed model with pseudoresponse , fixed effects , random effects , and .
Now define
as the marginal variance in the linear mixed pseudomodel, where is the parameter vector containing all unknowns in and . Based on this linearized model, an objective function can be defined, assuming that the distribution of is known. The GLIMMIX procedure assumes that has a normal distribution. The maximum log pseudolikelihood (MxPL) and restricted log pseudolikelihood (RxPL) for are then




with . f denotes the sum of the frequencies used in the analysis, and k denotes the rank of . The fixedeffects parameters are profiled from these expressions. The parameters in are estimated by the optimization techniques specified in the NLOPTIONS statement. The objective function for minimization is or . At convergence, the profiled parameters are estimated and the random effects are predicted as




With these statistics, the pseudoresponse and error weights of the linearized model are recomputed and the objective function is minimized again. The predictors are the estimated BLUPs in the approximated linear model. This process continues until the relative change between parameter estimates at two successive (outer) iterations is sufficiently small. See the PCONV= option in the PROC GLIMMIX statement for the computational details about how the GLIMMIX procedure compares parameter estimates across optimizations.
If the conditional distribution contains a scale parameter (Table 41.20), the GLIMMIX procedure profiles this parameter in GLMMs from the log pseudolikelihoods as well. To this end define
where is the covariance parameter vector with q – 1 elements. The matrices and are appropriately reparameterized versions of and . For example, if has a variance component structure and , then contains ratios of the variance components and , and . The solution for is
where m = f for MxPL and m = f – k for RxPL. Substitution into the previous functions yields the profiled log pseudolikelihoods,






Profiling of can be suppressed with the NOPROFILE option in the PROC GLIMMIX statement.
Where possible, the objective function, its gradient, and its Hessian employ the sweepbased Wtransformation ( Hemmerle and Hartley 1973; Goodnight 1979; Goodnight and Hemmerle 1979). Further details about the minimization process in the general linear mixed model can be found in Wolfinger, Tobias, and Sall (1994).
The GLIMMIX procedure produces estimates of the variability of , , and estimates of the prediction variability for , . Denote as the matrix
where all components on the right side are evaluated at the converged estimates. The mixed model equations (Henderson, 1984) in the linear mixed (pseudo)model are then
and




is the approximate estimated variancecovariance matrix of . Here, and .
The square roots of the diagonal elements of are reported in the Standard Error column of the “Parameter Estimates” table. This table is produced with the SOLUTION option in the MODEL statement. The prediction standard errors of the randomeffects solutions are reported in the Std Err Pred column of the “Solution for Random Effects” table. This table is produced with the SOLUTION option in the RANDOM statement.
As a cautionary note, tends to underestimate the true sampling variability of [, because no account is made for the uncertainty in estimating and . Although inflation factors have been proposed (Kackar and Harville, 1984; Kass and Steffey, 1989; Prasad and Rao, 1990), they tend to be small for data sets that are fairly well balanced. PROC GLIMMIX does not compute any inflation factors by default. The DDFM=KENWARDROGER option in the MODEL statement prompts PROC GLIMMIX to compute a specific inflation factor (Kenward and Roger, 1997), along with Satterthwaitebased degrees of freedom.
If is singular, or if you use the CHOL option of the PROC GLIMMIX statement, the mixed model equations are modified as follows. Let denote the lower triangular matrix so that . PROC GLIMMIX then solves the equations
and transforms and a generalized inverse of the leftside coefficient matrix by using .
The asymptotic covariance matrix of the covariance parameter estimator is computed based on the observed or expected Hessian matrix of the optimization procedure. Consider first the case where the scale parameter is not present or not profiled. Because is profiled from the pseudolikelihood, the objective function for minimization is for METHOD=MSPL and METHOD=MMPL and for METHOD=RSPL and METHOD=RMPL. Denote the observed Hessian (second derivative) matrix as
The GLIMMIX procedure computes the variance of by default as . If the Hessian is not positive definite, a sweepbased generalized inverse is used instead. When the EXPHESSIAN option of the PROC GLIMMIX statement is used, or when the procedure is in scoring mode at convergence (see the SCORING option in the PROC GLIMMIX statement), the observed Hessian is replaced with an approximated expected Hessian matrix in these calculations.
Following Wolfinger, Tobias, and Sall (1994), define the following components of the gradient and Hessian in the optimization:








Table 41.23 gives expressions for the Hessian matrix depending on estimation method, profiling, and scoring.
Table 41.23: Hessian Computation in GLIMMIX
Profiling 
Scoring 
MxPL 
RxPL 

No 
No 


No 
Yes 


No 
Modified 


Yes 
No 


Yes 
Yes 


Yes 
Modified 


The “Modified” expressions for the Hessian under scoring in RxPL estimation refer to a modified scoring method. In some cases, the modification leads to faster convergence than the standard scoring algorithm. The modification is requested with the SCOREMOD option in the PROC GLIMMIX statement.
Finally, in the case of a profiled scale parameter , the Hessian for the parameterization is converted into that for the parameterization as
where
There are two basic choices for the expansion locus of the linearization. A subjectspecific (SS) expansion uses
which are the current estimates of the fixed effects and estimated BLUPs. The populationaveraged (PA) expansion expands about the same fixed effects and the expected value of the random effects
To recompute the pseudoresponse and weights in the SS expansion, the BLUPs must be computed every time the objective function in the linear mixed model is maximized. The PA expansion does not require any BLUPs. The four pseudolikelihood methods implemented in the GLIMMIX procedure are the factorial combination between two expansion loci and residual versus maximum pseudolikelihood estimation. The following table shows the combination and the corresponding values of the METHOD= option (PROC GLIMMIX statement); METHOD=RSPL is the default.
Type of 
Expansion Locus 


PL 

E 
residual 
RSPL 
RMPL 
maximum 
MSPL 
MMPL 