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

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 Hotelling 𝑇² (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.

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, having 𝑘 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 𝐴 and 𝐵, each having 𝑎 and 𝑏 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 𝑇² (AHT) test. A Wald-type test statistic is used. Its null distribution is approximated by a Hotelling 𝑇²-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 Section 4 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 factors 𝐴 and 𝐵. Suppose that the nesting effect corresponds to the factor 𝐴 having 𝑎 factor levels, and the nested effect corresponds to the factor 𝐵 having a total of 𝑏=∑^𝑎_𝑖=1𝑏_𝑖 levels with 𝑏_𝑖 levels of 𝐵 nested within the 𝑖th level of the factor 𝐴  (𝑖=1,…,𝑎). Suppose a 𝑟-variate random sample of size 𝑛_𝑖𝑗 is available under (𝑖,𝑗)th level, 𝑖=1,…,𝑎; 𝑗=1,…,𝑏_𝑖. Let 𝐘_𝑖𝑗𝑘, 𝑖=1,…,𝑎; 𝑗=1,…,𝑏_𝑖; 𝑘=1,…,𝑛_𝑖𝑗 represent these random vectors and 𝐲_𝑖𝑗𝑘 represent their observed (sample) values. Assume that 𝑛_𝑖𝑗>𝑟 so that positive definite sample covariance matrices can be computed for each cell of the design. Suppose that 𝐘_𝑖𝑗𝑘 satisfy the following model: 𝐘_𝑖𝑗𝑘=𝝁_𝑖𝑗+𝐞_𝑖𝑗𝑘,𝐞_𝑖𝑗𝑘∼𝑁_𝑟𝟎,𝚺_𝑖𝑗,𝑖=1,…,𝑎;𝑗=1,…,𝑏_𝑖;𝑘=1,…,𝑛_𝑖𝑗, (1) where 𝝁_𝑖𝑗∶𝑟×1 and 𝚺_𝑖𝑗∶𝑟×𝑟 are the cell mean vector and cell covariance matrix of the random sample at the (𝑖,𝑗)th cell and 𝐞_𝑖𝑗𝑘 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 𝝁_𝑖𝑗 are usually decomposed into the form 𝝁_𝑖𝑗=𝝁₀+𝜶_𝑖+𝜷_𝑖𝑗,𝑖=1,2,…,𝑎;𝑗=1,2,…,𝑏_𝑖, (2) where 𝝁₀ is the grand mean vector, 𝜶_𝑖 is the effect vector of the 𝑖th level of 𝐴 on each of the 𝑟 variables in 𝐘_𝑖𝑗𝑘, and 𝜷_𝑖𝑗 is the effect due to the 𝑗th level of the factor 𝐵 nested within the 𝑖th level of the factor 𝐴 so that (1) can be further written as the following two-factor multivariate nested model with unequal error covariance matrices: 𝐘_𝑖𝑗𝑘=𝝁₀+𝜶_𝑖+𝜷_𝑖𝑗+𝐞_𝑖𝑗𝑘,𝐞_𝑖𝑗𝑘∼𝑁_𝑟𝟎,𝚺_𝑖𝑗,𝑖=1,…,𝑎;𝑗=1,…,𝑏_𝑖;𝑘=1,…,𝑛_𝑖𝑗. (3) In order to have 𝝁₀, 𝜶_𝑖, and 𝜷_𝑖𝑗 uniquely defined, we need to have additional constraints. Let 𝑢₁,…,𝑢_𝑎 and 𝑣_𝑖1,…,𝑣_𝑖𝑏𝑖 be nonnegative weights such that ∑^𝑎_𝑖=1𝑢_𝑖>0 and ∑^𝑏𝑖_𝑗=1𝑣_𝑖𝑗>0. We consider the following constraints: ^𝑎_𝑖=1𝑢_𝑖𝜶_𝑖=𝟎,^𝑏𝑖_𝑗=1𝑣_𝑖𝑗𝜷_𝑖𝑗=𝟎, (4) where 𝑢₁,…,𝑢_𝑎 and 𝑣_𝑖𝑗,…,𝑣_𝑖𝑏𝑖 are nonnegative weights such that ∑^𝑎_𝑖=1𝑢_𝑖>0 and ∑^𝑏𝑖_𝑗=1𝑣_𝑖𝑗>0 for each 𝑖.

In this section, we are interested in testing the null hypothesis, 𝐻_0𝜷∶𝜷_𝑖𝑗=𝟎;𝑖=1,…,𝑎,𝑗=1,…,𝑏_𝑖, (5) against its natural alternative hypothesis. The null hypothesis in (5) aims to test if the nested effect corresponding to the factor 𝐵 is statistically significant.

The sample mean vector and the sample covariance matrix of the (𝑖,𝑗)th cell are denoted by 𝐘_𝑖𝑗 and 𝐒_𝑖𝑗, respectively, where 𝐘_𝑖𝑗=1𝑛^𝑛_𝑖𝑗^𝑖𝑗_𝑘=1𝐘_𝑖𝑗𝑘,𝐒_𝑖𝑗=1𝑛_𝑖𝑗−1^𝑛𝑖𝑗_𝑘=1𝐘_𝑖𝑗𝑘−𝐘_𝑖𝑗𝐘_𝑖𝑗𝑘−𝐘_𝑖𝑗^𝑇, (6)𝑖=1,…,𝑎; 𝑗=1,…,𝑏_𝑖, and 𝐌^𝑇 denotes the transpose of the matrix 𝐌. The observed values of these random variables are denoted as 𝐲_𝑖𝑗 and 𝐬_𝑖𝑗, respectively, 𝑖=1,…,𝑎, 𝑗=1,…,𝑏_𝑖. We note that 𝐘_𝑖𝑗’s and 𝐒_𝑖𝑗’s are mutually independent with (𝑛_𝑖𝑗−1)𝐒_𝑖𝑗∼𝑊_𝑟(𝑛_𝑖𝑗−1,𝚺_𝑖𝑗) and the model for 𝐘_𝑖𝑗 as follows: 𝐘_𝑖𝑗=𝝁₀+𝜶_𝑖+𝜷_𝑖𝑗+𝐞_𝑖𝑗,𝐞_𝑖𝑗∼𝑁_𝑟10,𝑛_𝑖𝑗𝚺_𝑖𝑗,𝑖=1,…,𝑎;𝑗=1,…,𝑏_𝑖, (7) where 𝐞_𝑖𝑗=(1/𝑛_𝑖𝑗)∑^𝑛𝑖𝑗_𝑘=1𝐞_𝑖𝑗𝑘 and 𝑊_𝑟(𝑚,𝚪) denotes the 𝑟-dimensional Wishart distribution with degrees of freedom (df)=𝑚 and scale parameter matrix 𝚪. Writing 𝐘=(𝐘^𝑇₁₁,…,𝐘^𝑇_1𝑏1,…,𝐘^𝑇_𝑎1,…,𝐘^𝑇_𝑎𝑏𝑎)^𝑇 and 𝐞=(𝐞^𝑇₁₁,…,𝐞^𝑇_1𝑏1,…,𝐞^𝑇_𝑎1,…,𝐞^𝑇_𝑎𝑏𝑎)^𝑇, the model (7) can be written as 𝐘=𝟏_𝑏⊗𝐈_𝑟𝝁₀+diag𝟏_𝑏1,…,𝟏_𝑏𝑎⊗𝐈_𝑟𝜶+𝜷+𝐞, (8) where 𝜶=(𝜶^𝑇₁,…,𝜶^𝑇_𝑎)^𝑇, 𝜷=(𝜷^𝑇₁₁,…,𝜷^𝑇_1𝑏1,…,𝜷^𝑇_𝑎1,…,𝜷^𝑇_𝑎𝑏𝑎)^𝑇, 𝐞∼𝑁_𝑏𝑟(𝟎,𝚺), 𝚺=diag((1/𝑛₁₁)𝚺₁₁,…,(1/𝑛_1𝑏1)𝚺_1𝑏1,…,(1/𝑛_𝑎1)𝚺_𝑎1,…,(1/𝑛_𝑎𝑏𝑎)𝚺_𝑎𝑏𝑎), 𝟏_𝑚 denotes the 𝑚×1 vector of ones, 𝐈_𝑚 is an identity matrix with order 𝑚, 𝐕⊗𝐖 denotes the Kronecker product of matrices 𝐕 and 𝐖, and diag(𝐌₁,…,𝐌_𝑚) denotes a block-diagonal matrix with 𝐌₁,…,𝐌_𝑚 along the blocks.

Define the standardized sum of squares due to the factor 𝐵𝑆_𝜷𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎=^𝑎𝑏_𝑖=1^𝑖_𝑗=1𝐘_𝑖𝑗−𝝁₀−𝜶_𝑖^𝑇𝑛_𝑖𝑗𝚺⁻¹_𝑖𝑗𝐘_𝑖𝑗−𝝁₀−𝜶_𝑖, (9) where 𝝁₀ and 𝜶_𝑖 are solutions of 𝝁₀ and 𝜶_𝑖 that minimize the quadratic equation 𝑆𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎=^𝑎𝑏_𝑖=1^𝑖_𝑗=1𝐘_𝑖𝑗−𝝁₀−𝜶_𝑖^𝑇𝑛_𝑖𝑗𝚺⁻¹_𝑖𝑗𝐘_𝑖𝑗−𝝁₀−𝜶_𝑖, (10) subject to the constraints given in (4).

In fact, denoting 𝜽=(𝝁^𝑇₀,𝜶^𝑇₁,…,𝜶^𝑇_𝑎)^𝑇 and 𝜽=(𝝁^𝑇₀,𝜶^𝑇₁,…,𝜶^𝑇_𝑎)^𝑇, it follows from Theorem 5.2.5 in Wang and Chow [11] that 𝜽=𝝁^𝑇₀,𝜶^𝑇₁,…,𝜶^𝑇_𝑎^𝑇=𝐗^𝑇𝚺⁻¹𝐗+𝐋^𝑇𝐋⁻¹𝐗^𝑇𝚺⁻¹𝐘, (11) where 𝐗=(𝟏_𝑏,diag(𝟏_𝑏1,…,𝟏_𝑏𝑎))⊗𝐈_𝑟 and 𝐋=(0,𝑢₁,…,𝑢_𝑎)⊗𝐈_𝑟. Then 𝑆_𝜷𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎=𝐘−𝐗𝜽^𝑇𝚺⁻¹𝐘−𝐗𝜽. (12) Notice that 𝜽 indeed depends on 𝐋 of chosen weights. However, we can prove that the standardized sum of squares in (9) does not depend on 𝐋. This is the following Theorem.

Theorem 1. The standardized sum of squares 𝑆_𝜷(𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎) does not depend on 𝐋 of chosen weights, and 𝑆_𝜷𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎=𝐘^𝑇𝚺^−1/2𝐈_𝑏𝑟−𝚺^−1/2𝐗𝐗^𝑇𝚺⁻¹𝐗⁻𝐗^𝑇𝚺^−1/2𝚺^−1/2𝐘, (13) where 𝐌⁻ denotes any generalized inverse of matrix 𝐌.

Note from Theorem 1 that 𝑆_𝜷(𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎) does not depend on 𝐋. If 𝚺_𝑖𝑗 are known, then a natural statistic for testing (5) is 𝑆_𝜷(𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎). In fact, 𝚺^−1/2𝐘∼𝑁_𝑏𝑟(𝚺^−1/2𝝁,𝐈_𝑏𝑟), and 𝚫=𝐈_𝑏𝑟−𝚺^−1/2𝐗𝐗^𝑇𝚺⁻¹𝐗⁻𝐗^𝑇𝚺^−1/2 (14) is an idempotent matrix with rank 𝑟(𝑏−𝑎); we have 𝐘^𝑇𝚺^−1/2𝚫𝚺^−1/2𝐘∼𝜒²_{𝑟(𝑏−𝑎)}𝝁^𝑇𝚺^−1/2𝚫𝚺^−1/2𝝁, (15) where 𝝁=(𝝁^𝑇₁₁,…,𝝁^𝑇_1𝑏1,…,𝝁^𝑇_𝑎1,…,𝝁^𝑇_𝑎𝑏𝑎)^𝑇 and 𝜒²_𝑚(𝛿) denotes a noncentral chi-square random variable with degrees of freedom 𝑚 and noncentrality parameter 𝛿. The noncentrality parameter 𝝁^𝑇𝚺^−1/2𝚫𝚺^−1/2𝝁 is equal to zero when 𝜷₁₁=⋯=𝜷_𝑎𝑏𝑎. Let 𝐲=(𝐲^𝑇₁₁,…,𝐲^𝑇_1𝑏1,…,𝐲^𝑇_𝑎1,…,𝐲^𝑇_𝑎𝑏𝑎)^𝑇 be the observed value of 𝐘. Then, the test that rejects 𝐻_0𝜷 in (5) whenever 𝑆_𝜷𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎>𝜒²_{𝑟(𝑏−𝑎),𝜆} (16) is a size 𝜆 test, where 𝜒²_𝑚,𝜆 is the upper 𝜆th quantile of a chi-square distribution with 𝑑𝑓=𝑚.

In general, the covariance matrices 𝚺_𝑖𝑗 are unknown; in this case, a test statistic can be obtained by replacing 𝚺_𝑖𝑗 in (13) by 𝐒_𝑖𝑗, 𝑖=1,…,𝑎,𝑗=1,…,𝑏_𝑖 and is given by 𝑆_𝜷𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝐒₁₁,…,𝐒_𝑎𝑏𝑎=𝐘^𝑇𝐒^−1/2𝐈_𝑏𝑟−𝐒^−1/2𝐗𝐗^𝑇𝐒⁻¹𝐗⁻𝐗^𝑇𝐒^−1/2𝐒^−1/2𝐘, (17) where 𝐒=diag((1/𝑛₁₁)𝐒₁₁,…,(1/𝑛_1𝑏1)𝐒_1𝑏1,…,(1/𝑛_𝑎1)𝐒_𝑎1,…,(1/𝑛_𝑎𝑏𝑎)𝐒_𝑎𝑏𝑎).

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 𝐻_0𝜷∶𝐂𝝁=0, where 𝐂=𝐇_𝜷𝐀_𝜷⊗𝐈_𝑟, 𝐇_𝜷=𝐈_𝑏1−1,−𝟏_𝑏1−10𝐈000_𝑏2−1,−𝟏_𝑏2−1𝐈0000⋱0000_𝑏𝑎−1,−𝟏_𝑏𝑎−1,𝐀_𝜷=𝐈_𝑏1−𝟏_𝑏1𝑢₁𝐯^𝑇₁0000𝐈_𝑏2−𝟏_𝑏2𝑢₂𝐯^𝑇₂0000⋱0000𝐈_𝑏𝑎−𝟏_𝑏𝑎𝑢_𝑎𝐯^𝑇_𝑎, (18) and 𝐯_𝑖=(𝑣_𝑖1,…,𝑣_𝑖𝑏𝑖)^𝑇, 𝑖=1,2,…,𝑎. Then the associated Wald-type test statistic is given as 𝑇=(𝐂𝝁)^𝑇𝐂𝐒𝐂^𝑇−1(𝐂𝝁), (19) where 𝝁=𝐘, and 𝐒 is defined as in (17).

In the following, we shall describe two tests based on (17) and (19), respectively.

2.1. The Approximate Hotelling 𝑇² (AHT) Test

Similar to Zhang [8], we proposed a test referred as the AHT test, which is based on the test statistic 𝑑−𝑞+1𝑇𝑞𝑑, (20) where 𝑞=rank(𝐂)=𝑟(𝑏−𝑎), 𝑞𝑑=𝑞+1∑^𝑎_𝑖=1∑^𝑏𝑖_𝑗=1𝑛_𝑖𝑗−1⁻¹tr²Ω_𝑖𝑗+trΩ²_𝑖𝑗, (21) and Ω_𝑖𝑗=𝑛⁻¹_𝑖𝑗𝐇_𝑖𝑗𝐒_𝑖𝑗𝐇^𝑇_𝑖𝑗 and 𝐇_𝑖𝑗=(𝐂𝐒𝐂^𝑇)^−1/2𝐂_𝑖𝑗 with 𝐂=[𝐂₁₁,…𝐂_1𝑏1,…,𝐂_𝑎1,…,𝐂_𝑎𝑏𝑎], and 𝐂_𝑖𝑗∶𝑞×𝑟 is the (𝑖𝑗)th block matrix of 𝐂. Zhang [8] showed that, under 𝐻_0𝜷, (𝑑−𝑞+1)𝑇/𝑞𝑑 is approximately distributed as 𝐹_{𝑞,𝑑−𝑞+1} random variable. Thus, the AHT test rejects the null hypothesis in (5) when the critical value (𝑞𝑑/(𝑑−𝑞+1))𝐹_{𝑞,𝑑−𝑞+1}(1−𝜆) for the nominal significance level 𝜆 is exceeded by the observed test statistic 𝑇 in (20). The AHT test can also be conducted via computing the 𝑃-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 under 𝐻_0𝜷 the vector 𝐘 have the mean 𝐗𝜽, where 𝜽=(𝝁^𝑇₀,𝜶^𝑇₁,…,𝜶^𝑇_𝑎)^𝑇. As the test statistic in (17) is location invariant under the group of location transformations 𝒢={𝐘→𝐘+𝐗𝜼,𝜼∈ℝ^𝑟(𝑎+1)}, without loss of generality, we can take 𝐗𝜽=𝟎. Using these facts, the parametric bootstrap pivot variable can be developed as follows. For a given (𝐲₁₁,…,𝐲_𝑎𝑏𝑎;𝐬₁₁,…,𝐬_𝑎𝑏𝑎), let 𝐘_𝐵𝑖𝑗∼𝑁_𝑟(0,(1/𝑛_𝑖𝑗)𝐬_𝑖𝑗) and 𝐒_𝐵𝑖𝑗∼𝑊_𝑟(𝑛_𝑖𝑗−1,(1/(𝑛_𝑖𝑗−1))𝐬_𝑖𝑗), 𝑖=1,…,𝑎, 𝑗=1,…,𝑏_𝑖. Then the PB pivot variable based on the test statistic (17) is given by 𝑆_𝜷𝐁𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎}=𝐘^𝑇_𝐵𝐒_𝐵^−1/2𝐈_𝑏𝑟−𝐒_𝐵^−1/2𝐗𝐗^𝑇𝐒_𝐵⁻¹𝐗⁻𝐗^𝑇𝐒_𝐵^−1/2𝐒_𝐵^−1/2𝐘_𝐵, (22) where 𝐘_𝐵=(𝐘^𝑇_𝐵11,…,𝐘^𝑇_𝐵1𝑏1…,𝐘^𝑇_𝐵𝑎1,…,𝐘^𝑇_{𝐵𝑎𝑏𝑎})^𝑇 and 𝐒_𝐵=diag((1/𝑛₁₁)𝐒_𝐵11,…,(1/𝑛_1𝑏1)𝐒_𝐵1𝑏,…,(1/𝑛_𝑎1)𝐒_𝐵𝑎1,…,(1/𝑛_𝑎𝑏𝑎)𝐒_{𝐵𝑎𝑏𝑎}). For a given level 𝜆, the PB test rejects 𝐻_0𝜷 in (5) when 𝑃𝑆_𝜷𝐵𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎}>𝑠_𝜷̃<𝜆, (23) where 𝑠_𝜷̃=𝑆_𝜷𝐲₁₁,…,𝐲_𝑎𝑏𝑎;𝐬₁₁,…,𝐬_𝑎𝑏𝑎 (24) is an observed value of 𝑆_𝜷(𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝐒₁₁,…,𝐒_𝑎𝑏𝑎) in (17). For fixed (𝐲₁₁,…,𝐲_𝑎𝑏𝑎;𝐬₁₁,…,𝐬_𝑎𝑏𝑎), 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 (𝑛₁₁,…,𝑛_𝑎𝑏𝑎), (𝐲₁₁,…,𝐲_𝑎𝑏𝑎), and (𝐬₁₁,…,𝐬_𝑎𝑏𝑎), Compute 𝑆_𝜷(𝐲₁₁,…,𝐲_𝑎𝑏𝑎;𝐬₁₁,…,𝐬_𝑎𝑏𝑎) in (17) and call it 𝑠_𝜷̃. For ℎ=1,…,𝑚, Generate 𝐘_𝐵𝑖𝑗∼𝑁_𝑟(0,(1/𝑛_𝑖𝑗)𝐬_𝑖𝑗) and 𝐒_𝐵𝑖𝑗∼𝑊_𝑟(𝑛_𝑖𝑗−1,(1/(𝑛_𝑖𝑗−1))𝐬_𝑖𝑗),  𝑖=1,…,𝑎,𝑗=1,…,𝑏_𝑖, Compute 𝑆_𝜷𝐵(𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎}) using (22). If 𝑆_𝜷𝐵(𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎})>𝑠_𝜷̃, set 𝑄_ℎ=1. (End loop) (1/𝑚)∑^𝑚_ℎ=1𝑄_ℎ is a Monte Carlo estimate of the 𝑃 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: 𝐘^∗_𝑖𝑗𝑘=𝐃𝐘_𝑖𝑗𝑘+𝝃,𝑖=1,…,𝑎;𝑗=1,…,𝑏_𝑖;𝑘=1,…,𝑛_𝑖𝑗, (25) where 𝐃 is any nonsingular matrix and 𝝃 is any given vector. This property is desirable since in practice, the cell responses 𝐘_𝑖𝑗𝑘 (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 𝑠_𝜷̃ (24) and the distribution of PB pivot variable 𝑆_𝜷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) 𝑆_𝜷𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝐒₁₁,…,𝐒_𝑎𝑏𝑎=𝐘^𝑇𝐒^−1/2𝐈_𝑏𝑟−𝐒^−1/2𝐗𝐗^𝑇𝐒⁻¹𝐗⁻𝐗^𝑇𝐒^−1/2𝐒^−1/2𝐘=𝐘^𝑇𝐒^−1/2𝐈_𝑏𝑟−𝐒^−1/2𝐗𝐗^𝑇𝐒⁻¹𝐗+𝐋^𝑇𝐋⁻¹𝐗^𝑇𝐒^−1/2𝐒^−1/2𝐘. (26) 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 𝐋 such as 𝐋=(0,1/𝑎,…,1/𝑎)⊗𝐈_𝑟 in the simulation study.

3. Tests for the Nesting Effects

When the nested effects are present, the nesting effect 𝜶_𝑖 can not reflect the effect of 𝐴_𝑖 because it depends on which level of factor 𝐵 it is in. As pointed in Searle [12], Chap. 7, in the case of 𝑟=1, a popular solution to the problem is not quite a test for 𝜶_𝑖=0 in the presence of nested effects but rather to test the null hypothesis 𝐻_0𝜶∣𝜷∶𝜶_𝑖+𝜷_𝑖𝑗=𝟎,𝑖=1,…,𝑎,𝑗=1,…,𝑏_𝑖. (27) We will see that the problem of testing (27) is the same case of the problem of comparing 𝑏=∑^𝑎_𝑖=1𝑏_𝑖 normal mean vectors when the cell covariance matrices are unknown and arbitrary [5]. In fact, if Σ_𝑖𝑗 are known, then 𝝁₀=^𝑎𝑏_𝑖=1^𝑖_𝑗=1𝑛_𝑖𝑗𝚺⁻¹_𝑖𝑗^𝑎−1𝑏_𝑖=1^𝑖_𝑗=1𝑛_𝑖𝑗𝚺⁻¹_𝑖𝑗𝐘_𝑖𝑗 (28) is the best linear unbiased estimator of 𝝁₀, and a natural statistic for testing (27) is 𝑆_𝜶𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎=^𝑎𝑏_𝑖=1^𝑖_𝑗=1𝐘_𝑖𝑗−𝝁₀^𝑇𝑛_𝑖𝑗𝚺⁻¹_𝑖𝑗𝐘_𝑖𝑗−𝝁₀=𝐘^𝑇𝚺^−1/2𝐈_𝑏𝑟−𝚺^−1/2𝐗₁𝐗^𝑇₁𝚺⁻¹𝐗₁⁻𝐗^𝑇₁𝚺^−1/2𝚺^−1/2𝐘, (29) where 𝐗₁=𝟏_𝑏⊗𝐈_𝑟. Notice that 𝚺^−1/2𝐘∼𝑁_𝑏𝑟(𝚺^−1/2𝝁,𝐈_𝑏𝑟), and 𝚫₁=(𝐈_𝑏𝑟−𝚺^−1/2𝐗₁×(𝐗^𝑇₁𝚺⁻¹𝐗₁)⁻𝐗^𝑇₁𝚺^−1/2) is an idempotent matrix with rank 𝑟(𝑏−1); we have 𝐘^𝑇𝚺^−1/2𝚫₁𝚺^−1/2𝐘∼𝜒²_{𝑟(𝑏−1)}𝝁^𝑇𝚺^−1/2𝚫₁𝚺^−1/2𝝁. (30) The noncentrality parameter 𝝁^′𝚺^−1/2𝚫₁𝚺^−1/2𝝁 is equal to zero when 𝜶_𝑖+𝜷_𝑖𝑗=0, 𝑖=1,…,𝑎, 𝑗=1,…,𝑏_𝑖. Then, the test that rejects 𝐻_0𝜶∣𝜷 in (27) whenever 𝑆_𝜶𝐲₁₁,…,𝐲_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎>𝜒²_{𝑟(𝑏−1),𝜆} (31) is a size 𝜆 test. In general, the covariance matrices 𝚺_𝑖𝑗 are unknown; in this case, a test statistic can be obtained by replacing 𝚺_𝑖𝑗 in (29) by 𝐒_𝑖𝑗, 𝑖=1,…,𝑎, 𝑗=1,…,𝑏_𝑖 and is given by 𝑆_𝜶𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝐒₁₁,…,𝐒_𝑎𝑏𝑎=𝐘^𝑇𝐒^−1/2𝐈_𝑏𝑟−𝐒^−1/2𝐗₁𝐗^𝑇₁𝐒⁻¹𝐗₁⁻𝐗^𝑇₁𝐒^−1/2𝐒^−1/2𝐘. (32)

In the following, we describe the PB test for 𝐻_0𝜶∣𝜷 in (27). Recall that under 𝐻_0𝜶∣𝜷 the vector 𝐘 have the mean 𝐗₁𝝁₀. As the test statistic in (32) is location invariant under the group of location transformations 𝒢₁={𝐘+𝐗₁𝜼,𝜼∈ℝ^𝑟}, without loss of generality, we can take 𝐗₁𝝁₀=0. Using these facts, the parametric bootstrap pivot variable can be developed as follows. For a given (𝐲₁₁,…,𝐲_𝑎𝑏𝑎;𝐬₁₁,…,𝐬_𝑎𝑏𝑎), let 𝐘_𝐵𝑖𝑗∼𝑁_𝑟(0,(1/𝑛_𝑖𝑗)𝐬_𝑖𝑗) and 𝐒_𝐵𝑖𝑗∼𝑊_𝑟(𝑛_𝑖𝑗−1,(1/(𝑛_𝑖𝑗−1))𝐬_𝑖𝑗), 𝑖=1,…,𝑎,𝑗=1,…,𝑏_𝑖. Then the PB pivot variable based on the test statistic (32) is given by 𝑆_𝜶𝐵𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎}=𝐘^𝑇_𝐵𝐒_𝐵^−1/2𝐈_𝑏𝑟−𝐒_𝐵^−1/2𝐗₁𝐗^𝑇₁𝐒_𝐵⁻¹𝐗⁻𝐗^𝑇₁𝐒_𝐵^−1/2𝐒_𝐵^−1/2𝐘_𝐵, (33) where 𝐘_𝐵=(𝐘^𝑇_𝐵11,…,𝐘^𝑇_{𝐵𝑎𝑏𝑎})^𝑇 and 𝐒_𝐵=diag((1/𝑛₁₁)𝐒_𝐵11,…,(1/𝑛_𝑎𝑏𝑎)𝐒_{𝐵𝑎𝑏𝑎}). For a given level 𝜆, the PB test rejects 𝐻_0𝜶∣𝜷 in (27) when 𝑃𝑆_𝜶𝐵𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎}>𝑠_𝜶̃<𝜆, (34) where 𝑠_𝜶̃=𝑆_𝜶𝐲₁₁,…,𝐲_𝑎𝑏𝑎;𝐬₁₁,…,𝐬_𝑎𝑏𝑎 (35) is an observed value of 𝑆_𝜶(𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝐒₁₁,…,𝐒_𝑎𝑏𝑎) in (32). For fixed (𝐲₁₁,…,𝐲_𝑎𝑏𝑎; 𝐬₁₁,…,𝐬_𝑎𝑏𝑎), 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 (𝑛₁₁,…,𝑛_𝑎𝑏𝑎),(𝐲₁₁,…,𝐲_𝑎𝑏𝑎), and (𝐬₁₁,…,𝐬_𝑎𝑏𝑎), Compute 𝑆_𝜶(𝐲₁₁,…,𝐲_𝑎𝑏𝑎;𝐬₁₁,…,𝐬_𝑎𝑏𝑎) in (32) and call it 𝑠_𝜶̃. For ℎ=1,…,𝑚, Generate 𝐘_𝐵𝑖𝑗∼𝑁_𝑟(0,(1/𝑛_𝑖𝑗)𝐬_𝑖𝑗) and 𝐒_𝐵𝑖𝑗∼𝑊_𝑟(𝑛_𝑖𝑗−1,(1/(𝑛_𝑖𝑗−1))𝐬_𝑖𝑗), 𝑖=1,…,𝑎,𝑗=1,…,𝑏_𝑖, Compute 𝑆_𝜶𝐵(𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎}) using (33). If 𝑆_𝜶𝐵(𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎})>𝑠_𝜶̃, set 𝑄_ℎ=1. (End loop) (1/𝑚)∑^𝑚_ℎ=1𝑄_ℎ is a Monte Carlo estimate of the 𝑃 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 Section 3, the problem of testing (27) is the same case of the problem of comparing 𝑏=∑^𝑎_𝑖=1𝑏_𝑖 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 be 𝐴 and 𝐵, respectively. Suppose that the nesting effect of the factor 𝐴 has 𝑎 factor levels, and the nested effect of the factor 𝐵 has a total of 𝑏=∑^𝑎_𝑖=1𝑏_𝑖 levels with 𝑏_𝑖 levels of B nested within the 𝑖th level of the factor 𝐴. Let 𝐧=[𝑛₁₁,𝑛₁₂,…,𝑛_𝑎𝑏𝑎] denote the vector of cell sizes. For given 𝐧 and covariance matrices 𝚺_𝑖𝑗, 𝑖=1,2,…,𝑎, 𝑗=1,2,…,𝑏_𝑖, we first generate ∑^𝑎_𝑖=1∑^𝑏𝑖_𝑗=1𝑛_𝑖𝑗 multivariate samples as 𝐲_𝑖𝑗𝑘=𝝁_𝑖𝑗+𝚺^1/2_𝑖𝑗𝐞_𝑖𝑗𝑘,𝑘=1,2,…,𝑛_𝑖𝑗, (36) where the cell mean vectors 𝝁_𝑖𝑗=𝝁₁₁+𝑖𝑗𝛿𝐡/𝑏 with 𝝁₁₁ being the first cell mean vector, 𝐡 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 𝑟 entries of the error terms 𝐞_𝑖𝑗𝑘 from the 𝑁(0,1) distribution so that we always have E(𝐞_𝑖𝑗𝑘)=𝟎 and Cov(𝐞_𝑖𝑗𝑘)=𝐈_𝑟. This means that (36) will generate the (𝑖𝑗)th multivariate normal sample 𝐲_𝑖𝑗𝑘, 𝑘=1,2,…,𝑛_𝑖𝑗 with the given mean vector 𝝁_𝑖𝑗 and covariance matrix 𝚺_𝑖𝑗. Without loss of generality, we specify 𝝁₁₁ as 𝟎 and 𝐡 as 𝐡₀/‖𝐡₀‖ where 𝐡₀=[1,2,…,𝑟]^𝑇 for any given dimension 𝑟 and ‖𝐡₀‖ denotes the usual 𝐿²-norm of 𝐡₀. To estimate the sizes and powers of the AHT test, we used simulation consisting of 10,000 runs and recorded corresponding 𝑃 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 𝑃 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 𝑃 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 =𝑀⁻¹∑^𝑀_𝑚=1|𝜆_𝑚−𝜆|/𝜆×100 where 𝜆_𝑚 denotes the 𝑚th empirical size for 𝑚=1,2,…,𝑀, 𝜆=0.05, and 𝑀 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 (𝐮_𝑡) to denote “𝐮 repeats 𝑡 times.” Table 1 shows the empirical sizes and powers of the two tests for a bivariate case with 𝑎=2 and 𝑏₁=16, 𝑏₂=20. With 𝑏=𝑏₁+𝑏₂=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 𝑎=3 and 𝑏₁=8, 𝑏₂=10, 𝑏₃=12, and a 10-variate case with 𝑎=3 and 𝑏₁=3, 𝑏₂=5, 𝑏₃=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 [𝑎=2,𝑏₁=16,𝑏₂=20, (𝚺₁₁,𝚺₁₂,…,𝚺_𝑎𝑏𝑎)=𝟏^𝑇₁₈⊗(𝐈₂,diag(𝝀))].

Table 2 

Empirical sizes and powers of the two tests for nested effects for 3-variate two-factor multivariate nested designs [𝑎=3,𝑏₁=8,𝑏₂=10,𝑏₃=12, (𝚺₁₁,𝚺₁₂,…,𝚺_𝑎𝑏𝑎)=𝟏^𝑇₁₀⊗(𝐈₃,diag(𝝀),(1−𝜌)𝐈₃+𝜌𝟏₃𝟏^𝑇₃)].

Table 3 

Empirical sizes and powers of the two tests for nested effects for 10-variate two-factor multivariate nested designs [𝑎=3,𝑏₁=3,𝑏₂=5,𝑏₃=7, (𝚺₁₁,𝚺₁₂,…,𝚺_𝑎𝑏𝑎)=𝟏^𝑇₁₀⊗(𝐈₃,diag(𝝀),diag(𝜼))].

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 Tables 2 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.

Appendix

Proof of Theorem 1. When 𝐻_0𝜷 is true, the model (8) can be written as the following reduced model: 𝐘=𝐗𝜽+𝐞,𝐞∼𝑁_𝑏𝑟(𝟎,𝚺). (A.1) Premultiplying 𝚺^−1/2 in model (A.1), we have the following transformed model: 𝚺^−1/2𝐘=𝚺^−1/2𝐗𝜽+𝚺^−1/2𝐞,𝚺^−1/2𝐞∼𝑁_𝑏𝑟(𝟎,𝐈). (A.2) Under the constraint 𝐋𝜽=𝟎, the model (A.2) becomes 𝚺^−1/2𝐘=𝚺^−1/2𝐗𝜽+𝚺^−1/2𝐞,𝚺^−1/2𝐞∼𝑁_𝑏𝑟(𝟎,𝐈),𝐋𝜽=𝟎. (A.3) In view of Theorem 5.2.5 in Wang and Chow [11], P. 177, we shall prove that 𝐋𝜽=𝟎 is a side condition on the matrix 𝚺^−1/2𝐗; that is, 𝐋 satisfies conditions ℛ𝚺^−1/2𝐗^𝑇∩ℛ𝐋^𝑇={𝟎},𝚺rank^−1/2𝐗𝐋=(𝑎+1)𝑟, (A.4) where ℛ(𝐌) denotes the subspace spanned by the column of matrix 𝐌. Denote 𝐗₀=[𝟏_𝑏,diag(𝟏_𝑏1,…,𝟏_𝑏𝑎)] and 𝐋₀=(0,𝑢₁,⋯,𝑢_𝑎), respectively. Then 𝐗=𝐗₀⊗𝐈_𝑟, and 𝐋=𝐋₀⊗𝐈_𝑟. Note that 𝑢₁,…,𝑢_𝑎 are nonnegative weights such that ∑^𝑎_𝑖=1𝑢_𝑖>0. Then it can be easily verified that 𝐗₀ and 𝐋₀ satisfy the following conditions: ℛ𝐗^𝑇₀∩ℛ𝐋^𝑇₀={𝟎},𝐗rank⁰𝐋₀=𝑎+1. (A.5) Assume that 𝐝∈ℛ(𝐗^𝑇)∩ℛ(𝐋^𝑇). Then there exist two vectors 𝐠 and 𝐡 such that 𝐝=𝐗^𝑇𝐠=𝐋^𝑇𝐡. (A.6) Denote 𝐗^𝑇₀=𝑥₁₁⋯𝑥_1𝑏𝑥⋮⋯⋮_𝑎+1,1⋯𝑥_𝑎+1,𝑏,𝐋^𝑇₀=𝑙₁₁⋮𝑙_𝑎+1,1, (A.7) and 𝐠=(𝐠^𝑇₁,𝐠^𝑇₂,…,𝐠^𝑇_𝑏)^𝑇, where 𝐠_𝑖=(𝑔_𝑖1,…,𝑔_𝑖𝑟)^𝑇 and 𝐡=(ℎ₁,…,ℎ_𝑟)^𝑇 are 𝑟×1 vectors, respectively, 𝑖=1,2,…,𝑏. Then (A.6) may be rewritten as 𝐗^𝑇𝑥𝐠=¹¹𝐠₁+⋯+𝑥_1𝑏𝐠_𝑏⋮𝑥_𝑎+1,1𝐠₁+⋯+𝑥_𝑎+1,𝑏𝐠_𝑏=𝑙₁₁𝐡⋮𝑙_𝑎+1,1𝐡=𝐋^𝑇𝐡, (A.8) which is equivalent to 𝐗^𝑇₀𝐠^∗_𝑘=𝐋^𝑇₀ℎ_𝑘, where 𝐠^∗_𝑘=(𝑔_1𝑘,𝑔_2𝑘,…,𝑔_𝑏𝑘)^𝑇, 𝑘=1,…,𝑟. It follows from (A.5) that 𝐗^𝑇₀𝐠^∗_𝑘=𝐋^𝑇₀ℎ_𝑘,𝑘=1,…,𝑟, (A.9) and hence 𝐝=𝐗^𝑇𝑥𝐠=¹¹𝐠₁+⋯+𝑥_1𝑏𝐠_𝑏⋮𝑥_𝑎+1,1𝐠₁+⋯+𝑥_𝑎+1,𝑏𝐠_𝑏=𝑙₁₁𝐡⋮𝑙_𝑎+1,1𝐡=𝐋^𝑇𝐡=0. (A.10) This means that ℛ𝐗^𝑇∩ℛ𝐋^𝑇={𝟎}. (A.11) It follows from ℛ(𝐗^𝑇)=ℛ(𝐗^𝑇𝚺^−1/2) that ℛ𝚺^−1/2𝐗^𝑇∩ℛ𝐋^𝑇={𝟎}. (A.12) On the other hand, 𝐗𝐋𝐗rank=rank⁰𝐋₀⊗𝐈_𝑟𝐗=rank⁰𝐋₀×rank𝐈_𝑟=(𝑎+1)𝑟. (A.13) The second equation of (A.13) may be derived easily from Theorem 2.2.1 in Wang and Chow [11]. It follows from (A.13) that 𝚺rank^−1/2𝐗𝐋𝚺=rank^−1/2𝟎𝐗𝐋𝐗𝐋=()𝟎𝐈=rank𝑎+1𝑟. (A.14) Combining (A.12) with (A.14), we obtain that 𝐋𝜽=𝟎 is a side condition on the matrix 𝚺^−1/2𝐗. Thus it follows from Lemma 5.2.2 in Wang and Chow [11] that 𝐗^𝑇𝚺⁻¹𝐗𝐗^𝑇𝚺⁻¹𝐗+𝐋^𝑇𝐋⁻¹𝐗^𝑇𝚺^−1/2=𝐗^𝑇𝚺^−1/2. (A.15) Thus, we obtain that 𝐗^𝑇𝚺⁻¹𝐗(𝐗^𝑇𝚺⁻¹𝐗+𝐋^𝑇𝐋)⁻¹𝐗^𝑇𝚺⁻¹𝐗=𝐗^𝑇𝚺⁻¹𝐗. This imply that (𝐗^𝑇𝚺⁻¹𝐗+𝐋^𝑇𝐋)⁻¹ is a generalized inverse of 𝐗^𝑇𝚺⁻¹𝐗. Notice from Corollary 2.2.3 in Wang and Chow [11] that 𝚺^−1/2𝐗(𝐗^𝑇𝚺⁻¹𝐗)⁻𝐗^𝑇𝚺^−1/2 is invariant with respect to the choice of the involved generalized (𝐗^𝑇𝚺⁻¹𝐗)⁻, which leads to 𝚺^−1/2𝐗(𝐗^𝑇𝚺⁻¹𝐗)⁻𝐗^𝑇𝚺^−1/2=𝚺^−1/2𝐗(𝐗^𝑇𝚺⁻¹𝐗+𝐋^𝑇𝐋)⁻¹𝐗^𝑇𝚺^−1/2. Therefore, It follows that 𝐗𝜽=𝐗𝐗^𝑇𝚺⁻¹𝐗+𝐋^𝑇𝐋⁻¹𝐗^𝑇𝚺⁻¹𝐘=𝚺^1/2𝚺^−1/2𝐗𝐗^𝑇𝚺⁻¹𝐗+𝐋^𝑇𝐋⁻¹𝐗^𝑇𝚺^−1/2𝚺^−1/2𝐘=𝚺^1/2𝚺^−1/2𝐗𝐗^𝑇𝚺⁻¹𝐗⁻𝐗^𝑇𝚺^−1/2𝚺^−1/2𝐘=𝐗𝐗^𝑇𝚺⁻¹𝐗⁻𝐗^𝑇𝚺⁻¹𝐘. (A.16) Then we have 𝑆_𝜷𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝚺₁₁,…,𝚺_𝑎𝑏𝑎=𝐘−𝐗𝜽^𝑇𝚺⁻¹𝜽=𝐘−𝐗𝐘^𝑇𝚺^−1/2𝐈_𝑏𝑟−𝚺^−1/2𝐗𝐗^𝑇𝚺⁻¹𝐗⁻𝐗^𝑇𝚺^−1/2𝚺^−1/2𝐘, (A.17) as desired. Theorem 1 is proved.

Proof of Theorem 3. Since the observed values of 𝐘_𝑖𝑗𝑘, 𝑘=1,2,…,𝑛_𝑖𝑗, are denoted as 𝐲_𝑖𝑗𝑘, we let 𝐲^∗_𝑖𝑗𝑘 denote the observed values of the affine-transformed cell responses 𝐘^∗_𝑖𝑗𝑘, 𝑘=1,2,…,𝑛_𝑖𝑗, given by (25). Then we have 𝐲^∗_𝑖𝑗𝑘=𝐃𝐲^∗_𝑖𝑗𝑘+𝝃. On the other hand, we obtain that the sample mean vector and the sample covariance matrix of the (𝑖,𝑗)th affine-transformed cell are 𝐘^∗_𝑖𝑗=𝐃𝐘_𝑖𝑗+𝝃 and 𝐒^∗_𝑖𝑗=𝐃𝐒_𝑖𝑗𝐃^𝑇, respectively. The observed values of these random variables are 𝐲^∗_𝑖𝑗=𝐃𝐲_𝑖𝑗+𝝃 and 𝐬^∗_𝑖𝑗=𝐃𝐬_𝑖𝑗𝐃^𝑇, respectively, 𝑖=1,…,𝑎, 𝑗=1,…,𝑏_𝑖. Let 𝐘^∗=(𝐘^∗𝑇₁₁,…,𝐘^∗𝑇_1𝑏1,…,𝐘^∗𝑇_𝑎1,…,𝐘^∗𝑇_𝑎𝑏𝑎)^𝑇 and 𝐲^∗=(𝐲^∗𝑇₁₁,…,𝐲^∗𝑇_1𝑏1,…,𝐲^∗𝑇_𝑎1,…,𝐲^∗𝑇_𝑎𝑏𝑎)^𝑇; from the reduced model (A.1), we have the following affine-transformed reduced model: 𝐘^∗−𝝃_𝑏=𝐃_𝑏𝐗𝜽+𝐃_𝑏𝐞,𝐃_𝑏𝐞∼𝑁_𝑏𝑟𝟎,𝐃_𝑏𝚺𝐃^𝑇_𝑏, (A.18) where 𝐃_𝑏=𝐈_𝑏⊗𝐃, 𝝃_𝑏=𝟏_𝑏⊗𝝃.
Denote 𝐬^∗=diag((1/𝑛₁₁)𝐬^∗₁₁,…,(1/𝑛_1𝑏1)𝐬^∗_1𝑏1,…,(1/𝑛_𝑎1)𝐬^∗_𝑎1,…,(1/𝑛_𝑎𝑏𝑎)𝐬^∗_𝑎𝑏𝑎). Let 𝑠^∗_𝜷̃ be affine-transformed observed test statistic of 𝑠_𝜷̃ (24). Then, by direct calculations, we have 𝑠^∗_𝜷̃=𝑆_𝜷𝐲^∗₁₁−𝝃,…,𝐲^∗_𝑎𝑏𝑎−𝝃;𝐬^∗₁₁,…,𝐬^∗_𝑎𝑏𝑎=𝐲^∗−𝝃_𝑏^𝑇𝐬^∗−1/2𝐈_𝑏𝑟−𝐬^∗−1/2𝐃_𝑏𝐗𝐗^𝑇𝐃^𝑇_𝑏𝐬^∗−1𝐃_𝑏𝐗⁻×𝐗^𝑇𝐃^𝑇_𝑏𝐬^∗−1/2𝐬^∗−1/2𝐲^∗−𝝃_𝑏=𝐲^𝑇𝐬^−1/2𝐈_𝑏𝑟−𝐬^−1/2𝐗𝐗^𝑇𝐬⁻¹𝐗⁻𝐗^𝑇𝐬^−1/2𝐬^−1/2𝐲=𝑠_𝜷̃, (A.19) where 𝐬=diag((1/𝑛₁₁)𝐬₁₁,…,(1/𝑛_1𝑏1)𝐬_1𝑏1,…,(1/𝑛_𝑎1)𝐬_𝑎1,…,(1/𝑛_𝑎𝑏𝑎)𝐬_𝑎𝑏𝑎) is the observed values of 𝐒. The affine-invariance of the observed test statistic 𝑠_𝜷̃ (24) then follows.
Let 𝐘^∗_𝐵𝑖𝑗∼𝑁_𝑟10,𝑛_𝑖𝑗𝐬^∗_𝑖𝑗,𝐒^∗_𝐵𝑖𝑗∼𝑊_𝑟𝑛_𝑖𝑗1−1,𝑛_𝑖𝑗𝐬−1^∗𝑖𝑗 (A.20) be mutually independent, 𝑖=1,…,𝑎, 𝑗=1,…,𝑏_𝑖. Then 𝐘^∗𝑑_𝐵𝑖𝑗=𝐃𝐘_𝐵𝑖𝑗 and 𝐒^∗𝑑_𝐵𝑖𝑗=𝐃𝐒_𝐵𝑖𝑗𝐃^𝑇, where 𝐇=^𝑑𝐆 means that 𝐇 and 𝐆 have the same distribution. Using these facts, the affine-transformed parametric bootstrap pivot variable can be given by 𝑆^∗_𝜷𝐵𝐘^∗_𝐵11,…,𝐘^∗_{𝐵𝑎𝑏𝑎};𝐒^∗_𝐵11,…,𝐒^∗_{𝐵𝑎𝑏𝑎}=𝐘_𝐵^∗𝑇𝐒_𝐵^∗−1/2𝐈_𝑏𝑟−𝐒_𝐵^∗−1/2𝐃_𝑏𝐗𝐗^𝑇𝐃^𝑇_𝑏𝐒_𝐵^∗−1𝐃_𝑏𝐗⁻×𝐗^𝑇𝐃^𝑇_𝑏𝐒_𝐵^∗−1/2𝐒_𝐵^∗−1/2𝐘^∗_𝐵𝑑=𝐘^𝑇_𝐵𝐃^𝑇_𝑏𝐃_𝑏𝐒_𝐵𝐃^𝑇_𝑏^−1/2𝐈_𝑏𝑟−𝐃_𝑏𝐒_𝐵𝐃^𝑇_𝑏^−1/2𝐃_𝑏𝐗×𝐗^𝑇𝐃^𝑇_𝑏𝐃_𝑏𝐒_𝐵𝐃^𝑇_𝑏⁻¹𝐃_𝑏𝐗⁻×𝐗^𝑇𝐃^𝑇_𝑏𝐃_𝑏𝐒_𝐵𝐃^𝑇_𝑏^−1/2×𝐃_𝑏𝐒_𝐵𝐃^𝑇_𝑏^−1/2𝐃_𝑏𝐘_𝐵=𝐘^𝑇_𝐵𝐒_𝐵^−1/2𝐈_𝑏𝑟−𝐒_𝐵^−1/2𝐗𝐗^𝑇𝐒_𝐵⁻¹𝐗⁻𝐗^𝑇𝐒_𝐵^−1/2𝐒_𝐵^−1/2𝐘_𝐵=𝑆_𝜷𝐵𝐘_𝐵11,…,𝐘_{𝐵𝑎𝑏𝑎};𝐒_𝐵11,…,𝐒_{𝐵𝑎𝑏𝑎}, (A.21) where 𝐘^∗_𝐵=(𝐘^∗𝑇_𝐵11,…,𝐘^∗𝑇_𝐵1𝑏1…,𝐘^∗𝑇_𝐵𝑎1,…,𝐘^∗𝑇_{𝐵𝑎𝑏𝑎})^𝑇 and 𝐒^∗_𝐵=diag((1/𝑛₁₁)𝐒^∗_𝐵11,…,(1/𝑛_1𝑏1)𝐒^∗_𝐵1𝑏1,…,(1/𝑛_𝑎1)𝐒^∗_𝐵𝑎1,…,(1/𝑛_𝑎𝑏𝑎)𝐒^∗_{𝐵𝑎𝑏𝑎}). The affine-invariance of the distribution of PB pivot variable 𝑆_𝜷𝐵 (22) follows. Theorem 3 is then proved.

Proof of Theorem 4. The proof of Theorem 4 is obvious and hence omitted here.

Proof of Theorem 5. Using the same proof as that of Theorem 1, it follows that 𝐋𝜽=𝟎 is a side condition on the matrix 𝐒^−1/2𝐗; that is, 𝐋 satisfies conditions ℛ𝐒^−1/2𝐗^𝑇∩ℛ𝐋^𝑇={𝟎}𝐒,a.e.,rank^−1/2𝐗𝐋=(𝑎+1)𝑟,a.e., (A.22) where a.e. denotes “almost everywhere.” Thus it follows from the similar proof of Theorem 1 that 𝐗𝐗^𝑇𝐒⁻¹𝐗+𝐋^𝑇𝐋⁻¹𝐗^𝑇𝐒^−1/2=𝐗𝐗^𝑇𝐒⁻¹𝐗⁻𝐗^𝑇𝐒^−1/2,a.e. (A.23) Then we obtain from (A.23) that the test statistic (17) 𝑆_𝜷𝐘₁₁,…,𝐘_𝑎𝑏𝑎;𝐒₁₁,…,𝐒_𝑎𝑏𝑎=𝐘^𝑇𝐒^−1/2𝐈_𝑏𝑟−𝐒^−1/2𝐗𝐗^𝑇𝐒⁻¹𝐗⁻𝐗^𝑇𝐒^−1/2𝐒^−1/2𝐘=𝐘^𝑇𝐒^−1/2𝐈_𝑏𝑟−𝐒^−1/2𝐗𝐗^𝑇𝐒⁻¹𝐗+𝐋^𝑇𝐋⁻¹𝐗^𝑇𝐒^−1/2𝐒^−1/2𝐘 (A.24) as desired. Theorem 5 is proved.

Conflict of Interests

The author declares that there is not any conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (11171002) and the Importation and Development of High-Caliber Talents Project of Beijing Municipal Institutions (CIT&TCD201404002).