Parametric test based on the bootstrapping approach for the MANOVA under a Behrens-Fisher problem

Abstract

This article presents a comparison of multivariate normal mean vectors under covariance positive definite matrices. We introduce an improved parametric bootstrap (IPB) approach for addressing the multivariate Behrens-Fisher problem, specifically focusing on cases with unequal covariance matrices. Additionally, we evaluate the performance of the IPB test by comparing it with three existing tests: the parametric bootstrap (PB) test, the generalized variable (GV) test, and the Johansen test. Through Monte Carlo simulation, our results demonstrate that both the IPB test and the PB test exhibit superior control over Type I error rates compared to the GV and Johansen tests. Notably, the IPB test outperforms the PB test in terms of controlling Type I error rates. Consequently, our study concludes that the IPB test represents a preferred statistical method for testing the equality of mean vectors in the multivariate Behrens-Fisher problem.

1. Introduction

Multivariate analysis of variance (MANOVA) is used for comparing the mean vectors with multivariate normal populations. The observed vectors are $Y_{ij}$ where $i=1,2,\mathrm{\dots },k$ denotes the sample and $j=1,2,\mathrm{\dots },n_i$ denotes the observation within the sample. Each observation vector $Y_{ij}$ is a $p$ -multiple multivariate normal vectors with a mean vector ${\mu }_i$ and equal common covariance matrix $\sum$ . The model for MANOVA can be written as

$Y_{ij}={\mu }_i+{\epsilon }_{ij},(i=1,2,\mathrm{\dots },k;j=1,2,\mathrm{\dots },n_i),$(1)

where ${\mu }_i$ is the parameter mean vector and ${\epsilon }_{ij}$ represents error terms that are independently drawn from a multivariate normal distribution with mean vector 0 and equal common covariance matrix $\sum$ . Therefore, $E(Y_{ij})={\mu }_i$ and $\text{𝑉𝑎𝑟}(Y_{ij})=\sum$ . The hypotheses for testing the equality of the mean vectors in MANOVA are

$H_0:{\mu }_1={\mu }_2=\mathrm{\dots }={\mu }_k\text{ versus }{H}_1:{\mu }_i\ne {\mu }_j\text{ for some }i\ne j.$(2)

In general, a test statistic for testing the equality of the mean vectors in Eq. (2) is the Hotelling $T^2$ test (Zhang & Xu, 2009), which is a test statistic based on assumptions. Most tests available in the literature for testing the equality of the mean vectors for the hypotheses in Eq. (2) based on equal common covariance matrices $\sum$ have been presented by many researchers, such as Lawley (1938), Hotelling (1951), Wilks’ (1932) and Bartlett (1939).

Sometime, the assumption of the equal common covariance matrix in MANOVA is violated, which is the Behrens-Fisher problem (unequal common covariance matric ${\sum }_i$ ). With an unequal common covariance matric ${\sum }_i$ , the aforementioned statistics show that it cannot control the Type I error rate, which is reflected in the performance of the test statistics. Therefore, many researchers have presented test statistics for testing the equality of the mean vectors for the hypotheses in Eq. (1) based on an unequal common covariance matrix.

A natural test statistic for testing the equality of mean vectors in Eq. (2) is given by

$T({\overline{Y}}_1,\mathrm{\dots },{\overline{Y}}_k;{\stackrel{\sim }{S}}_1,\mathrm{\dots },{\stackrel{\sim }{S}}_k)={\displaystyle \sum _{i=1}^k}{({\overline{Y}}_i-{\widehat{\mu }}_0^{*})}^{\prime }{W}_i({\overline{Y}}_i-{\widehat{\mu }}_0^{*})={\displaystyle \sum _{i=1}^k}{\overline{Y}}_i^{\prime }W{\overline{Y}}_i-{\left({\displaystyle \sum _{i=1}^k}{W}_i{\overline{Y}}_i\right)}^{\prime }{W}^{-1}\left({\displaystyle \sum _{i=1}^k}{W}_i{\overline{Y}}_i\right)$(3)

Where ${\overline{Y}}_i=\sum _{j=1}^{n_i}{Y}_{ij}$ , $W=\sum _{i=1}^k{W}_i$ and $W_i={\stackrel{\sim }{S}}^{-1}$ with ${\stackrel{\sim }{S}}_i=\frac{1}{n_i}{S}_i$ and $S_i=\frac{1}{n_i-1}\sum _{j=1}^{n_i}(Y_{ij}-{\overline{Y}}_i){(Y_{ij}-{\overline{Y}}_i)}^{\prime }$ . The null hypothesis $H_0$ in 2 is rejected when $T({\overline{Y}}_1,\mathrm{\dots },{\overline{Y}}_k;{\stackrel{\sim }{S}}_1,\mathrm{\dots },{\stackrel{\sim }{S}}_k)>{\chi }_{p(k-1)}^2$ .

Johansen (1980) developed the test statistic in Eq. (3) for testing the equality of the mean vectors for the hypotheses in Eq. (2), based on an unequal common covariance matrix which is

$J_{OH}={\displaystyle \frac{T({\overline{Y}}_1,\mathrm{\dots },{\overline{Y}}_k;{\stackrel{\sim }{S}}_1,\mathrm{\dots },{\stackrel{\sim }{S}}_k)}{c}}$(4)

where $c$ is

$c=p(k-1)+2A-{\displaystyle \frac{6A}{p(k-1)+2}}$(5)

and

$A={\displaystyle \frac{\sum _{i=1}^{k}tr{(I-W^{-1}{W}_i)}^2+{[tr(I-W^{-1}{W}_i)]}^2}{2(n_i-1)}}$(6)

The Johansen test is distributed as an $F$ -distribution with degrees of freedom $f_1=p(k-1)$ and $f_2=p(k-1)[p(k-1)+2]/(3A)$ . Therefore, Johansen’s test rejects the null hypothesis $H_0$ in Eq. (2) when $J_{OH}>F_{f_1,f_2,1-\alpha }$ , where $F_{m,r;q}$ denotes the $q$ th quantile of an $F$ distribution with degrees of freedom $m$ and $r$ .

Gamage et al. (2004) used the concept of a generalized $p$ -value supported by Tsui and Weerahandi and Tsui (1989) for developing the generalized variable (GV) test. The GV test is described as follows. Let $({\stackrel{\sim }{y}}_i,{\stackrel{\sim }{s}}_i)$ be an observed value of $({\overline{Y}}_i,{\stackrel{\sim }{S}}_i)$ , and let

$R_i^{*}={[{\stackrel{\sim }{s}}_i^{-1/2}{\stackrel{\sim }{\mathrm{\Sigma }}}_i{\stackrel{\sim }{s}}_i^{-1/2}]}^{-1/2}[{\stackrel{\sim }{s}}_i^{-1/2}{\stackrel{\sim }{S}}_i{\stackrel{\sim }{s}}_i^{-1/2}]{[{\stackrel{\sim }{s}}_i^{-1/2}{\stackrel{\sim }{\mathrm{\Sigma }}}_i{\stackrel{\sim }{s}}_i^{-1/2}]}^{-1/2}\hspace{1em},i=1,\mathrm{\dots },k$(7)

where ${\stackrel{\sim }{s}}_i$ and $R_i^{*}$ are independent, and ${\stackrel{\sim }{s}}_i^{-1/2}{\stackrel{\sim }{S}}_i{\stackrel{\sim }{s}}_i^{-1/2}\sim W_p(n_i-1,{\stackrel{\sim }{s}}_i^{-1/2}{\stackrel{\sim }{\mathrm{\Sigma }}}_i{\stackrel{\sim }{s}}_i^{-1/2})$ , $R_i^{*}\sim W_p(n_i-1,1/(n_i-1)I_p),i=1,\mathrm{\dots },k$ .

The GV test is presented as

$G={\displaystyle \frac{T({\overline{Y}}_1,\mathrm{\dots },{\overline{Y}}_k;{\stackrel{\sim }{S}}_1,\mathrm{\dots },{\stackrel{\sim }{S}}_k)}{T({\stackrel{\sim }{y}}_1,\mathrm{\dots },{\stackrel{\sim }{y}}_k;{\stackrel{\sim }{s}}_1^{-1/2}{R}_1^{*-1}{\stackrel{\sim }{s}}_1^{-1/2},\mathrm{\dots },{\stackrel{\sim }{s}}_k^{-1/2}{R}_k^{*-1}{\stackrel{\sim }{s}}_k^{-1/2})}}$(8)

since $T({\overline{Y}}_1,\mathrm{\dots },{\overline{Y}}_k;{\stackrel{\sim }{S}}_1,\mathrm{\dots },{\stackrel{\sim }{S}}_k)$ is distributed as chi-squared with degree of freedom $p(k-1)$ . The generalized $p$ -value is given by

$P_{{\chi }_{p(k-1)}^2,R_1^{*},\mathrm{\dots },R_k^{*}}({\displaystyle \frac{T({\overline{Y}}_1,\mathrm{\dots },{\overline{Y}}_k;{\stackrel{\sim }{S}}_1,\mathrm{\dots },{\stackrel{\sim }{S}}_k)}{T({\stackrel{\sim }{y}}_1,\mathrm{\dots },{\stackrel{\sim }{y}}_k;{\stackrel{\sim }{s}}_1^{-1/2}{R}_1^{*-1}{\stackrel{\sim }{s}}_1^{-1/2},\mathrm{\dots },{\stackrel{\sim }{s}}_k^{-1/2}{R}_k^{*-1}{\stackrel{\sim }{s}}_k^{-1/2})}}>1)$(9)

where the GV test rejects the null hypothesis in Eq. (2) when the GV test in Eq. (9) is less than a given nominal level $\alpha$ .

Krishnamoorthy and Yu (2010) proposed the parametric bootstrap (PB) test by improving the sample mean ${\overline{Y}}_i$ and sample variance ${\stackrel{\sim }{S}}_i$ by using the bootstrap approach in which the PB pivotal quantity can be obtained as follows. Let ${\overline{Y}}_{Bi}\sim N_p(0,{\stackrel{\sim }{s}}_i)$ and ${\stackrel{\sim }{S}}_i\sim W_p(n_i-1,(1/(n_i-1)){\stackrel{\sim }{s}}_i)$ . The PB pivotal quantity can be written as

$T({\overline{Y}}_{B1},\mathrm{\dots },{\overline{Y}}_{Bk};{\stackrel{\sim }{S}}_{B1},\mathrm{\dots },{\stackrel{\sim }{S}}_{Bk})={\displaystyle \sum _{i=1}^k}{({\overline{Y}}_{Bi}-{\widehat{\mu }}_B^{*})}^{\prime }{\stackrel{\sim }{S}}_{Bi}^{-1}({\overline{Y}}_{Bi}-{\widehat{\mu }}_B^{*})$(10)

where ${\widehat{\mu }}_B^{*}={\left(\sum _{i=1}^k{\stackrel{\sim }{S}}_{Bi}^{-1}\right)}^{-1}\sum _{i=1}^k{\stackrel{\sim }{S}}_{Bi}^{-1}{\overline{Y}}_{Bi}$ or equivalently,

$T({\overline{Y}}_{B1},\mathrm{\dots },{\overline{Y}}_{Bk};{\stackrel{\sim }{S}}_{B1},\mathrm{\dots },{\stackrel{\sim }{S}}_{Bk})={\displaystyle \sum _{i=1}^k}{\overline{Y}}_{Bi}^{\prime }{\stackrel{\sim }{S}}_{Bi}^{-1}{\overline{Y}}_{Bi}-\left({\displaystyle \sum _{i=1}^k}{\overline{Y}}_{Bi}^{\prime }{\stackrel{\sim }{S}}^{-1}\right){\left({\displaystyle \sum _{i=1}^k}{\stackrel{\sim }{S}}_{Bi}^{-1}\right)}^{-1}\left({\displaystyle \sum _{i=1}^k}{\stackrel{\sim }{S}}_{Bi}^{-1}{\overline{Y}}_{Bi}\right)$(11)

The PB $p$ -value is defined as

$P(T({\overline{Y}}_{B1},\mathrm{\dots },{\overline{Y}}_{Bk};{\stackrel{\sim }{S}}_{B1},\mathrm{\dots },{\stackrel{\sim }{S}}_{Bk})>T_0),$(12)

where $T_0$ is the observed value in Eq. (3).

The PB pivotal quantity can be estimated using the Monte Carlo simulation as described below. Let $t_i$ be the Cholesky factor of ${\stackrel{\sim }{s}}_i$ , so that ${\stackrel{\sim }{s}}_i=t_i{t}_i^{\prime }$ , $i=1,\mathrm{\dots },k$ . Then ${\overline{Y}}_{Bi}\sim t_i{Z}_i$ and ${\stackrel{\sim }{S}}_{Bi}^{-1}\sim t_i{V}_i{t}_i^{\prime }/(n_i-1)$ , where $Z_i$ and $V_i$ are independent with $Z_i\sim N_p(0,I_p)$ and $V_i\sim W_p(n_i-1,I_p)$ . The PB test in Eq. (11) is distributed as

$T_B(Z_i,V_i)={\displaystyle \sum _{i=1}^k}{f}_{i}Z^{\prime }{}_i{V}_i^{-1}{Z}_i-\left({\displaystyle \sum _{i=1}^k}{f}_i{Z}_i^{\prime }{V}_i^{-1}{t}_i^{-1}\right){\left({\displaystyle \sum _{i=1}^k}{f}_i{(t_i{V}_i{t}_i^{\prime })}^{-1}\right)}^{-1}\left({\displaystyle \sum _{i=1}^k}{f}_i{t}_i^{\prime -1}{V}_i^{-1}{Z}_i\right),$(13)

where $f_i=n_i-1$ . Hence the $p$ -value of the PB test is

$P(T_B(Z_i,V_i)>T_0).$(14)

The PB test rejects the null hypothesis in Eq. (2) when the $p$ -value of the PB test in Eq. (14) is less than the nominal level $\alpha$ . According to the results presented by Krishnamoorthy and Yu (2010) the PB test could control type I error rates better than the Johansen test, and the GV test performed very satisfactorily when the sample sizes are small.

The Bootstrap approach is a technical statistic for reducing the error of the test statistics for hypothesis testing, where it uses the process of the resampling technique. The Bootstrap approach is employed to reduce the error for estimating the test statistics supported by Jiang and Simon. Hence, this research uses the Bootstrap approach for developing the PB test for testing the mean vector ${\mu }_i$ with multivariate normal populations under unequal variance, where the Bootstrap approach is used to improve the $p$ -value of the PB test, and it is called the Improved Parametric Bootstrap (IPB) test.

The objective of this research is to compare the performance of the Johansen test, PB test, GV test and IPB test for testing the equality of the vector means in MANOVA when the unequal common covariance matrix is based on Type I error rates, which reflect the performance of the test statistic. The paper is organized as follows. The IPB test in MANOVA is described in Section 2. Section 3 show the results of simulation of the performance of the tests based on Type I error rates.and Section 4 presents illustrative examples. Finally, Section 5 contains conclusions.

2. Methodology

We proposed IPB test for comparing mean vectors under a Behrens–Fisher problem as follows:

Step 1: Compute the test statistic in Eq. (3) as

$T({\overline{Y}}_1,\mathrm{\dots },{\overline{Y}}_k;{\stackrel{\sim }{S}}_1,\mathrm{\dots },{\stackrel{\sim }{S}}_k)={\displaystyle \sum _{i=1}^k}{\overline{Y}}_i^{\prime }W{\overline{Y}}_i-{\left({\displaystyle \sum _{i=1}^k}{W}_i{\overline{Y}}_i\right)}^{\prime }{W}^{-1}\left({\displaystyle \sum _{i=1}^k}{W}_i{\overline{Y}}_i\right).$(15)

Step 2: Calculate PB tests in Eq. (13) with m values.

Step 3: From step 2, obtain the PB test of m values which is

$X=(P{B}_1,\mathrm{\dots },P{B}_{i-1},P{B}_i,P{B}_{i+1},\mathrm{\dots },P{B}_m).$(16)

Step 4: From step 3, the sample of the Bootstrap approach is resampled in Eq. (16). The sample of the Bootstrap approach is

$X^{\ast }=(P{B}_1^{\ast },\mathrm{\dots },P{B}_{i-1}^{\ast },P{B}_i^{\ast },P{B}_{i+1}^{\ast },\mathrm{\dots },P{B}_m^{\ast }).$(17)

Step 5: Compare the values in $X^{\ast }$ with the test statistic in step 1 when the values in $X^{\ast }$ are more than the test statistic in step 1, set $Q_j=1$ .

Step 6: From step 5, we obtain the $p$ -value as

${\mathit{\text{p-value}}}_i={\displaystyle \frac{\sum _{j=1}^m{Q}_j}{m}}.$(18)

Step 7: Repeat the steps 4, 5 and 6 n times. We obtain the $p$ -value of IPB is

${\mathit{\text{p-value}}}_{\text{𝐼𝑃𝐵}}={\displaystyle \frac{\sum _{i=1}^n{\mathit{\text{p-value}}}_i}{n}}$(19)

where IPB test rejects the null hypothesis in Eq. (2) when the $p$ -value in Eq. (19) is less than the nominal level $\alpha$ .

3. Results

For simulation studies, we used Monte Carlo simulation with the R statistical package. Monte Carlo simulation is presented for calculating Type I error rates of the Johansen test, GV test, PB test and IPB test. We denote, without loss of generality, that $\mu =({\mu }_1,{\mu }_2,\mathrm{\dots },{\mu }_k)$ to be a vector of zeroes to examine Type I error rates of the Johansen test, GV test, PB test and IPB test. For comparing the population mean vectors, we can assume ${\sum }_1=I_p$ , ${\sum }_2=\text{diag}({\lambda }_1,\mathrm{\dots },{\lambda }_p)$ and other matrices are arbitrarily positive definite. For $k=$ 3 and 5 under various values of $n_i$ in our simulation studies, we take the covariance matrices to range from $-1<\rho <1$ .

For Type I error rates of the Johansen test, GV test, PB test and IPB test, the sample mean ${\overline{y}}_i$ and the sample covariance matrix ${\stackrel{\sim }{s}}_i$ based on the $i$ th sample are generated independently as ${\overline{y}}_i\sim N_p(0,\frac{1}{n_i}{\mathrm{\Sigma }}_i)$ and ${\stackrel{\sim }{s}}_i\sim W_p(n_i-1,\frac{1}{n_i-1}{\stackrel{\sim }{\mathrm{\Sigma }}}_i)$ with ${\stackrel{\sim }{\mathrm{\Sigma }}}_i=\frac{1}{n_i}{\mathrm{\Sigma }}_i$ . We used 10,000 observed vectors $({\overline{y}}_1,\mathrm{\dots },{\overline{y}}_k,{\stackrel{\sim }{s}}_1,\mathrm{\dots },{\stackrel{\sim }{s}}_k)$ to compute the observed value in Eq. (3). To calculate the Type I error rates of the Johansen test, they are determined by the proportions of Johansen test that exceed the critical value. For estimating Type I error rates of the GV test, PB test and IPB test, we used 10,000 observed vectors $({\overline{y}}_1,\mathrm{\dots },{\overline{y}}_k,{\stackrel{\sim }{s}}_1,\mathrm{\dots },{\stackrel{\sim }{s}}_k)$ to compute the observed value in Eq. (3). We use 10,000 runs to estimate the $p$ -value of the GV test, PB test and IPB test. Finally, the Type I error rates are estimated by the proportion of the 10,000 $p$ -values that are less than the nominal level 0.05. The values of the Type I error rates are shown in Tables 1–3.

Table 1

Monte Carlo estimates of Type I error rates for comparing bivariate normal mean vectors

$k=$ 3, $p=$ 2, ${\sum }_1=I_2$ , ${\sum }_2=\text{diag}({\lambda }_1,{\lambda }_2)$ , ${\sum }_3=\left(\begin{array}{cccccccccccccccccccc}\hfill 1\hfill & \hfill {\rho }_3\hfill & & & & & & & & & & & & & & & & & & \\ \hfill {\rho }_3\hfill & \hfill 1\hfill & & & & & & & & & & & & & & & & & & \end{array}\right)$
$(n_1,n_2,n_3)$	$({\lambda }_1,{\lambda }_2,{\rho }_3)$	PB	IPB	Johansen	GV
(7, 7, 7)	(1, 1, 0)	0.044	0.046	0.057	0.054
	(1, 0.9, 0.1)	0.048	0.047	0.054	0.052
	(1, 0.5, 0.2)	0.043	0.049	0.058	0.047
	(1, 0.1, 0.3)	0.047	0.045	0.068	0.061
	(0.2, 0.6, 0.5)	0.056	0.053	0.067	0.073
	(0.9, 0.9, 0.6)	0.049	0.049	0.055	0.056
	(0.7, 0.8, $-$ 0.2)	0.046	0.045	0.057	0.058
(7, 10, 20)	(1, 1, 0)	0.052	0.051	0.060	0.095
	(1, 0.9, 0.1)	0.054	0.051	0.059	0.072
	(1, 0.5, 0.2)	0.052	0.051	0.062	0.089
	(1, 0.1, 0.3)	0.054	0.049	0.070	0.086
	(0.2, 0.6, 0.5)	0.053	0.051	0.067	0.078
	(0.9, 0.9, 0.6)	0.054	0.052	0.061	0.096
	(0.7, 0.8, $-$ 0.2)	0.053	0.051	0.064	0.079
(10, 10, 40)	(1, 1, 0)	0.058	0.051	0.055	0.100
	(1, 0.9, 0.1)	0.049	0.048	0.055	0.090
	(1, 0.5, 0.2)	0.043	0.042	0.054	0.093
	(1, 0.1, 0.3)	0.052	0.049	0.055	0.096
	(0.2, 0.6, 0.5)	0.052	0.049	0.054	0.110
	(0.9, 0.9, 0.6)	0.057	0.051	0.055	0.111
	(0.7, 0.8, $-$ 0.2)	0.045	0.043	0.055	0.099
(25, 20, 20)	(1, 1, 0)	0.049	0.051	0.049	0.043
	(1, 0.9, 0.1)	0.050	0.048	0.050	0.059
	(1, 0.5, 0.2)	0.049	0.046	0.051	0.048
	(1, 0.1, 0.3)	0.052	0.049	0.049	0.054
	(0.2, 0.6, 0.5)	0.053	0.049	0.052	0.054
	(0.9, 0.9, 0.6)	0.050	0.048	0.053	0.059
	(0.7, 0.8, $-$ 0.2)	0.050	0.048	0.050	0.058
$k=$ 5, $p=$ 2, ${\sum }_1=I_2$ , ${\sum }_2=\text{diag}({\lambda }_1,{\lambda }_2)$ , ${\sum }_i=\left(\begin{array}{cccccccccccccccccccc}\hfill 1\hfill & \hfill {\rho }_i\hfill & & & & & & & & & & & & & & & & & & \\ \hfill {\rho }_i\hfill & \hfill 1\hfill & & & & & & & & & & & & & & & & & & \end{array}\right)$ , $i=3,4,5$
$(n_1,\mathrm{\dots },n_5)$	$({\lambda }_1,{\lambda }_2,{\rho }_3,{\rho }_4,{\rho }_5)$	PB	IPB	Johansen	GV
(7, 7, 7, 7, 7)	(1, 1, 1, 0.5, 0.5, 0.5)	0.050	0.049	0.071	0.104
	(0.1, 0.1, 0.1, 0.3, 0.3, 0.3)	0.050	0.048	0.072	0.114
	(0.1, 0.4, 0.7, 0, 0, 0)	0.048	0.048	0.072	0.113
	(0.1, 0.3, 0.9, 0.1, 0.4, 0.9)	0.048	0.048	0.074	0.123
	(0.1, 0.2, 0.3, $-$ 0.1, 0.1, 0.9)	0.051	0.050	0.076	0.124
	(0.4, 0.4, 0.5, $-$ 0.3, 0.4, 0.3)	0.053	0.051	0.072	0.118
	(0.9, 0.9, 0.9, $-$ 0.4, 0.6, 0.9)	0.052	0.050	0.077	0.133
(12, 12, 12, 12, 12)	(1, 1, 1, 0.5, 0.5, 0.5)	0.050	0.045	0.055	0.075
	(0.1, 0.1, 0.1, 0.3, 0.3, 0.3)	0.053	0.046	0.056	0.078
	(0.1, 0.4, 0.7, 0, 0, 0)	0.052	0.045	0.056	0.085
	(0.1, 0.3, 0.9, 0.1, 0.4, 0.9)	0.051	0.046	0.056	0.083
	(0.1, 0.2, 0.3, $-$ 0.1, 0.1, 0.9)	0.050	0.047	0.057	0.086
	(0.4, 0.4, 0.5, $-$ 0.3, 0.4, 0.3)	0.050	0.046	0.056	0.082
	(0.9, 0.9, 0.9, $-$ 0.4, 0.6, 0.9)	0.048	0.048	0.057	0.084
(20, 20, 20, 20, 20)	(1, 1, 1, 0.5, 0.5, 0.5)	0.054	0.049	0.053	0.054
	(0.1, 0.1, 0.1, 0.3, 0.3, 0.3)	0.051	0.045	0.051	0.057
	(0.1, 0.4, 0.7, 0, 0, 0)	0.052	0.046	0.052	0.065
	(0.1, 0.3, 0.9, 0.1, 0.4, 0.9)	0.047	0.047	0.052	0.061
	(0.1, 0.2, 0.3, $-$ 0.1, 0.1, 0.9)	0.051	0.049	0.052	0.067
	(0.4, 0.4, 0.5, $-$ 0.3, 0.4, 0.3)	0.053	0.050	0.052	0.054
	(0.9, 0.9, 0.9, $-$ 0.4, 0.6, 0.9)	0.048	0.047	0.053	0.065

In Table 1, for $k=$ 3 and $p=$ 2, in the case of the balanced sample sizes, the four tests control Type I error rates quite well. Meanwhile, in the case of the unbalanced sample sizes $n=(7,10,20)$ , the PB test and IBP test could control Type I error rates more satisfactorily than the Johansen test and GV test. Again, observe that the GV test is extremely poor in the case of unbalanced sample sizes $n=(10,10,40)$ , but the PB test, IPB test and Johansen test are quite good in this case. For $k=$ 5 and $p=$ 2, in the case of the balanced sample sizes, the GV test is poor at controlling Type I error rates, however the PB test and IBP test are very good in the case of balanced sample sizes. In unbalanced sample sizes $n=(15,20,10,32,7)$ , the PB test and IPB test could control Type I error rates more than the Johansen test and GV test. Considering the PB test and IPB test it was found that the Type I error rates of the IPB test do not exceed the nominal level 0.05 better than for the PB test in case $k=$ 5 and $p=$ 2.

Table 2

Monte Carlo estimates of Type I error rates for comparing trivariate normal mean vectors

$k=$ 3, $p=$ 3, ${\sum }_1=I_3$ , ${\sum }_2=\text{diag}({\lambda }_1,{\lambda }_2,{\lambda }_3)$ , ${\sum }_3=\left(\begin{array}{cccccccccccccccccccc}\hfill 1\hfill & \hfill \rho \hfill & \hfill \rho \hfill & & & & & & & & & & & & & & & & & \\ \hfill \rho \hfill & \hfill 1\hfill & \hfill \rho \hfill & & & & & & & & & & & & & & & & & \\ \hfill \rho \hfill & \hfill \rho \hfill & \hfill 1\hfill & & & & & & & & & & & & & & & & & \end{array}\right)$
$(n_1,n_2,n_3)$	$({\lambda }_1,{\lambda }_2,{\lambda }_3,\rho )$	PB	IPB	Johansen	GV
(7, 7, 7)	(1, 1, 1, 0)	0.039	0.041	0.079	0.077
	(1, 1, 0.1, 0.1)	0.043	0.042	0.090	0.081
	(1, 0.1, 0.1, 0.5)	0.049	0.047	0.105	0.096
	(0.2, 0.6, 0.9, $-$ 0.3)	0.048	0.048	0.089	0.086
	(0.6, 0.6, 0.6, 0)	0.045	0.046	0.086	0.076
	(0.3, 0.9, 0.1, $-$ 0.1)	0.045	0.048	0.102	0.081
	(0.8, 0.5, 0.5, 0.1)	0.047	0.046	0.087	0.079
(7, 10, 20)	(1, 1, 1, 0)	0.054	0.050	0.080	0.101
	(1, 1, 0.1, 0.1)	0.058	0.051	0.091	0.114
	(1, 0.1, 0.1, 0.5)	0.064	0.058	0.111	0.117
	(0.2, 0.6, 0.9, $-$ 0.3)	0.061	0.053	0.094	0.108
	(0.6, 0.6, 0.6, 0)	0.052	0.045	0.086	0.109
	(0.3, 0.9, 0.1, $-$ 0.1)	0.053	0.048	0.098	0.111
	(0.8, 0.5, 0.5, 0.1)	0.052	0.049	0.087	0.108
(10, 10, 10)	(1, 1, 1, 0)	0.045	0.041	0.057	0.063
	(1, 1, 0.1, 0.1)	0.058	0.049	0.062	0.069
	(1, 0.1, 0.1, 0.5)	0.055	0.047	0.064	0.069
	(0.2, 0.6, 0.9, $-$ 0.3)	0.057	0.052	0.065	0.070
	(0.6, 0.6, 0.6, 0)	0.056	0.053	0.068	0.072
	(0.3, 0.9, 0.1, $-$ 0.1)	0.053	0.046	0.064	0.073
	(0.8, 0.5, 0.5, 0.1)	0.052	0.049	0.060	0.072
(10, 10, 40)	(1, 1, 1, 0)	0.060	0.052	0.063	0.131
	(1, 1, 0.1, 0.1)	0.056	0.051	0.064	0.133
	(1, 0.1, 0.1, 0.5)	0.053	0.049	0.065	0.136
	(0.2, 0.6, 0.9, $-$ 0.3)	0.054	0.053	0.063	0.123
	(0.6, 0.6, 0.6, 0)	0.052	0.048	0.062	0.121
	(0.3, 0.9, 0.1, $-$ 0.1)	0.053	0.051	0.061	0.126
	(0.8, 0.5, 0.5, 0.1)	0.052	0.052	0.063	0.131
(25, 20, 20)	(1, 1, 1, 0)	0.056	0.052	0.050	0.057
	(1, 1, 0.1, 0.1)	0.052	0.049	0.053	0.054
	(1, 0.1, 0.1, 0.5)	0.047	0.047	0.055	0.063
	(0.2, 0.6, 0.9, $-$ 0.3)	0.054	0.053	0.051	0.057
	(0.6, 0.6, 0.6, 0)	0.052	0.052	0.052	0.045
	(0.3, 0.9, 0.1, $-$ 0.1)	0.059	0.055	0.052	0.053
	(0.8, 0.5, 0.5, 0.1)	0.048	0.042	0.051	0.042
$k=$ 5, $p=$ 3, ${\sum }_1=I_3$ , ${\sum }_2=\text{diag}({\lambda }_1,{\lambda }_2,{\lambda }_3)$ , ${\sum }_3=\left(\begin{array}{cccccccccccccccccccc}\hfill 1\hfill & \hfill {\rho }_i\hfill & \hfill {\rho }_i\hfill & & & & & & & & & & & & & & & & & \\ \hfill {\rho }_i\hfill & \hfill 1\hfill & \hfill {\rho }_i\hfill & & & & & & & & & & & & & & & & & \\ \hfill {\rho }_i\hfill & \hfill {\rho }_i\hfill & \hfill 1\hfill & & & & & & & & & & & & & & & & & \end{array}\right)$ , $i=3,4,5$
$(n_1,n_2,n_3,n_4,n_5)$	$({\lambda }_1,{\lambda }_2,{\lambda }_3,{\rho }_3,{\rho }_4,{\rho }_5)$	PB	IPB	Johansen	GV
(7, 7, 7, 7, 7)	(1, 1, 1, 0.5, 0.5, 0.5)	0.047	0.045	0.112	0.208
	(0.1, 0.1, 0.1, 0.3, 0.3, 0.3)	0.048	0.044	0.172	0.211
	(0.1, 0.4, 0.7, 0, 0, 0)	0.050	0.046	0.132	0.222
	(0.1, 0.3, 0.9, 0.1, 0.4, 0.9)	0.051	0.046	0.145	0.232
	(0.1, 0.2, 0.3, $-$ 0.1, 0.1, 0.9)	0.049	0.045	0.134	0.225
	(0.4, 0.4, 0.5, $-$ 0.3, 0.4, 0.3)	0.048	0.043	0.135	0.226
	(0.9, 0.9, 0.9, $-$ 0.4, 0.6, 0.9)	0.049	0.046	0.130	0.221
(12, 12, 12, 12, 12)	(1, 1, 1, 0.5, 0.5, 0.5)	0.047	0.045	0.073	0.133
	(0.1, 0.1, 0.1, 0.3, 0.3, 0.3)	0.051	0.050	0.083	0.143
	(0.1, 0.4, 0.7, 0, 0, 0)	0.049	0.046	0.074	0.124
	(0.1, 0.3, 0.9, 0.1, 0.4, 0.9)	0.050	0.048	0.081	0.134
	(0.1, 0.2, 0.3, $-$ 0.1, 0.1, 0.9)	0.051	0.049	0.082	0.121
	(0.4, 0.4, 0.5, $-$ 0.3, 0.4, 0.3)	0.050	0.050	0.067	0.125
	(0.9, 0.9, 0.9, $-$ 0.4, 0.6, 0.9)	0.050	0.050	0.073	0.141
(20, 20, 20, 20, 20)	(1, 1, 1, 0.5, 0.5, 0.5)	0.047	0.045	0.061	0.086
	(0.1, 0.1, 0.1, 0.3, 0.3, 0.3)	0.052	0.049	0.065	0.105

Table 2, continued
$k=$ 5, $p=$ 3, ${\sum }_1=I_3$ , ${\sum }_2=\text{diag}({\lambda }_1,{\lambda }_2,{\lambda }_3)$ , ${\sum }_3=\left(\begin{array}{cccccccccccccccccccc}\hfill 1\hfill & \hfill {\rho }_i\hfill & \hfill {\rho }_i\hfill & & & & & & & & & & & & & & & & & \\ \hfill {\rho }_i\hfill & \hfill 1\hfill & \hfill {\rho }_i\hfill & & & & & & & & & & & & & & & & & \\ \hfill {\rho }_i\hfill & \hfill {\rho }_i\hfill & \hfill 1\hfill & & & & & & & & & & & & & & & & & \end{array}\right)$ , $i=3,4,5$
$(n_1,n_2,n_3,n_4,n_5)$	$({\lambda }_1,{\lambda }_2,{\lambda }_3,{\rho }_3,{\rho }_4,{\rho }_5)$	PB	IPB	Johansen	GV
(20, 20, 20, 20, 20)	(0.1, 0.4, 0.7, 0, 0, 0)	0.048	0.048	0.065	0.076
	(0.1, 0.3, 0.9, 0.1, 0.4, 0.9)	0.051	0.049	0.058	0.089
	(0.1, 0.2, 0.3, $-$ 0.1, 0.1, 0.9)	0.051	0.050	0.061	0.088
	(0.4, 0.4, 0.5, $-$ 0.3, 0.4, 0.3)	0.051	0.052	0.062	0.089
	(0.9, 0.9, 0.9, $-$ 0.4, 0.6, 0.9)	0.050	0.051	0.063	0.089
(15, 20, 10, 32, 7)	(1, 1, 1, 0.5, 0.5, 0.5)	0.049	0.045	0.092	0.158
	(0.1, 0.1, 0.1, 0.3, 0.3, 0.3)	0.051	0.049	0.126	0.234
	(0.1, 0.4, 0.7, 0, 0, 0)	0.050	0.051	0.116	0.156
	(0.1, 0.3, 0.9, 0.1, 0.4, 0.9)	0.050	0.050	0.086	0.165
	(0.1, 0.2, 0.3, $-$ 0.1, 0.1, 0.9)	0.052	0.050	0.101	0.145
	(0.4, 0.4, 0.5, $-$ 0.3, 0.4, 0.3)	0.051	0.051	0.111	0.153
	(0.9, 0.9, 0.9, $-$ 0.4, 0.6, 0.9)	0.050	0.051	0.076	0.152

In Table 2, the results indicate that the Type I error rates of the PB test and IPB test, in the case of balanced sample sizes and unbalanced sample sizes, control Type I error rates very well under $k=$ 5 and $p=$ 3, but some Type I error rates of the Johansen test and GV test exceed 0.1. However, the four tests control Type I error rates satisfactorily in the case of the unbalanced sample sizes $n=(25,20,20)$ , $k=$ 5 and $p=$ 3. As observed for balanced and unbalanced sample sizes, for $k=$ 5 and $p=$ 3, the results show that the PB test and IPB test control quite well, but the Type I error rates of the Johansen test and GV test exceed the nominal level 0.05. Among the four tests, both the PB test and the IPB test are the best test statistics in terms of Type I error rates, however the IPB test inhibited Type I error rates quite well.

Table 3

Monte Carlo estimates of power of the test for four tests

	${\mu }_1=0,({\lambda }_1,{\lambda }_2,{\lambda }_3,\rho )=(1,0.1,0.1,0.5),k=3,p=3,{\sum }_1=I_3$
${\sum }_2=\text{diag}({\lambda }_1,{\lambda }_2,{\lambda }_3)$ , ${\sum }_3=\left(\begin{array}{cccccccccccccccccccc}\hfill 1\hfill & \hfill \rho \hfill & \hfill \rho \hfill & & & & & & & & & & & & & & & & & \\ \hfill \rho \hfill & \hfill 1\hfill & \hfill \rho \hfill & & & & & & & & & & & & & & & & & \\ \hfill \rho \hfill & \hfill \rho \hfill & \hfill 1\hfill & & & & & & & & & & & & & & & & & \end{array}\right)$	${\mu }_2=0$ , ${\mu }_3=0$	${\mu }_2=0$ , ${\mu }_3=0.2$	${\mu }_2=0$ ,
${\mu }_3=0.5$	${\mu }_2=0$ , ${\mu }_3=0.7$	${\mu }_2=0.5$ , ${\mu }_3=1$	${\mu }_2=1.5$ , ${\mu }_3=1$
$n=(7,7,7)$
PB	0.049	0.252	0.355	0.524	0.742	1
IPB	0.047	0.284	0.405	0.611	0.802	1
Johansen	0.105	0.304	0.444	0.702	0.823	1
GV	0.096	0.301	0.424	0.700	0.820	1
$n=(7,10,20)$
PB	0.058	0.342	0.355	0.644	0.901	1
IPB	0.051	0.374	0.405	0.721	0.912	1
Johansen	0.091	0.424	0.444	0.742	0.914	1
GV	0.114	0.421	0.424	0.743	0.917	1
$n=(20,20,20)$
PB	0.049	0.352	0.412	0.727	0.812	1
IPB	0.047	0.374	0.423	0.742	0.825	1
Johansen	0.051	0.367	0.434	0.745	0.827	1
GV	0.091	0.364	0.435	0.747	0.824	1

In Table 3, we present the power of four tests for $k=$ 3. Upon examination, we find that the following tests as PB test, IPB test, PB test, Johansen test, and GV test demonstrate significant power in the case of sample sizes $n=$ (7, 7, 7), $n=$ (7, 10, 20) and $n=$ (20, 20, 20). Notably, the power of the IPB test in Table 3 surpasses that of the PB test, especially when sample sizes are equal or unequal under unequal common covariance matric.

4. Illustrative examples

We used a data set for explaining the four tests by Bryan and Jorge (2016). The data set consisted of fives samples of 30 skulls for each sample, where there are the early predynastic period (circa 4000 BC), the late predynastic period (circa 3300 BC), the 12th and 13th dynasties (circa 1850 BC), the ptolemaic period (circa 200 BC), and the Roman period (circa AD 150). For each time period, 30 skulls are measured with four variables, being maximal breadth, basibregmatic height, basialveolar length and nasal height. We conducted using a situation by Krishnamoorthy and Yu (2010), which is $(n_1,n_2,n_3,n_4)=(15,15,15,15)$ , the number of groups is $k=$ 4, and the number of variables is $p=$ 4. Hence, the hypothesis testing is

$H_0:{\mu }_1={\mu }_2={\mu }_3={\mu }_4\text{ versus }{H}_1:{\mu }_i\ne {\mu }_j\text{ for some }i\ne j.$

Where ${\mu }_1=\left(\begin{array}{cccccccccccccccccccc}\hfill {\mu }_{11}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{12}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{13}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{14}\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right)$ , ${\mu }_2=\left(\begin{array}{cccccccccccccccccccc}\hfill {\mu }_{21}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{22}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{23}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{24}\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right)$ , ${\mu }_3=\left(\begin{array}{cccccccccccccccccccc}\hfill {\mu }_{31}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{32}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{33}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{34}\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right)$ and ${\mu }_4=\left(\begin{array}{cccccccccccccccccccc}\hfill {\mu }_{41}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{42}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{43}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {\mu }_{44}\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right)$ . The Johansen’s test was calculated by Krishnamoorthy and Yu (2010) where the $p$ -value is computed as 0.0304. Therefore, we will calculate the GV test, PB test and IPB test as follows.

4.1 Generalized variable (GV) test

In the first step, we calculate statistics as ${\stackrel{\sim }{y}}_i$ and ${\stackrel{\sim }{s}}_i$ . After that the GV test in Eq. (9) has 10,000 values generated, and then the p value of the GV test is estimated by the proportion of these 10,000 generated values that are greater than 1. We have the $p$ value of the GV test as 0.001, which is the null hypothesis in (20). Hence, the Egyptian skulls have statistically significant changes over those four periods.

4.2 Parametric bootstrap (PB) test

The Cholesky factors $t_i$ are calculated where ${\stackrel{\sim }{s}}_i=t_i{t}_i^{\prime }(i=1,2,3,4)$ . The Cholesky factors $t_i$ are as follows.

$t_1=\left(\begin{array}{cccccccccccccccccccc}\hfill 1.542\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill 0.320\hfill & \hfill 1.330\hfill & \hfill 0\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill 0.264\hfill & \hfill -0.374\hfill & \hfill 1.465\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill 0.486\hfill & \hfill 0.004\hfill & \hfill -0.129\hfill & \hfill 0.534\hfill & & & & & & & & & & & & & & & & \end{array}\right),t_2=\left(\begin{array}{cccccccccccccccccccc}\hfill 1.433\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill -0.226\hfill & \hfill 1.361\hfill & \hfill 0\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill 0.216\hfill & \hfill 0.309\hfill & \hfill 0.717\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill 0.050\hfill & \hfill 0.561\hfill & \hfill 0.229\hfill & \hfill 0.607\hfill & & & & & & & & & & & & & & & & \end{array}\right)$ $t_3=\left(\begin{array}{cccccccccccccccccccc}\hfill 1.053\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill -0.167\hfill & \hfill 1.305\hfill & \hfill 0\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill -0.128\hfill & \hfill 0.005\hfill & \hfill 1.299\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill -0.048\hfill & \hfill 0.159\hfill & \hfill 0.185\hfill & \hfill 0.794\hfill & & & & & & & & & & & & & & & & \end{array}\right),t_4=\left(\begin{array}{cccccccccccccccccccc}\hfill 0.871\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill 0.036\hfill & \hfill 1.206\hfill & \hfill 0\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill -0.062\hfill & \hfill 0.428\hfill & \hfill 1.314\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill 0.173\hfill & \hfill 0.441\hfill & \hfill 0.219\hfill & \hfill 0.678\hfill & & & & & & & & & & & & & & & & \end{array}\right)$

when the Cholesky factors $t_i(i=1,2,3,4)$ are calculated. After that the $p$ -value of the PB test is estimated with a simulation consisting of 10,000 runs, then the $p$ -value of the PB test is equal to 0.044. Therefore, we reject the null hypothesis in (20). Hence, the Egyptian skulls have statistically significant changes over those four periods.

4.3 Parametric bootstrap (PB) test

The PB test is estimated from a simulation consisting of 10,000 runs. We estimated the $p$ value with the bootstrap approach consisting of 10,000 times. We have the sample bootstrap approach as $X^{\ast }=(P{B}_1^{\ast },\mathrm{\dots },P{B}_{i-1}^{\ast },P{B}_i^{\ast },P{B}_{i+1}^{\ast },\mathrm{\dots },P{B}_{10000}^{\ast })$ . After that the values in $X^{\ast }$ are compared with the test statistic in Eq. (15), when the values in $X^{\ast }$ are more than the test statistic in Eq. (15), set $Q_j=1$ . The $p$ -value is ${\mathit{\text{p-value}}}_i=\frac{\sum _{j=1}^{10000}{Q}_j}{10000}$ . The $p$ -value repeated with 10,000 times as ${\mathit{\text{p-value}}}_{\text{𝐼𝑃𝐵}}=\frac{\sum _{i=1}^{10000}{\mathit{\text{p-value}}}_i}{10,000}$ , then PB $p$ -value is obtained as 0.041. Therefore, we reject the null hypothesis in Eq. (20). Hence, the Egyptian skulls have statistically significant changes over those four periods.

5. Conclusion

The Hotelling T² test for MANOVA is a test statistic based on equal common covariance matrices $\sum$ , and it could not control Type I error rates under the Behrens-Fisher problem (unequal common covariance matric ${\sum }_i)$ . In this research, we propose the IBP test and compare it with three tests (PB test, GV test and Johansen’s test) based on Type I error rates. For $p=$ 2, the PB test and IPB test control Type I error rates very well, which is better than the GV test and Johansen’s test at the nominal level 0.05 under the Behrens-Fisher problem. When $p=$ 3, the PB test and IPB test control Type I error rates quite well when the nominal level is 0.05. For $k=$ 10, the PB test and IPB test control Type I error rates quite well at the nominal level 0.05 under the Behrens-Fisher problem. Moreover, power of the four tests appear that IPB test, Johansen test, and GV test are very well power of the test. From the results of the Type I error rates and power of the tests, we suggest that the IBF test be used as an alternative approach for comparing the mean vectors with multivariate normal populations under the Behrens-Fisher problem in MANOVA, since it can control Type I error rates and power of the test very well.