On the Sub-D’s mean square error on “one-way” random designs

Abstract

Through simulations, it was shown that either in balanced “one-way” random designs or in unbalanced “one-way” random designs, Sub-D estimates are in general more accurate than those provided by Anova-based estimators. Moreover, the Sub-D estimates exhibit less dispersion magnitude. Such estimates reveal to be also slightly more accurate than those provide by REML-based estimator, although this latter one presented dispersion with slightly less magnitude, which, indeed, is not bigger than 0.0261. In order to have somehow a robust tool which will allow us to infer over Sub-D’s efficiency in “one-way” random designs, this paper aims to deduce and discuss its MSE.

Keywords

Variance components Sub-D MSE efficiency

1. Introduction

Searching for new contributions concerning the inference in Mixed Linear Models (MLM), Silva (2017) developed the new estimators for variance components called Sub-D and it’s corresponding improvement called Sub-DI. It’s relevant to emphasize that both estimators apply to MLM with an arbitrary number of variance components. It was proved that in MLM with two variance components, Sub-D and Sub-DI are equivalent, having better performance than the Anova-based method estimators and nearly the same results as the Likelihood-based estimators. More over, since Sub-D depends only on information retained by the eigenvalues of the design matrix and quadratic errors of sub-models (see Silva (2017), Silva and Fonseca (2018) and Ferreira et al. (2017)) it provides unbiased estimates whether the data are balanced or unbalanced in both crossed and nested designs, even having empty cells, which does not hold in case of Anova and Likelihood-based estimators.

Through simulations carried out in balanced “one-way” random designs it was shown that the Sub-D and Anova estimates as well as the corresponding dispersion magnitude are exactly the same (one may find the one-way random designs characterization in Silva (2017) or, for instance, Anderson (1975) and Anderson and Crump (1967)). Such estimates reveal more accurate than those provide by REML-based one, although the latter one presented dispersion with slightly less magnitude, more precisely, not bigger than 0.0261. By carrying out simulations in unbalanced “one-way” random designs it was shown that Sub-D still producing accurate estimates while those based on Anova and REML produce estimates with low accuracy (see Silva (2017)); also, it is worth highlight that the Sub-D estimates presents a slightly more dispersion than the REML ones.

Since the simulations studies carried out by Silva (2017) suggested that Sub-D has a good performance in “one-way” random designs (MLM with two variance components), better than those based on Anova methods, and somehow comparable to those based on REML methods, for the point of view of efficiency, it seems to be interesting to study its MSE. In fact, this problem has already been addressed by Klotz et al. (2012), although for likelihood-based estimators and only in balanced “one-way” random designs. Here we only set down the MSE of Sub-D, leaving the discussion regarding its efficiency versus the efficiency of Anova-based and likelihood-based estimators for the next issue.

Section 2 is dedicated to the Sub-D’s deduction (the background). In Section 3 the Sub-D’s MSE is introduced and discussed, while in Section 4 a practical example is explored. Finally, in Section 5 a final remark is made.

2. Background

We approach the random “one-way” design

$\displaystyle z_{ij}=\mu+\alpha_{i}+e_{ij},∼{}i=1,\ldots,k,∼{}j=1,\ldots,n_{i},$ (1)

where $\mu$ denotes the general mean, $\alpha_{i}$ the observed random effect due to the $i$ th group of treatments, with mean zero and variance $\gamma_{1}\geqslant 0$ , and $e_{ij}$ the random residual errors with mean zero and variance $\gamma_{2}>0$ . It is assumed that $\alpha_{i}$ and $e_{ij}$ are independently observed.

Let ${\bm{1}}_{s}$ and ${\bm{0}}_{s}$ denote the $s$ dimension vectors of 1’s and 0’s, respectively. Let also $\mathcal{M}^{s,t}$ denotes the $s\times t$ matrices (matrices with $s$ rows and $t$ columns) and $I_{s}$ the identity matrices of dimension $s$ . Mor over, let $w\sim(u,A)$ denotes a random vector $w$ with mean $u$ and variance-covariance matrix $A$ .

With

$\displaystyle m=\sum_{i=1}^{k}n_{i},∼{}z=[z_{11}\ldots z_{kn_{k}}]^{\rm T},∼{}% X={\bm{1}}_{m},∼{}X_{1}=\left[\begin{array}[]{@{}cccc@{}}{1}_{n_{1}}&{0}_{n_{1% }}&\ldots&{0}_{n_{1}}\\ {0}_{n_{2}}&{1}_{n_{2}}&\ldots&{0}_{n_{2}}\\ \vdots&\vdots&\ddots&\vdots\\ {0}_{n_{k}}&{0}_{n_{k}}&\ldots&{1}_{n_{k}}\end{array}\right]\in\mathcal{M}^{m,% k},\text{∼{}and∼{}}\epsilon=[e_{11}∼{}\ldots∼{}e_{kn_{k}}]^{\rm T},$

a matrix formulation for model Eq. (1) is

$\displaystyle z=X\mu+X_{1}\alpha+\epsilon,$ (2)

where $\alpha\sim({0}_{k},\gamma_{1}I_{k})$ and $\epsilon\sim({0}_{m},\gamma_{2}I_{m})$ ; consequently it follows that

$\displaystyle z\sim(X\mu,\gamma_{1}X_{1}X_{1}^{\rm T}+\gamma_{2}I_{m}).$

We are also assuming that $m\gg k$ which is indeed an usual assumption, once the number of groups of treatment is largely small then the total number of sum of observations for all groups.

It must be remarked that when there is empty cells (which means no observation for certain group or groups of treatments) $X_{1}$ will not be a full rank matrix. Indeed, in such case it holds $r(X_{1})\leqslant k$ .

Since the parameters to be estimated do not depend on the fixed effect, Ferreira et al. (2017) and Silva (2017) thought it was convenient to remove the model dependence on the fixed effect, remarking that such an action cause no loss of information needed to estimate those parameters and consequently the reduction on the complexity of the model for the algebraic manipulation, as well as the bias in estimation process. The strategy followed is in all similar to the first phase of REML: the observations vector is projected on the orthogonal complement of $R(X)$ , the vectorial subspace spanned by the columns of the main vector $X$ . Thus, following such a strategy, we found the following restricted model:

$\displaystyle y=B^{\rm T}z\sim({0}_{n},\gamma_{1}M+\gamma_{2}I_{n}),$ (3)

with $M=B^{\rm T}X_{1}X_{1}^{\rm T}B$ and $B$ the matrix whose columns are the eigenvectors associated to the null eigenvalues of $P_{R(X)}$ , that is the projection matrix onto $R(X)$ , holding therefore

$\displaystyle B^{\rm T}B=I_{m-r(P_{R(X)})}\text{ and }BB^{\rm T}=I_{m}-P_{R(X)% }=P_{R(X)^{\bot}},$

noting that $n=m-r(P_{R(X)})=m-1$ and $r(B)=n$ .

2.1 Discussing the rank of

M

Recall $M=B^{\rm T}X_{1}X_{1}^{\rm T}B=B^{\rm T}NB$ , with $N=X_{1}X_{1}^{\rm T}$ , and let $h-1=r(M)$ .

With $X_{1}$ being a full rank matrix and so $N$ (that is $r(N)=k$ ), according with Theorem 2.10 of Schott (1997), we will have

$\displaystyle h-1=r(B^{\rm T}NB)\leqslant\min\{r(B^{\rm T}N),r(B)\}\leqslant% \min\{\min\{r(B^{\rm T}),r(N)\},m-1\}\leqslant\min\{\min\{m-1,k\},m-1\}% \leqslant k,$ (4)

holding $h\leqslant k+1$ . On the other hand

$\displaystyle h-1=r(B^{\rm T}NB)\geqslant r(B^{\rm T}N)+r(NB)-r(N)\geqslant(m-% 1+k-m)+(k+m-1-m)-k\geqslant k-2,$

holding $h\geqslant k-1$ . Then, $k-1\leqslant h\leqslant k+1$ .

Now, with $N$ not necessairelly a full rank matrix, that is $r(N)\leqslant k$ , we have the following:

$\displaystyle h-1=r(B^{\rm T}NB)\leqslant\min\{r(B^{\rm T}N),r(B)\}\leqslant% \min\{\min\{r(B^{\rm T}),r(N)\},m-1\}\leqslant\min\{\min\{m-1,k\},m-1\}% \leqslant k,$ (5)

holding $h\leqslant k+1$ .

Thus, in both case we have that $h\leqslant k+1$ .

2.2 Estimator Sub-D

With $Q=[A_{1}^{\rm T}\ldots A_{h}^{\rm T}]^{\rm T}$ , where $A_{i}^{\rm T}$ , $i=1,\ldots,h$ , is the matrix whose columns are the $g_{i}=r(A_{i}^{\rm T})$ orthonormal eigenvectors associated to the $i$ th eigenvalue of $M$ , say $\theta_{i}$ , holding therefore $Q^{\rm T}Q=QQ^{\rm T}=I_{n}$ , we may produce the new restricted model:

$\displaystyle y^{*}=Q^{\rm T}y\sim({0}_{\sum_{i=1}^{h}g_{i}},\gamma_{1}QMQ^{% \rm T}+\gamma_{2}QQ^{\rm T}),$ $\displaystyle\gamma_{1}QMQ^{\rm T}+\gamma_{2}QQ^{\rm T}=\gamma_{1}\left[\begin% {array}[]{@{}cccc@{}}\theta_{1}I_{g_{1}}&{0}_{g_{1},g_{2}}&\ldots&{0}_{g_{1},g% _{h}}\\ {0}_{g_{2},g_{1}}&\theta_{2}I_{g_{2}}&\ldots&{0}_{g_{2},g_{h}}\\ \vdots&\vdots&\ddots&\vdots\\ {0}_{g_{h},g_{1}}&{0}_{g_{h},g_{2}}&\ldots&\theta_{h}I_{g_{h}}\end{array}% \right]+\gamma_{2}I_{\sum_{i=1}^{h}g_{i}}=\left[\begin{array}[]{@{}cccc@{}}% \lambda_{1}I_{g_{1}}&{0}_{g_{1},g_{2}}&\ldots&{0}_{g_{1},g_{h}}\\ {0}_{g_{2},g_{1}}&\lambda_{2}I_{g_{2}}&\ldots&{0}_{g_{2},g_{h}}\\ \vdots&\vdots&\ddots&\vdots\\ {0}_{g_{h},g_{1}}&{0}_{g_{h},g_{2}}&\ldots&\lambda_{h}I_{g_{h}}\end{array}% \right],$ (6)

where $\lambda_{i}=\gamma_{1}\theta_{i}+\gamma_{2}$ and $\theta_{h}=0$ , since $h-1=r(M)\leqslant k<n$ (notice $M$ has $h-1$ nonnull eigenvalues, say $\theta_{1},\ldots,\theta_{h-1}$ , with respective multiplicity $g_{1},\ldots,g_{h-1}$ , and one null eigenvalues with multiplicity $g_{h}=n-\sum_{i=1}^{h}g_{i}$ ).

Now, with $\gamma=[\gamma_{1}\gamma_{2}]^{\rm T}$ , the estimator Sub-D for $\gamma$ in model Eq. (3) is given as

$\displaystyle\hat{\gamma}=\left[\begin{array}[]{@{}c@{}}\hat{\gamma_{1}}\\ \hat{\gamma_{2}}\end{array}\right]=(\Theta^{\rm T}\Theta)^{-1}\Theta^{\rm T}S,$ (7)

where $\Theta=\left[\begin{array}[]{@{}cc@{}}\theta_{1}&1\\ \vdots&\vdots\\ \theta_{h}&1\end{array}\right]$ , $S=\left[\begin{array}[]{@{}c@{}}\frac{y_{1}^{\rm T}y_{1}}{g_{1}}\\ \vdots\\ \frac{y_{h}^{\rm T}y_{h}}{g_{h}}\end{array}\right]$ and $y_{i}=A_{i}y$ , $i=1,\ldots,h$ , the sub-models. See Silva (2017) for additional explanation.

3. Mean square error

For the several aspects of the efficiency of an estimator - namely Asymptotic efficiency, relative efficiency, etc - discussing the efficiency of an estimator goes through a strait analysis of its variance-covariance and/or MSE.

Let $E(x)$ and $\Sigma(x)$ , respectively, denote the expectation and variance-covariance matrix of a random variable (or vector) $x$ , and $\operatorname{tr}(A)$ the trace of a real matrix $A$ (see Schott (1997)). The cross-covariance between the random variable (or vector) $x$ and $s$ is denoted as $\Sigma(x,s)$ . Let also $\text{MSE}(\hat{\gamma})$ denotes the MSE of $\hat{\gamma}$ .

$\displaystyle E(\hat{\gamma})=(\Theta^{\rm T}\Theta)^{-1}\Theta^{\rm T}E(S)=(% \Theta^{\rm T}\Theta)^{-1}\Theta^{\rm T}[\lambda_{1}\ldots\lambda_{h}]^{\rm T}% =(\Theta^{\rm T}\Theta)^{-1}\Theta^{\rm T}\Theta\gamma=\gamma,$ (8)

since $y_{i}\sim({0}_{g_{i}},\lambda_{i}I_{g_{i}})$ and $E\left(\frac{y_{1}^{\rm T}y_{1}}{g_{1}}\right)=\frac{\lambda_{i}}{g_{i}}% \operatorname{tr}(I_{g_{i}})=\lambda_{i}$ .

Since $\hat{\gamma}$ is an unbiased estimator according with Eq. (8), it follows that

$\displaystyle\text{MSE}(\hat{\gamma})=E[(\hat{\gamma}-\gamma)^{2}]=% \operatorname{tr}(V(\hat{\gamma}))=\operatorname{tr}[(\Theta^{\rm T}\Theta)^{-% 1}\Theta^{\rm T}\Sigma(S)\Theta(\Theta^{\rm T}\Theta)^{-1}],=f(\gamma_{1},% \gamma_{2})$ (9)

with $V(\hat{\gamma})$ denoting the variance-covariance matrix of $\hat{\gamma}$ and $f(\gamma_{1},\gamma_{2})$ a function fo $\gamma_{1}$ and $\gamma_{2}$ . We will handle it later in this section. It is clear that the unique unknown matrix in $f(\gamma_{1},\gamma_{2})$ is $\Sigma(S)$ . In what follows we prove that such a matrix is a diagonal one.

Proposition 1. Recall $\lambda_{i}=\gamma_{1}\theta_{i}+\gamma_{2}$ , $i=1,\ldots,h$ . Then:

$\displaystyle\Sigma(S)=2\left[\begin{array}[]{@{}cccc@{}}\frac{\lambda_{1}^{2}% }{g_{1}}&0&\ldots&0\\ 0&\frac{\lambda_{2}^{2}}{g_{2}}&\ldots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\ldots&\frac{\lambda_{h}^{2}}{g_{h}}\end{array}\right].$ (10)

Proof Let $s_{i}^{2}=\frac{y_{i}^{\rm T}y_{i}}{g_{i}}$ , $i=1,\ldots,h$ . The proof follows if we notice the variance of $s_{i}^{2}$ , $i=1,\ldots,h$ , is given as

$\displaystyle\Sigma(s_{i}^{2})=\Sigma\left(\frac{y^{\rm T}A_{i}^{\rm T}A_{i}y}% {g_{i}}\right)=\frac{2}{g_{i}^{2}}\operatorname{tr}(\lambda_{i}^{2}I_{g_{i}})=% \frac{2\lambda^{2}}{g_{i}},$ (11)

and the cross-covariance between $s_{i}^{2}$ and $s_{j}^{2}$ , for $i\neq j$ , given as

$\displaystyle\Sigma(s_{i}^{2},s_{j}^{2})=\Sigma\left(\frac{y^{\rm T}A_{i}^{\rm T% }A_{i}y}{g_{i}},\frac{y^{\rm T}A_{j}^{\rm T}A_{j}}{g_{j}}\right)=2% \operatorname{tr}\left(\frac{A_{i}^{\rm T}A_{i}}{g_{i}}(\gamma_{1}M+\gamma_{2}% )\frac{A_{j}^{\rm T}A_{j}}{g_{j}}(\gamma_{1}M+\gamma_{2})\right)=0.$ (12)

Since the columns of $A_{i}^{\rm T}$ are the orthonormal eigenvectors of $M$ , the latest result is due to the following facts: $A_{i}MA_{j}^{\rm T}={0}_{g_{i},g_{j}}$ , $A_{i}^{\rm T}A_{j}={0}_{g_{i},g_{j}}$ and the commutativity of the trace operator. For the variance and cross-covariance calculation we used the Theorem 9.18 of Schott (1997).

Now we resume this section main result: the MSE of $\hat{\gamma}$ stated at Eq. (9).

Proposition 2.Let

$\displaystyle A=2h_{*}^{2}\left[h^{2}+\left(\sum_{i=1}^{h}\theta_{i}\right)^{2% }\right],B=4h_{*}^{2}\left[\sum_{i=1}^{h}h\theta_{i}+\left(\sum_{i=1}^{h}% \theta_{i}\right)^{2}\right]\text{ and }$ $\displaystyle C=4h_{*}^{2}\left(\sum_{i=1}^{h}\theta_{i}\right)^{2},\text{ % with }h_{*}=\frac{1}{h\sum_{i=1}^{h}\theta_{i}^{2}-\left(\sum_{i=1}^{h}\theta_% {i}\right)^{2}}.$

Then

$\displaystyle\text{MSE}(\hat{\gamma})=\operatorname{tr}[(\Theta^{\rm T}\Theta)% ^{-1}\Theta^{\rm T}\Sigma(S)\Theta(\Theta^{\rm T}\Theta)^{-1}].=\gamma_{1}^{2}% \left(A\sum_{i=1}^{h}\frac{\theta_{i}^{4}}{g_{i}}-B\sum_{i=1}^{h}\frac{\theta_% {i}^{3}}{g_{i}}+C\sum_{i=1}^{h}\frac{\theta_{i}^{2}}{g_{i}}\right)+2\gamma_{1}% \gamma_{2}\left(A\sum_{i=1}^{h}\frac{\theta_{i}^{3}}{g_{i}}-B\sum_{i=1}^{h}% \frac{\theta_{i}^{2}}{g_{i}}+C\sum_{i=1}^{h}\frac{\theta_{i}}{g_{i}}\right)+% \gamma_{2}^{2}\left(A\sum_{i=1}^{h}\frac{\theta_{i}^{2}}{g_{i}}-B\sum_{i=1}^{h% }\frac{\theta_{i}^{2}}{g_{i}}+C\sum_{i=1}^{h}\frac{1}{g_{i}}\right),$

(recall $h\sum_{i=1}^{h}\theta_{i}^{2}-\left(\sum_{i=1}^{h}\theta_{i}\right)^{2}\neq 0$ . See Silva and Fonseca (2018) for the proof).

Proof One should notice

$\displaystyle(\Theta^{\rm T}\Theta)^{-1}=h_{*}\left[\begin{array}[]{@{}cc@{}}h% &\displaystyle-\sum_{i=1}^{h}\theta_{i}\\ \displaystyle-\sum_{i=1}^{h}\theta_{i}&\displaystyle\sum_{i=1}^{h}\theta_{i}% \end{array}\right],$ $\displaystyle\Theta^{\rm T}\Sigma(S)\Theta=2\left[\begin{array}[]{@{}cc@{}}% \displaystyle\sum_{i=1}^{h}\frac{\lambda_{i}^{2}\theta_{i}^{2}}{g_{i}}&% \displaystyle\sum_{i=1}^{h}\frac{\lambda_{i}^{2}\theta_{i}}{g_{i}}\\ \displaystyle\sum_{i=1}^{h}\frac{\lambda_{i}^{2}\theta_{i}}{g_{i}}&% \displaystyle\sum_{i=1}^{h}\frac{\lambda_{i}^{2}}{g_{i}}\end{array}\right],$ (13)

with $h_{*}=\frac{1}{h\sum_{i=1}^{h}\theta_{i}^{2}-\left(\sum_{i=1}^{h}\theta_{i}% \right)^{2}}$ , and consequently

$\displaystyle(\Theta^{\rm T}\Theta)^{-1}\Theta^{\rm T}\Sigma(S)\Theta(\Theta^{% \rm T}\Theta)^{-1}=2h_{*}^{2}\left[\begin{array}[]{@{}cc@{}}b_{1}&b\\ b&b_{2}\end{array}\right]$

where

$\displaystyle b_{1}=h^{2}\sum_{i=1}^{h}\frac{\lambda_{i}^{2}\theta_{i}^{2}}{g_% {i}}-2h\sum_{i=1}^{h}\theta_{i}\sum_{i=1}^{h}\frac{\lambda_{i}^{2}\theta_{i}}{% g_{i}}+\left(\sum_{i=1}^{h}\theta_{i}\right)^{2}\sum_{i=1}^{h}\frac{\lambda_{i% }^{2}}{g_{i}},$ $\displaystyle b_{2}=\left(\sum_{i=1}^{h}\theta_{i}\right)^{2}\left[\sum_{i=1}^% {h}\frac{\lambda_{i}^{2}\theta_{i}^{2}}{g_{i}}-2\sum_{i=1}^{h}\frac{\lambda_{i% }^{2}\theta_{i}}{g_{i}}+\sum_{i=1}^{h}\frac{\lambda_{i}^{2}}{g_{i}}\right]% \text{ and }$ $\displaystyle b=-h\sum_{i=1}^{h}\theta\sum_{i=1}^{h}\frac{\lambda_{i}^{2}% \theta_{i}^{2}}{g_{i}}+\left(\left(\sum_{i=1}^{h}\theta_{i}\right)^{2}+h\sum_{% i=1}^{h}\theta_{i}\right)\sum_{i=1}^{h}\frac{\lambda_{i}^{2}\theta_{i}}{g_{i}}% -\left(\sum_{i=1}^{h}\theta_{i}\right)^{2}\sum_{i=1}^{h}\frac{\lambda_{i}^{2}}% {g_{i}}.$

According with Eq. (9) and results in Eq. (3), and notice $\lambda_{i}=\gamma_{1}\theta_{i}+\gamma_{2}$ , it follows that

$\displaystyle\text{MSE}(\hat{\gamma})=2h_{*}^{2}(b_{1}+b_{2}).$

Thus, after a straightforward calculation the proposed result follows.

If we denote $a_{t}=A\sum_{i=1}^{h}\frac{\theta_{i}^{t}}{g_{i}}-B\sum_{i=1}^{h}\frac{\theta_% {i}^{t-1}}{g_{i}}+C\sum_{i=1}^{h}\frac{\theta_{i}^{t-2}}{g_{i}}$ , for a natural number $t>1$ , with $\theta_{i}^{0}=0$ , the MSE of $\hat{\gamma}$ can be written as follows:

$\displaystyle\text{MSE}(\hat{\gamma})=a_{4}\gamma_{1}^{2}+2a_{3}\gamma_{1}% \gamma_{2}+a_{2}\gamma_{2}^{2}=\gamma^{\rm T}\left[\begin{array}[]{@{}cc@{}}a_% {4}&a_{3}\\ a_{3}&a_{2}\end{array}\right]\gamma=f(\gamma_{1},\gamma_{2}).$ (14)

Due to its non-negativity, for fixed $\gamma_{1}$ and $\gamma_{2}$ , the magnitude of $f(\gamma_{1},\gamma_{2})$ is directly connected to the parameter $a_{t}$ ; that is, the smaller the value of $a_{t}$ the smaller the value of the MSE of $\gamma$ .

4. Numerical example

Recall the numerical example in Silva (2018) in which its is considered $n_{i}$ , $i=1,\ldots,20$ , object produced on the ith of the 20 random chosen machines, with a continuous score recorded for each objects. The data (160 observations) was considered to be from the following random “one-way” designs:

$\displaystyle y_{ij}=\mu+\alpha_{i}+\epsilon_{ij},i=1,\ldots,20;j=1,\ldots,n_{% i},$ (15)

where $\mu$ is the overall mean, each $\alpha_{i}$ is the independent random effect due to machine $A_{i}$ and independent of the error $\epsilon_{ij}$ (the error for $j$ th observation fo the group $A_{i}$ ). The data includes $\sum_{i=1}^{20}n_{i}=160$ observations (sum of the observed objects for all the 20 machines).

The method Sub-D was used to infer about the variance for the machine’s effect as well as the one for the error effect; that is, the variance components for the considered design.

For such a model we have $h=18$ different eigenvalues:

$\displaystyle\theta_{1}=17.3568,\theta_{2}=15.1238,\theta_{3}=12.5421,$ $\displaystyle\theta_{4}=11.0000,\theta_{5}=10.5664,\theta_{6}=10.0000,$ $\displaystyle\theta_{7}=9.47630,\theta_{8}=9.00000,\theta_{9}=8.4603,$ $\displaystyle\theta_{10}=8.000\text{ (multiplicity 2)},\theta_{11}=7.1431,$ $\displaystyle\theta_{12}=5.2635,\theta_{13}=4.4000,$ $\displaystyle\theta_{14}=4.0000\text{ (multiplicity 2)},\theta_{15}=3.0883,$ $\displaystyle\theta_{16}=1.0925,\theta_{17}=1.0000,$ $\displaystyle\theta_{18}=0.0000\text{ (multiplicity 140)}.$

It was found that $g_{10}=2$ , $g_{18}=140$ and $g_{i}=r(A_{i})=1$ for all $i\neq 10$ and $i\neq 18$ . It is also found that $a2=$ 0.8475, $a3=$ 10,5947 and $a4=$ 143.5984. Thus, the MSE for Sub-D in the corresponding design is:

$\displaystyle\text{MSE}(\hat{\gamma})=143.5984\gamma_{1}+21.1894\gamma_{1}% \gamma_{2}+08{,}475\gamma_{2}.$

5. Final remark

Sub-D have been vastly tested in mixed linear models, specifically thought simulations in balanced and unbalanced “one-way” random designs, nested and crossed “two-way” designs and nested “three-way” designs, as may be checked on Silva (2017), Silva (2017) and Silva and Fonseca (2018). Although the simulations suggest Sub-D has in general better performance than the Anova-based and likelihood-based estimators, we have yet no tools wherewith infer over its efficiency, so that the deduction of its variance-covariance matrix as well as its mean square errors seems to be very timely. In fact, discussing the efficiency of an estimator (on its several aspects) goes through its variance-covariance and/or MSE.

As it may be seen, the Sub-D MSE structure that its greatness is intrinsically linked to the number of the different eigenvalues ( $h$ ) and the values of eigenvalues itself; that is the smaller the $h$ eigenvalues are, the smaller the values of $a^{t}$ ( $a_{2},a_{3},a_{4}$ ), and consequently the smaller the value of $\text{MSE}(\hat{\gamma})$ (see Eq. (14)).

The next issue will addressee the efficiency of Sub-D related to the ones of Anova-based and likelihood-based estimators.

Footnotes

Acknowledgments

This work is funded by National Funds through the FCT – Fundação para a Ciência e a Tecnologia, I.P., under the scope of the project UIDB/00297/2020 (Center for Mathematics and Applications).

References

Anderson

(1975). Designs and estimators for variance components. Statistical Design and linear Model, 1-30.

Anderson

, & Crump

(1967). Comparisons of designs and estimation procedures for estimating parameters in a two-stages nested process. Technometrics, 499-516.

Ferreira

Nunes

Fonseca

Silva

, & Mexia

(2017). Estimation and incommutativity in mixed models. Journal of Multivariate Analysis, 161, 58-67.

Klotz

Milton

, & Zacks

(2012). Mean square efficiency of estimators of variance components. Journal of the American Statistical Association, 64, 1383-1402.

Schott

(1997). Matrix analysis for statistics. John Wiley & Sons, New York.

Silva

(2017). Variance components estimation in mixed linear models. Ph.D Thesis.

Silva

, & Fonseca

(2018). The method Sub-D for variance componentes estimation in random one-way designs. International Journal of Multivariate Data Analysis, 1(4), 327-336.