Multivariate analysis of covariance when standard assumptions are violated

TL;DR

This paper develops a flexible multivariate analysis of covariance method that remains valid under violations of classical assumptions, including heteroskedasticity and singular covariance matrices, using bootstrap techniques and simulations.

Contribution

It introduces a modified MANCOVA ANOVA-type statistic (MANCATS) combined with bootstrap methods, extending applicability to non-normal, heteroskedastic, and singular covariance scenarios.

Findings

Parametric bootstrap MANCATS controls type I error well.

Bootstrap MANCATS shows higher power in simulations.

Method applied successfully to real achievement test data.

Abstract

In applied research, it is often sensible to account for one or several covariates when testing for differences between multivariate means of several groups. However, the "classical" parametric multivariate analysis of covariance (MANCOVA) tests (e.g., Wilks' Lambda) are based on quite restrictive assumptions (homoscedasticity and normality of the errors), which might be difficult to justify in small sample size settings. Furthermore, existing potential remedies (e.g., heteroskedasticity-robust approaches) become inappropriate in cases where the covariance matrices are singular. Nevertheless, such scenarios are frequently encountered in the life sciences and other fields, when for example, in the context of standardized assessments, a summary performance measure as well as its corresponding subscales are analyzed. In the present manuscript, we consider a general MANCOVA model, allowing for potentially heteroskedastic and even singular covariance matrices as well as non-normal errors. We combine heteroskedasticity-consistent covariance matrix estimation methods with our proposed modified MANCOVA ANOVA-type statistic (MANCATS) and apply two different bootstrap approaches. We provide the proofs of the asymptotic validity of the respective testing procedures as well as the results from an extensive simulation study, which indicate that especially the parametric bootstrap version of the MANCATS outperforms its competitors in most scenarios, both in terms of type I error rates and power. These considerations are further illustrated and substantiated by examining real-life data from standardized achievement tests.

Tables

Table 1. Table 1: Empirical type I error rates (in %) of the Wilks’ Lambda test (WI), the Wald-type test (WT), the wild bootstrap MANCATS test (MW), and the parametric bootstrap MANCATS test (MP) for a = 2 𝑎 2 a=2 groups, p = 2 𝑝 2 p=2 dimensions, c = 2 𝑐 2 c=2 fixed covariates, and different covariance scenarios, namely I: 𝚺 i = I 2 + 0.5 ​ ( J 2 − I 2 ) subscript 𝚺 𝑖 subscript 𝐼 2 0.5 subscript 𝐽 2 subscript 𝐼 2 {\bf\Sigma}_{i}=I_{2}+0.5(J_{2}-I_{2}) , II: 𝚺 i = i ⋅ I 2 + 0.5 ​ ( J 2 − I 2 ) subscript 𝚺 𝑖 ⋅ 𝑖 subscript 𝐼 2 0.5 subscript 𝐽 2 subscript 𝐼 2 {\bf\Sigma}_{i}=i\cdot I_{2}+0.5(J_{2}-I_{2}) , III: 𝚺 i = d ​ i ​ a ​ g ​ ( 1 , 0.25 ) + 0.5 ​ ( J 2 − I 2 ) subscript 𝚺 𝑖 𝑑 𝑖 𝑎 𝑔 1 0.25 0.5 subscript 𝐽 2 subscript 𝐼 2 {\bf\Sigma}_{i}=diag(1,0.25)+0.5(J_{2}-I_{2}) , i ∈ { 1 , 2 } 𝑖 1 2 i\in\{1,2\} . Errors were drawn either from standard normal, χ 5 2 superscript subscript 𝜒 5 2 \chi_{5}^{2} , standard lognormal, or double exponential distribution.

		Normal				Chi-squared(5)				Lognormal				Double exp.
$𝚺$	$(n_1,n_2)$	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP
I	$(20,20)$	$5.0$	$8.4$	$6.6$	$5.3$	$5.0$	$8.2$	$6.4$	$4.8$	$5.0$	$6.1$	$5.3$	$4.5$	$5.0$	$7.7$	$5.8$	$4.8$
	$(10,10)$	$5.3$	$12.0$	$7.2$	$5.0$	$5.0$	$11.8$	$6.7$	$4.8$	$4.9$	$8.1$	$4.6$	$4.0$	$4.7$	$11.2$	$6.5$	$4.8$
	$(10,20)$	$5.4$	$12.4$	$7.0$	$5.0$	$4.9$	$11.6$	$7.0$	$5.3$	$5.4$	$8.2$	$5.4$	$5.0$	$4.8$	$11.2$	$6.9$	$5.1$
	$(20,10)$	$5.2$	$10.8$	$7.5$	$5.1$	$5.2$	$10.3$	$7.6$	$5.1$	$5.9$	$6.9$	$5.8$	$4.4$	$5.1$	$10.2$	$7.0$	$4.8$
II	$(20,20)$	$2.7$	$8.1$	$6.1$	$5.0$	$3.0$	$7.3$	$6.0$	$4.6$	$3.6$	$5.1$	$4.8$	$4.5$	$2.7$	$7.1$	$5.6$	$4.7$
	$(10,10)$	$3.4$	$11.3$	$6.7$	$4.9$	$3.7$	$11.0$	$6.3$	$4.6$	$4.1$	$7.6$	$4.3$	$3.9$	$3.5$	$10.6$	$6.2$	$4.5$
	$(10,20)$	$3.7$	$12.9$	$7.5$	$5.2$	$3.6$	$12.3$	$7.4$	$5.2$	$4.4$	$9.5$	$5.7$	$4.9$	$3.5$	$11.8$	$7.3$	$5.2$
	$(20,10)$	$11.5$	$12.3$	$7.6$	$5.0$	$10.9$	$12.6$	$7.4$	$5.0$	$10.1$	$10.1$	$6.6$	$4.6$	$10.9$	$11.9$	$6.9$	$4.5$
III	$(20,20)$	$-$	$2.5$	$6.4$	$5.4$	$-$	$2.4$	$6.3$	$5.2$	$-$	$1.6$	$5.7$	$4.7$	$-$	$2.1$	$5.7$	$4.8$
	$(10,10)$	$-$	$3.3$	$7.4$	$5.0$	$-$	$3.0$	$6.8$	$4.9$	$-$	$2.2$	$5.3$	$4.2$	$-$	$3.2$	$6.9$	$5.0$
	$(10,20)$	$-$	$3.4$	$7.3$	$4.9$	$-$	$3.0$	$7.4$	$5.1$	$-$	$2.2$	$5.5$	$4.9$	$-$	$3.1$	$6.9$	$5.0$
	$(20,10)$	$-$	$2.8$	$7.6$	$5.0$	$-$	$2.7$	$7.6$	$5.0$	$-$	$1.7$	$5.9$	$4.6$	$-$	$2.6$	$6.7$	$4.5$

Table 2. Table 2: Empirical type I error rates (in %) of the Wilks’ Lambda test (WI), the Wald-type test (WT), the wild bootstrap MANCATS test (MW), and the parametric bootstrap MANCATS test (MP) for a = 4 𝑎 4 a=4 groups, p = 2 𝑝 2 p=2 dimensions and c = 2 𝑐 2 c=2 fixed covariates, and different covariance scenarios, namely I: 𝚺 i = I 2 + 0.5 ​ ( J 2 − I 2 ) subscript 𝚺 𝑖 subscript 𝐼 2 0.5 subscript 𝐽 2 subscript 𝐼 2 {\bf\Sigma}_{i}=I_{2}+0.5(J_{2}-I_{2}) , II: 𝚺 i = i ⋅ I 2 + 0.5 ​ ( J 2 − I 2 ) subscript 𝚺 𝑖 ⋅ 𝑖 subscript 𝐼 2 0.5 subscript 𝐽 2 subscript 𝐼 2 {\bf\Sigma}_{i}=i\cdot I_{2}+0.5(J_{2}-I_{2}) , III: 𝚺 i = d ​ i ​ a ​ g ​ ( 1 , 0.25 ) + 0.5 ​ ( J 2 − I 2 ) subscript 𝚺 𝑖 𝑑 𝑖 𝑎 𝑔 1 0.25 0.5 subscript 𝐽 2 subscript 𝐼 2 {\bf\Sigma}_{i}=diag(1,0.25)+0.5(J_{2}-I_{2}) , i ∈ { 1 , 2 , 3 , 4 } 𝑖 1 2 3 4 i\in\{1,2,3,4\} . The group sizes were 𝐧 1 = ( 30 , 30 , 30 , 30 ) subscript 𝐧 1 30 30 30 30 {\bf n}_{1}=(30,30,30,30) , 𝐧 2 = ( 15 , 15 , 15 , 15 ) subscript 𝐧 2 15 15 15 15 {\bf n}_{2}=(15,15,15,15) , 𝐧 3 = ( 5 , 10 , 20 , 25 ) subscript 𝐧 3 5 10 20 25 {\bf n}_{3}=(5,10,20,25) , and 𝐧 4 = ( 25 , 20 , 10 , 5 ) subscript 𝐧 4 25 20 10 5 {\bf n}_{4}=(25,20,10,5) . Errors were drawn either from standard normal, χ 5 2 subscript superscript 𝜒 2 5 \chi^{2}_{5} , standard lognormal, or double exponential distribution.

		Normal				Chi-squared(5)				Lognormal				Double exp.
$𝚺$	$N$	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP
I	${𝐧}_1$	$4.2$	$8.5$	$5.1$	$4.3$	$5.2$	$8.8$	$5.2$	$4.2$	$4.7$	$6.5$	$5.3$	$4.6$	$5.0$	$9.1$	$5.8$	$4.9$
	${𝐧}_2$	$5.3$	$13.6$	$6.0$	$4.6$	$5.1$	$12.0$	$5.4$	$4.1$	$4.8$	$8.6$	$5.7$	$4.8$	$4.9$	$14.0$	$6.0$	$4.6$
	${𝐧}_3$	$5.0$	$13.6$	$5.7$	$5.6$	$4.9$	$13.8$	$5.6$	$5.8$	$5.1$	$9.4$	$5.2$	$5.7$	$5.1$	$13.0$	$5.6$	$5.6$
	${𝐧}_4$	$5.0$	$13.7$	$5.3$	$5.6$	$5.2$	$12.4$	$5.1$	$5.9$	$4.9$	$8.1$	$4.8$	$6.3$	$4.9$	$13.0$	$5.2$	$5.4$
II	${𝐧}_1$	$4.2$	$8.9$	$5.0$	$4.1$	$4.7$	$9.3$	$5.5$	$4.4$	$4.4$	$7.2$	$5.6$	$4.6$	$4.5$	$9.1$	$5.7$	$4.7$
	${𝐧}_2$	$4.9$	$13.9$	$5.8$	$4.4$	$5.0$	$12.3$	$5.2$	$3.9$	$4.8$	$9.6$	$5.7$	$5.0$	$4.6$	$14.3$	$6.0$	$4.4$
	${𝐧}_3$	$2.7$	$13.2$	$5.4$	$6.0$	$2.5$	$13.2$	$5.5$	$6.2$	$2.6$	$10.1$	$5.0$	$5.8$	$2.5$	$12.5$	$5.4$	$5.8$
	${𝐧}_4$	$13.5$	$14.2$	$4.8$	$4.8$	$12.8$	$12.6$	$4.9$	$5.4$	$10.9$	$8.5$	$5.4$	$7.1$	$12.6$	$13.6$	$4.7$	$4.8$
III	${𝐧}_1$	$-$	$1.1$	$5.2$	$4.5$	$-$	$1.0$	$5.1$	$4.3$	$-$	$0.6$	$5.0$	$4.6$	$-$	$1.1$	$5.7$	$5.1$
	${𝐧}_2$	$-$	$2.0$	$5.9$	$4.5$	$-$	$1.7$	$5.3$	$4.3$	$-$	$1.1$	$5.5$	$4.6$	$-$	$2.0$	$5.7$	$4.6$
	${𝐧}_3$	$-$	$2.0$	$5.9$	$5.6$	$-$	$1.9$	$5.7$	$5.4$	$-$	$1.0$	$5.2$	$5.8$	$-$	$1.8$	$5.5$	$5.6$
	${𝐧}_4$	$-$	$2.1$	$5.3$	$5.4$	$-$	$1.8$	$5.1$	$5.7$	$-$	$0.9$	$4.8$	$5.9$	$-$	$1.7$	$5.2$	$5.3$

Table 3. Table 3: Empirical type I error rates (in %) of the Wilks’ Lambda test (WI), the Wald-type test (WT), the wild bootstrap MANCATS test (MW), and the parametric bootstrap MANCATS test (MP) for a = 2 𝑎 2 a=2 groups, p = 4 𝑝 4 p=4 dimensions, c = 2 𝑐 2 c=2 fixed covariates, and different covariance scenarios, namely I: 𝚺 i = I 4 + 0.5 ​ ( J 4 − I 4 ) subscript 𝚺 𝑖 subscript 𝐼 4 0.5 subscript 𝐽 4 subscript 𝐼 4 {\bf\Sigma}_{i}=I_{4}+0.5(J_{4}-I_{4}) , II: 𝚺 i = i ⋅ I 4 + 0.5 ​ ( J 4 − I 4 ) subscript 𝚺 𝑖 ⋅ 𝑖 subscript 𝐼 4 0.5 subscript 𝐽 4 subscript 𝐼 4 {\bf\Sigma}_{i}=i\cdot I_{4}+0.5(J_{4}-I_{4}) , i ∈ { 1 , 2 , 3 , 4 } 𝑖 1 2 3 4 i\in\{1,2,3,4\} . III: 𝚺 i = ( 𝐬 1 , 𝐬 2 , 𝐬 3 , 𝐬 4 ) subscript 𝚺 𝑖 subscript 𝐬 1 subscript 𝐬 2 subscript 𝐬 3 subscript 𝐬 4 {\bf\Sigma}_{i}=({\bf s}_{1},{\bf s}_{2},{\bf s}_{3},{\bf s}_{4}) , where 𝐬 1 = 𝐬 2 = 𝐬 3 = ( 1 , 1 , 1 , 1 / 2 ) ′ subscript 𝐬 1 subscript 𝐬 2 subscript 𝐬 3 superscript 1 1 1 1 2 ′ {\bf s}_{1}={\bf s}_{2}={\bf s}_{3}=(1,1,1,1/2)^{\prime} , 𝐬 4 = ( 1 / 2 , 1 / 2 , 1 / 2 , 1 ) ′ subscript 𝐬 4 superscript 1 2 1 2 1 2 1 ′ {\bf s}_{4}=(1/2,1/2,1/2,1)^{\prime} . Errors were drawn either from standard normal, χ 5 2 subscript superscript 𝜒 2 5 \chi^{2}_{5} , standard lognormal, or double exponential distribution.

		Normal				Chi-squared(5)				Lognormal				Double exp.
$𝚺$	$(n_1,n_2)$	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP
I	$(20,20)$	$5.1$	$15.8$	$6.6$	$5.2$	$4.8$	$14.3$	$6.4$	$5.2$	$4.8$	$11.3$	$5.7$	$4.6$	$5.2$	$14.7$	$6.5$	$5.2$
	$(10,10)$	$5.3$	$28.7$	$7.2$	$5.0$	$5.0$	$27.8$	$6.5$	$4.6$	$4.7$	$23.2$	$5.1$	$4.2$	$5.1$	$27.9$	$6.9$	$4.7$
	$(10,20)$	$5.0$	$27.7$	$7.3$	$5.2$	$5.1$	$26.7$	$7.2$	$5.4$	$5.5$	$21.4$	$6.3$	$5.3$	$5.2$	$27.0$	$7.2$	$5.1$
	$(20,10)$	$5.0$	$27.0$	$7.6$	$5.0$	$5.3$	$24.3$	$7.2$	$4.7$	$6.0$	$22.2$	$6.8$	$4.8$	$5.2$	$25.1$	$7.3$	$4.4$
II	$(20,20)$	$2.3$	$13.7$	$6.4$	$4.9$	$2.4$	$12.3$	$6.3$	$5.0$	$2.9$	$9.5$	$5.2$	$4.6$	$2.5$	$12.9$	$6.0$	$5.0$
	$(10,10)$	$3.7$	$27.0$	$6.6$	$4.8$	$3.6$	$25.9$	$6.0$	$4.3$	$3.9$	$20.7$	$4.5$	$4.1$	$3.5$	$25.7$	$6.3$	$4.5$
	$(10,20)$	$3.2$	$28.7$	$7.2$	$5.0$	$3.2$	$28.2$	$7.4$	$5.3$	$4.2$	$21.8$	$5.6$	$4.9$	$3.4$	$28.5$	$7.3$	$5.2$
	$(20,10)$	$15.1$	$33.0$	$7.4$	$4.5$	$14.4$	$30.9$	$6.9$	$4.0$	$12.9$	$26.6$	$6.5$	$4.3$	$14.9$	$15.1$	$6.7$	$3.9$
III	$(20,20)$	$-$	$2.5$	$6.3$	$5.1$	$-$	$2.5$	$6.4$	$5.0$	$-$	$1.6$	$5.7$	$5.0$	$-$	$2.8$	$6.3$	$5.2$
	$(10,10)$	$-$	$5.7$	$7.4$	$5.1$	$-$	$5.2$	$7.1$	$4.7$	$-$	$3.7$	$5.5$	$4.3$	$-$	$5.0$	$6.8$	$4.8$
	$(10,20)$	$-$	$5.2$	$7.4$	$4.9$	$-$	$5.2$	$7.4$	$5.3$	$-$	$3.0$	$5.8$	$5.0$	$-$	$5.0$	$7.3$	$5.0$
	$(20,10)$	$-$	$4.6$	$7.6$	$5.0$	$-$	$4.0$	$7.2$	$4.7$	$-$	$2.8$	$6.4$	$4.5$	$-$	$4.0$	$7.0$	$4.8$

Table 4. Table 4: Empirical type I error rates (in %) of the Wilks’ Lambda test (WI), the Wald-type test (WT), the wild bootstrap MANCATS test (MW), and the parametric bootstrap MANCATS test (MP) for a = 2 𝑎 2 a=2 groups (low vs. high), c = 1 𝑐 1 c=1 covariate (sum of the learning task results), and the following covariance matrix settings: Regular scenario (R), where 𝚺 1 = 𝚺 ^ h ( r ) subscript 𝚺 1 superscript subscript ^ 𝚺 ℎ 𝑟 {\bf\Sigma}_{1}=\hat{{\bf\Sigma}}_{h}^{(r)} and 𝚺 2 = 𝚺 ^ ℓ ( r ) subscript 𝚺 2 superscript subscript ^ 𝚺 ℓ 𝑟 {\bf\Sigma}_{2}=\hat{{\bf\Sigma}}_{\ell}^{(r)} ; singular scenario (S), with 𝚺 1 = 𝚺 ^ h ( s ) subscript 𝚺 1 superscript subscript ^ 𝚺 ℎ 𝑠 {\bf\Sigma}_{1}=\hat{{\bf\Sigma}}_{h}^{(s)} and 𝚺 2 = 𝚺 ^ ℓ ( s ) subscript 𝚺 2 superscript subscript ^ 𝚺 ℓ 𝑠 {\bf\Sigma}_{2}=\hat{{\bf\Sigma}}_{\ell}^{(s)} . Errors were drawn either from standard normal, χ 5 2 subscript superscript 𝜒 2 5 \chi^{2}_{5} , standard lognormal, or double exponential distribution.

		Normal				Chi-square(5)				Lognormal				Double exp.
$𝚺$	$(n_1,n_2)$	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP
R	$(32,37)$	$5.9$	$5.9$	$5.1$	$4.9$	$6.0$	$6.1$	$5.1$	$5.0$	$5.6$	$6.2$	$6.0$	$5.0$	$5.7$	$6.2$	$5.5$	$5.3$
	$(23,46)$	$8.4$	$6.7$	$5.5$	$5.3$	$9.0$	$7.0$	$5.6$	$5.2$	$7.5$	$9.8$	$9.2$	$8.0$	$8.2$	$7.0$	$5.7$	$5.6$
	$(46,23)$	$2.8$	$5.7$	$5.1$	$5.0$	$2.6$	$5.7$	$5.1$	$4.9$	$3.4$	$4.4$	$4.9$	$4.0$	$3.2$	$6.0$	$5.1$	$5.1$
S	$(32,37)$	$-$	$2.5$	$5.1$	$4.9$	$-$	$2.7$	$5.5$	$5.4$	$-$	$2.5$	$5.4$	$4.7$	$-$	$2.8$	$5.2$	$5.0$
	$(23,46)$	$-$	$3.0$	$5.1$	$5.1$	$-$	$3.7$	$5.6$	$5.3$	$-$	$4.2$	$7.4$	$6.7$	$-$	$3.1$	$5.7$	$5.6$
	$(46,23)$	$-$	$2.6$	$5.3$	$5.1$	$-$	$2.5$	$4.8$	$4.7$	$-$	$1.9$	$5.0$	$4.3$	$-$	$2.3$	$5.1$	$4.9$

Table 5. Table S1: Empirical type I error rates (in %) of the Wald-type test (WT), the wild bootstrap MANCATS test (MW), and the parametric bootstrap MANCATS test (MP) for a = 2 𝑎 2 a=2 groups, p = 8 𝑝 8 p=8 dimensions, c = 2 𝑐 2 c=2 fixed covariates, and the following homoskedastic singular covariance scenario: 𝚺 i = ( 𝐬 1 , … , 𝐬 8 ) subscript 𝚺 𝑖 subscript 𝐬 1 … subscript 𝐬 8 {\bf\Sigma}_{i}=({\bf s}_{1},\ldots,{\bf s}_{8}) , where 𝐬 1 = … = 𝐬 5 = ( 1 , 1 , 1 , 1 , 1 , 1 / 2 , 1 / 2 , 1 / 2 ) ′ subscript 𝐬 1 … subscript 𝐬 5 superscript 1 1 1 1 1 1 2 1 2 1 2 ′ {\bf s}_{1}=\ldots={\bf s}_{5}=(1,1,1,1,1,1/2,1/2,1/2)^{\prime} , 𝐬 6 = ( 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 , 1 / 2 , 1 / 2 ) ′ subscript 𝐬 6 superscript 1 2 1 2 1 2 1 2 1 2 1 1 2 1 2 ′ {\bf s}_{6}=(1/2,1/2,1/2,1/2,1/2,1,1/2,1/2)^{\prime} , 𝐬 7 = ( 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 , 1 / 2 ) ′ subscript 𝐬 7 superscript 1 2 1 2 1 2 1 2 1 2 1 2 1 1 2 ′ {\bf s}_{7}=(1/2,1/2,1/2,1/2,1/2,1/2,1,1/2)^{\prime} , 𝐬 8 = ( 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 / 2 , 1 ) ′ subscript 𝐬 8 superscript 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 ′ {\bf s}_{8}=(1/2,1/2,1/2,1/2,1/2,1/2,1/2,1)^{\prime} , i ∈ { 1 , 2 } 𝑖 1 2 i\in\{1,2\} . Errors were drawn either from a standard normal, a Chi-square(5), a standard lognormal, or a double exponential distribution.

	Normal			Chi-square(5)			Lognormal			Double exp.
$(n_1,n_2)$	WT	MW	MP	WT	MW	MP	WT	MW	MP	WT	MW	MP
$(20,20)$	$4.2$	$6.1$	$5.0$	$4.1$	$6.4$	$5.1$	$2.7$	$6.6$	$5.3$	$4.2$	$6.1$	$5.1$
$(10,10)$	$13.9$	$7.6$	$5.2$	$12.6$	$7.3$	$5.3$	$9.8$	$5.9$	$4.5$	$13.0$	$7.0$	$4.9$
$(10,20)$	$13.2$	$8.2$	$5.4$	$11.9$	$7.7$	$5.3$	$8.2$	$6.9$	$5.8$	$11.7$	$7.2$	$4.5$
$(20,10)$	$12.5$	$8.0$	$5.1$	$11.0$	$7.9$	$5.5$	$7.6$	$7.3$	$4.7$	$11.2$	$7.1$	$4.9$

Table 6. Table S2: Empirical type I error rates (in %) of the Wald-type test (WT), the wild bootstrap MANCATS test (MW), and the parametric bootstrap MANCATS test (MP) for a = 2 𝑎 2 a=2 groups, p = 2 𝑝 2 p=2 dimensions, and the following homoskedastic singular covariance scenario: 𝚺 i = d ​ i ​ a ​ g ​ ( 1 , 0.25 ) + 0.5 ​ ( J 2 − I 2 ) subscript 𝚺 𝑖 𝑑 𝑖 𝑎 𝑔 1 0.25 0.5 subscript 𝐽 2 subscript 𝐼 2 {\bf\Sigma}_{i}=diag(1,0.25)+0.5(J_{2}-I_{2}) , i ∈ { 1 , 2 } 𝑖 1 2 i\in\{1,2\} . The values of the first covariate were equally spaced between − 10 10 -10 and 10 10 10 , whereas the second covariate was a sample from a standard normal or lognormal distribution. Errors were drawn either from a standard normal, a Chi-square(5), a standard lognormal, or a double exponential distribution.

		Normal			Chi-square(5)			Lognormal			Double exp.
Covariate	$(n_1,n_2)$	WT	MW	MP	WT	MW	MP	WT	MW	MP	WT	MW	MP
Normal	$(20,20)$	$2.2$	$5.8$	$5.0$	$1.8$	$5.6$	$4.8$	$1.3$	$4.8$	$4.6$	$2.0$	$5.7$	$4.8$
	$(10,10)$	$2.9$	$7.0$	$4.7$	$2.7$	$6.8$	$5.2$	$1.7$	$4.8$	$4.3$	$2.8$	$6.7$	$4.9$
	$(10,20)$	$2.4$	$6.2$	$5.2$	$2.5$	$6.1$	$5.4$	$1.3$	$4.6$	$4.6$	$2.3$	$6.1$	$5.3$
	$(20,10)$	$2.2$	$5.6$	$5.1$	$2.6$	$5.9$	$5.3$	$1.5$	$4.6$	$4.7$	$2.2$	$5.5$	$4.8$
Lognormal	$(20,20)$	$1.9$	$5.9$	$5.0$	$1.6$	$5.3$	$4.8$	$1.1$	$4.5$	$4.7$	$1.8$	$5.3$	$4.8$
	$(10,10)$	$2.4$	$6.4$	$4.8$	$2.3$	$6.5$	$5.1$	$1.6$	$4.9$	$4.8$	$2.5$	$6.5$	$5.2$
	$(10,20)$	$2.5$	$5.9$	$5.7$	$2.4$	$6.0$	$5.8$	$1.4$	$4.6$	$5.4$	$2.4$	$6.1$	$5.9$
	$(20,10)$	$2.0$	$5.4$	$4.8$	$2.1$	$6.0$	$4.9$	$1.2$	$4.6$	$4.6$	$1.9$	$5.4$	$4.7$

Table 7. Table S3: Empirical type I error rates (in %) of the Wilks’ Lambda test (WI), the Wald-type test (WT), the wild bootstrap MANCATS test (MW), and the parametric bootstrap MANCATS test (MP) for a = 2 𝑎 2 a=2 groups, p = 2 𝑝 2 p=2 dimensions, c = 2 𝑐 2 c=2 fixed covariates, and the following heteroskedastic covariance scenario: 𝚺 i = i ⋅ I 2 + 0.5 ​ ( J 2 − I 2 ) subscript 𝚺 𝑖 ⋅ 𝑖 subscript 𝐼 2 0.5 subscript 𝐽 2 subscript 𝐼 2 {\bf\Sigma}_{i}=i\cdot I_{2}+0.5(J_{2}-I_{2}) , i ∈ { 1 , 2 } 𝑖 1 2 i\in\{1,2\} . Errors were drawn either from a standard normal, a Chi-square(5), a standard lognormal, or a double exponential distribution. Within each of the n s ​ i ​ m = 1000 subscript 𝑛 𝑠 𝑖 𝑚 1000 n_{sim}=1000 simulation runs, n b ​ o ​ o ​ t = 500 subscript 𝑛 𝑏 𝑜 𝑜 𝑡 500 n_{boot}=500 bootstrap samples were drawn.

	Normal				Chi-square(5)				Lognormal				Double exp.
$(n_1,n_2)$	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP	WI	WT	MW	MP
$(200,200)$	$2.9$	$5.8$	$5.3$	$5.4$	$1.8$	$5.1$	$4.5$	$4.5$	$2.7$	$4.8$	$6.0$	$6.2$	$1.8$	$6.0$	$5.4$	$5.0$
$(100,100)$	$2.4$	$5.9$	$5.2$	$4.1$	$3.3$	$6.6$	$7.7$	$6.9$	$2.6$	$5.0$	$6.1$	$4.5$	$2.2$	$6.8$	$5.6$	$5.5$
$(150,50)$	$9.4$	$6.2$	$4.1$	$3.8$	$9.6$	$5.3$	$5.4$	$4.9$	$9.3$	$4.9$	$5.1$	$6.2$	$9.7$	$6.1$	$5.7$	$4.1$
$(50,150)$	$4.5$	$5.4$	$5.4$	$5.6$	$4.0$	$6.9$	$7.0$	$6.5$	$3.7$	$6.4$	$5.2$	$5.2$	$3.7$	$6.4$	$5.7$	$5.7$

Equations

$${\mathbf{Y}}_{ij}={\mathbf{\mu }}_i+({\mathbf{z}}_{ij}^{\prime }\otimes {\mathbf{I}}_p)\mathbf{\nu }+{\mathbf{ϵ}}_{ij}={\mathbf{\mu }}_i+\sum _{w=1}^c{z}_{ij}^{(w)}{\mathbf{\nu }}^{(w)}+{\mathbf{ϵ}}_{ij}.$$

$$\mathbf{Y}=\tilde{\mathbf{M}}\mathbf{\mu }+\tilde{\mathbf{Z}}\mathbf{\nu }+\mathbf{ϵ},$$

$${\hat{\mathbf{\mu }}}_i={\bar{\mathbf{Y}}}_{i.}-\sum _{w=1}^c{\bar{z}}_{i.}^{(w)}{\hat{\mathbf{\nu }}}^{(w)},$$

$$W(\mathbf{T})={\hat{\mathbf{\mu }}}^{\prime }{\mathbf{T}}^{\prime }(\mathbf{T}\hat{\mathbf{\Sigma }}\mathbf{T}{)}^{+}\mathbf{T}\hat{\mathbf{\mu }},$$

$$A(\mathbf{T}):={\hat{\mathbf{\mu }}}^{\prime }{\mathbf{T}}^{\prime }(\mathbf{T}\hat{\mathbf{D}}\mathbf{T}{)}^{+}\mathbf{T}\hat{\mathbf{\mu }}.$$

$$\hat{\mathbf{D}}=\underset{i=1}{\overset{a}{⨁}}\frac{1}{n_i}\underset{k=1}{\overset{p}{⨁}}{\hat{\sigma }}_{ik}^2=diag(n_1^{-1}{\hat{\sigma }}_{11}^2,\dots ,n_1^{-1}{\hat{\sigma }}_{1p}^2,\dots ,n_a^{-1}{\hat{\sigma }}_{ap}^2),$$

$${\hat{\sigma }}_{ik}^2=\frac{1}{n_i-c-1}\sum _{j=1}^{n_i}{u}_{ijk}^2,$$

$$U=\sum _{i=1}^a\sum _{k=1}^p{\lambda }_{ik}{U}_{ik},$$

$${\mathbf{Y}}_{ij}^{\ast }:={\mathbf{u}}_{ij}\frac{T_{ij}}{\sqrt{1-p_{ij,ij}}},1\le i\le a,1\le j\le n_i,$$

$$\underset{x\in \mathbb{R}}{\sup }|P_{\mathbf{\beta }}(A^{\ast }(\mathbf{T})\le x|\mathbf{Y})-P_{{\mathbf{\beta }}_0}(A(\mathbf{T})\le x)|\stackrel{P}{⟶}0$$

$${\mathbf{Y}}_{\mathbf{ij}}^{\star }\stackrel{i.i.d.}{\sim }\mathcal{N}(0,{\hat{\mathbf{\Sigma }}}_{\mathbf{i}}),1\le i\le a,1\le j\le n_i,$$

$${\hat{\mathbf{\Sigma }}}_{\mathbf{i}}:=\frac{1}{n_i-c-1}\sum _{j=1}^{n_i}{\mathbf{u}}_{\mathbf{ij}}{\mathbf{u}}_{\mathbf{ij}}^{\prime },1\le i\le a.$$

$$\underset{x\in \mathbb{R}}{\sup }|P_{\mathbf{\beta }}(A^{\star }(\mathbf{T})\le x|\mathbf{Y})-P_{{\mathbf{\mu }}_0}(A(\mathbf{T})\le x)|\stackrel{P}{⟶}0$$

$${\mathbf{ϵ}}_{ij}={\mathbf{\Sigma }}_i^{1/2}{\tilde{\mathbf{\xi }}}_{ij},$$

$${\hat{\mathbf{\Sigma }}}_h^{(r)}=\left(\begin{array}{cc}145.88& 113.18\\ 113.18& 1073.21\end{array}\right)\text{and}{\hat{\mathbf{\Sigma }}}_{\ell }^{(r)}=\left(\begin{array}{cc}99.07& 60.85\\ 60.85& 458.41\end{array}\right),$$

$${\hat{\mathbf{\Sigma }}}_h^{(s)}=\left(\begin{array}{ccc}145.88& 113.18& 259.06\\ 113.18& 1073.21& 1186.39\\ 259.06& 1186.39& 1445.45\end{array}\right)\text{and}{\hat{\mathbf{\Sigma }}}_{\ell }^{(s)}=\left(\begin{array}{ccc}99.07& 60.85& 159.92\\ 60.85& 458.41& 519.26\\ 159.92& 519.26& 679.19\end{array}\right),$$

$$\sqrt{N}(\hat{\mathbf{\beta }}-\mathbf{\beta })=\sqrt{N}(({\mathbf{X}}^{\prime }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\prime }\otimes {\mathbf{I}}_p)\mathbf{ϵ}=\sqrt{N}(({\mathbf{X}}^{\prime }\mathbf{X}{)}^{-1}\otimes {\mathbf{I}}_p)\sum _{i=1}^a\sum _{j=1}^{n_i}({\mathbf{x}}_{ij}\otimes {\mathbf{I}}_p){\mathbf{ϵ}}_{ij}.$$

$$\sum _{i=1}^a\sum _{j=1}^{n_i}Cov({\mathbf{\eta }}_{ij})$$

$$=N\sum _{i=1}^a\sum _{j=1}^{n_i}(({\mathbf{X}}^{\prime }\mathbf{X}{)}^{-1}{\mathbf{x}}_{ij}\otimes {\mathbf{I}}_p){\mathbf{\Sigma }}_i({\mathbf{x}}_{ij}^{\prime }({\mathbf{X}}^{\prime }\mathbf{X}{)}^{-1}\otimes {\mathbf{I}}_p)$$

$$=({\mathbf{X}}^{\prime }\mathbf{X}/N{)}^{-1}\otimes {\mathbf{I}}_p)(N^{-1}({\mathbf{X}}^{\prime }\otimes {\mathbf{I}}_p)\mathbf{S}(\mathbf{X}\otimes {\mathbf{I}}_p))({\mathbf{X}}^{\prime }\mathbf{X}/N{)}^{-1}\otimes {\mathbf{I}}_p)$$

$$\to \mathbf{\Lambda }.$$

$$\sum _{i=1}^a\sum _{j=1}^{n_i}E(||{\mathbf{\eta }}_{ij}|{|}^2\text{mathbbm}1\{||{\mathbf{\eta }}_{ij}|{|}^2>\delta \})\to 0\forall \delta >0.$$

$$||{\mathbf{\eta }}_{ij}|{|}^2$$

$$=\frac{1}{N}{\mathbf{ϵ}}_{ij}^{\prime }({\mathbf{x}}_{ij}^{\prime }({\mathbf{X}}^{\prime }\mathbf{X}/N{)}^{-1}\otimes {\mathbf{I}}_p)(({\mathbf{X}}^{\prime }\mathbf{X}/N{)}^{-1}{\mathbf{x}}_{ij}\otimes {\mathbf{I}}_p){\mathbf{ϵ}}_{ij}.$$

$$E(||{\tilde{\mathbf{\eta }}}_{ij}|{|}^2)$$

$$=tr\left(({\mathbf{x}}_{ij}\otimes {\mathbf{I}}_p){\mathbf{\Sigma }}_i({\mathbf{x}}_{ij}^{\prime }\otimes {\mathbf{I}}_p){\left[\frac{1}{N}({\mathbf{X}}^{\prime }\mathbf{X}\otimes {\mathbf{I}}_p)\right]}^{-1}{\left[\frac{1}{N}({\mathbf{X}}^{\prime }\mathbf{X}\otimes {\mathbf{I}}_p)\right]}^{-1}\right).$$

$$\underset{i,j,N}{\sup }E(||{\tilde{\mathbf{\eta }}}_{ij}|{|}^2)\le C\text{for some }C>0.$$

$$||{\tilde{\mathbf{\eta }}}_{ij}|{|}^2={\mathbf{ϵ}}_{ij}^{\prime }{\mathbf{A}}_{ij}{\mathbf{ϵ}}_{ij}=\sum _{k=1}^p\sum _{\ell =1}^p{a}_{ij,k\ell }{ϵ}_{ijk}{ϵ}_{ij\ell },$$

$$\sum _{i=1}^a\sum _{j=1}^{n_i}E[||{\mathbf{\eta }}_{ij}|{|}^2\text{mathbbm}1\{||{\mathbf{\eta }}_{ij}|{|}^2>\delta \}]$$

$$\le \frac{1}{N}\sum _{i=1}^a\sum _{j=1}^{n_i}\sqrt{E[||{\tilde{\mathbf{\eta }}}_{ij}|{|}^4]}\sqrt{E[\text{mathbbm}1\{||{\tilde{\mathbf{\eta }}}_{ij}|{|}^2>\delta N\}]}$$

$$\le \frac{\tilde{C}}{N}\sum _{i=1}^a\sum _{j=1}^{n_i}\sqrt{E[\text{mathbbm}1\{||{\tilde{\mathbf{\eta }}}_{ij}|{|}^2>\delta N\}]}$$

$$\le \tilde{C}{\left(\frac{C}{\delta N}\right)}^{1/2},$$

$$E[\text{mathbbm}1\{||{\tilde{\mathbf{\eta }}}_{ij}|{|}^2>\delta N\}]$$

$$\frac{1}{N}\sum _{\ell =1}^N{\mathbf{x}}_{\ell }{\mathbf{x}}_{\ell }^{\prime }\otimes (\widehat{{\mathbf{\Sigma }}_{\mathbf{\ell }}}-{\mathbf{\Sigma }}_{\mathbf{\ell }})\stackrel{P}{⟶}0$$

Full text

Keywords: Covariate adjustment, MANCOVA, Small sample, Resampling, Heteroskedasticity-consistent covariance matrix

1 Introduction

Consider a randomized clinical trial, where each subject is randomized to one out of $a$ groups. Frequently, multiple outcomes (e.g., one primary and one co-primary outcome, several secondary outcomes) need to be compared between the groups. One straightforward way of accomplishing this goal would be to conduct univariate tests (e.g., $t$ tests) and adjust for multiplicity. Nevertheless, this approach has the drawback of not taking potential correlations between the outcomes (e.g., due to physiological reasons, there might be correlations between laboratory measurements) into account. Moreover, when analyzing global assessment variables, further safety and efficacy outcomes, or objective variables related to the global assessment, should be analyzed [International Council for Harmonisation of Technical Requirements for Pharmaceuticals for Human Use, 1998]. In these instances, employing a multivariate analysis of variance model would be an attractive option. Still, however, the estimators of the treatment effects might be biased, and the power of the test might be low, if covariates (i.e., variables that are thought to influence the outcomes) are not included in the model. Therefore, a multivariate analysis of covariance (MANCOVA) model might be a better choice. Indeed, MANCOVA models have been used in several recent publications in medical research [Roldan-Valadez et al., 2013, Jackson et al., 2018, Setyowibowo et al., 2018, Tournikioti et al., 2018] and psychology [Memarmoghaddam et al., 2016, Lyndon et al., 2017, Freire et al., 2018, Hyams et al., 2018].

In those papers, the “classical” multivariate tests, which are based on the determinant or the eigenvalues of certain matrices (e.g., Wilks’ Lambda, Roy’s largest root test, see Rencher and Christensen [2002], Anderson [2003]), were used. Analogously to their univariate counterparts, they rely on the assumption that the error vectors are i.i.d. multivariate normal, with expectation ${\bf 0}$ and identical covariance matrix ${\bf\Sigma}$ across groups (i.e., homoskedasticity). However, assessing multivariate normality and homoskedasticity is difficult, in particular when the sample sizes are small. Moreover, even with reasonable sample sizes, there might be good reasons for assuming non-normal errors or heterogeneous covariance matrices. In contrast to the classical multivariate tests discussed in standard textbooks [Rencher and Christensen, 2002, Huitema, 2011], there are only few methods available that are potentially applicable in a more general setting. Recently, Fan and Zhang [2017] have proposed a rank repeated measures analysis of covariance method; however, only limited simulation evidence regarding the performance in small and possibly unbalanced sample size settings is provided. Apart from that, despite the very broad applicability of the approach of Fan and Zhang [2017], which renders this method useful for analyzing ordinal data, too, it cannot be regarded as generalizing the scope of classical mean-based MANCOVA because adjusted mean vectors are not used as effect measures. As an alternative, the sandwich variance estimation technique that was proposed in the context of regression models with clustered errors [Arellano, 1987] and later examined in combination with various bootstrap procedures [Cameron et al., 2008] stays within a semiparametric framework, but allows for heteroskedasticity. Although this considerably extends the scope of applications beyond the classical assumptions that were mentioned before, some requirements might still be crucial in a considerable number of practically relevant cases:

Firstly, the performance in MANCOVA models with possibly small and unequal group sizes combined with heteroskedasticity and/or non-normality has not been sufficiently examined so far. Secondly, the model under consideration in the cluster-robust estimation approach [Cameron et al., 2008] assumes that the regression parameter vector is one and the same across all subjects and clusters. In the context of regression with clustered errors, this is perfectly fine. However, when being translated to the MANCOVA setting, this means that the model does not allow for regression coefficients, which might vary between the respective coordinates of the outcome vector. But, most likely, the association between a particular covariate and the outcomes might not be the same across all components of the outcome. Thirdly, from a theoretical perspective, it should be noted that the existing methods are based on Wald-type statistics and, therefore, require the covariance matrices to be positive definite. Nevertheless, this assumption might be violated in several ways:

Especially in psychology, but also in medical research (e.g., quality of life assessments [Hyams et al., 2018, Setyowibowo et al., 2018]), groups of subjects are often compared with respect to a score. Provided that it is sensible to interpret the score as a metric variable – which requires careful thoughts case-by-case – employing an ANCOVA model would be appropriate. Most likely, a multivariate ANCOVA would be used, in order to allow for between-group comparisons of, for example, the overall score (i.e., a linear combination of the sub-scores) as well as the subscales on which the score is based. However, the covariance matrix would be singular, then. A similar problem might arise in epilepsy research: Consider the situation of examining whether a particular antiepileptic drug is more efficient than placebo with respect to reducing the number of seizures compared to baseline. In this setting, the change from baseline represents the outcome, and the baseline seizure frequency is regarded as a covariate. However, since it might be sensible to consider not only the overall seizure frequency, but also the reduction in the number of focal and generalized seizures, respectively, we would have a trivariate outcome, and again, the covariance matrix would be singular, then. Apart from singularity due to linear dependencies between the components of the outcome vector, computational issues can also lead to singular covariance structures.

In order to overcome the potential problems mentioned above, we propose a modified ANCOVA ANOVA-type statistic (MANCATS), which is inspired by the modified ANOVA-type statistic proposed for heteroskedastic multivariate analysis of variance without covariates [Friedrich and Pauly, 2018]. The main feature of this test statistic is that the “full” covariance matrix estimator, which is used in the Wald-type statistic, is replaced by a diagonal matrix that contains only the group-specific variances on the diagonal. So, the “full” covariance matrix estimator is not used, yet not being removed completely, as in the ANOVA-type approach (see, for example, Brunner et al. [1997]). Consequently, our proposed method shares the advantages of both the Wald-type and the ANOVA-type statistic, since the MANCATS is invariant with respect to the scales of the outcome variables, yet not requiring the assumption of positive definite covariance matrices within the groups. Moreover, it should be noted that, in contrast to repeated measures ANCOVA models, the outcomes are allowed to be measured on different scales. In order to approximate the distribution of the MANCATS and improve its finite-sample performance, we apply two different bootstrap procedures.

The present manuscript is organized as follows: In Section 2, we introduce some notations and set up the model as well as the assumptions that are required to ensure the asymptotic validity of the multivariate Wald-type statistic and the two bootstrap MANCATS approaches, respectively. Actually, it turns out that those assumptions are quite weak, which renders our proposed method applicable to a broad range of practically relevant settings. Then, in Section 3, we describe how the two proposed bootstrap procedures work. The finite-sample performance in terms of type I error rates and power is investigated in an extensive simulation study discussed in Section 4. In order to further substantiate the findings of the previous section and to illustrate the aforementioned practical applicability, simulations that are based on real-life data from standardized student achievement tests are discussed in Section 5. Finally, Section 6 contains some closing remarks and ideas for future research. The proofs of the theorems stated in the main body of the manuscript as well as additional simulation results are included in the supplementary material.

2 The general multivariate analysis of covariance (MANCOVA) model

In the sequel, let ${\bf I}_{m}$ denote the $m$-dimensional identity matrix, and let $\oplus$ and $\otimes$ denote the direct sum and the Kronecker product of matrices, respectively. Let ${\bf Y}_{ij}=(Y_{ij1},Y_{ij2},\dots,Y_{ijp})^{\prime}$ and ${\bf z}_{ij}=(z_{ij}^{(1)},z_{ij}^{(2)},\dots,z_{ij}^{(c)})^{\prime}$ denote the $p$- and $c$-dimensional outcome and covariate vectors of subject $j$ in group $i$, $i\in\{1,2,\dots,a\}$, $j\in\{1,2,\dots,n_{i}\}$. We assume that the covariates are fixed, whereas the outcomes are random variables, satisfying the model equation

[TABLE]

Thereby, $\boldsymbol{\mu}_{i}=(\mu_{i1},\dots,\mu_{ip})^{\prime}$ denotes the vector of adjusted means in group $i$, and $\boldsymbol{\nu}=(\boldsymbol{\nu}^{(1)\prime},\dots,\boldsymbol{\nu}^{(c)\prime})^{\prime}$, where $\boldsymbol{\nu}^{(w)}$ contains the $p$ regression coefficients modelling the association between the $w$-th covariate and the components $1,\dots,p$ of the outcome, $w\in\{1,2,\dots,c\}$. It should be noted that model (1) allows for unequal regression coefficients for different components of the outcome. However, the coefficients are assumed to be the same across groups, which corresponds to the classical assumption of equal regression slopes in the simple univariate ANCOVA setting. Regarding the errors, we assume that the $\boldsymbol{\epsilon}_{ij}$ are independent, with $E[\boldsymbol{\epsilon}_{ij}]={\bf 0}$ and group-specific covariance matrix ${\bf\Sigma}_{i}:=Cov[\boldsymbol{\epsilon}_{ij}]$, $i\in\{1,\dots,a\},j\in\{1,\dots n_{i}\}$.

Model (1) can be expressed in a more compact form by using matrix notation: Let ${\bf Y}=({\bf Y}_{11}^{\prime},{\bf Y}_{12}^{\prime},\dots,{\bf Y}_{an_{a}}^{\prime})^{\prime}$, ${\bf Z}=({\bf z}_{11},{\bf z}_{12},\dots,{\bf z}_{an_{a}})^{\prime}$, $\boldsymbol{\epsilon}=(\boldsymbol{\epsilon}_{11}^{\prime},\boldsymbol{\epsilon}_{12}^{\prime},\dots,\boldsymbol{\epsilon}_{an_{a}}^{\prime})^{\prime}$, and $\boldsymbol{\mu}=(\boldsymbol{\mu}_{1}^{\prime},\dots,\boldsymbol{\mu}_{a}^{\prime})^{\prime}$. Moreover, let ${\bf 1}_{m}$ denote the $m$-dimensional vector containing all $1$’s. Then, model (1) is equivalent to

[TABLE]

where ${\bf\tilde{M}}=\bigoplus_{i=1}^{a}({\bf 1}_{n_{i}}\otimes{\bf I}_{p})$, ${\bf\tilde{Z}}={\bf Z}\otimes{\bf I}_{p}$ and $Cov[\boldsymbol{\epsilon}]=\bigoplus_{i=1}^{a}({\bf I}_{n_{i}}\otimes{\bf\Sigma}_{i})$. So, in particular, $\boldsymbol{\mu}$ can be estimated by the ordinary least squares estimator $\boldsymbol{\hat{\mu}}=(\boldsymbol{\hat{\mu}}_{1}^{\prime},\dots,\boldsymbol{\hat{\mu}}_{a}^{\prime})^{\prime}$, that is,

[TABLE]

where the dot notation indicates averaging over all subjects in the respective group, and $\boldsymbol{\hat{\nu}}^{(w)}$ denotes the OLS estimator of $\boldsymbol{\nu}^{(w)}$, $w=1,\dots,c$. We are interested in testing hypotheses about the adjusted mean vectors $\boldsymbol{\mu}_{1},\dots,\boldsymbol{\mu}_{a}$ of the form $H_{0}:{\bf H}\boldsymbol{\mu}={\bf 0}$, where ${\bf H}$ denotes a contrast matrix (i.e.,${\bf H}{\bf 1}_{ap}={\bf 0}$) of full row rank $(a-1)p$. In the present manuscript, however, we specify the hypothesis matrix in its projector form, that is, we use ${\bf T}:={\bf H}^{\prime}({\bf HH}^{\prime})^{+}{\bf H}$. Thereby, $({\bf HH}^{\prime})^{+}$ denotes the Moore-Penrose inverse of ${\bf HH}^{\prime}$. Note that ${\bf T}$ is uniquely defined, and ${\bf T}\boldsymbol{\mu}={\bf 0}\Leftrightarrow{\bf H}\boldsymbol{\mu}={\bf 0}$. For example, the hypothesis $H_{0}:\boldsymbol{\mu}_{1}=\ldots=\boldsymbol{\mu}_{a}$ is equivalent to $H_{0}:{\bf T}\boldsymbol{\mu}={\bf 0}$, with ${\bf T}={\bf P}_{a}\otimes{\bf I}_{p}$. Thereby, ${\bf P}_{a}={\bf I}_{a}-\frac{1}{a}{\bf J}_{a}$, where ${\bf J}_{a}$ denotes the quadratic $a$-dimensional matrix containing all $1$’s. A more detailed discussion of various multi-factorial and hierarchically nested designs, which are frequently used in practice, can be found in Konietschke et al. [2015].

Now, for testing the hypothesis $H_{0}:{\bf T}\boldsymbol{\mu}={\bf 0}$, one may consider the Wald-type statistic (WTS)

[TABLE]

where ${\bf\hat{\Sigma}}$ denotes the upper-left block of the $2\times 2$ block diagonal matrix $({\bf\tilde{X}}^{\prime}{\bf\tilde{X}})^{-1}{\bf\tilde{X}}^{\prime}{\bf\hat{S}}{\bf\tilde{X}}({\bf\tilde{X}}^{\prime}{\bf\tilde{X}})^{-1}$. Thereby, ${\bf\tilde{X}}=({\bf\tilde{M}},{\bf\tilde{Z}})$, and ${\bf\hat{S}}:=\bigoplus_{i=1}^{a}\bigoplus_{j=1}^{n_{i}}{\bf\hat{\Sigma}}_{ij}$ (i.e., ${\bf\hat{S}}$ is a block-diagonal matrix with matrices ${\bf\hat{\Sigma}}_{11}$,${\bf\hat{\Sigma}}_{12}$,$\dots,{\bf\hat{\Sigma}}_{an_{a}}$ on the diagonal), and ${\bf\hat{\Sigma}}_{ij}={\bf u}_{ij}{\bf u}_{ij}^{\prime}$, ${\bf u}_{ij}={\bf Y}_{ij}-\boldsymbol{\hat{\mu}}_{i}-({\bf z}_{ij}^{\prime}\otimes{\bf I}_{p})\boldsymbol{\hat{\nu}}$, $i\in\{1,2,\dots,a\}$, $j\in\{1,2,\dots,n_{i}\}$. Observe that this covariance matrix estimator basically represents the multivariate generalization of the approach proposed in White [1980] for univariate regression models. Since the publication of this seminal work, several re-scaled estimators have been proposed, with the aim of improving the performance in moderate and small samples [MacKinnon and White, 1985, Cribari-Neto, 2004]. Regarding the theorems stated in the present manuscript, modifications of this type are covered as well, because the scaling factors converge uniformly to $1$. For example, in our simulation studies discussed in Section 4, we propose a multivariate generalization of the so-called HC4 version of ${\bf\hat{\Sigma}}_{ij}$, by multiplying the latter with $1/(1-p_{ij,ij})^{\delta_{ij}}$, $\delta_{ij}:=\min\left(4,p_{ij,ij}/(N^{-1}\sum_{r=1}^{a}\sum_{s=1}^{n_{r}}p_{rs,rs})\right)$, where $p_{rs,rs}$ denotes the diagonal element of the hat matrix ${\bf X}({\bf X}^{\prime}{\bf X})^{-1}{\bf X}^{\prime}$ corresponding to individual $s$ in group $r$, $r\in\{1,2,\dots,a\},s\in\{1,2,\dots,n_{r}\}$. Thereby, ${\bf X}=({\bf M},{\bf Z})=(\bigoplus_{i=1}^{a}{\bf 1}_{n_{i}},{\bf Z})$. This adjustment can be interpreted as a natural generalization of the univariate HC4 estimator, which was proposed by Cribari-Neto [Cribari-Neto, 2004] because for $p=1$, the covariance matrix estimator ${\bf\hat{\Sigma}}_{ij}$ actually reduces to one single value (i.e., the estimator of the subject-specific variance in the univariate case), and the scaling factor defined above does not depend on the dimension of the outcome.

In order to derive the asymptotic distribution of ${\bf W}({\bf T})$ under $H_{0}$, the following assumptions are required, where all convergences are understood for $N\to\infty$:

1. (M1)

The errors $\boldsymbol{\epsilon}_{ij}$ are independent, with $E[\boldsymbol{\epsilon}_{ij}]={\bf 0}$, $Cov[\boldsymbol{\epsilon}_{ij}]={\bf\Sigma}_{i}$, and $E(\epsilon_{ijk}^{4})\leq C_{1}<\infty$, uniformly for all $i\in\{1,2,\dots,a\},j\in\{1,2,\dots,n_{i}\},k\in\{1,2,\dots,p\}$. 2. (M2)

${\bf\Sigma}_{i}>0$ for all $i\in\{1,2,\dots,a\}$. 3. (M3)

For all $i\in\{1,2,\dots,a\}:\frac{N}{n_{i}}\rightarrow\kappa_{i}>0$. 4. (M4)

The columns of ${\bf Z}$ are linearly independent of each other and of the columns of $\bigoplus_{i=1}^{a}{\bf 1}_{n_{i}}$. 5. (M5)

For all $i\in\{1,2,\dots,a\},w\in\{1,\dots,c\}:\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}z_{ij}^{(w)}\longrightarrow c_{i}^{(w)}\in\mathbb{R}$. 6. (M6)

For all $i\in\{1,2,\dots,a\}:\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}{\bf z}_{ij}{\bf z}_{ij}^{\prime}\longrightarrow{\bf\Xi}_{i}\in\mathbb{R}^{c\times c}$.

The following theorem establishes the basic result concerning the asymptotic distribution of the Wald-type statistic $W({\bf T})$:

Theorem 1**.**

Assuming that model (2) as well as (M1) – (M6) hold, the test statistic $W({\bf T})$ defined in (3) has an asymptotic ($N\to\infty$) $\chi_{f}^{2}$ distribution under $H_{0}:{\bf T}\boldsymbol{\mu}={\bf 0}$, with $f=rank({\bf T})$.

The assumptions (M1) and (M3) represent standard requirements in an asymptotic framework. Likewise, a violation of (M4) would mean that collinearity was present; in such a case, the validity of the results would be questionable anyway. Conditions (M5) and (M6) are most likely met in virtually any practical application, too, because an unstable average and/or varying spread of the covariates “on the long run” (i.e., as more and more subjects are enrolled) would indicate a serious flaw in the design or the conduct of the study. However, (M2) might be violated in a considerable number of practically relevant settings, as outlined in Section 1. This restricts the scope of potential applications of the asymptotic Wald-type test statistic, since it requires positive definiteness. Therefore, we propose a different approach, allowing for possibly singular covariance matrices: Analogously to the modified ANOVA-type statistic (MATS) [Friedrich and Pauly, 2018], we consider a modified ANCOVA ANOVA-type statistic (MANCATS), which is defined as

[TABLE]

Thereby,

[TABLE]

where

[TABLE]

and $u_{ijk}=Y_{ijk}-\hat{\mu}_{ik}-\sum_{w=1}^{c}z_{ij}^{(w)}\hat{\nu}_{k}^{(w)}$ is the residual of outcome $k\in\{1,2,\ldots,p\}$ in group $i\in\{1,2,\ldots,a\}$. It should be noted that if no covariates were included in the model, we would have $c=0$, resulting in $u_{ijk}=Y_{ijk}-\bar{Y}_{i.k}$. Consequently, $\hat{\sigma}_{ik}^{2}$ would be equal to the empirical variance estimator, which was used in Friedrich and Pauly [2018]. So, (5) is a natural generalization of the MATS approach to the MANCOVA setting. As already mentioned before, when using $A({\bf T})$ instead of $W({\bf T})$, assumption (M2) can be replaced by the weaker requirement

1. (M2)∗

${\bf\Sigma}_{i}\geq 0$ and $\sigma_{ik}^{2}>0$, for all $i\in\{1,2,\dots,a\}$ and $k\in\{1,2,\dots,p\}$.

This assumption is supposed to be met in a very broad range of practically relevant settings, excluding only cases where, for example, at least one component of the outcome vector is a discrete variable with very few distinct values.

Since the “full” covariance matrix estimator is replaced by $\hat{{\bf D}}$ in the MANCATS, the asymptotic chi-squared limit distribution from Theorem 1 will not hold any more. Nevertheless, there is still at least a formal result regarding the asymptotic distribution of $A({\bf T})$.

Theorem 2**.**

Under conditions (M1), (M(M2)∗)∗, (M3)–(M6), and under $H_{0}:{\bf T}\boldsymbol{\mu}={\bf 0}$, the MANCATS test statistic $A({\bf T})$ defined in (4) has, asymptotically ($N\to\infty$), the same distribution as the weighted sum

[TABLE]

where $U_{ik}\overset{i.i.d.}{\sim}\chi_{1}^{2}$, and the weights $\lambda_{ik}$ are the eigenvalues of ${\bf T}({\bf TDT})^{+}{\bf T}{\bf\Lambda}_{11}$. Thereby, ${\bf D}=\bigoplus_{i=1}^{a}\bigoplus_{k=1}^{p}\kappa_{i}\sigma_{ik}^{2}$, and ${\bf\Lambda}_{11}$ denotes the upper-left block of the $2\times 2$ block matrix ${\bf\Lambda}=({\bf\Xi}^{-1}\otimes{\bf I}_{p}){\bf\Psi}({\bf\Xi}^{-1}\otimes{\bf I}_{p})=\lim_{N\to\infty}[(N^{-1}{\bf\tilde{X}}^{\prime}{\bf\tilde{X}})^{-1}(N^{-1}{\bf\tilde{X}}^{\prime}Cov(\boldsymbol{\epsilon}){\bf\tilde{X}})(N^{-1}{\bf\tilde{X}}^{\prime}{\bf\tilde{X}})^{-1}]$.

However, the weights $\lambda_{ik}$ in Theorem 2 cannot be calculated, because they represent the eigenvalues of a matrix that contains unknown quantities. Therefore, we will discuss two bootstrap-based approximations in the next section.

3 Bootstrapping the MANCATS

Firstly, we consider a so-called wild (or multiplier) bootstrap approach, which has been developed and put forward by the work of Wu [Wu, 1986], Liu [Liu, 1988] and Mammen [Mammen, 1993]. Wild bootstrap techniques have already been applied to Wald-type test statistics that are based on White-type covariance matrix estimation techniques in heteroskedastic univariate settings [Cribari-Neto, 2004, Zimmermann et al., 2019], in regression models with clustered errors [Cameron et al., 2008] and in multivariate factorial designs [Friedrich et al., 2017, Friedrich and Pauly, 2018], in order to improve the small-sample performance of the respective tests. In the present work, we use the following procedure: We generate $N=\sum_{i=1}^{a}n_{i}$ i.i.d. random variables $T_{ij}$ independently from the data, with $E(T_{11})=0$, $Var(T_{11})=1$, and $\sup_{i,j}E(T_{ij}^{4})<\infty$, $1\leq i\leq a,1\leq j\leq n_{i}$. Then, the wild bootstrap observations are defined as

[TABLE]

where $p_{ij,ij}:={\bf x}_{ij}^{\prime}({\bf X}^{\prime}{\bf X})^{-1}{\bf x}_{ij}$ denotes the diagonal element of the hat matrix corresponding to individual $j$ in group $i$, and ${\bf u}_{ij}$ is the $p$-dimensional residual vector of individual $j$ in group $i$, $i\in\{1,\dots,a\}$, $j\in\{1,\dots,n_{i}\}$. The scaling factor $(1-p_{ij,ij})^{-1/2}$ has been introduced by Wu [1986] in case of $p=1$. Observe that we only generate one single random variable $T_{ij}$ per subject, because coordinate-wise bootstrapping would potentially destroy the dependence structure within subjects. Moreover, it should be mentioned that our proof of Theorem 3 works for any particular choice of $T_{ij}$, as long as the fundamental moment conditions mentioned above are met.

Once the bootstrap observations have been generated, the bootstrap least squares estimator $\hat{\boldsymbol{\beta}}^{\ast}:=({\bf\tilde{X}}^{\prime}{\bf\tilde{X}})^{-1}{\bf\tilde{X}}^{\prime}{\bf Y}^{\ast}$, the bootstrap covariance matrix estimator $\hat{{\bf D}}^{\ast}:=diag(n_{1}^{-1}\hat{\sigma}_{11}^{\ast 2},\dots,n_{a}^{-1}\hat{\sigma}_{ap}^{\ast 2})$ and, by plugging in the bootstrap versions instead of the original estimators in (4), the bootstrap analogon $A^{\ast}({\bf T})$ of the MANCATS test statistic $A({\bf T})$ are calculated. The following theorem states that using the conditional $(1-\alpha)$-quantile of the empirical distribution of the wild bootstrap test statistic $A^{\ast}({\bf T})$ as the critical value indeed yields an asymptotically valid test.

Theorem 3**.**

For any parameter vector $\boldsymbol{\beta}:=(\boldsymbol{\mu}^{\prime},\boldsymbol{\nu}^{\prime})^{\prime}\in\mathbb{R}^{(a+c)p}$ and any $\boldsymbol{\beta}_{0}=(\boldsymbol{\mu}_{0}^{\prime},\boldsymbol{\nu}_{0}^{\prime})^{\prime}\in\mathbb{R}^{(a+c)p}$ with ${\bf T}\boldsymbol{\mu}_{0}={\bf 0}$, we have that given the data, as $N\to\infty$,

[TABLE]

holds, provided that (M1), (M(M2)∗)∗ and (M3)–(M6) are fulfilled. Thereby, $\overset{P}{\longrightarrow}$ denotes convergence in probability.

Secondly, we propose a parametric bootstrap approach, which follows an idea that is similar to existing methods for multivariate analysis of variance models [Konietschke et al., 2015, Friedrich and Pauly, 2018]. Given the observed data, we impose a group-wise i.i.d. structure in the bootstrap world by drawing

[TABLE]

where

[TABLE]

Once the bootstrap observations have been generated, the parametric bootstrap version $A^{\star}({\bf T})$ of the MANCATS is obtained analogously to the procedure for the wild bootstrap that was outlined above. Again, the bootstrap test is asymptotically valid:

Theorem 4**.**

For any parameter vector $\boldsymbol{\beta}:=(\boldsymbol{\mu}^{\prime},\boldsymbol{\nu}^{\prime})^{\prime}\in\mathbb{R}^{(a+c)p}$ and any $\boldsymbol{\beta}_{0}=(\boldsymbol{\mu}_{0}^{\prime},\boldsymbol{\nu}_{0}^{\prime})^{\prime}\in\mathbb{R}^{(a+c)p}$ with ${\bf T}\boldsymbol{\mu}_{0}={\bf 0}$, we have that given the data, as $N\to\infty$,

[TABLE]

holds, provided that (M1), (M(M2)∗)∗ and (M3)–(M6) are fulfilled.

4 Simulations

In order to investigate the finite-sample performance of the three methods, which have been discussed in Sections 2 and 3, we have conducted simulations for a broad range of settings, using R version 3.5.1 [R Development Core Team, 2018]. In addition to the three test statistics under consideration, we also included the Wilks’ Lambda test for comparison. We chose Rademacher variables as wild bootstrap weights, which means that independently from the data, we generated i.i.d. random variables $T_{11},\ldots,T_{an_{a}}$ with $P(T_{11}=-1)=P(T_{11}=1)=1/2$. These weights have already been applied successfully with respect to the performance of the corresponding tests in the univariate ANCOVA setting [Zimmermann et al., 2019]. For each scenario, the data generation process was repeated $n_{sim}=10000$ times, and within each simulation run, $n_{boot}=5000$ bootstrap iterations were performed. We considered data from $a=2$ groups, with a bivariate outcome (i.e., $p=2$) and $c=2$ fixed covariates for each subject. More precisely, the values of the first covariate were equally spaced between $-10$ and $10$, whereas the first and second half of the components of the second covariate vector were equally spaced in $[0,5]$ and $[-2,-1]$, respectively, sorted in descending order. In order to generate the error vectors $\boldsymbol{\epsilon}_{ij}$, at first we drew a random sample of $Np$ independent observations $\xi_{111},\xi_{112},\ldots,\xi_{an_{a}p}$ from one out of several distributions (normal, $\chi_{5}^{2}$, lognormal, double exponential). Then, we transformed them to standardized random variables $\tilde{\xi}_{ijk}:=(\xi_{ijk}-E[\xi_{ijk}])/(Var[\xi_{ijk}])^{1/2}$ and subsequently calculated the products

[TABLE]

where ${\bf\Sigma}_{i}^{1/2}$ denotes the matrix square root of the covariance matrix of group $i$, and $\boldsymbol{\tilde{\xi}}_{ij}=(\tilde{\xi}_{ij1},\dots,\tilde{\xi}_{ijp})^{\prime}$, $i\in\{1,2,\dots,a\},j\in\{1,2,\dots,n_{i}\}$. Doing so, it is ensured that $E[\boldsymbol{\tilde{\xi}}_{ij}]={\bf 0}$ and $Cov[\boldsymbol{\tilde{\xi}}_{ij}]={\bf\Sigma}_{i}$. For each of the aforementioned error distributions, we considered three different choices of ${\bf\Sigma}_{1}$ and ${\bf\Sigma}_{2}$, namely ${\bf\Sigma}_{1}={\bf\Sigma}_{2}={\bf I}_{2}+0.5({\bf J}_{2}-{\bf I}_{2})$ (I), ${\bf\Sigma}_{i}=i\cdot{\bf I}_{2}+0.5({\bf J}_{2}-{\bf I}_{2})$ (II), and ${\bf\Sigma}_{1}={\bf\Sigma}_{2}=diag(1,0.25)+0.5({\bf J}_{2}-{\bf I}_{2})$ (III). Observe that the latter covariance matrix is singular. Finally, the observations were obtained as ${\bf Y}_{ij}=\boldsymbol{\mu}_{i}+({\bf z}_{ij}^{\prime}\otimes{\bf I}_{p})\boldsymbol{\nu}+\boldsymbol{\epsilon}_{ij}$, where $\boldsymbol{\nu}=(-0.5,-1.0,1.5,1.0)^{\prime}$. In order to check whether these specifications are sensible, for some scenarios considered in this section, we exemplarily investigated the significance of the conditional associations between the covariates and each of the two outcomes in univariate linear regression models. Exactly speaking, we simulated data for two groups of size $n=40$ each, assuming lognormal or normal errors, and covariance matrix scenarios I and II. In each of the 4 settings, all $10000$ simulation runs yielded univariate conditional associations between the components of the outcome and the covariates that were significant at the 5 percent level. Hence, when comparing the means between the groups, an adjustment for covariates is definitely warranted. For the singular scenario III, we set $\boldsymbol{\nu}=(-0.5,-1.0,1.5,3.0)^{\prime}$, in order to reflect the positive linear dependence between the two components of the outcome vector (observe that $(-1.0,3.0)=2(-0.5,1.5)$). For each of the 12 combinations of distribution and covariance-structure, 4 different sample size scenarios $(n_{1},n_{2})\in\{(20,20),(10,10),(10,20),(20,10)\}$ were simulated. As mentioned before, for estimating the respective variances and covariances, we multiplied the residual vector of subject $j$ in group $i$ with $(1-p_{ij,ij})^{-\delta_{ij}/2}$. Thereby, $\delta_{ij}:=\min\left(4,p_{ij,ij}/(N^{-1}\sum_{r=1}^{a}\sum_{s=1}^{n_{r}}p_{rs,rs})\right)$, where $p_{rs,rs}$ denotes the diagonal element of the hat matrix corresponding to individual $s$ in group $r$, $r\in\{1,2,\dots,a\},s\in\{1,2,\dots,n_{r}\}$. The adjustment was carried out for all variance and covariance estimators that were used in the three test statistics under consideration.

At first, we investigated the empirical type I error rates. Without loss of generality, we set $\boldsymbol{\mu}_{i}={\bf 0}$, $i\in\{1,2\}$. The results are displayed in Table LABEL:Table:TypeIbasic.

The Wilks’ Lambda test showed its well-known conservative behavior in case of positive pairing (i.e., the smallest group has the “smallest” variance) and, on the other hand, exceeded the target 5 percent level if negative pairing (i.e., the smallest group has the “largest” variance) was present. Although the asymptotic Wald-type test might be considered as a remedy, it tended to be either liberal (covariance settings I and II) or conservative (covariance setting III), even under standard assumptions (i.e., normality and homoskedasticity). As indicated by the results that are displayed in Table S3 in the supplementary material, quite large sample sizes are required in order to achieve accurate type I error control for the Wald-type test. By contrast, the wild bootstrap version of the MANCATS showed a good performance also in small samples, yet being more liberal than the parametric bootstrap MANCATS, which maintained the target level very well in all cases. For neither the asymptotic nor the bootstrap approaches, positive and negative pairing (i.e., setting II with $(n_{1},n_{2})=(10,20)$ and $(n_{1},n_{2})=(20,10)$, respectively) affected the results substantially, which could be expected due to the fact that the White estimator and its modified versions have been intended to be “heteroskedasticity-consistent” estimators [White, 1980, MacKinnon, 2012]. However, the results might be dependent on the error distribution, since we see that all approaches yielded slightly smaller type I error rates in the lognormal case compared to the other scenarios.

Secondly, we simulated the empirical power of the three tests under consideration by staying with the setup that was described at the beginning of this section, but specifying $\boldsymbol{\mu}_{1}=(0,0)^{\prime}$ and $\boldsymbol{\mu}_{2}=(\delta,0)^{\prime}$, with $\delta\in\{0.5,1.0,\dots,3.0\}$. For the singular setting III, we set $\boldsymbol{\mu}=(\delta,\delta)^{\prime}$, in order to reflect the positive linear dependence between the components of the outcome vector. Although it was difficult to find scenarios where all approaches provided sufficient type I error level control, we tried to cover different distributions and covariance structures. The results are displayed in Figure 1. Additionally, in order to account for differences in the type I error rates, we provide a plot of the “achieved power” (i.e., the difference between empirical power and the scenario-specific type I error rate) in the supplementary material (Figure S1). The two bootstrap MANCATS approaches showed a similar performance, yet having lower power than the Wald-type test. It should be noted, however, that in case of an underlying normal distribution, the Wald-type approach tended to be liberal (Table LABEL:Table:TypeIbasic). Therefore, the power comparison must be interpreted with caution. Still, in the lognormal case, there was a remarkable difference especially for the homoskedastic setting. However, for the two lognormal scenarios and commonly used target power values around $80$ and $90$ percent (i.e., $\delta\in\{2,2.5\}$), the average power loss of the bootstrap methods compared to the asymptotic Wald-type approach was about $9$ percentage points. The average decrease in achieved power was even lower (about $5$ percentages points). So, summing up, we admit that some caution is needed when applying our proposed bootstrap methods in cases with potentially skewed errors; nevertheless, the price to pay in terms of power might still be acceptable, given the generality of the MANCATS approach. In case of singular covariance matrices, the Wald-type approach should not be used at all, because it had considerably lower power (Figure 2). In these scenarios, the bootstrap-based tests showed a similar performance, with the wild bootstrap approach being somewhat more powerful, which might be due to its slightly liberal behavior (see Figure S2 in the supplementary material).

Furthermore, we conducted several additional simulations, with the aim of investigating whether or not our main findings change substantially under an even broader range of settings. Firstly, all simulation parameters were specified as described above, except from increasing the number of groups to $a=4$, with group size settings ${\bf n}_{1}=(30,30,30,30)$, ${\bf n}_{2}=(15,15,15,15)$, ${\bf n}_{3}=(5,10,20,25)$, and ${\bf n}_{4}=(25,20,10,5)$. The results are reported in Table LABEL:Table:TypeI4gr2dim. Compared to the results for the two-group setting, the Wald-type test was even slightly more liberal for covariance matrix scenarios I and II, whereas the wild bootstrap MANCATS yielded type I error rates that were closer to the target 5 percent level. The parametric bootstrap MANCATS again performed well in all scenarios. The most prominent finding was, however, that the conservative behavior of the Wald-type test was even more severe than for $a=2$ groups, with very low empirical type I error rates for some scenarios. By contrast, the MANCATS-based approaches were close to the target 5 percent level.

Secondly, we stayed with the two-group setting, but increased the number of dimensions to $p=4$. We modified the covariance matrices accordingly, so, I: $\Sigma_{i}=I_{4}+0.5(J_{4}-I_{4})$, II: $\Sigma_{i}=i\cdot I_{4}+0.5(J_{4}-I_{4})$, and III: ${\bf\Sigma}_{i}=({\bf s}_{1},{\bf s}_{2},{\bf s}_{3},{\bf s}_{4})$, where ${\bf s}_{1}={\bf s}_{2}={\bf s}_{3}=(1,1,1,1/2)^{\prime}$, ${\bf s}_{4}=(1/2,1/2,1/2,1)^{\prime}$, $i\in\{1,2,3,4\}$. The vectors modelling the associations between a particular covariate and the components of the outcome were set to $\boldsymbol{\nu}^{(1)}=(-0.5,-1.0,-2.0,-0.02)^{\prime}$ and $\boldsymbol{\nu}^{(2)}=(1.5,1.0,1.0,0.2)^{\prime}$, respectively. For the singular setting III, we set $\boldsymbol{\nu}^{(1)}=(-0.5,-1.0,-1.5,-0.02)^{\prime}$ and $\boldsymbol{\nu}^{(2)}=(1.5,3.0,4.5,0.2)^{\prime}$. Apart from that, all other specifications were the same as before. From Table LABEL:Table:TypeI2gr4dim, it is obvious that in particular the liberality of Wilks’ Lambda and the asymptotic Wald-type test was more pronounced than for $p=2$ dimensions. By contrast, the parametric bootstrap MANCATS again yielded type I error rates that were very close to the pre-specified 5 percent level. Moreover, interestingly, the conservative behavior of the asymptotic Wald-type test in case of singular covariance matrices (scenario III) was slightly ameliorated. When further increasing the dimension to $p=8$, the test even became liberal in small and unbalanced group size settings, yet being quite close to the target level for $n_{1}=n_{2}=20$ (see Table S1 in the supplementary material). These results might indicate that for very small sample sizes, the asymptotic Wald-type test is getting more and more liberal for higher dimensions, but performs reasonably well for moderate sample sizes. Anyway, the parametric bootstrap MANCATS showed very good type I error control regardless of the dimension $p$.

Finally, in order to examine whether the results are sensitive to the choice of the covariates, we considered the aforementioned settings for dimension $p=2$ and $a=2$ groups once more, but replaced the values of the second covariate by a random sample from either the standard normal or the standard lognormal distribution, respectively. The results of the simulations for the homoskedastic singular setting III are displayed in Table S2 in the supplementary material. Compared to the corresponding part of Table LABEL:Table:TypeIbasic, the results were similar. The parametric bootstrap MANCATS appeared to perform even slightly better, whereas the asymptotic Wald-type test was somewhat more conservative.

5 Real-life data example

In this section, we consider standardized test data, which has been collected by William D. Rohwer (University of California, Berkeley). The dataset was taken from Timm [2002]. For each of the $N=69$ kindergarten children that were included in the study, the socioeconomic status (low/high), the number of points on three different achievement tests, and the performance on 5 paired associate learning tasks were recorded. For example, it might be of interest to compare the results from the achievement tests between the two socioeconomic groups, while adjusting for the performance on the learning tasks. For ease of illustration, we shall consider a two-group bivariate MANCOVA model with the performance on the Peabody Picture Vocabulary Test (PPVT) and the Student Achievement Test (SAT) as the components of the outcome vector, and the sum of the 5 paired associate learning task results as a covariate. At first, calculating the estimated residual covariance matrices for the high- $(n=32)$ and low-socioeconomic-status groups ($n=37$) yielded

[TABLE]

respectively. Secondly, we considered the following practically relevant setting: In general, when comparing test results between groups, the primary research question is often focused on differences in a particular overall score. Of course, however, the comparison of subscales might also be of interest. Ideally, both the overall score as well as the subscales, which make up the former overall value, should be analyzed jointly (i.e., by merging them into one single vector and subsequently applying a multivariate method), since the correlation would then be accounted for. By contrast, if two separate tests were applied, the correlation would be ignored. It should be noted that such situations might also arise quite frequently when analyzing, for example, global assessment variables in medical research [International Council for Harmonisation of Technical Requirements for Pharmaceuticals for Human Use, 1998]. However, in cases where the total score is a linear combination (e.g., the sum) of the subscores, the group-specific covariance matrices will be singular. We shall consider such a scenario in the sequel, by adding the sum of the two achievement test results as a third component to the outcome vector, while leaving all other specifications unchanged. The corresponding group-specific estimated residual covariance matrices were

[TABLE]

respectively. Observe that heteroskedasticity is present both in the regular and in the singular setting. In order to investigate whether or not the performance of the approaches that are considered in the present manuscript was similar to the results from Section 4, we examined the empirical type I error rates, adopting the same specificiations for the data generation process as in the previous section, but with ${\bf\Sigma}_{1}$ and ${\bf\Sigma}_{2}$ replaced by $\hat{{\bf\Sigma}}_{h}^{(r)}$ and $\hat{{\bf\Sigma}}_{l}^{(r)}$, or $\hat{{\bf\Sigma}}_{h}^{(s)}$ and $\hat{{\bf\Sigma}}_{l}^{(s)}$, respectively. Regarding the group sizes, apart from the original values (i.e., $(n_{1},n_{2})=(32,37)$), we also considered two additional scenarios, namely $(n_{1},n_{2})=(23,46)$ and $(n_{1},n_{2})=(46,23)$ (note that if we were in a randomized clinical trial setting, these group sizes would correspond to 1:2 and 2:1 allocation ratios). These specifications allow for investigating the effects of positive and negative pairing, while leaving the total sample size $N=69$ unchanged. As in Section 4, the target alpha level was set to 5 percent, and the number of simulation and bootstrap runs were $n_{sim}=10000$ and $n_{boot}=5000$, respectively. The results are displayed in Table LABEL:Table:TypeITimm. Obviously, the empirical type I error rates of Wilks’ Lambda deviated from the target level in presence of negative and positive pairing. By contrast, especially the parametric bootstrap MANCATS performed very well. Interestingly, the asymptotic Wald-type test showed a less liberal behavior compared to the settings in Section 4, which might be explained by the fact that the total sample size $N=69$ was larger. In the singular case, however, the Wald-type test was conservative again, whereas the two bootstrap-based approaches were close to the target level. Only in case of positive pairing and lognormal errors, our proposed MANCATS tests tended to yield liberal results. This might be explained by the somewhat suboptimal behavior of the White estimator in terms of bias and variance, which has been observed in the univariate case [Zimmermann et al., 2019]. Summing up, apart from some issues in that particular scenario, the parametric bootstrap MANCATS was close to the target type I error level in all settings and might therefore be considered as the method of choice for conducting MANCOVA analyses of small to moderately large samples.

6 Discussion and conclusions

In the present manuscript, we have proposed an alternative to the MANCOVA Wald-type test statistic. Although the latter already allowed for heteroskedasticity by incorporating a multivariate generalization of the White sandwich estimator, it still required assumptions which are not met in a considerable number of practically relevant settings. In particular, if the covariance matrix is singular due to computational reasons or linear dependencies between the components of the outcome vector, it is not appropriate to use the Wald-type statistic. Therefore, we have proposed and examined a modified ANOVA-type statistic for the MANCOVA setting, where the full covariance matrix estimator is replaced by a diagonal matrix, which contains basically the re-scaled group-wise variance estimators of the outcomes on the diagonal. As a consequence, the positive definiteness of the group covariance matrices is not required any more. We have proposed two different bootstrap versions of the MANCATS, which both yield an asymptotically valid test. We would like to emphasize that the respective proofs do not require any additional assumptions compared to the Wald-type setting. Moreover, the results of extensive simulations show that the MANCATS-based approaches outperform the asymptotic Wald-type test statistic in most scenarios under consideration, for singular as well as non-singular covariance matrix specifications. On top of that, the Wald-type test statistic is clearly inferior in terms of power in singular settings and, therefore, should not be used. So, summing up, both MANCATS tests can be recommended and are preferable compared to the classical asymptotic Wald-type test statistic, except perhaps for the case of lognormal (i.e., skewed) errors, where all three methods have a somewhat suboptimal performance.

There is a wide variety of different bootstrap methods available; the same applies to heteroskedasticity-consistent covariance matrix estimation techniques. Hence, not all combinations of methods from these fields can be considered. Nevertheless, we have carefully thought about some approaches, which are not covered in the present manuscript. Firstly, various variants of White’s sandwich estimator [White, 1980] for uni- and multivariate linear models have been proposed in the literature (see, for example, MacKinnon [2012] for an overview). We decided to use a multivariate generalization of the HC4 estimator [Cribari-Neto, 2004], since its performance turned out to be superior to other variants (HC0, HC2) in the univariate case [Zimmermann et al., 2019]. Secondly, one might wonder why a nonparametric bootstrap procedure was not considered, given that it is straightforward to implement, and it has been used in multivariate factorial designs [Konietschke et al., 2015, Friedrich and Pauly, 2018]. The reason is that in the context of heteroskedastic regression models, the nonparametric bootstrap has suboptimal properties with respect to variance estimation [Wu, 1986]. Moreover, in multivariate factorial designs, the nonparametic bootstrap performed somewhat worse than the respective competitors [Konietschke et al., 2015, Friedrich and Pauly, 2018].

We would like to emphasize that all three methods, which were discussed in our manuscript, are applicable to a broad range of practically relevant settings not only due to the generality of the underlying assumptions (e.g., allowing for group-specific covariance matrices, and not assuming any particular distribution of the errors), but also with respect to the model we have considered. In particular, the associations between the covariates and a certain component of the outcome vector are allowed to vary from coordinate to coordinate. Nevertheless, still, the scope of applications is restricted to metric outcome variables. Some evidence indicates that with respect to type I errors and power, ANCOVA methods might be valid even for ordinal outcomes [Sullivan and D’Agostino, 2003]. However, this cannot be regarded as a satisfactory solution, both from a theoretical point of view, and with respect to the interpretation of the estimates. Since the case of ordinal outcomes is highly relevant in practice, as briefly outlined in Section 1, efforts towards developing methods which allow for MANCOVA-like analyses of multivariate ordinal outcomes should be made in future research.

Supplementary Material

Section 1 contains the proofs of Theorems 1 to 3, which are stated in the manuscript. Some additional simulation results concerning the type I error rates in settings with higher dimension, different choices of the covariates, and large sample sizes are provided in Section 2. Finally, in order to account for differences in type I error rates when comparing the MANCATS and the White approaches with respect to empirical power, we display achieved power values (i.e., empirical power minus the corresponding type I error rate) in some plots in Section 3.

1 Proofs

1.1 Notations and assumptions

Let ${\bf I}_{m}$ denote the $m$-dimensional identity matrix, and let $\oplus$ and $\otimes$ denote the direct sum and the Kronecker product of matrices, respectively. In order to simplify notations, we consider ${\bf Y}=({\bf X}\otimes{\bf I}_{p})\boldsymbol{\beta}+\boldsymbol{\epsilon}$, where $E(\boldsymbol{\epsilon})={\bf 0},{\bf S}:=Cov(\boldsymbol{\epsilon})=\bigoplus_{i=1}^{a}({\bf I_{n_{i}}}\otimes{\bf\Sigma}_{i})$, ${\bf X}=(\bigoplus_{i=1}^{a}{\bf 1}_{n_{i}},{\bf Z})$, $\boldsymbol{\beta}=(\boldsymbol{\mu}^{\prime},\boldsymbol{\nu}^{\prime})^{\prime}$. Moreover, let $\boldsymbol{\hat{\beta}}=(\boldsymbol{\hat{\mu}}^{\prime},\boldsymbol{\hat{\nu}}^{\prime})^{\prime}$ denote the ordinary least squares estimator of $\boldsymbol{\beta}$, ${\bf\Sigma}_{N}=(({\bf X}^{\prime}{\bf X}/N)^{-1}\otimes{\bf I}_{p})(N^{-1}({\bf X}^{\prime}\otimes{\bf I}_{p}){\bf S}({\bf X}\otimes{\bf I}_{p}))(({\bf X}^{\prime}{\bf X}/N)^{-1}\otimes{\bf I}_{p})$, ${\bf\hat{\Sigma}}_{N}=(({\bf X}^{\prime}{\bf X}/N)^{-1}\otimes{\bf I}_{p})(N^{-1}({\bf X}^{\prime}\otimes{\bf I}_{p}){\bf\hat{S}}({\bf X}\otimes{\bf I}_{p}))(({\bf X}^{\prime}{\bf X}/N)^{-1}\otimes{\bf I}_{p})$, ${\bf\hat{S}}=\bigoplus_{i=1}^{a}\bigoplus_{j=1}^{n_{i}}{\bf\hat{\Sigma}}_{ij}$, ${\bf\hat{\Sigma}}_{ij}={\bf u}_{ij}{\bf u}_{ij}^{\prime}$, $N=\sum_{i=1}^{a}n_{i}$. Thereby, ${\bf u}_{ij}$ denotes the $p$-dimensional residual vector corresponding to subject $j$ in group $i$, $1\leq i\leq a$, $1\leq j\leq n_{i}$.

Observe that assumption (M5) implies that the elements of ${\bf X}$ are uniformly bounded (at least for sufficiently large $N$). Moreover, $(\frac{1}{N}{\bf X}^{\prime}{\bf X})$ converges to, say, ${\bf\Xi}$, due to (M3), (M5), and (M6). According to condition (M4), ${\bf X}$ is of full column rank, which in turn implies that $\frac{1}{N}{\bf X}^{\prime}{\bf X}$ is positive definite. Consequently, $(\frac{1}{N}{\bf X}^{\prime}{\bf X})^{-1}$ converges to ${\bf\Xi}^{-1}$. Analogously, we get that $\frac{1}{N}({\bf X}^{\prime}\otimes{\bf I}_{p}){\bf S}({\bf X}\otimes{\bf I}_{p})$ converges to a matrix $\Psi$, which is positive definite if we additionally assume (M2). In the following proofs, we omit all scaling factors that are based on elements of the hat matrix, for sake of notational simplicity. Observe that the derivations are still valid, though, because $|p_{ij}|=|{\bf x}_{ij}^{\prime}({\bf X}^{\prime}{\bf X})^{-1}{\bf x}_{ij}|\leq C/N$.

1.2 Proof of Theorem 1

The proof comprises two main steps:

1. (A1)

$\sqrt{N}(\boldsymbol{\hat{\beta}}-\boldsymbol{\beta})\overset{d}{\longrightarrow}\mathcal{N}({\bf 0},{\bf\Lambda})$, where ${\bf\Lambda}:=({\bf\Xi}^{-1}\otimes{\bf I}_{p}){\bf\Psi}({\bf\Xi}^{-1}\otimes{\bf I}_{p})$, ${\bf\Psi}:=\lim_{N\to\infty}(N^{-1}({\bf X}^{\prime}\otimes{\bf I}_{p}){\bf S}({\bf X}\otimes{\bf I}_{p}))$, ${\bf\Xi}^{-1}:=\lim_{N\to\infty}(N^{-1}{\bf X}^{\prime}{\bf X})^{-1}$. Note that both limits are well-defined (cf. Section 1.1). 2. (A2)

${\bf\hat{\Sigma}}_{N}-{\bf\Sigma}_{N}\overset{P}{\longrightarrow}{\bf 0}$.

*Proof of (A1):

Let ${\bf X}=({\bf x}_{11},{\bf x}_{12},\dots,{\bf x}_{an_{a}})^{\prime}$. Doing some algebra yields

[TABLE]

Hence, we define the triangular array $\boldsymbol{\eta}_{ij}:=\sqrt{N}(({\bf X}^{\prime}{\bf X})^{-1}\otimes{\bf I}_{p})({\bf x}_{ij}\otimes{\bf I}_{p})\boldsymbol{\epsilon}_{ij},1\leq i\leq a$, $1\leq j\leq n_{i}$, and apply the multivariate Lindeberg-Feller theorem. Obviously, $E(\boldsymbol{\eta}_{ij})={\bf 0}$, and $(\boldsymbol{\eta}_{ij})_{i,j}$ are independent, due to (M1). Secondly, we get

[TABLE]

Finally, we have to show that the Lindeberg condition holds. That is,

[TABLE]

Firstly,

[TABLE]

Let $\boldsymbol{\tilde{\eta}}_{ij}:=(({\bf X}^{\prime}{\bf X}/N)^{-1}\otimes{\bf I}_{p})({\bf x}_{ij}\otimes{\bf I}_{p})\boldsymbol{\epsilon}_{ij}$ (consequently, $||\boldsymbol{\eta}_{ij}||^{2}=N^{-1}||\boldsymbol{\tilde{\eta}}_{ij}||^{2}$). Using Lancaster’s theorem, we get

[TABLE]

So, we have to calculate the trace of a product of matrices, which are uniformly bounded element-wise and have dimensions independent of $N$. Therefore,

[TABLE]

Next, let ${\bf A}_{ij}:=({\bf x}_{ij}^{\prime}\otimes{\bf I}_{p})(({\bf X}^{\prime}{\bf X}/N)^{-1}\otimes{\bf I}_{p})(({\bf X}^{\prime}{\bf X}/N)^{-1}\otimes{\bf I}_{p})({\bf x}_{ij}\otimes{\bf I}_{p})$. Then, $||\boldsymbol{\tilde{\eta}}_{ij}||^{2}$ simplifies to $\boldsymbol{\epsilon}_{ij}^{\prime}{\bf A}_{ij}\boldsymbol{\epsilon}_{ij}$, which can be explicitely written as

[TABLE]

where $a_{ij,k\ell}$ denotes the element in the $k$-th row and the $\ell$-th column of ${\bf A}_{ij}$. Consequently, the expectation $E[||\boldsymbol{\tilde{\eta}}_{ij}||^{4}]$ is equal to a sum of $p^{4}$ expectations of products $\epsilon_{ij{k_{1}}}\epsilon_{ij{\ell_{1}}}\epsilon_{ij{k_{2}}}\epsilon_{ij{\ell_{2}}}$, scaled with elements of the matrix ${\bf A_{ij}}$. The latter quantities are uniformly bounded from above, and using the Cauchy-Schwarz inequality iteratively, the expectation of the product of the error terms can be bounded from above by $E[\epsilon_{ij{k_{1}}}^{4}]E[\epsilon_{ij{\ell_{1}}}^{4}]E[\epsilon_{ij{k_{2}}}^{4}]E[\epsilon_{ij{\ell_{2}}}^{4}]$, which is uniformly bounded, due to assumption (M1). So, all in all, the fourth moments $E[||\boldsymbol{\tilde{\eta}}_{ij}||^{4}]$ can be uniformly bounded from above by, say, $\tilde{C}$. Consequently,

[TABLE]

where the last step was due to

[TABLE]

according to (2). This completes the proof of the Lindeberg condition (1) and, thus, the proof of (A1).

*Proof of (A2):

At first, in order to keep the notation simple (otherwise, we would need double-indices for the diagonal elements of the hat matrix, see below) we merge the group and subject indices into one single index $\ell$ and show that

[TABLE]

holds element-wise, where ${\bf\Sigma}_{\ell}$ and ${\bf\hat{\Sigma}}_{\ell}$ denote the covariance matrix and its estimator for subject $\ell$, respectively, $1\leq\ell\leq N$. Note that the MANCOVA setting (i.e., identical covariance matrices ${\bf\Sigma}_{\ell}$ for all subjects within a particular group) is just a special case thereof. Due to the uniform boundedness of ${\bf x}_{\ell}$, it suffices to show that

[TABLE]

holds element-wise. First, recall that ${\bf u}_{\ell}={\bf Y}_{\ell}-({\bf x}_{\ell}^{\prime}\otimes{\bf I}_{p})\boldsymbol{\hat{\beta}}$. Consequently, we have

[TABLE]

Therefore, we get (recalling that $p_{\ell m}={\bf x}_{\ell}^{\prime}({\bf X}^{\prime}{\bf X})^{-1}{\bf x}_{m}$)

[TABLE]

where we have used that we can think of $X$ as a column vector containing ${\bf x}_{i}^{\prime}\otimes{\bf I}_{p}$ in the $i$-th row, and $\boldsymbol{\epsilon}=(\boldsymbol{\epsilon}_{1}^{\prime},\dots,\boldsymbol{\epsilon}_{N}^{\prime})^{\prime}$. Now, since $\boldsymbol{\epsilon}_{i}$ and $\boldsymbol{\epsilon}_{j}$ are independent for $i\neq j$, and $E[\boldsymbol{\epsilon}_{i}]={\bf 0}$, the expressions from above can be further simplified to

[TABLE]

where the latter step was due to $p_{\ell\ell}=\sum_{m=1}^{N}p_{\ell m}^{2}$. Consequently, the $k$-th diagonal element of the difference $E[{\bf u}_{\ell}{\bf u}_{\ell}^{\prime}]-{\bf\Sigma}_{\ell}$ is

[TABLE]

uniformly in $\ell$, due to the uniform boundedness assumption on the variances, and because the uniform boundedness of the covariates and the elements of $(N^{-1}{\bf X}^{\prime}{\bf X})^{-1}$ implies $p_{\ell m}=O(N^{-1})$.

Next, we consider $Var[({\bf u}_{\ell}{\bf u}_{\ell}^{\prime}-{\bf\Sigma}_{\ell})_{k,k}]=Var[({\bf u}_{\ell}{\bf u}_{\ell}^{\prime})_{k,k}]=Var[u_{\ell k}^{2}]$. Since ${\bf u}_{\ell}=\boldsymbol{\epsilon}_{\ell}-\sum_{m=1}^{N}p_{\ell m}\boldsymbol{\epsilon}_{m}$, we have

[TABLE]

Using Bienayme’s equality yields that (7) can be further simplified to

[TABLE]

By using some algebra and the uniform boundedness conditions on $E[\epsilon_{\ell k}^{4}]$ and the covariates, one can show that (8) - (13) are at least uniformly bounded for all $\ell$ and $k$ (if not converging to [math]). Analogously, we can show that for any $\ell_{1}\neq\ell_{2}$, we have

[TABLE]

uniformly in $\ell_{1},\ell_{2}$ and $k$. So, summing up, we get that

[TABLE]

Combining (6) and (14) and using Chebyshev’s inequality, we get that

[TABLE]

for all $k\in\{1,2,\dots,p\}$.

Finally, we prove that an analogous result holds for the off-diagonal elements. Define $\sigma_{\ell(jk)}:=Cov(\epsilon_{\ell j},\epsilon_{\ell k})$ for $\ell\in\{1,2,\ldots,N\}$, $j,k\in\{1,\ldots,p\}$, with $j\neq k$. Firstly, from (5), one finds that

[TABLE]

uniformly in $\ell,j,k$, using the same arguments as before and the Cauchy-Schwarz inequality in the last step. Secondly, applying the Cauchy-Schwarz inequality again yields

[TABLE]

According to the results from above, the latter expression is bounded from above, uniformly in $\ell$, for all $j,k$.

Moreover, one can show by tedious calculations that

[TABLE]

uniformly in $\ell_{1},\ell_{2}$. Therefore, summing up, we get

[TABLE]

So, all in all, using Chebyshev’s inequality, we now have convergence in probability for the off-diagonal elements, too, which completes the proof of (4). As already mentioned, (3) immediately follows, then, since the arguments used in the derivations above are still valid if the respective quantities are multiplied with the corresponding elements of the design matrix, due the uniform boundedness of the latter. This in turn implies that ${\bf\hat{\Sigma}}_{N}$ is weakly consistent for ${\bf\Sigma}_{N}$, because the deterministic matrix $(N^{-1}{\bf X}^{\prime}{\bf X})^{-1}$ converges element-wise and has dimensions independent of $N$. This completes the proof of (A2).

Finally, the statements (A1) and (A2) are combined as follows: According to the assumptions from Section 2 of the manuscript, we have ${\bf\Sigma}_{N}\rightarrow({\bf\Xi}^{-1}\otimes{\bf I}_{p}){\bf\Psi}({\bf\Xi}^{-1}\otimes{\bf I}_{p})$. Thus, (A2) implies that ${\bf\hat{\Sigma}}_{N}\rightarrow({\bf\Xi}^{-1}\otimes{\bf I}_{p}){\bf\Psi}({\bf\Xi}^{-1}\otimes{\bf I}_{p})$. Since ${\bf\hat{\Sigma}}_{N}$ is of full rank for sufficiently large $N$, there is no rank jump in the convergence ${\bf\tilde{T}}{\bf\hat{\Sigma}}_{N}{\bf\tilde{T}}\overset{P}{\longrightarrow}{\bf\tilde{T}}({\bf\Xi}^{-1}\otimes{\bf I}_{p}){\bf\Psi}({\bf\Xi}^{-1}\otimes{\bf I}_{p}){\bf\tilde{T}}$, which implies that the convergence result also holds for the Moore-Penrose inverse. Hence, applying Slutzky’s theorem to (A1) yields the final result, namely that the quadratic form

[TABLE]

has, asymptotically, a central $\chi_{f}^{2}$ distribution, where $f=rank({\bf\tilde{T}})$. Eventually, the statement of Theorem 1 follows by setting ${\bf\tilde{T}}=diag({\bf T},{\bf 0})$ and applying $H_{0}:{\bf\tilde{T}}\boldsymbol{\beta}={\bf T}\boldsymbol{\mu}={\bf 0}$.

∎

1.3 Proof of Theorem 2

Again, the proof comprises two main steps:

1. (B1)

$\sqrt{N}(\boldsymbol{\hat{\beta}}-\boldsymbol{\beta})\overset{d}{\longrightarrow}\mathcal{N}({\bf 0},{\bf\Lambda})$, where ${\bf\Lambda}:=({\bf\Xi}^{-1}\otimes{\bf I}_{p}){\bf\Psi}({\bf\Xi}^{-1}\otimes{\bf I}_{p})$, ${\bf\Psi}:=\lim_{N\to\infty}(N^{-1}({\bf X}^{\prime}\otimes{\bf I}_{p}){\bf S}({\bf X}\otimes{\bf I}_{p}))$, ${\bf\Xi}^{-1}:=\lim_{N\to\infty}(N^{-1}{\bf X}^{\prime}{\bf X})^{-1}$. 2. (B2)

$({\bf T}N{\bf\hat{D}T})^{+}\overset{P}{\longrightarrow}({\bf TDT})^{+}.$

Statement (B1) has been proven in the previous section; the proof of (B2) works analogously to the derivations in the second part of that section, because

[TABLE]

Therefore, this proof is omitted. Finally, the results are combined as follows: Let ${\bf\tilde{T}}=diag({\bf T},{\bf 0})$ and ${\bf\tilde{D}}=diag({\bf D},{\bf I}_{cp})$. Note that actually, the particular choice of the bottom-right block of the $2\times 2$ block matrix ${\bf\tilde{D}}$ does not matter, because this part is merely inserted for a technical reason, namely in order to allow for the following calculations:

[TABLE]

Now, according to the representation theorem and the continuous mapping theorem, statement (B1) implies that the quadratic form

[TABLE]

has, asymptotically, the same distribution as a weighted sum $\sum_{i=1}^{a+c}\sum_{k=1}^{p}\lambda_{ik}U_{ik}$, where $U_{ik}\overset{i.i.d.}{\sim}\chi_{1}^{2}$, and the weights are the eigenvalues of the matrix of the quadratic form times the asymptotic covariance matrix of $\boldsymbol{\hat{\beta}}$, that is, ${\bf\tilde{T}}({\bf\tilde{T}\tilde{D}\tilde{T}})^{+}{\bf\tilde{T}}\cdot{\bf\Lambda}$. However, since ${\bf\tilde{T}}=diag({\bf T},{\bf 0})$ and $({\bf\tilde{T}\tilde{D}\tilde{T}})^{+}=diag(({\bf TDT})^{+},{\bf 0})$, that product is equal to $diag({\bf T}({\bf TDT})^{+}{\bf T}{\bf\Lambda}_{11},{\bf 0})$, where ${\bf\Lambda}_{11}$ denotes the upper-left part of the matrix ${\bf\Lambda}$. Thus, the weighted some from above eventually reduces to

[TABLE]

where the $\lambda_{ik}$ are the eigenvalues of ${\bf T}({\bf TDT})^{+}{\bf T}{\bf\Lambda}_{11}$. Note that per definition, ${\bf\Lambda}_{11}$ is just the asymptotic covariance matrix of the vector of adjusted means. Now, using (B2) and Slutzky’s theorem yields the desired result, namely that under $H_{0}:{\bf T}\boldsymbol{\mu}={\bf 0}$,

[TABLE]

∎

1.4 Proof of Theorem 3

We have to show that given the data, the wild bootstrap version of the MANCATS weakly converges to the null distribution of the original MANCATS in probability, i.e., that given the data,

1. (C1)

$\sqrt{N}\boldsymbol{\hat{\beta}}^{\ast}\overset{d}{\longrightarrow}\mathcal{N}({\bf 0},{\bf\Lambda})$ in probability, for any general linear model with coefficient vector $\boldsymbol{\beta}$,where ${\bf\Lambda}:=({\bf\Xi}^{-1}\otimes{\bf I}_{p}){\bf\Psi}({\bf\Xi}^{-1}\otimes{\bf I}_{p})$. 2. (C2)

$({\bf T}N{\bf\hat{D}^{\ast}T})^{+}\overset{P}{\longrightarrow}({\bf TDT})^{+}$ in probability.

Proof of (C1):

Firstly,

[TABLE]

So, let $\boldsymbol{\zeta}_{ij}:=\sqrt{N}(({\bf X}^{\prime}{\bf X})^{-1}\otimes{\bf I}_{p})({\bf x}_{ij}\otimes{\bf I}_{p}){\bf u}_{ij}T_{ij}.$ Again, we use the Lindeberg-Feller theorem: $E(\boldsymbol{\zeta}_{ij}|{\bf Y})=0$ is obviously fulfilled, and conditionally on the data, the $\boldsymbol{\zeta}_{ij}$’s are independent, due to construction of the $T_{ij}$’s. Next,

[TABLE]

since it has been shown in the proof of Theorem 1 that

[TABLE]

Thirdly, we check the Lindeberg condition: Doing some algebra yields

[TABLE]

Now, analogously to the proof of Theorem 1, we let ${\bf A}_{ij}:=({\bf x}_{ij}^{\prime}({\bf X^{\prime}X}/N)^{-1}({\bf X^{\prime}X}/N)^{-1}{\bf x}_{ij}\otimes{\bf I}_{p})$, and we denote the $(k,l)$-th element of $A_{ij}$ by $a_{ij,k\ell}$. Consequently, we get

[TABLE]

Now, using

[TABLE]

yields

[TABLE]

Recall that the fourth moments of the wild bootstrap variables $T_{ij}$ were assumed to be uniformly bounded (see Section 3 of the manuscript). Furthermore, using analogous arguments as in the proof of Theorem 1, one can show that $E[({\bf u}_{ij}^{\prime}{\bf A}_{ij}{\bf u}_{ij})^{2}]$ is equal to the sum

[TABLE]

Since the fourth moments of the residuals are uniformly bounded, it follows by applying the Cauchy-Schwarz inequality repeatedly that the sum is uniformly bounded, too. Consequently, the expectation of the expression on the right handside of inequality (1.4) converges to [math]. Hence,

[TABLE]

So, from the conditional version of the Markov inequality, it follows that given the data, the Lindeberg condition holds in probability.

Proof of (C2):

Similarly to the proof in Section 1.2, the result follows if one can show that given the data,

1. (a)

$\lim_{n_{i}\to\infty}E[(n_{i}-c-1)^{-1}\sum_{j=1}^{n_{i}}(u_{ijk}^{\ast 2}-u_{ijk}^{2})|{\bf Y}]\longrightarrow 0$ in probability, and 2. (b)

$\lim_{n_{i}\to\infty}Var[(n_{i}-c-1)^{-1}\sum_{j=1}^{n_{i}}(u_{ijk}^{\ast 2}-u_{ijk}^{2})|{\bf Y}]\longrightarrow 0$ in probability,

for all $i\in\{1,2,\dots,a\}$ and $k\in\{1,2,\dots,p\}$. Thereby, $u_{ijk}^{\ast}$ denotes the $k$-th component of the wild bootstrap residual vector ${\bf u}_{ij}^{\ast}:={\bf Y}_{ij}^{\ast}-({\bf x}_{ij}^{\prime}\otimes{\bf I}_{p})\boldsymbol{\hat{\beta}}^{\ast}$ corresponding to subject $j$ in group $i$. The proof works works completely analogously as in the univariate case [Zimmermann et al., 2019] and is, therefore, omitted.

∎

1.5 Proof of Theorem 4

Recall that the bootstrap observations are generated according to

[TABLE]

where

[TABLE]

Therefore, we have that given the data,

[TABLE]

where

[TABLE]

Using similar arguments as in the proof of Theorem 1 as well as in the consistency proof below, we can show that

[TABLE]

which implies

[TABLE]

Consequently, given the data in probability, the asymptotic distribution of the quadratic form

[TABLE]

is the same as the distribution of $U$ as defined in Theorem 2.

What is left to show is that given the data,

[TABLE]

Again, it suffices to show that given the data,

1. (a)

$\lim_{n_{i}\to\infty}E[(n_{i}-c-1)^{-1}\sum_{j=1}^{n_{i}}(u_{ijk}^{\star 2}-u_{ijk}^{2})|{\bf Y}]\longrightarrow 0$ in probability, and 2. (b)

$\lim_{n_{i}\to\infty}Var[(n_{i}-c-1)^{-1}\sum_{j=1}^{n_{i}}(u_{ijk}^{\star 2}-u_{ijk}^{2})|{\bf Y}]\longrightarrow 0$ in probability,

for all $i\in\{1,2,\dots,a\}$ and $k\in\{1,2,\dots,p\}$. Thereby, $u_{ijk}^{\star}$ denotes the $k$-th component of the parametric bootstrap residual vector ${\bf u}_{ij}^{\star}=(u_{ijk}^{\star})_{k=1}^{p}:={\bf Y}_{ij}^{\star}-({\bf x}_{ij}^{\prime}\otimes{\bf I}_{p})\boldsymbol{\hat{\beta}}^{\star}$ corresponding to subject $j$ in group $i$. Firstly, observe that

[TABLE]

where $p_{ij,i_{1}j_{1}}:={\bf x}_{ij}^{\prime}({\bf X}^{\prime}{\bf X})^{-1}{\bf x}_{i_{1}j_{1}}$. Consequently, we have

[TABLE]

Since $E[Y_{ijk}^{\star 2}|{\bf Y}]=\hat{\sigma}_{ik}^{2}$, and the parametric bootstrap vectors are independent, we get

[TABLE]

Due to the fact that the elements of the hat matrix are bounded by $C/N$, uniformly in $i,j$, and $\hat{\sigma}_{ik}^{2}\overset{P}{\longrightarrow}\sigma_{ik}^{2}$, the expression above converges to [math] in probability. So, the proof of Statement (a) is complete.

Next, we rewrite the conditional variance as follows:

[TABLE]

Now, let

[TABLE]

If both $A$ and $B$ converge to [math] in probability, we are done. Firstly, using (19), Bienayme’s equality and the fact that $Y_{ij}$ and $Y_{\ell m}$ are independent for $(i,j)\neq(l,m)$, we get

[TABLE]

Since given the data, $Y_{ijk}^{\star}\overset{ind.}{\sim}\mathcal{N}(0,\hat{\sigma}_{ik}^{2})$ for fixed $k$, we have $Var[Y_{ijk}^{\star}|{\bf Y}]=\hat{\sigma}_{ik}^{2}$ and $E[Y_{ijk}^{\star 4}|{\bf Y}]=3\hat{\sigma}_{ik}^{4}$, for $i\in\{1,2,\ldots,a\}$. Moreover, recall that the elements of the hat matrix are uniformly bounded by $C/N$. Therefore, each term of the sum converges to [math] in probability. For example,

[TABLE]

and consequently,

[TABLE]

since $\frac{n_{i_{2}}}{N}\hat{\sigma}_{i_{2}k}^{2}\overset{P}{\longrightarrow}\kappa_{i_{2}}^{-1}\sigma_{i_{2}k}^{2}$, $\frac{n_{i_{3}}}{N}\hat{\sigma}_{i_{3}k}^{2}\overset{P}{\longrightarrow}\kappa_{i_{3}}^{-1}\sigma_{i_{3}k}^{2}$. The proof for the other summands of $A$ works analogously and is, therefore, omitted. Likewise, by using the bilinearity of the covariance and considering each part of the resulting sum separately, it follows that $B\overset{P}{\longrightarrow}0$. So, summing up, the proof of Statement (b) and, thus, of Theorem 4 is complete.

∎

2 Further simulation results: Type I error rates

3 Further simulation results: Achieved power