Method-of-Moments Estimation for Variance Components

Why Use Quadratic-Form Moments?

A covariance-component model represents variation as a sum of structured contributions. For example, repeated measurements on the same subject may vary because of both subject-to-subject differences and measurement noise. Method-of-moments estimation turns the expected values of selected data summaries into equations for the unknown component variances. Quadratic forms are especially useful because their expectations depend linearly on the covariance matrix.

The resulting equations are often simple to solve, and their construction makes the source of each estimate visible. But a solution to the moment equations is not automatically a valid covariance-component estimate: it can be noisy, nonunique, or outside the allowed parameter space. These are separate questions from whether it is unbiased.

Quadratic-Form Moments

Definition: Suppose $Y$ is a random vector with mean $\mu$ and covariance $V$, and $B$ is a fixed symmetric matrix of matching dimension. The quadratic form $Q=Y'BY$ has expectation $$ E[Q]=\operatorname{tr}(BV)+\mu'B\mu, $$ provided the second moments exist. If $Y$ is centered, this reduces to $E[Q]=\operatorname{tr}(BV)$.

The identity follows by expanding the entries: $Y'BY=\sum_{i,j}B_{ij}Y_iY_j$, so its expectation is $\sum_{i,j}B_{ij}E[Y_iY_j]$. Since $E[YY']=V+\mu\mu'$, this equals $\operatorname{tr}(B(V+\mu\mu'))=\operatorname{tr}(BV)+\mu'B\mu$. No normality assumption is needed.

Consider a linear covariance-component model $$ V(\theta)=\sum_{r=1}^{k}\theta_r A_r, $$ where the matrices $A_r$ are known symmetric matrices and $\theta=(\theta_1,\ldots,\theta_k)'$ contains unknown component parameters. In the usual variance-component case, the parameters are variances and must be nonnegative. If the mean is known to be zero, or has been removed by suitable contrasts, the expectation of each quadratic form is linear in $\theta$.

More precisely, let $B_1,\ldots,B_m$ be fixed symmetric matrices and $Q_j=Y'B_jY$. Define $M_{jr}=\operatorname{tr}(B_jA_r)$, and write $Q=(Q_1,\ldots,Q_m)'$. Then $E[Q]=M\theta$. When $m=k$ and $M$ is nonsingular, solving the moment equations gives $\hat{\theta}_{\mathrm{MM}}=M^{-1}Q$.

Theorem: Under the centered model above, if $M$ is nonsingular and the quadratic forms have finite expectations, the estimator $\hat{\theta}_{\mathrm{MM}}=M^{-1}Q$ is unbiased: $E[\hat{\theta}_{\mathrm{MM}}]=\theta$. More generally, any fixed matrix $C$ satisfying $CM=I_k$ gives the unbiased estimator $CQ$, even when more than $k$ moments are used.

Proof. By linearity of expectation and the quadratic-form identity, $E[Q]=M\theta$. Therefore $E[M^{-1}Q]=M^{-1}E[Q]=M^{-1}M\theta=\theta$. In the more general case, $E[CQ]=CE[Q]=CM\theta=\theta$. The rank condition is essential: without enough independent moment equations, this inversion argument cannot identify all components. $\square$

When the mean is not known, using the raw quadratic form generally adds $\mu'B_j\mu$ to its expectation. One standard way to remove a linear fixed mean is to use contrasts. If $E[Y]=X\beta$, let the columns of $K$ form an orthonormal basis for the null space of $X'$, so $K'X=0$. Then $Z=K'Y$ has mean zero and covariance $K'V(\theta)K$. Quadratic forms in $Z$ can be used in the same construction, with component matrices $K'A_rK$. This requires that the chosen contrasts retain enough information to distinguish the components.

Worked Example: Two Components From Two Quadratic Forms

Suppose $Y=(Y_1,Y_2)'$ has known mean zero and covariance $V=\theta_1 I_2+\theta_2 J$, where $I_2$ is the $2$-by-$2$ identity matrix and $J$ is the $2$-by-$2$ matrix of ones. Use $Q_1=Y'Y$ and $Q_2=Y'JY$.

The moment coefficients are $M_{jr}=\operatorname{tr}(B_jA_r)$, with $B_1=I_2$, $B_2=J$, $A_1=I_2$, and $A_2=J$. Thus $$ M= \begin{pmatrix} \operatorname{tr}(I_2I_2)&\operatorname{tr}(I_2J)\\ \operatorname{tr}(JI_2)&\operatorname{tr}(JJ) \end{pmatrix} = \begin{pmatrix}2&2\\2&4\end{pmatrix}. $$

Its determinant is $2(4)-2(2)=4$, so the equations identify both components. For observed $Y=(1,2)'$, $Q_1=1^2+2^2=5$, while $Q_2=(1+2)^2=9$. The equations are $2\hat\theta_1+2\hat\theta_2=5$ and $2\hat\theta_1+4\hat\theta_2=9$. Subtracting the first from the second gives $2\hat\theta_2=4$, so $\hat\theta_2=2$; then $2\hat\theta_1+4=5$, giving $\hat\theta_1=0.5$. The estimates solve the moments, but this single calculation does not establish their precision.

Balanced Random-Intercept Model

A familiar special case shows how quadratic-form moments lead to analysis-of-variance estimators. Suppose there are $g$ groups, each with $m$ observations, and $$ Y_{ij}=\mu+a_i+e_{ij},\qquad i=1,\ldots,g,\quad j=1,\ldots,m. $$ Assume the random intercepts $a_i$ are independent with mean zero and variance $\tau^2$; the errors $e_{ij}$ are independent with mean zero and variance $\sigma^2$; and all intercepts and errors are mutually independent. The component parameters are $\tau^2$ and $\sigma^2$. Normality is not needed for the expectation calculations below.

Let $\bar Y_{i\cdot}$ be group $i$'s mean and $\bar Y_{\cdot\cdot}$ the grand mean. Define the within-group and between-group sums of squares by $$ SS_W=\sum_{i=1}^{g}\sum_{j=1}^{m}(Y_{ij}-\bar Y_{i\cdot})^2,\qquad SS_B=m\sum_{i=1}^{g}(\bar Y_{i\cdot}-\bar Y_{\cdot\cdot})^2. $$

Within each group, subtracting its sample mean removes the group-constant terms $\mu+a_i$. The within-group residuals are therefore $e_{ij}-\bar e_{i\cdot}$. Their squared sum has expectation $(m-1)\sigma^2$: expanding it gives $\sum_j e_{ij}^2-m\bar e_{i\cdot}^{\,2}$, whose expectation is $m\sigma^2-m(\sigma^2/m)=(m-1)\sigma^2$. Summing over groups yields $E[SS_W]=g(m-1)\sigma^2$.

For the between-group term, $\bar Y_{i\cdot}=\mu+a_i+\bar e_{i\cdot}$. These group means are independent, each with variance $\tau^2+\sigma^2/m$. For $g$ independent, identically distributed variables with variance $v$, the expected sum of squared deviations from their sample mean is $(g-1)v$, by the same expansion used above. Multiplying by $m$ gives $E[SS_B]=m(g-1)\tau^2+(g-1)\sigma^2$.

Consequently, the mean squares $MS_W=SS_W/[g(m-1)]$ and $MS_B=SS_B/(g-1)$ have expectations $\sigma^2$ and $\sigma^2+m\tau^2$, respectively. Solving these two moment equations gives the following estimates.

Result: For the balanced random-intercept model, $\hat\sigma^2_{\mathrm{MM}}=MS_W$ and $\hat\tau^2_{\mathrm{MM}}=(MS_B-MS_W)/m$ are unbiased under the stated assumptions. Their unbiasedness does not ensure that both realized estimates are nonnegative.

Worked Example: Computing the Balanced-Design Estimates

Suppose there are $g=3$ groups of size $m=2$, with observations $(2,4)$, $(5,7)$, and $(8,10)$. The group means are $3,6,9$, and the grand mean is $6$.

Within each group, the two deviations from its mean are $-1$ and $1$, contributing $2$ to $SS_W$. Thus $SS_W=2+2+2=6$, and $MS_W=6/[3(2-1)]=2$.

The between-group sum of squares is $SS_B=2[(3-6)^2+(6-6)^2+(9-6)^2]=2(9+0+9)=36$. Since $g-1=2$, $MS_B=36/2=18$. Therefore $\hat\sigma^2_{\mathrm{MM}}=2$ and $\hat\tau^2_{\mathrm{MM}}=(18-2)/2=8$. Both happen to be nonnegative for these data; that is a feature of this example, not a guarantee of the method.

Unbiasedness, Bias, and Feasibility

An estimator is unbiased if its expectation equals the target parameter under the specified model. Feasibility is a different property: a variance-component estimate is feasible if it satisfies the parameter constraints and produces an allowed covariance matrix. In a model with nonnegative variance components, feasibility requires $\hat\theta_r\geq 0$ for every component. More generally, the defining requirement is $V(\hat\theta)\succeq 0$, meaning that the estimated covariance matrix is positive semidefinite.

An unbiased moment estimator can fail this check because it is a linear combination of random quadratic forms, and such a combination need not be nonnegative for every data vector. The following calculation also shows why simply truncating a negative estimate has a cost.

Worked Example: A Negative Estimate and the Bias of Truncation

Return to the two-dimensional model $V=\theta_1I_2+\theta_2J$. The moment equations give $\hat\theta_2=Y_1Y_2$. For observed $Y=(1,-1)'$, $Q_1=2$ and $Q_2=0$, so $\hat\theta_2=(Q_2-Q_1)/2=-1$. This estimate violates the usual constraint $\theta_2\geq0$, even though it solves the unbiased moment equations.

Now suppose the true $\theta_2=0$ and $\theta_1>0$. Then $Y_1$ and $Y_2$ are independent $N(0,\theta_1)$ variables, and $E[Y_1Y_2]=0$. If we replace the estimate by $\max(0,Y_1Y_2)$, its expectation is $$ E[\max(0,Y_1Y_2)] =\tfrac12 E[|Y_1Y_2|] =\tfrac12 E[|Y_1|]E[|Y_2|] =\tfrac12\left(\sqrt{\frac{2\theta_1}{\pi}}\right)^2 =\frac{\theta_1}{\pi}>0. $$ Here the first equality uses the symmetry of $Y_1Y_2$ about zero, and the second uses independence. The truncated estimate is feasible but biased upward at this parameter value.

Key distinction: “Admissible” can refer to a decision-theoretic property, which requires specifying a loss function and a class of competing estimators. In covariance-component calculations, it is important to state whether one means that an estimate is feasible for the covariance model or that it is decision-theoretically admissible. Feasibility alone does not establish decision-theoretic admissibility, and unbiasedness alone does not establish feasibility.

Common Pitfalls and Practical Checks

Ignoring the mean term. For an uncentered vector, $E[Y'BY]$ includes $\mu'B\mu$. Use justified contrasts or account for the mean before treating the moment as a covariance equation.
Solving moments without checking rank. A numerical solution requires enough independent equations. Check the rank of $M$; if it is deficient, the selected moments do not identify every component.
Calling a negative variance estimate impossible under the model. The true variance component may be constrained to be nonnegative while an unconstrained unbiased estimator can still be negative in a particular sample.
Truncating and retaining the unbiasedness claim. Truncation changes the estimator's sampling distribution. Its bias must be assessed under the model; it cannot be assumed to remain zero.

Moment estimates are useful as transparent summaries and as initial estimates for further analysis. Their central checks are distinct: verify that the moment equations match the stated mean and covariance model, that the equations identify the components, and that the resulting estimates are compatible with the covariance constraints. A failure of feasibility is a warning about the realized estimate, not by itself evidence that the moment equations were derived incorrectly.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

Why does $E[Y'BY]$ contain a mean term when $E[Y]\ne0$, and what condition on contrasts removes it for a linear mean model?
What rank condition allows the square system of quadratic-form moments to identify all covariance components?
In the balanced random-intercept model, what are $E[MS_W]$ and $E[MS_B]$, and how do they yield estimates of $\sigma^2$ and $\tau^2$?
How can an estimator be unbiased and yet produce a negative variance-component estimate in a particular sample?
Why does replacing a negative estimate by zero generally require a new bias analysis?