# MANOVA for Nested Designs with Unequal Cell Sizes and Unequal Cell Covariance Matrices

Authors: Li-Wen Xu
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We propose and study parametric bootstrap (PB) tests for heteroscedastic two-factor MANOVA with nested designs. For the problem of testing “main effects” of both factors, we develop a flexible test based on a parametric bootstrap approach. The PB test is shown to be invariant under affine-transformations. Moreover, the PB test does not depend on the chosen weights used to define the parameters uniquely. The proposed test is compared with the approximate HotellingT2 (AHT) test by the simulations. Simulation results indicate that the PB test performs
satisfactorily for various cell sizes and parameter configurations and generally outperforms the AHT test in terms of controlling the nominal size. For the heteroscedastic cases, the PB test outperforms the AHT test in terms of power. In addition, the PB test does not lose too much power when the homogeneity assumption is actually valid.


---

## Body

## 1. Introduction

There are many situations where we record more than one response variable from each sampling or experimental unit and where these units are allocated to or occur in treatment groups. Ecologists often record the abundances of many species from each sampling or experimental unit and physiologists commonly measure more than one variable (e.g., blood pressure, heart rate, etc.) on experimental animals. With multiple response variables, we might be more interested in whether there are group differences on all the response variables considered simultaneously. This is the aim of multivariate analysis of variance (MANOVA), the analogue of univariate ANOVA when we have multiple response variables for each experimental or sampling unit.The simplest design is single factor MANOVA (or one-way MANOVA) which tests the effect of a factor, havingk levels. If the cell covariance matrices are assumed to be equal, then there are some popular tests available to test the equality of the mean vectors. The tests that are commonly used are the Wilks likelihood ratio (WLR), Lawley-Hotelling trace (LHT), Bartlett-Nanda-Pillai (BNP), and Roys largest root tests ([1], chap. 8, sec. 6). When there are some departures from the standard assumption, that is, unequal cell covariance matrices, these solutions were proposed by James [2], Johansen [3], Gamage et al. [4], and Krishnamoorthy and Lu [5], among others. The more complex design is multifactor MANOVA especially when the homogeneity of the cell covariance matrices assumption is seriously violated. There has been a continuous interest in the heteroscedastic two-factor MANOVA which tests in checking the significance of the effects of two factors A and B, each having a and b levels, respectively, in a multivariate factorial experiment with crossed designs. Recently, Harrar and Bathke [6] attacked this problem via modifying the WLR, LHT, and BNP tests. Their main ideas focus on modifying the degrees of freedom of the random matrices involved in the test statistics so that the heteroscedasticity of the cell covariance matrices is taken into account and the WLR, LHT, and BNP tests can still be used but with the degrees of freedom estimated from the data via matching the first two moments; see some details in Section  2.2 of Zhang and Xiao [7]. Zhang [8] proposed an approximate Hotelling T2 (AHT) test. A Wald-type test statistic is used. Its null distribution is approximated by a Hotelling T2-distribution with one parameter estimated from the data. Some simulation studies conducted in Zhang [8] showed that the AHT test outperforms the modified LHT test of Harrar and Bathke [6].Another useful two-factor design is the nested MANOVA [9]. Notice that each of the testing problems associated with the three null hypotheses in the heteroscedastic two-way MANOVA can be equivalently expressed in the form of the general linear hypothesis testing (GLHT) problem as described in Zhang [8]. The GLHT problem is very general. It includes not only the main and interaction effect tests but also various post hoc and contrast tests as special cases. To construct the test for two-factor nested MANOVA, we may use the same Wald-type test statistic as in Zhang [8]. However, according to our simulation studies, the empirical sizes of the AHT test for the two-factor nested MANOVA model may far exceed the nominal level. On the other hand, the AHT test needs to consider the the selected weights when the cell sizes are unequal. Therefore, it is important to develop a test procedure for the nesting effects and the nested effects with satisfactory size and power regardless of number of factorial effects and the sample sizes.Accordingly, the present paper will develop a parametric bootstrap (PB) test for heteroscedastic two-factor nested MANOVA. We use standardized effects sum of squares and a natural test statistic obtained by replacing cell covariance matrices by the corresponding sample cell covariance matrices. The PB test admits several nice properties: (1) it can be simply conducted by a routine Monte Carlo algorithm; (2) it is shown to be invariant under affine-transformations; (3) The PB test does not depend on choices of the weights used to identify the parameters; and (4) it works well. Simulation results reported in Section4 indicate that the PB test performs satisfactorily for various cell sizes and parameter configurations when the homogeneity assumption is seriously violated and generally outperforms the AHT test in terms of power and controlling size. The main ideas of the proposed PB test are closely related to the work by Xu et al. [10]. The methodologies for the PB test are presented in Sections 2 and 3. Simulation results are presented in Section 4. Some concluding remarks are given in Section 5. Finally, some technical proofs of the main results are given in the Appendix.


## 2. Tests for Nested Effects

Consider a two-factor nested design model with factorsA and B. Suppose that the nesting effect corresponds to the factor A having a factor levels, and the nested effect corresponds to the factor B having a total of b=∑i=1abi levels with bi levels of B nested within the ith level of the factor A  (i=1,…,a). Suppose a r-variate random sample of size nij is available under (i,j)th level, i=1,…,a; j=1,…,bi. Let Yijk, i=1,…,a; j=1,…,bi; k=1,…,nij represent these random vectors and yijk represent their observed (sample) values. Assume that nij>r so that positive definite sample covariance matrices can be computed for each cell of the design. Suppose that Yijk satisfy the following model:
(1)Yijk=μij+eijk,eijk~Nr(0,Σij),i=1,…,a;aaaaaaaaaaaaij=1,…,bi;k=1,…,nij,
where μij:r×1 and Σij:r×r are the cell mean vector and cell covariance matrix of the random sample at the (i,j)th cell and eijk is an experimental error vector term. All these samples are assumed to be independent with each other. In the two-factor nested model, the cell mean vectors μij are usually decomposed into the form
(2)μij=μ0+αi+βij,i=1,2,…,a;j=1,2,…,bi,
where μ0 is the grand mean vector, αi is the effect vector of the ith level of A on each of the r variables in Yijk, and βij is the effect due to the jth level of the factor B nested within the ith level of the factor A so that (1) can be further written as the following two-factor multivariate nested model with unequal error covariance matrices:
(3)Yijk=μ0+αi+βij+eijk,eijk~Nr(0,Σij),aii=1,…,a;j=1,…,bi;k=1,…,nij.
In order to have μ0, αi, and βij uniquely defined, we need to have additional constraints. Let u1,…,ua and vi1,…,vibi be nonnegative weights such that ∑i=1a‍ui>0 and ∑j=1bi‍vij>0. We consider the following constraints:
(4)∑i=1auiαi=0,∑j=1bivijβij=0,
where u1,…,ua and vij,…,vibi are nonnegative weights such that ∑i=1a‍ui>0 and ∑j=1bi‍vij>0 for each i.In this section, we are interested in testing the null hypothesis,(5)H0β:βij=0;i=1,…,a,j=1,…,bi,
against its natural alternative hypothesis. The null hypothesis in (5) aims to test if the nested effect corresponding to the factor B is statistically significant.The sample mean vector and the sample covariance matrix of the(i,j)th cell are denoted by Y¯ij and Sij, respectively, where
(6)Y¯ij=1nij∑k=1nijYijk,Sij=1nij-1∑k=1nij(Yijk-Y¯ij)(Yijk-Y¯ij)T,i=1,…,a; j=1,…,bi, and MT denotes the transpose of the matrix M. The observed values of these random variables are denoted as y¯ij and sij, respectively, i=1,…,a, j=1,…,bi. We note that Y¯ij’s and Sij’s are mutually independent with (nij-1)Sij~Wr(nij-1,Σij) and the model for Y¯ij as follows:
(7)Y¯ij=μ0+αi+βij+e¯ij,e¯ij~Nr(0,1nijΣij),aaaaaaaaaaaaaai=1,…,a;j=1,…,bi,
where e¯ij=(1/nij)∑k=1nij‍eijk and Wr(m,Γ) denotes the r-dimensional Wishart distribution with degrees of freedom (df)=m and scale parameter matrix Γ. Writing Y¯=(Y¯11T,…,Y¯1b1T,…,Y¯a1T,…,Y¯abaT)T and e¯=(e¯11T,…,e¯1b1T,…,e¯a1T,…,e¯abaT)T, the model (7) can be written as
(8)Y¯=(1b⊗Ir)μ0+(diag⁡(1b1,…,1ba)⊗Ir)α+β+e¯,
where α=(α1T,…,αaT)T, β=(β11T,…,β1b1T,…,βa1T,…,βabaT)T, e¯~Nbr(0,Σ), Σ=diag⁡((1/n11)Σ11,…,(1/n1b1)Σ1b1,…,(1/na1)Σa1,…,(1/naba)Σaba), 1m denotes the m×1 vector of ones, Im is an identity matrix with order m, V⊗W denotes the Kronecker product of matrices V and W, and diag⁡(M1,…,Mm) denotes a block-diagonal matrix with M1,…,Mm along the blocks.Define the standardized sum of squares due to the factorB(9)S~β(Y¯11,…,Y¯aba;Σ11,…,Σaba)=∑i=1a‍∑j=1bi‍(Y¯ij-μ^0-α^i)T(nijΣij-1)(Y¯ij-μ^0-α^i),
where μ^0 and α^i are solutions of μ0 and αi that minimize the quadratic equation
(10)S~(Y¯11,…,Y¯aba;Σ11,…,Σaba)=∑i=1a‍∑j=1bi‍(Y¯ij-μ0-αi)T(nijΣij-1)(Y¯ij-μ0-αi),
subject to the constraints given in (4).In fact, denotingθ=(μ0T,α1T,…,αaT)T and θ^=(μ^0T,α^1T,…,α^aT)T, it follows from Theorem 5.2.5 in Wang and Chow [11] that
(11)θ^=(μ^0T,α^1T,…,α^aT)T=(XTΣ-1X+LTL)-1XTΣ-1Y¯,
where X=(1b,diag⁡(1b1,…,1ba))⊗Ir and L=(0,u1,…,ua)⊗Ir. Then
(12)S~β(Y¯11,…,Y¯aba;Σ11,…,Σaba)=(Y¯-Xθ^)TΣ-1(Y¯-Xθ^).
Notice that θ^ indeed depends on L of chosen weights. However, we can prove that the standardized sum of squares in (9) does not depend on L. This is the following Theorem.Theorem 1.
The standardized sum of squaresS~β(Y¯11,…,Y¯aba;Σ11,…,Σaba) does not depend on L of chosen weights, and
(13)S~β(Y¯11,…,Y¯aba;Σ11,…,Σaba)=Y¯TΣ-1/2(Ibr-Σ-1/2X(XTΣ-1X)-XTΣ-1/2)Σ-1/2Y¯,
where M- denotes any generalized inverse of matrix M.Note from Theorem1 that S~β(Y¯11,…,Y¯aba;Σ11,…,Σaba) does not depend on L. If Σij are known, then a natural statistic for testing (5) is S~β(Y¯11,…,Y¯aba;Σ11,…,Σaba). In fact, Σ-1/2Y¯~Nbr(Σ-1/2μ,Ibr), and
(14)Δ=(Ibr-Σ-1/2X(XTΣ-1X)-XTΣ-1/2)
is an idempotent matrix with rank r(b-a); we have
(15)Y¯TΣ-1/2ΔΣ-1/2Y¯~χr(b-a)2(μTΣ-1/2ΔΣ-1/2μ),
where μ=(μ11T,…,μ1b1T,…,μa1T,…,μabaT)T and χm2(δ) denotes a noncentral chi-square random variable with degrees of freedom m and noncentrality parameter δ. The noncentrality parameter μTΣ-1/2ΔΣ-1/2μ is equal to zero when β11=⋯=βaba. Let y¯=(y¯11T,…,y¯1b1T,…,y¯a1T,…,y¯abaT)T be the observed value of Y¯. Then, the test that rejects H0β in (5) whenever
(16)S~β(Y¯11,…,Y¯aba;Σ11,…,Σaba)>χr(b-a),λ2
is a size λ test, where χm,λ2 is the upper λth quantile of a chi-square distribution with df=m.In general, the covariance matricesΣij are unknown; in this case, a test statistic can be obtained by replacing Σij in (13) by Sij, i=1,…,a,j=1,…,bi and is given by
(17)S~β(Y¯11,…,Y¯aba;S11,…,Saba)=Y¯TS-1/2(Ibr-S-1/2X(XTS-1X)-XTS-1/2)S-1/2Y¯,
where S=diag⁡((1/n11)S11,…,(1/n1b1)S1b1,…,(1/na1)Sa1,…,(1/naba)Saba).As shown in Zhang [8], the testing problem associated with the null hypothesis (5) can also be equivalently expressed in the form of the general linear hypothesis testing (GLHT) problem as H0β:Cμ=0, where C=HβAβ⊗Ir,
(18)Hβ=((Ib1-1,-1b1-1)0000(Ib2-1,-1b2-1)0000⋱0000(Iba-1,-1ba-1)),Aβ=(Ib1-1b1u1v1T0000Ib2-1b2u2v2T0000⋱0000Iba-1bauavaT),
and vi=(vi1,…,vibi)T, i=1,2,…,a. Then the associated Wald-type test statistic is given as
(19)T=(Cμ^)T(CSCT)-1(Cμ^),
where μ^=Y¯, and S is defined as in (17).In the following, we shall describe two tests based on (17) and (19), respectively.


### 2.1. The Approximate HotellingT2 (AHT) Test

Similar to Zhang [8], we proposed a test referred as the AHT test, which is based on the test statistic
(20)(d^-q+1)Tqd^,
where q=rank⁡(C)=r(b-a),
(21)d^=q(q+1)∑i=1a∑j=1bi(nij-1)-1[tr2(Ω^ij)+tr⁡(Ω^ij2)],
and Ω^ij=nij-1H^ijSijH^ijT and H^ij=(CSCT)-1/2Cij with C=[C11,…C1b1,…,Ca1,…,Caba], and Cij:q×r is the (ij)th block matrix of C. Zhang [8] showed that, under H0β, (d^-q+1)T/qd^ is approximately distributed as Fq,d^-q+1 random variable. Thus, the AHT test rejects the null hypothesis in (5) when the critical value (qd^/(d^-q+1))Fq,d^-q+1(1-λ) for the nominal significance level λ is exceeded by the observed test statistic T in (20). The AHT test can also be conducted via computing the P-value based on the approximate null distribution.


### 2.2. The Parametric Bootstrap (PB) Test

The parametric bootstrap involves sampling from the estimated models. That is, samples or sample statistics are generated from parametric models with the parameters replaced by their estimates. Recall that underH0β the vector Y¯ have the mean Xθ, where θ=(μ0T,α1T,…,αaT)T. As the test statistic in (17) is location invariant under the group of location transformations G={Y¯→Y¯+Xη,η∈Rr(a+1)}, without loss of generality, we can take Xθ=0. Using these facts, the parametric bootstrap pivot variable can be developed as follows. For a given (y¯11,…,y¯aba;s11,…,saba), let Y¯Bij~Nr(0,(1/nij)sij) and SBij~Wr(nij-1,(1/(nij-1))sij), i=1,…,a, j=1,…,bi. Then the PB pivot variable based on the test statistic (17) is given by
(22)S~βB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)=Y¯BTSB-1/2(Ibr-SB-1/2X(XTSB-1X)-XTSB-1/2)SB-1/2Y¯B,
where Y¯B=(Y¯B11T,…,Y¯B1b1T…,Y¯Ba1T,…,Y¯BabaT)T and SB=diag⁡((1/n11)SB11,…,(1/n1b1)SB1b,…,(1/na1)SBa1,…,(1/naba)SBaba). For a given level λ, the PB test rejects H0β in (5) when
(23)P(S~βB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)>s~β)<λ,
where
(24)s~β=S~β(y¯11,…,y¯aba;s11,…,saba)
is an observed value of S~β(Y¯11,…,Y¯aba;S11,…,Saba) in (17). For fixed (y¯11,…,y¯aba;s11,…,saba), the above probability does not depend on any unknown parameters, and so it can be estimated using Monte Carlo simulation given in the following Algorithm 2.Algorithm 2.
For a given(n11,…,naba), (y¯11,…,y¯aba), and (s11,…,saba),ComputeS~β(y¯11,…,y¯aba;s11,…,saba) in (17) and call it s~β.Forh=1,…,m,GenerateY¯Bij~Nr(0,(1/nij)sij) and SBij~Wr(nij-1,(1/(nij-1))sij),  i=1,…,a,j=1,…,bi,ComputeS~βB(Y¯B11,…,Y¯Baba;SB11,…,SBaba) using (22).IfS~βB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)>s~β, set Qh=1.(End loop)(1/m)∑h=1m‍Qh is a Monte Carlo estimate of the P value in (23).


### 2.3. Some Desirable Properties of the PB Test

The PB test (23) has several desirable invariance properties. First of all, the PB test is affine-invariant. That is, it is invariant under the following affine-transformation:
(25)Yijk*=DYijk+ξ,i=1,…,a;aaj=1,…,bi;k=1,…,nij,
where D is any nonsingular matrix and ξ is any given vector. This property is desirable since in practice, the cell responses Yijk (1) are often recentered or rescaled before an inference is conducted. The recentering and rescaling transformations are special cases of (25).Theorem 3.
The PB test (23) is affine-invariant in the sense that both the observed test statistic s~β (24) and the distribution of PB pivot variable S~βB (22) are affine-invariant.The model (3) is identifiable when the constraints (4) are imposed. In many situations there are no natural weights to justify a particular test procedure. For this, we have the following result.Theorem 4.
The PB test (23) is invariant on the different weight choices in constraints (4).Notice that the test statistic (17) does not depend on the chosen weights. The conclusion in Theorem 4 is obvious. However, we have the following property.Theorem 5.
The test statistic (17)
(26)S~β(Y¯11,…,Y¯aba;S11,…,Saba)=Y¯TS-1/2(Ibr-S-1/2X(XTS-1X)-XTS-1/2)S-1/2Y¯=Y¯TS-1/2(Ibr-S-1/2X(XTS-1X+LTL)-1XTS-1/2)S-1/2Y¯.
This property is desirable since in practice, in order to avoid computation of generalized inverse of matrices and enhance the accuracy of computation, we have taken a particular L such as L=(0,1/a,…,1/a)⊗Ir in the simulation study.


## 2.1. The Approximate HotellingT2 (AHT) Test

Similar to Zhang [8], we proposed a test referred as the AHT test, which is based on the test statistic
(20)(d^-q+1)Tqd^,
where q=rank⁡(C)=r(b-a),
(21)d^=q(q+1)∑i=1a∑j=1bi(nij-1)-1[tr2(Ω^ij)+tr⁡(Ω^ij2)],
and Ω^ij=nij-1H^ijSijH^ijT and H^ij=(CSCT)-1/2Cij with C=[C11,…C1b1,…,Ca1,…,Caba], and Cij:q×r is the (ij)th block matrix of C. Zhang [8] showed that, under H0β, (d^-q+1)T/qd^ is approximately distributed as Fq,d^-q+1 random variable. Thus, the AHT test rejects the null hypothesis in (5) when the critical value (qd^/(d^-q+1))Fq,d^-q+1(1-λ) for the nominal significance level λ is exceeded by the observed test statistic T in (20). The AHT test can also be conducted via computing the P-value based on the approximate null distribution.


## 2.2. The Parametric Bootstrap (PB) Test

The parametric bootstrap involves sampling from the estimated models. That is, samples or sample statistics are generated from parametric models with the parameters replaced by their estimates. Recall that underH0β the vector Y¯ have the mean Xθ, where θ=(μ0T,α1T,…,αaT)T. As the test statistic in (17) is location invariant under the group of location transformations G={Y¯→Y¯+Xη,η∈Rr(a+1)}, without loss of generality, we can take Xθ=0. Using these facts, the parametric bootstrap pivot variable can be developed as follows. For a given (y¯11,…,y¯aba;s11,…,saba), let Y¯Bij~Nr(0,(1/nij)sij) and SBij~Wr(nij-1,(1/(nij-1))sij), i=1,…,a, j=1,…,bi. Then the PB pivot variable based on the test statistic (17) is given by
(22)S~βB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)=Y¯BTSB-1/2(Ibr-SB-1/2X(XTSB-1X)-XTSB-1/2)SB-1/2Y¯B,
where Y¯B=(Y¯B11T,…,Y¯B1b1T…,Y¯Ba1T,…,Y¯BabaT)T and SB=diag⁡((1/n11)SB11,…,(1/n1b1)SB1b,…,(1/na1)SBa1,…,(1/naba)SBaba). For a given level λ, the PB test rejects H0β in (5) when
(23)P(S~βB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)>s~β)<λ,
where
(24)s~β=S~β(y¯11,…,y¯aba;s11,…,saba)
is an observed value of S~β(Y¯11,…,Y¯aba;S11,…,Saba) in (17). For fixed (y¯11,…,y¯aba;s11,…,saba), the above probability does not depend on any unknown parameters, and so it can be estimated using Monte Carlo simulation given in the following Algorithm 2.Algorithm 2.
For a given(n11,…,naba), (y¯11,…,y¯aba), and (s11,…,saba),ComputeS~β(y¯11,…,y¯aba;s11,…,saba) in (17) and call it s~β.Forh=1,…,m,GenerateY¯Bij~Nr(0,(1/nij)sij) and SBij~Wr(nij-1,(1/(nij-1))sij),  i=1,…,a,j=1,…,bi,ComputeS~βB(Y¯B11,…,Y¯Baba;SB11,…,SBaba) using (22).IfS~βB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)>s~β, set Qh=1.(End loop)(1/m)∑h=1m‍Qh is a Monte Carlo estimate of the P value in (23).


## 2.3. Some Desirable Properties of the PB Test

The PB test (23) has several desirable invariance properties. First of all, the PB test is affine-invariant. That is, it is invariant under the following affine-transformation:
(25)Yijk*=DYijk+ξ,i=1,…,a;aaj=1,…,bi;k=1,…,nij,
where D is any nonsingular matrix and ξ is any given vector. This property is desirable since in practice, the cell responses Yijk (1) are often recentered or rescaled before an inference is conducted. The recentering and rescaling transformations are special cases of (25).Theorem 3.
The PB test (23) is affine-invariant in the sense that both the observed test statistic s~β (24) and the distribution of PB pivot variable S~βB (22) are affine-invariant.The model (3) is identifiable when the constraints (4) are imposed. In many situations there are no natural weights to justify a particular test procedure. For this, we have the following result.Theorem 4.
The PB test (23) is invariant on the different weight choices in constraints (4).Notice that the test statistic (17) does not depend on the chosen weights. The conclusion in Theorem 4 is obvious. However, we have the following property.Theorem 5.
The test statistic (17)
(26)S~β(Y¯11,…,Y¯aba;S11,…,Saba)=Y¯TS-1/2(Ibr-S-1/2X(XTS-1X)-XTS-1/2)S-1/2Y¯=Y¯TS-1/2(Ibr-S-1/2X(XTS-1X+LTL)-1XTS-1/2)S-1/2Y¯.
This property is desirable since in practice, in order to avoid computation of generalized inverse of matrices and enhance the accuracy of computation, we have taken a particular L such as L=(0,1/a,…,1/a)⊗Ir in the simulation study.


## 3. Tests for the Nesting Effects

When the nested effects are present, the nesting effectαi can not reflect the effect of Ai because it depends on which level of factor B it is in. As pointed in Searle [12], Chap. 7, in the case of r=1, a popular solution to the problem is not quite a test for αi=0 in the presence of nested effects but rather to test the null hypothesis
(27)H0α∣β:αi+βij=0,i=1,…,a,j=1,…,bi.
We will see that the problem of testing (27) is the same case of the problem of comparing b=∑i=1a‍bi normal mean vectors when the cell covariance matrices are unknown and arbitrary [5]. In fact, if Σij are known, then
(28)μ^0=(∑i=1a‍∑j=1bi‍nijΣij-1)-1∑i=1a‍∑j=1binijΣij-1Y¯ij
is the best linear unbiased estimator of μ0, and a natural statistic for testing (27) is
(29)S~α(Y¯11,…,Y¯aba;Σ11,…,Σaba)=∑i=1a‍∑j=1bi(Y¯ij-μ^0)T(nijΣij-1)(Y¯ij-μ^0)=Y¯TΣ-1/2(Ibr-Σ-1/2X1(X1TΣ-1X1)-X1TΣ-1/2)Σ-1/2Y¯,
where X1=1b⊗Ir. Notice that Σ-1/2Y¯~Nbr(Σ-1/2μ,Ibr), and Δ1=(Ibr-Σ-1/2X1×(X1TΣ-1X1)-X1TΣ-1/2) is an idempotent matrix with rank r(b-1); we have
(30)Y¯TΣ-1/2Δ1Σ-1/2Y¯~χr(b-1)2(μTΣ-1/2Δ1Σ-1/2μ).
The noncentrality parameter μ′Σ-1/2Δ1Σ-1/2μ is equal to zero when αi+βij=0, i=1,…,a, j=1,…,bi. Then, the test that rejects H0α∣β in (27) whenever
(31)S~α(y¯11,…,y¯aba;Σ11,…,Σaba)>χr(b-1),λ2
is a size λ test. In general, the covariance matrices Σij are unknown; in this case, a test statistic can be obtained by replacing Σij in (29) by Sij, i=1,…,a, j=1,…,bi and is given by
(32)S~α(Y¯11,…,Y¯aba;S11,…,Saba)=Y¯TS-1/2(Ibr-S-1/2X1(X1TS-1X1)-X1TS-1/2)S-1/2Y¯.In the following, we describe the PB test forH0α∣β in (27). Recall that under H0α∣β the vector Y¯ have the mean X1μ0. As the test statistic in (32) is location invariant under the group of location transformations G1={Y¯+X1η,η∈Rr}, without loss of generality, we can take X1μ0=0. Using these facts, the parametric bootstrap pivot variable can be developed as follows. For a given (y¯11,…,y¯aba;s11,…,saba), let Y¯Bij~Nr(0,(1/nij)sij) and SBij~Wr(nij-1,(1/(nij-1))sij), i=1,…,a,j=1,…,bi. Then the PB pivot variable based on the test statistic (32) is given by
(33)S~αB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)=Y¯BTSB-1/2(Ibr-SB-1/2X1(X1TSB-1X)-X1TSB-1/2)SB-1/2Y¯B,
where Y¯B=(Y¯B11T,…,Y¯BabaT)T and SB=diag⁡((1/n11)SB11,…,(1/naba)SBaba). For a given level λ, the PB test rejects H0α∣β in (27) when
(34)P(S~αB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)>s~α)<λ,
where
(35)s~α=S~α(y¯11,…,y¯aba;s11,…,saba)
is an observed value of S~α(Y¯11,…,Y¯aba;S11,…,Saba) in (32). For fixed (y¯11,…,y¯aba;s11,…,saba), the above probability does not depend on any unknown parameters, and so it can be estimated using Monte Carlo simulation given in the following Algorithm 6.Algorithm 6.
For a given(n11,…,naba),(y¯11,…,y¯aba), and (s11,…,saba),ComputeS~α(y¯11,…,y¯aba;s11,…,saba) in (32) and call it s~α.Forh=1,…,m,GenerateY¯Bij~Nr(0,(1/nij)sij) and SBij~Wr(nij-1,(1/(nij-1))sij),i=1,…,a,j=1,…,bi,ComputeS~αB(Y¯B11,…,Y¯Baba;SB11,…,SBaba) using (33).IfS~αB(Y¯B11,…,Y¯Baba;SB11,…,SBaba)>s~α, set Qh=1.(End loop)(1/m)∑h=1m‍Qh is a Monte Carlo estimate of the P value in (34). Note that the PB test provided in this section has the same invariance properties as in Section 2.


## 4. Monte Carlo Studies

In this section, we need to compare the PB test with the AHT test via comparing their empirical sizes (Type I error rates) and powers for two-factor nested MANOVA models via simulations. As pointed out in Section3, the problem of testing (27) is the same case of the problem of comparing b=∑i=1a‍bi normal mean vectors when the cell covariance matrices are unknown and arbitrary. Due to this reason, we will look only at the nested effects for comparisons.Let the two factors beA and B, respectively. Suppose that the nesting effect of the factor A has a factor levels, and the nested effect of the factor B has a total of b=∑i=1a‍bi levels with bi levels of B nested within the ith level of the factor A. Let n=[n11,n12,…,naba] denote the vector of cell sizes. For given n and covariance matrices Σij, i=1,2,…,a, j=1,2,…,bi, we first generate ∑i=1a‍∑j=1bi‍nij multivariate samples as
(36)yijk=μij+Σij1/2eijk,k=1,2,…,nij,
where the cell mean vectors μij=μ11+ijδh/b with μ11 being the first cell mean vector, h a constant unit vector specifying the direction of the cell mean differences, and δ a tuning parameter controlling the amount of the cell mean differences. We independently generate the r entries of the error terms eijk from the N(0,1) distribution so that we always have E(eijk)=0 and Cov(eijk)=Ir. This means that (36) will generate the (ij)th multivariate normal sample yijk, k=1,2,…,nij with the given mean vector μij and covariance matrix Σij. Without loss of generality, we specify μ11 as 0 and h as h0/∥h0∥ where h0=[1,2,…,r]T for any given dimension r and ∥h0∥ denotes the usual L2-norm of h0. To estimate the sizes and powers of the AHT test, we used simulation consisting of 10,000 runs and recorded corresponding P values. Notice that two nested “do loops” are required to estimate the sizes and powers of the PB test, we used 2500 runs for outer “do loops” (for generating the data) and 5000 runs for inner “do loops” for estimating the bootstrap P values. The empirical sizes (when δ=0) and powers (when δ>0) of the two tests are the proportions of rejecting the null hypothesis, that is, when their P values are less than the nominal significance level λ. In all the simulations conducted, we used λ=5% for simplicity.Notice that the PB test does not depend on chosen weights. Here we report only the comparative studies for the equal-weight method to specify the weights of the AHT test so that the AHT test in Zhang [8] may be used for showing properties of the PB test. The empirical sizes and powers of the two tests for nested effect tests, together with the associated tuning parameters, are presented in Tables 1–3, in the columns labeled with AHT and PB under “δ=0” and “δ>0”, respectively. As seen from the three tables, three sets of the tuning parameters for the cell covariance matrices are examined, with the first set specifying the homogeneous cases; four sets of the cell sizes are specified, with the first two sets specifying the balanced cell size cases. To measure the overall performance of a test in terms of maintaining the nominal size λ, we use the average relative error defined in Zhang [8] as ARE =M-1∑m=1M‍|λ^m-λ|/λ×100 where λ^m denotes the mth empirical size for m=1,2,…,M, λ=0.05, and M is the number of empirical sizes under consideration. The smaller ARE value indicates the better overall performance of the associated test. Usually, when ARE≤10, the test performs very well; when 10<ARE≤20, the test performs reasonably well; and when ARE>20, the test does not perform well since its empirical sizes are either too liberal or too conservative and hence may be unacceptable. Notice that for a good test, the larger the cell sizes, the smaller the ARE values. The ARE values of the two tests under the two error schemes are also presented in these three tables. Notice that for simplicity, in the specification of the covariance and size tuning parameters, we often use (ut) to denote “u repeats t times.” Table 1 shows the empirical sizes and powers of the two tests for a bivariate case with a=2 and b1=16, b2=20. With b=b1+b2=36, one may be able to check how the two tests behave when one of the factors has a large number of levels. Tables 2 and 3 show the empirical sizes and powers of the two tests for a three-variate case with a=3 and b1=8, b2=10, b3=12, and a 10-variate case with a=3 and b1=3, b2=5, b3=7, respectively. These two tables allow us to compare the two tests for higher-dimensional data.Table 1
Empirical sizes and powers of the two tests for nested effects for bivariate two-factor multivariate nested designs [a=2,b1=16,b2=20, (Σ11,Σ12,…,Σaba)=118T⊗(I2,diag⁡(λ))].
λ
n
δ
=
0
δ
=
1
δ
=
2
δ
=
3
AHT
PB
AHT
PB
AHT
PB
AHT
PB
λ
1
n
1
0.056
0.046
0.210
0.162
0.831
0.793
0.999
0.999
n
2
0.057
0.050
0.394
0.372
0.995
0.994
1.000
1.000
n
3
0.070
0.052
0.295
0.232
0.957
0.930
1.000
1.000
n
4
0.057
0.054
0.926
0.923
1.000
1.000
1.000
1.000
λ
2
n
1
0.058
0.055
0.108
0.124
0.379
0.551
0.826
0.964
n
2
0.068
0.044
0.173
0.232
0.718
0.930
0.995
1.000
n
3
0.071
0.053
0.152
0.150
0.549
0.717
0.958
0.995
n
4
0.059
0.051
0.475
0.776
0.999
1.000
1.000
1.000
λ
3
n
1
0.056
0.045
0.093
0.114
0.301
0.507
0.722
0.943
n
2
0.058
0.046
0.151
0.204
0.614
0.878
0.977
0.999
n
3
0.070
0.048
0.131
0.145
0.458
0.654
0.903
0.990
n
4
0.061
0.055
0.397
0.755
0.994
1.000
1.000
1.000
ARE
24.03
6.83
λ
1
=
(
1,1
), λ2=(1,5), and λ3=(1,10). n1=(7,7)20, n2=(10,10)20, n3=(7,10)20, and n4=(30,15)20.Table 2
Empirical sizes and powers of the two tests for nested effects for 3-variate two-factor multivariate nested designs [a=3,b1=8,b2=10,b3=12, (Σ11,Σ12,…,Σaba)=110T⊗(I3,diag⁡(λ),(1-ρ)I3+ρ1313T)].
(
λ
,
ρ
)
n
δ
=
0
δ
=
0.6
δ
=
1.2
δ
=
1.8
AHT
PB
AHT
PB
AHT
PB
AHT
PB
(
λ
1
,
ρ
1
)
n
1
0.052
0.054
0.099
0.103
0.369
0.353
0.840
0.815
n
2
0.056
0.041
0.164
0.152
0.705
0.691
0.994
0.991
n
3
0.092
0.041
0.302
0.194
0.944
0.891
1.000
1.000
n
4
0.096
0.044
0.288
0.185
0.927
0.847
1.000
1.000
(
λ
2
,
ρ
2
)
n
1
0.054
0.049
0.181
0.213
0.756
0.860
0.997
1.000
n
2
0.055
0.048
0.335
0.422
0.985
0.997
1.000
1.000
n
3
0.090
0.050
0.629
0.519
1.000
1.000
1.000
1.000
n
4
0.089
0.061
0.522
0.545
0.999
1.000
1.000
1.000
(
λ
3
,
ρ
3
)
n
1
0.049
0.050
0.271
0.194
0.727
0.852
0.995
1.000
n
2
0.052
0.056
0.308
0.427
0.978
0.998
1.000
1.000
n
3
0.094
0.060
0.593
0.519
1.000
0.999
1.000
1.000
n
4
0.098
0.057
0.486
0.551
0.999
1.000
1.000
1.000
ARE
46.67
10.83
(
λ
1
,
ρ
1
)
=
(
1
3
,
0
), (λ2,ρ2)=(1,15,0.1,0.1), and (λ3,ρ3)=(1,10,0.1,0.5). n1=(10,10,10)10, n2=(15,15,15)10, n3=(10,20,40)10, and n4=(40,20,10)10.Table 3
Empirical sizes and powers of the two tests for nested effects for 10-variate two-factor multivariate nested designs [a=3,b1=3,b2=5,b3=7, (Σ11,Σ12,…,Σaba)=110T⊗(I3,diag⁡(λ),diag⁡(η))].
(
λ
,
η
)
n
δ
=
0
δ
=
1
δ
=
2
δ
=
3
AHT
PB
AHT
PB
AHT
PB
AHT
PB
(
λ
1
,
η
1
)
n
1
0.046
0.042
0.438
0.415
0.999
0.998
1.000
1.000
n
2
0.052
0.050
0.937
0.932
1.000
1.000
1.000
1.000
n
3
0.062
0.046
0.788
0.762
1.000
1.000
1.000
1.000
n
4
0.058
0.056
0.871
0.846
1.000
1.000
1.000
1.000
(
λ
2
,
η
2
)
n
1
0.052
0.044
0.842
0.635
1.000
1.000
1.000
1.000
n
2
0.052
0.047
1.000
0.996
1.000
1.000
1.000
1.000
n
3
0.067
0.049
0.995
0.965
1.000
1.000
1.000
1.000
n
4
0.059
0.046
0.999
0.972
1.000
1.000
1.000
1.000
(
λ
3
,
η
3
)
n
1
0.049
0.042
0.066
0.101
0.115
0.401
0.235
0.877
n
2
0.050
0.054
0.090
0.236
0.255
0.930
0.660
1.000
n
3
0.071
0.045
0.089
0.145
0.183
0.622
0.442
0.982
n
4
0.055
0.047
0.095
0.236
0.260
0.961
0.662
1.000
ARE
13.78
8.67
λ
1
=
(
1
10
)
5, η1=(110)5; λ2=(123,13,243,1)5, η2=(13,0.13,22,24,21)5; and λ3=(13,33,93,20)5, η3=(53,153,453,100)5. n1=(253)5, n2=(503)5, n3=(25,35,50)5, and n4=(70,40,35)5.In the following, let us compare the AHT and PB tests via examining their empirical sizes and powers. It is seen from the three tables that for all the cases, the PB test generally outperforms the AHT test for all the cases under consideration as shown by their empirical sizes and the associated ARE values presented in the three tables. The AHT test performs reasonably well only for 10-dimensional data. In the homogeneous cases, the AHT test appears to be more powerful than the PB tests because of its inflated empirical sizes exceeding the nominal level considerably. For the heteroscedastic cases, we once again observe from Tables2 and 3 that the PB test performs superior to the AHT test in terms of controlling nominal level and powers.We conclude this simulation section via summarizing all the simulation studies conducted. In terms of size, the overall conclusion is that the PB test is a flexible procedure that performs satisfactorily, regardless of the cell sizes, values of the cell covariance matrices, and the number of effects being compared. In terms of power, the PB test generally outperforms the AHT test for most heteroscedastic cases under consideration. Thus, one may recommend to use the PB test as a good alternative in practical applications because of its simplicity and accuracy.


## 5. Concluding Remarks

The available classical tests for the two-factor MANOVA model with crossed designs and heteroscedastic cell covariance matrices have serious Type I error problems that have been overlooked. To address this serious problem, Zhang [8] proposed the AHT test. In this paper, we proposed and studied the PB test and the AHT test for the two-factor MANOVA model with nested designs and heteroscedastic cell covariance matrices. We showed that the PB test is invariant under affine-transformations and different choices of weights used to define the parameters uniquely. We demonstrated via intensive simulations that the PB test generally performs well and outperforms the AHT test in terms of size and power for most cell sizes and parameter configurations.



---

*Source: 649202-2014-07-16.xml*