The VARCOMP Procedure

Computational Methods

Subsections:

Four methods of estimation can be specified in the PROC VARCOMP statement by using the METHOD= option. They are described in the following sections.

The Type I Method

This method (METHOD= TYPE1) computes the Type I sum of squares for each effect, equates each mean square involving only random effects to its expected value, and solves the resulting system of equations (Gaylor, Lucas, and Anderson 1970). The $\bX ’\bX ~ |~ \bX ’\bY $ matrix is computed and adjusted in segments whenever memory is not sufficient to hold the entire matrix.

The MIVQUE0 Method

Based on the technique suggested by Hartley, Rao, and LaMotte (1978), the MIVQUE0 method (METHOD=MIVQUE0) produces unbiased estimates that are invariant with respect to the fixed effects of the model and that are locally best quadratic unbiased estimates given that the true ratio of each component to the residual error component is zero. The technique is similar to TYPE1 except that the random effects are adjusted only for the fixed effects. This affords a considerable timing advantage over the TYPE1 method; thus, MIVQUE0 is the default method used in PROC VARCOMP. The $\bX ’\bX | \bX ’\bY $ matrix is computed and adjusted in segments whenever memory is not sufficient to hold the entire matrix. Each element $(i,j)$ of the form

\[ \mbox{SSQ}(\bX _ i’\bM \bX _ j) \]

is computed, where

\[ \bM =\bI -\bX _0(\bX _0’\bX _0)^{-}\bX _0’ \]

and where $\bX _0$ is part of the design matrix for the fixed effects, $\bX _ i$ is part of the design matrix for one of the random effects, and SSQ is an operator that takes the sum of squares of the elements. For more information see Rao (1971, 1972) and Goodnight (1978).

The Maximum Likelihood Method

The maximum likelihood method ( METHOD= ML) computes maximum likelihood estimates of the variance components; see Searle, Casella, and McCulloch (1992). The computing algorithm makes use of the W-transformation developed by Hemmerle and Hartley (1973) and Goodnight and Hemmerle (1979). The procedure uses a Newton-Raphson algorithm, iterating until the log-likelihood objective function converges.

The objective function for METHOD=ML is $\ln (|\bV |) + \mb{r}’ \bV ^{-1} \mb{r}$, where

\[ \bV =\sigma _0^2\bI +\sum _{i=1}^{n_ r}\sigma _ i^2\bX _ i\bX _ i’ \]

and where $\sigma _0^2$ is the residual variance, $n_ r$ is the number of random effects in the model, $\sigma _ i^2$ represents the variance components, $\bX _ i$ is part of the design matrix for one of the random effects, and

\[ \mb{r} = \mb{y} - \bX _0 (\bX _0’ \bV ^{-1} \bX _0)^{-} \bX _0’ \bV ^{-1} \mb{y} \]

is the vector of residuals.

The Restricted Maximum Likelihood Method

The restricted maximum likelihood method ( METHOD= REML) is similar to the maximum likelihood method, but it first separates the likelihood into two parts: one that contains the fixed effects and one that does not (Patterson and Thompson 1971). The procedure uses a Newton-Raphson algorithm, iterating until convergence is reached for the log-likelihood objective function of the portion of the likelihood that does not contain the fixed effects. Using notation from earlier methods, the objective function for METHOD=REML is $\ln (|\bV |) + \mb{r}’ \bV ^{-1} \mb{r} + \ln (|\bX _0’ \bV ^{-1} \bX _0|)$. See Searle, Casella, and McCulloch (1992) for additional details.

The GRR Method

Based on the technique suggested by Burdick, Borror, and Montgomery (2005), the GRR method (METHOD= GRR) produces minimum variance unbiased estimators.