The GLIMMIX Procedure

Quasi-likelihood for Independent Data

Quasi-likelihood estimation uses only the first and second moment of the response. In the case of independent data, this requires only a specification of the mean and variance of your data. The GLIMMIX procedure estimates parameters by quasi-likelihood, if the following conditions are met:

  • The response distribution is unknown, because of a user-specified variance function.

  • There are no G-side random effects.

  • There are no R-side covariance structures or at most an overdispersion parameter.

Under some mild regularity conditions, the function

\[ Q(\mu _ i,y_ i) = \int _{y_ i}^{\mu _ i} \frac{y_ i-t}{\phi a(t)} \, dt \]

known as the log quasi-likelihood of the ith observation, has some properties of a log-likelihood function (McCullagh and Nelder, 1989, p. 325). For example, the expected value of its derivative is zero, and the variance of its derivative equals the negative of the expected value of the second derivative. Consequently,

\[ QL(\bmu ,\phi ,\mb{y}) = \sum _{i=1}^{n} f_ i w_ i \frac{Y_ i - \mu _ i}{\phi a(\mu _ i)} \]

can serve as the score function for estimation. Quasi-likelihood estimation takes as the gradient and "Hessian" matrix—with respect to the fixed-effects parameters $\bbeta $—the quantities

\begin{align*} \mb{g}_{ql} & = \left[g_{ql,j}\right] = \left[ \frac{\partial QL(\bmu ,\phi ,\mb{y})}{\partial \beta _ j} \right] = \mb{D}’\mb{V}^{-1}(\mb{Y}-\bmu )/\phi \\ \mb{H}_{ql} & = \left[ h_{ql,jk} \right] = \left[ \frac{\partial ^2 QL(\bmu ,\phi ,\mb{y})}{\partial \beta _ j \partial \beta _ k} \right] = \mb{D}’\mb{V}^{-1}\mb{D}/\phi \end{align*}

In this expression, $\bD $ is a matrix of derivatives of $\bmu $ with respect to the elements in $\bbeta $, and $\mb{V}$ is a diagonal matrix containing variance functions, $\mb{V} = [a(\mu _1), \cdots , a(\mu _ n)]$. Notice that $\bH _{ql}$ is not the second derivative matrix of $Q(\bmu ,\mb{y})$. Rather, it is the negative of the expected value of $\partial \mb{g}_{ql} / \partial \bbeta $. $\bH _{ql}$ thus has the form of a "scoring Hessian."

The GLIMMIX procedure fixes the scale parameter $\phi $ at 1.0 by default. To estimate the parameter, add the statement 

random _residual_;

The resulting estimator (McCullagh and Nelder, 1989, p. 328) is

\[ \widehat{\phi } = \frac{1}{m}\sum _{i=1}^ n f_ i w_ i \frac{\left(y_ i - \widehat{\mu }_ i\right)^2}{a(\widehat{\mu }_ i)} \]

where $m = f-\mr{rank}\{ \mb{X}\} $ if the NOREML option is in effect, m = f otherwise, and f is the sum of the frequencies.

See Example 44.4 for an application of quasi-likelihood estimation with PROC GLIMMIX.