MATS: Inference for potentially Singular and Heteroscedastic MANOVA

TL;DR

This paper develops a new inference method for multivariate analysis of variance that handles singular, heteroscedastic, and high-dimensional data without strict distributional assumptions, using resampling techniques.

Contribution

It extends existing MANOVA tests to accommodate singular and heteroscedastic covariance matrices in factorial designs, with inference based on resampling quantiles.

Findings

Method performs well in simulations under various conditions.

Successfully applied to real election data with complex covariance structures.

Provides confidence regions and ellipsoids for multivariate parameters.

Abstract

In many experiments in the life sciences, several endpoints are recorded per subject. The analysis of such multivariate data is usually based on MANOVA models assuming multivariate normality and covariance homogeneity. These assumptions, however, are often not met in practice. Furthermore, test statistics should be invariant under scale transformations of the data, since the endpoints may be measured on different scales. In the context of high-dimensional data, Srivastava and Kubokawa (2013) proposed such a test statistic for a specific one-way model, which, however, relies on the assumption of a common non-singular covariance matrix. We modify and extend this test statistic to factorial MANOVA designs, incorporating general heteroscedastic models. In particular, our only distributional assumption is the existence of the group-wise covariance matrices, which may even be singular. We base inference on quantiles of resampling distributions, and derive confidence regions and ellipsoids based on these quantiles. In a simulation study, we extensively analyze the behavior of these procedures. Finally, the methods are applied to a data set containing information on the 2016 presidential elections in the USA with unequal and singular empirical covariance matrices.

Figures

Tables

Table 1. Table 1: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the WTS with χ 2 superscript 𝜒 2 \chi^{2} -approximation and parametric bootstrap (PBS) and the MATS with wild bootstrap (wild), parametric bootstrap (PBS) and nonparametric bootstrap (NPBS) in the one-way layout for the normal distribution.

$d$	Cov	$𝒏$	WTS (${\chi }^2$)	WTS (PBS)	MATS (wild)	MATS (PBS)	MATS (NPBS)
$d=4$	S1	(10, 10)	15.2	4.5	6.9	5.2	4.4
		(10, 20)	14.5	5.9	6.9	5.1	4.5
		(20, 10)	14	5.6	7.2	5.3	4.8
		(20, 20)	9.5	5.3	6	5.1	5.1
	S2	(10, 10)	15.2	4.5	7	5	4.5
		(10, 20)	14.5	5.8	6.9	5	4.5
		(20, 10)	14	5.6	7.3	5.5	4.9
		(20, 20)	9.5	5.4	6.3	5.2	5
	S3	(10, 10)	18.3	5.5	7.3	4.8	3.6
		(10, 20)	10.9	4.7	6.4	4.8	4.4
		(20, 10)	21.4	6.6	7.8	4.8	3.4
		(20, 20)	11.2	5.7	6.3	5.1	4.6
	S4	(10, 10)	18.3	5.6	7.5	4.8	3.9
		(10, 20)	11	5.2	6.1	4.6	4.3
		(20, 10)	21	6.7	7.9	4.7	3.2
		(20, 20)	10.9	5.7	6.2	5.0	4.7
$d=8$	S1	(10, 10)	38.6	4.7	7.7	5.1	4.3
		(10, 20)	31	6.2	6.9	5	4.2
		(20, 10)	32.1	6.1	6.6	4.6	4
		(20, 20)	17.0	4.9	5.8	4.8	4.8
	S2	(10, 10)	38.6	4.5	7.9	4.3	3.4
		(10, 20)	31	6.3	7.4	4.3	3.6
		(20, 10)	32.1	6.1	7	4.1	3.4
		(20, 20)	17.0	4.7	6.2	4.8	4.5
	S3	(10, 10)	50.1	6.6	7.9	4.2	2.8
		(10, 20)	21.8	4.1	6.3	4.4	4.1
		(20, 10)	55	10.3	8.5	3.6	2.2
		(20, 20)	21.9	5.4	6.1	4.0	3.6
	S4	(10, 10)	48.9	6.3	7.8	3.6	2.4
		(10, 20)	21.9	4.2	6.3	3.8	3.4
		(20, 10)	54.1	10.4	8.4	3.5	2
		(20, 20)	21.8	5.2	6.0	3.9	3.6

Table 2. Table 2: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the WTS with χ 2 superscript 𝜒 2 \chi^{2} -approximation and parametric bootstrap (PBS) and the MATS with wild bootstrap (wild), parametric bootstrap (PBS) and nonparametric bootstrap (NPBS) in the one-way layout for the χ 3 2 subscript superscript 𝜒 2 3 \chi^{2}_{3} -distribution.

$d$	Cov	$𝒏$	WTS (${\chi }^2$)	WTS (PBS)	MATS (wild)	MATS (PBS)	MATS (NPBS)
$d=4$	S1	(10, 10)	15.3	4	7.1	4.8	3.4
		(10, 20)	13.9	5.5	7.3	5.6	4.6
		(20, 10)	14.6	5.7	7.7	5.9	4.6
		(20, 20)	8.9	4.7	6.3	5.5	5
	S2	(10, 10)	15.3	4.1	7.1	4.5	3.2
		(10, 20)	13.9	5.6	7.5	5.5	4.5
		(20, 10)	14.6	5.8	7.7	5.5	4.5
		(20, 20)	8.9	4.7	6.3	5.3	4.7
	S3	(10, 10)	20.6	7.1	9.5	6.1	3.8
		(10, 20)	11.2	4.8	6.9	4.8	3.7
		(20, 10)	26.2	10.9	12.3	8.9	5.6
		(20, 20)	12.8	6.6	7.6	6	4.7
	S4	(10, 10)	21.2	7.2	9.6	6.2	3.8
		(10, 20)	11.1	5	6.9	4.7	3.4
		(20, 10)	26.5	10.7	12.7	8.9	5.6
		(20, 20)	12.9	6.7	7.7	6.2	4.7
$d=8$	S1	(10, 10)	39.3	3.8	7.7	4.9	3.4
		(10, 20)	32.3	5.5	7.6	5.9	4.7
		(20, 10)	33.4	6.3	7.2	5.1	4.2
		(20, 20)	16.9	4.5	5.9	4.9	4.6
	S2	(10, 10)	39.3	3.8	8.1	4.3	2.7
		(10, 20)	32.3	5.5	8.6	5.2	4
		(20, 10)	33.4	6.3	7.6	4.9	3.9
		(20, 20)	16.9	4.5	6.2	4.5	4.0
	S3	(10, 10)	53.1	6.8	10.2	5.5	3.1
		(10, 20)	23.4	4.8	6.8	4.6	3.5
		(20, 10)	59.9	13.9	13.7	8.1	4.6
		(20, 20)	24.8	6.9	7.9	5.7	4.1
	S4	(10, 10)	52.5	6.3	11	5.3	2.6
		(10, 20)	24.3	4.5	7.1	4.1	2.8
		(20, 10)	59	13.6	14.8	8.4	4.5
		(20, 20)	24.3	6.9	7.7	5.7	4.0

Table 3. Table 3: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the WTS with χ 2 superscript 𝜒 2 \chi^{2} -approximation and parametric bootstrap (PBS) and the MATS with wild bootstrap (wild), parametric bootstrap (PBS) and nonparametric bootstrap (NPBS) in the one-way layout with singular covariance matrices for the normal distribution.

$d$	Cov	$𝒏$	WTS (${\chi }^2$)	WTS (PBS)	MATS (wild)	MATS (PBS)	MATS (NPBS)
$d=4$	S5	(10, 10)	3.1	5.1	5.9	4.8	4.4
		(10, 20)	2.6	5.1	5.7	4.9	4.6
		(20, 10)	2.7	5.3	6.6	5.3	4.7
		(20, 20)	1.7	4.8	5.6	5.2	5.0
	S6	(10, 10)	16.7	5.3	7.1	5	4.6
		(10, 20)	12.4	5.1	6	4.7	4.3
		(20, 10)	17.1	6.1	7.1	5.3	4.6
		(20, 20)	10.2	5.5	6.0	5.2	5.1
	S7	(10, 10)	16.6	5.2	7.3	4.7	4
		(10, 20)	12.3	5.8	6.6	4.6	4.2
		(20, 10)	16.3	5.7	6.9	4.5	4
		(20, 20)	9.4	4.8	5.9	4.8	4.8
$d=8$	S5	(10, 10)	2.8	4.5	6.2	5	4.7
		(10, 20)	2.3	4.9	5.5	4.9	4.7
		(20, 10)	2.6	4.4	5.6	4.7	4.3
		(20, 20)	1.5	4.6	5.5	4.9	4.8
	S6	(10, 10)	39.5	4.4	8.2	5	4.2
		(10, 20)	28.8	5.4	6.6	4.6	4.2
		(20, 10)	35.7	7.0	7.5	4.8	4.0
		(20, 20)	17.3	4.7	6.1	4.8	4.5
	S7	(10, 10)	38.8	4.2	7.4	4	2.9
		(10, 20)	27.4	5.2	6.6	3.6	3
		(20, 10)	36.3	6.3	7.4	3.8	3.2
		(20, 20)	17.3	5.1	5.9	4.2	4.0

Table 4. Table 4: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the WTS with χ 2 superscript 𝜒 2 \chi^{2} -approximation and parametric bootstrap (PBS) and the MATS with wild bootstrap (wild), parametric bootstrap (PBS) and nonparametric bootstrap (NPBS) in the one-way layout with singular covariance matrices for the χ 3 2 subscript superscript 𝜒 2 3 \chi^{2}_{3} -distribution.

$d$	Cov	$𝒏$	WTS (${\chi }^2$)	WTS (PBS)	MATS (wild)	MATS (PBS)	MATS (NPBS)
$d=4$	S5	(10, 10)	2.7	4.2	6.7	5.4	4.5
		(10, 20)	1.9	4.8	5.9	4.9	4.5
		(20, 10)	3.5	6.5	7.3	5.9	5.2
		(20, 20)	1.7	4.9	6.2	5.7	5.5
	S6	(10, 10)	19.7	7.1	7.4	5.2	4.1
		(10, 20)	14.6	7.1	6.4	5	4.3
		(20, 10)	20.1	8.5	8.1	6.4	5.4
		(20, 20)	11.4	6.3	6.5	5.7	5.3
	S7	(10, 10)	19.4	7	7.1	4.1	3.1
		(10, 20)	14.5	6.7	6.4	4.2	3.4
		(20, 10)	20.3	8.7	8.3	6.1	4.5
		(20, 20)	11.7	6.4	6.1	5.1	4.5
$d=8$	S5	(10, 10)	2.4	4.7	6.1	5.1	4.6
		(10, 20)	2.6	5.3	6.1	5.3	5
		(20, 10)	3	5.6	6	5.1	4.6
		(20, 20)	1.2	4.5	5.9	5.3	5.1
	S6	(10, 10)	43.1	5.4	8.2	5.1	3.9
		(10, 20)	30.7	6.6	7.3	5.2	4.2
		(20, 10)	39.2	8.7	8.3	5.6	4.5
		(20, 20)	19.3	5.6	6.8	5.1	4.7
	S7	(10, 10)	42.4	5.5	7.5	3.3	1.7
		(10, 20)	31.1	6.3	7.1	4	2.4
		(20, 10)	39.5	9	9.2	5.1	3.4
		(20, 20)	18.7	5.3	5.1	3.2	2.6

Table 5. Table 5: Descriptive statistics of the data example: Reported are the sample sizes and the 7-dimensional mean vectors for each of the 43 states.

State	n	PST045214	SEX255214	RHI125214	RHI225214	RHI325214	RHI425214	RHI525214
AK	29	25404.55	45.73	52.51	1.92	31.89	5.68	0.55
AL	67	72378.76	51.26	68.26	28.66	0.80	0.73	0.11
AR	75	39551.59	50.54	80.41	16.13	0.89	0.78	0.10
AZ	15	448765.60	49.63	78.99	2.33	14.76	1.45	0.20
CA	58	669008.62	49.52	81.59	3.57	3.16	7.50	0.39
CO	64	83685.41	48.05	92.62	1.81	2.06	1.30	0.12
FL	67	296914.88	48.60	80.61	15.00	0.73	1.70	0.10
GA	159	63505.30	50.31	68.19	28.36	0.49	1.30	0.12
IA	99	31385.11	50.15	95.69	1.43	0.46	1.13	0.08
ID	44	37146.91	49.26	94.36	0.58	2.03	0.82	0.16
IL	102	126280.20	49.93	91.67	5.30	0.35	1.23	0.03
IN	92	71704.95	50.20	94.42	2.85	0.36	0.99	0.04
KS	105	27657.34	49.74	93.64	2.11	1.23	0.91	0.08
KY	120	36778.81	50.25	93.87	3.86	0.29	0.57	0.06
LA	64	72651.19	49.95	64.63	32.08	0.84	0.94	0.05
MD	24	249016.96	50.88	73.14	20.58	0.44	3.41	0.10
ME	16	83130.56	50.89	95.67	0.98	0.87	0.89	0.02
MI	83	119396.11	49.67	91.18	4.05	1.69	1.03	0.03
MN	87	62726.13	49.88	92.93	1.59	2.24	1.47	0.06
MO	115	52726.86	50.11	93.08	3.71	0.62	0.74	0.12
MS	82	36513.16	51.01	56.48	41.21	0.68	0.57	0.04
MT	56	18278.20	49.16	88.81	0.41	8.00	0.48	0.05
NC	100	99439.64	50.82	74.15	20.80	1.93	1.26	0.10
ND	53	13952.49	48.69	90.26	0.79	6.79	0.59	0.04
NE	93	20231.22	49.82	95.43	0.93	1.73	0.59	0.07
NJ	21	425627.38	51.06	76.97	13.19	0.56	7.16	0.10
NM	33	63199.15	49.37	86.02	1.83	8.70	1.11	0.14
NV	17	167005.82	47.41	87.34	2.88	4.27	2.26	0.31
NY	62	318487.53	50.22	87.43	6.98	0.70	2.90	0.06
OH	88	131751.85	50.38	92.70	4.33	0.28	0.97	0.02
OK	77	50364.30	49.84	78.41	3.78	11.18	0.86	0.14
OR	36	110284.42	50.04	91.49	0.89	2.48	1.77	0.30
PA	67	190853.87	50.03	91.99	4.86	0.27	1.46	0.03
SC	46	105053.96	50.93	60.75	36.10	0.64	0.93	0.09
SD	66	12926.89	49.48	82.51	0.78	14.02	0.63	0.04
TN	95	68940.55	50.46	89.75	7.50	0.45	0.73	0.05
TX	254	106129.76	49.20	89.24	6.82	1.17	1.16	0.08
UT	29	101479.38	49.12	92.88	0.66	3.30	1.04	0.36
VA	134	62136.49	50.25	75.67	18.84	0.51	2.09	0.07
WA	39	181064.87	49.85	88.87	1.61	3.01	2.79	0.34
WI	72	79966.17	49.64	92.55	1.71	2.96	1.31	0.04
WV	55	33642.29	50.12	95.57	2.38	0.25	0.49	0.01
WY	23	25397.96	49.04	94.04	1.17	2.13	0.83	0.10

Table 6. Table 7: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the parametric bootstrap (PBS) of the WTS and the MATS in the one-way layout for the log-normal distribution.

$d$	Cov	$n$	WTS (PBS)	MATS (PBS)
$d=4$	S1	(10, 10)	2.0	3.5
		(10, 20)	4.4	5.3
		(20, 10)	4.3	5.5
		(20, 20)	3.5	4.0
	S2	(10, 10)	1.9	3.2
		(10, 20)	4.5	4.9
		(20, 10)	4.4	5.4
		(20, 20)	3.4	3.9
	S3	(10, 10)	6.8	6.1
		(10, 20)	4.2	3.6
		(20, 10)	14.9	13.8
		(20, 20)	8.4	6.7
	S4	(10, 10)	7.0	6.2
		(10, 20)	4.2	3.6
		(20, 10)	15.4	13.5
		(20, 20)	8.7	6.8
$d=8$	S1	(10, 10)	2.9	4.1
		(10, 20)	3.9	5.6
		(20, 10)	4.2	5.3
		(20, 20)	3.0	4.0
	S2	(10, 10)	2.8	2.5
		(10, 20)	4.0	4.8
		(20, 10)	4.1	4.1
		(20, 20)	3.0	3.1
	S3	(10, 10)	6.3	6.0
		(10, 20)	3.5	3.6
		(20, 10)	15.0	12.7
		(20, 20)	7.7	6.7
	S4	(10, 10)	7.0	6.0
		(10, 20)	3.6	2.8
		(20, 10)	15.5	13.6
		(20, 20)	8.0	7.0

Table 7. Table 8: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the parametric bootstrap (PBS) of the WTS and the MATS in the one-way layout for the t 3 subscript 𝑡 3 t_{3} distribution.

$d$	Cov	$n$	WTS (PBS)	MATS (PBS)
$d=4$	S1	(10, 10)	3.6	4.4
		(10, 20)	4.8	4.8
		(20, 10)	3.9	4.3
		(20, 20)	4.0	4.5
	S2	(10, 10)	3.7	4.1
		(10, 20)	4.7	4.6
		(20, 10)	4.1	4.2
		(20, 20)	4.0	4.3
	S3	(10, 10)	4.1	3.1
		(10, 20)	4.3	4.2
		(20, 10)	4.9	3.3
		(20, 20)	3.9	3.8
	S4	(10, 10)	4.2	3.1
		(10, 20)	4.4	3.9
		(20, 10)	4.9	3.2
		(20, 20)	4.0	4.0
$d=8$	S1	(10, 10)	3.5	4.7
		(10, 20)	4.8	4.7
		(20, 10)	5.0	4.2
		(20, 20)	3.9	4.7
	S2	(10, 10)	3.4	3.8
		(10, 20)	4.8	4.1
		(20, 10)	5.0	3.5
		(20, 20)	3.9	3.9
	S3	(10, 10)	5.2	3.2
		(10, 20)	3.5	4.2
		(20, 10)	8.6	2.7
		(20, 20)	4.0	3.6
	S4	(10, 10)	5.0	2.6
		(10, 20)	3.6	3.6
		(20, 10)	8.3	2.5
		(20, 20)	3.8	3.2

Table 8. Table 9: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the parametric bootstrap (PBS) of the WTS and the MATS in the one-way layout for the double-exponential distribution.

$d$	Cov	$n$	WTS (PBS)	MATS (PBS)
$d=4$	S1	(10, 10)	4.1	4.5
		(10, 20)	5.0	4.8
		(20, 10)	4.1	4.2
		(20, 20)	5.1	5.6
	S2	(10, 10)	4.1	4.5
		(10, 20)	4.9	4.8
		(20, 10)	4.0	4.2
		(20, 20)	5.1	5.4
	S3	(10, 10)	4.5	3.3
		(10, 20)	4.3	4.4
		(20, 10)	5.1	3.2
		(20, 20)	4.7	4.9
	S4	(10, 10)	4.6	3.7
		(10, 20)	4.6	4.3
		(20, 10)	5.0	3.3
		(20, 20)	4.7	4.7
$d=8$	S1	(10, 10)	3.4	4.6
		(10, 20)	5.2	4.8
		(20, 10)	4.7	4.9
		(20, 20)	4.3	4.7
	S2	(10, 10)	3.6	3.4
		(10, 20)	5.3	4.3
		(20, 10)	4.6	4.3
		(20, 20)	4.3	4.3
	S3	(10, 10)	4.8	3.0
		(10, 20)	4.1	4.2
		(20, 10)	8.1	3.0
		(20, 20)	5.0	3.6
	S4	(10, 10)	4.7	2.5
		(10, 20)	4.2	3.7
		(20, 10)	7.9	2.5
		(20, 20)	4.9	3.5

Table 9. Table 10: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the WTS with χ 2 superscript 𝜒 2 \chi^{2} -approximation and parametric bootstrap (PBS) and the MATS with wild bootstrap (wild), parametric bootstrap (PBS) and nonparametric bootstrap (NPBS) when testing for an effect of factor A 𝐴 A in a two-way layout with d = 4 𝑑 4 d=4 dimensional observations.

distr	Cov	n	WTS (${\chi }^2$)	WTS (PBS)	MATS (wild)	MATS (PBS)	MATS (NPBS)
normal	S8	${𝒏}^{\left(1\right)}$	10.4	4.7	6.1	4.6	4.4
		${𝒏}^{\left(2\right)}$	9.5	4.5	6.4	5.1	4.9
		${𝒏}^{\left(3\right)}$	11.2	5.4	6.6	5.3	4.8
		${𝒏}^{\left(4\right)}$	6.8	5.0	5.6	5.3	5.2
	S9	${𝒏}^{\left(1\right)}$	10.4	4.8	6.2	4.7	4.4
		${𝒏}^{\left(2\right)}$	9.5	4.4	6.6	5.3	4.9
		${𝒏}^{\left(3\right)}$	11.2	5.5	6.5	4.9	4.6
		${𝒏}^{\left(4\right)}$	6.8	4.9	5.7	5.1	5.1
	S10	${𝒏}^{\left(1\right)}$	9	4.7	5.8	4.5	4.2
		${𝒏}^{\left(2\right)}$	11.1	5.2	6.4	4.8	4.3
		${𝒏}^{\left(3\right)}$	14.8	6.1	7.3	4.9	4
		${𝒏}^{\left(4\right)}$	7.4	5	5.5	4.7	4.7
	S11	${𝒏}^{\left(1\right)}$	9.2	4.7	5.8	4.4	4.1
		${𝒏}^{\left(2\right)}$	10.2	4.9	6	4.4	4.1
		${𝒏}^{\left(3\right)}$	13.3	6	7	4.8	4.2
		${𝒏}^{\left(4\right)}$	7	5.1	5.4	4.6	4.5
${\chi }_3^2$	S8	${𝒏}^{\left(1\right)}$	9.7	4.6	6.7	5.1	4.5
		${𝒏}^{\left(2\right)}$	9.2	3.9	6.7	5.1	4.5
		${𝒏}^{\left(3\right)}$	10.7	4.7	6.7	5	4.4
		${𝒏}^{\left(4\right)}$	6.6	4.3	4.7	4.3	4.1
	S9	${𝒏}^{\left(1\right)}$	9.7	4.6	6.5	4.7	4.1
		${𝒏}^{\left(2\right)}$	9.2	3.9	7	5.1	4.3
		${𝒏}^{\left(3\right)}$	10.7	4.8	6.6	5	4.2
		${𝒏}^{\left(4\right)}$	6.6	4.5	4.9	4.4	4.2
	S10	${𝒏}^{\left(1\right)}$	8.8	4.5	6.1	4.3	3.6
		${𝒏}^{\left(2\right)}$	11.2	5	7.6	5.2	4.1
		${𝒏}^{\left(3\right)}$	16	7	9.4	6.5	5.3
		${𝒏}^{\left(4\right)}$	7.7	5.5	5.6	4.9	4.5
	S11	${𝒏}^{\left(1\right)}$	8.5	4.2	5.6	3.7	3
		${𝒏}^{\left(2\right)}$	9.9	4.3	6.7	4.7	3.8
		${𝒏}^{\left(3\right)}$	13.9	6.4	8.7	5.9	4.4
		${𝒏}^{\left(4\right)}$	7.2	5	5.6	4.9	4.3

Table 10. Table 11: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the WTS with χ 2 superscript 𝜒 2 \chi^{2} -approximation and parametric bootstrap (PBS) and the MATS with wild bootstrap (wild), parametric bootstrap (PBS) and nonparametric bootstrap (NPBS) when testing for an effect of factor B 𝐵 B in a two-way layout with d = 4 𝑑 4 d=4 dimensional observations.

distr	Cov	n	WTS(${\chi }^2$)	WTS (PBS)	MATS (wild)	MATS (PBS)	MATS (NPBS)
normal	S8	${𝒏}^{\left(1\right)}$	10.7	5.1	6.9	5.3	4.8
		${𝒏}^{\left(2\right)}$	9.2	4.6	6	4.8	4.5
		${𝒏}^{\left(3\right)}$	10.2	4.6	6.3	5	4.7
		${𝒏}^{\left(4\right)}$	6.5	4.4	5.0	4.7	4.7
	S9	${𝒏}^{\left(1\right)}$	10.7	5.2	6.6	5.1	4.8
		${𝒏}^{\left(2\right)}$	9.2	4.7	6.1	4.9	4.8
		${𝒏}^{\left(3\right)}$	10.2	4.5	6.2	4.5	4.3
		${𝒏}^{\left(4\right)}$	6.5	4.5	5	4.5	4.4
	S10	${𝒏}^{\left(1\right)}$	8.7	4.9	5.8	4.7	4.4
		${𝒏}^{\left(2\right)}$	10.3	4.6	5.9	4.5	4.1
		${𝒏}^{\left(3\right)}$	13.4	5.3	6.4	4	3.3
		${𝒏}^{\left(4\right)}$	6.9	4.6	5.2	4.7	4.6
	S11	${𝒏}^{\left(1\right)}$	9	4.9	5.7	4.4	4.3
		${𝒏}^{\left(2\right)}$	9.4	4.7	6.1	4.4	4.1
		${𝒏}^{\left(3\right)}$	12.4	5.1	6.4	4.1	3.4
		${𝒏}^{\left(4\right)}$	6.4	4.6	5.2	4.4	4.3
${\chi }_3^2$	S8	${𝒏}^{\left(1\right)}$	10.3	4.9	6.4	4.9	4.3
		${𝒏}^{\left(2\right)}$	9.4	4	6.2	4.8	4.3
		${𝒏}^{\left(3\right)}$	10.7	4.5	5.9	4.7	4.1
		${𝒏}^{\left(4\right)}$	7.2	4.7	5.1	4.7	4.5
	S9	${𝒏}^{\left(1\right)}$	10.3	4.9	6.8	4.8	4.2
		${𝒏}^{\left(2\right)}$	9.4	4	6.3	4.6	4
		${𝒏}^{\left(3\right)}$	10.7	4.5	6	4.4	4
		${𝒏}^{\left(4\right)}$	7.2	4.9	5.4	5	4.8
	S10	${𝒏}^{\left(1\right)}$	9	4.5	6.1	4.4	3.8
		${𝒏}^{\left(2\right)}$	10.6	4.5	6.5	4.4	3.2
		${𝒏}^{\left(3\right)}$	13.7	5.4	7	4.4	3.3
		${𝒏}^{\left(4\right)}$	7.7	5.1	6	5.1	4.6
	S11	${𝒏}^{\left(1\right)}$	9.2	4.6	6.2	4.4	3.7
		${𝒏}^{\left(2\right)}$	9.7	4.2	6.3	4.1	2.9
		${𝒏}^{\left(3\right)}$	12.4	5.1	6.4	4.2	3.1
		${𝒏}^{\left(4\right)}$	7.3	5.2	6	5.1	4.3

Table 11. Table 12: Type-I error rates in % (nominal level α = 5 % 𝛼 percent 5 \alpha=5\% ) for the WTS with χ 2 superscript 𝜒 2 \chi^{2} -approximation and parametric bootstrap (PBS) and the MATS with wild bootstrap (wild), parametric bootstrap (PBS) and nonparametric bootstrap (NPBS) when testing the interaction hypothesis in a two-way layout with d = 4 𝑑 4 d=4 dimensional observations.

distr	Cov	n	WTS(${\chi }^2$)	WTS (PBS)	MATS (wild)	MATS (PBS)	MATS(NPBS)
normal	S8	${𝒏}^{\left(1\right)}$	10.8	4.8	6.9	5.4	5
		${𝒏}^{\left(2\right)}$	10	4.7	6.5	5.4	5.2
		${𝒏}^{\left(3\right)}$	10	4.7	6.8	5.3	5
		${𝒏}^{\left(4\right)}$	6.3	4.3	5.0	4.6	4.6
	S9	${𝒏}^{\left(1\right)}$	10.8	4.9	6.9	5.1	4.9
		${𝒏}^{\left(2\right)}$	10	4.8	6.3	5.1	4.8
		${𝒏}^{\left(3\right)}$	10	4.8	6.4	5	4.8
		${𝒏}^{\left(4\right)}$	6.3	4.3	5	4.7	4.4
	S10	${𝒏}^{\left(1\right)}$	8.8	5.1	6.2	4.8	4.5
		${𝒏}^{\left(2\right)}$	10.9	4.8	6.6	4.9	4.5
		${𝒏}^{\left(3\right)}$	13.8	5.7	7.5	4.9	4.1
		${𝒏}^{\left(4\right)}$	6.9	4.5	5	4.3	4.2
	S11	${𝒏}^{\left(1\right)}$	9.1	4.9	6.4	4.8	4.4
		${𝒏}^{\left(2\right)}$	10.1	4.8	6.6	4.7	4.5
		${𝒏}^{\left(3\right)}$	12.9	5.4	7.3	4.7	4.1
		${𝒏}^{\left(4\right)}$	6.5	4.5	4.9	4.4	4.2
${\chi }_3^2$	S8	${𝒏}^{\left(1\right)}$	9.5	4.3	6.3	4.8	4.3
		${𝒏}^{\left(2\right)}$	9.4	4.3	6.3	4.8	4.4
		${𝒏}^{\left(3\right)}$	9.8	4.1	5.9	4.6	4
		${𝒏}^{\left(4\right)}$	6.6	4.6	4.8	4.4	4.2
	S9	${𝒏}^{\left(1\right)}$	9.5	4.3	6.3	4.6	4.1
		${𝒏}^{\left(2\right)}$	9.4	4.1	6.5	4.8	4.2
		${𝒏}^{\left(3\right)}$	9.8	4.2	5.9	4.6	4
		${𝒏}^{\left(4\right)}$	6.6	4.6	5.1	4.6	4.4
	S10	${𝒏}^{\left(1\right)}$	7.9	4.1	5.9	4.3	3.6
		${𝒏}^{\left(2\right)}$	10.3	4.2	5.9	3.9	2.9
		${𝒏}^{\left(3\right)}$	13	4.3	6.4	3.7	2.4
		${𝒏}^{\left(4\right)}$	6.8	4.7	5.5	4.6	4
	S11	${𝒏}^{\left(1\right)}$	8.2	4	5.8	4	3.4
		${𝒏}^{\left(2\right)}$	9.7	4.2	6	3.9	2.9
		${𝒏}^{\left(3\right)}$	11.7	4.3	6.2	3.5	2.4
		${𝒏}^{\left(4\right)}$	6.7	4.8	5.4	4.4	4.2

Table 12. Table 13: Descriptive statistics of the cardiology data. Volume measurements (EDV, ESV and SV) are in m ​ ℓ 𝑚 ℓ m\ell , peak strain rate measurements are in 1 / s ​ e ​ c . 1 𝑠 𝑒 𝑐 1/sec.

Gender	$n$	Mean					Sd
		EDV	ESV	SV	PSSR	PDSR	EDV	ESV	SV	PSSR	PDSR
female	92	124.37	41.39	82.98	-1.16	1.07	25.25	13.38	15.77	0.32	0.30
male	96	157.02	54.98	102.04	-1.07	0.94	31.02	15.98	18.70	0.39	0.35

Equations

$${\mathbf{X}}_{ik}$$

$$n_i/N\to {\kappa }_i>0,i=1,\dots ,a$$

$${\mathbf{\Sigma }}_N=\mathrm{Cov}(\sqrt{N}{\overline{\mathbf{X}}}_{\bullet })=\mathrm{diag}\left(\frac{N}{n_i}{\mathbf{V}}_i:1\le i\le a\right),$$

$${\hat{\mathbf{V}}}_i=\frac{1}{n_i-1}\sum _{k=1}^{n_i}({\mathbf{X}}_{ik}-{\overline{\mathbf{X}}}_{i\cdot })({\mathbf{X}}_{ik}-{\overline{\mathbf{X}}}_{i\cdot }{)}^{\top }$$

$${\hat{\mathbf{\Sigma }}}_N=\mathrm{diag}\left(\frac{N}{n_i}{\hat{\mathbf{V}}}_i:1\le i\le a\right).$$

$$T_N=N{\overline{\mathbf{X}}}_{\bullet }^{\top }\mathbf{T}(\mathbf{T}{\hat{\mathbf{\Sigma }}}_N\mathbf{T}{)}^{+}\mathbf{T}{\overline{\mathbf{X}}}_{\bullet },$$

$$\mathrm{Pr}(T_N>{\chi }_{4;0.95}^2)=1-\mathrm{Pr}(T_N\le {\chi }_{4;0.95}^2)\approx 0.0087,$$

$${\tilde{Q}}_N=N{\overline{\mathbf{X}}}_{\bullet }^{\top }\mathbf{T}{\overline{\mathbf{X}}}_{\bullet }.$$

$$Q_N=N{\overline{\mathbf{X}}}_{\bullet }^{\top }\mathbf{T}(\mathbf{T}{\hat{\mathbf{D}}}_N\mathbf{T}{)}^{+}\mathbf{T}{\overline{\mathbf{X}}}_{\bullet }.$$

$$Q_N=N{\overline{\mathbf{X}}}_{\bullet }^{\top }\mathbf{T}(\mathbf{T}{\hat{\mathbf{D}}}_N\mathbf{T}{)}^{+}\mathbf{T}{\overline{\mathbf{X}}}_{\bullet }\stackrel{\mathcal{D}}{\to }Z=\sum _{i=1}^a\sum _{s=1}^d{\lambda }_{is}{Z}_{is},$$

$${\mathbf{X}}_{i1}^{\ast },\dots ,{\mathbf{X}}_{i{n}_i}^{\ast }\stackrel{\text{i.i.d.}}{\sim }\mathcal{N}(0,{\hat{\mathbf{V}}}_i),i=1,\dots ,a.$$

$$Q_N^{\ast }=N({\overline{\mathbf{X}}}_{\bullet }^{\ast }{)}^{\top }\mathbf{T}(\mathbf{T}{\hat{\mathbf{D}}}_N^{\ast }\mathbf{T}{)}^{+}\mathbf{T}{\overline{\mathbf{X}}}_{\bullet }^{\ast }.$$

$$\underset{x\in \mathbb{R}}{\sup }|\underset{\mathbf{\mu }}{\mathrm{Pr}}(Q_N^{\ast }\le x|\mathbf{X})-\underset{{\mathbf{\mu }}_0}{\mathrm{Pr}}(Q_N\le x)|\stackrel{\mathrm{Pr}}{\to }0,$$

$${\mathbf{X}}_{ik}^{\star }=W_{ik}({\mathbf{X}}_{ik}-{\overline{\mathbf{X}}}_{i\cdot }),i=1,\dots ,a,k=1,\dots ,n_i.$$

$$\underset{x\in \mathbb{R}}{\sup }|\underset{\mathbf{\mu }}{\mathrm{Pr}}(Q_N^{\star }\le x|\mathbf{X})-\underset{{\mathbf{\mu }}_0}{\mathrm{Pr}}(Q_N\le x)|\stackrel{\mathrm{Pr}}{\to }0.$$

$$\underset{x\in \mathbb{R}}{\sup }|\underset{\mathbf{\mu }}{\mathrm{Pr}}(Q_N^{†}\le x|\mathbf{X})-\underset{{\mathbf{\mu }}_0}{\mathrm{Pr}}(Q_N\le x)|\stackrel{\mathrm{Pr}}{\to }0.$$

$$\text{p-value}=\frac{1}{B}\sum _{b=1}^B\text{mathds}1\{Q_N\le Q_N^{\ast ,b}\}.$$

$$N(\mathbf{H}{\overline{\mathbf{X}}}_{\bullet }-\mathbf{H}\mathbf{\mu }{)}^{\top }(\mathbf{H}{\hat{\mathbf{D}}}_N{\mathbf{H}}^{\top }{)}^{+}(\mathbf{H}{\overline{\mathbf{X}}}_{\bullet }-\mathbf{H}\mathbf{\mu })\le c_{1-\alpha }^{\ast }.$$

$$\mathbf{H}{\overline{\mathbf{X}}}_{\bullet }\pm \sqrt{{\hat{\lambda }}_s\cdot c_{1-\alpha }^{\ast }/N}\cdot {\hat{\mathbf{e}}}_s,s=1,\dots ,d.$$

$$\mathbf{\Omega }=\{H_0:{\mathbf{h}}_{\ell }^{\top }\mathbf{\mu }=0\text{ with }{\mathbf{h}}_{\ell }\ne 0,\ell =1,\dots ,q\}.$$

$$Q_N^{\ell }=N({\mathbf{h}}_{\ell }^{\top }{\overline{\mathbf{X}}}_{\bullet }{)}^{\top }({\mathbf{h}}_{\ell }^{\top }{\hat{\mathbf{D}}}_N{\mathbf{h}}_{\ell }{)}^{-1}({\mathbf{h}}_{\ell }^{\top }{\overline{\mathbf{X}}}_{\bullet })=N\frac{{\left({\sum }_{i=1}^a{\sum }_{s=1}^d{h}_{\ell ,is}{\overline{X}}_{i\cdot s}\right)}^2}{{\sum }_{i=1}^a{\sum }_{s=1}^d{h}_{\ell ,is}^2{\hat{\sigma }}_{is}^2}$$

$$S_N\equiv N(\mathbf{H}{\overline{\mathbf{X}}}_{\bullet }{)}^{\top }\mathrm{diag}\left({\left({\mathbf{h}}_{\ell }^{\top }{\hat{\mathbf{D}}}_N{\mathbf{h}}_{\ell }\right)}^{-1}:\ell =1,\dots ,q\right)\mathbf{H}{\overline{\mathbf{X}}}_{\bullet }=\sum _{\ell =1}^q{Q}_N^{\ell }$$

$$\mathrm{Pr}\left(\bigcap _{\ell =1}^q\{Q_N^{\ell }\le q_{1-\alpha }^{\ast }\}\right)\le \mathrm{Pr}\left(\sum _{\ell =1}^q{Q}_N^{\ell }\le q_{1-\alpha }^{\ast }\right)\to 1-\alpha ,N\to \infty .$$

$${\mathbf{h}}_{\ell }^{\top }{\overline{\mathbf{X}}}_{\bullet }\pm \sqrt{q_{1-\alpha }^{\ast }\cdot {\mathbf{h}}_{\ell }^{\top }{\hat{\mathbf{D}}}_N{\mathbf{h}}_{\ell }/N}.$$

$${\mathbf{V}}_1={\mathbf{I}}_d+0.5({\mathbf{J}}_d-{\mathbf{I}}_d)={\mathbf{V}}_2,$$

$${\mathbf{V}}_1={\left((0.6{)}^{|r-s|}\right)}_{r,s=1}^d={\mathbf{V}}_2,$$

$${\mathbf{V}}_1={\mathbf{I}}_d+0.5({\mathbf{J}}_d-{\mathbf{I}}_d)\text{ and }{\mathbf{V}}_2={\mathbf{I}}_d\cdot 3+0.5({\mathbf{J}}_d-{\mathbf{I}}_d),$$

$${\mathbf{V}}_1={\left((0.6{)}^{|r-s|}\right)}_{r,s=1}^d\text{ and }{\mathbf{V}}_2={\left((0.6{)}^{|r-s|}\right)}_{r,s=1}^d+{\mathbf{I}}_d\cdot 2.$$

$${\mathbf{X}}_{ik}={\mathbf{\mu }}_i+{\mathbf{V}}_i^{1/2}{\mathbf{ϵ}}_{ik},i=1,\dots ,a;k=1,\dots ,n_i,$$

$${ϵ}_{iks}=\frac{Y_{iks}-\mathrm{E}(Y_{iks})}{\sqrt{\mathrm{Var}(Y_{iks})}}$$

$$\displaystyle\boldsymbol{V}_{1}=\left(\begin{array}[]{cccc}1&1/2&1&1\\ 1/2&1&1/2&1/2\\ 1&1/2&1&1\\ 1&1/2&1&1\\ \end{array}\right),\boldsymbol{V}_{2}=\boldsymbol{V}_{1}+0.5\cdot\boldsymbol{J}_{d}$$

$$\displaystyle\boldsymbol{V}_{1}=\left(\begin{array}[]{cccc}1&0.6&0.36&0.18\\ 0.6&1&0.6&0.3\\ 0.36&0.6&1&0.5\\ 0.18&0.3&0.5&0.25\\ \end{array}\right),\boldsymbol{V}_{2}=\boldsymbol{V}_{1}+0.5\cdot\boldsymbol{J}_{d}$$

Full text

1 Motivation and Introduction

In many experiments in, e.g., biology, ecology and psychology several endpoints, potentially measured on different scales, are recorded per subject. As an example, we consider a data set on the 2016 presidential elections in the USA containing demographic data on counties from US census. For our exemplary analysis, we aim to investigate whether the states differ with respect to some demographic factors. In addition to unequal empirical covariance matrices between groups, analysis is further complicated by their singularity.

The analysis of such multivariate data is typically based on classical MANOVA models assuming multivariate normality and/or homogeneity of the covariance matrices, see, e.g., [1, 14, 15, 20, 25, 34, 44]. These assumptions, however, are often not met in practice (as in the motivating example) and it is well known that the methods perform poorly in case of heterogeneous data [23, 40]. Furthermore, the test statistic should be invariant under scale transformations of the components, since the endpoints may be measured on different scales. Thus, test statistics of multivariate ANOVA-type (ATS) as, e.g., proposed in [5] and studied in [16], are only applicable if all endpoints are measured on the same scale, i.e., for repeated measures designs. Assuming non-singular covariance matrices and certain moment assumptions the scale invariance is typically accomplished by utilizing test statistics of Wald-type (WTS). However, inference procedures based thereon require (extremely) large sample sizes for being accurate, see [23, 37, 41]. In particular, even the novel approaches of [23] and [37] showed a more or less liberal behavior in case of skewed distributions. Moreover, their procedures cannot be used to analyze the motivating data example with possibly singular covariance matrices. Therefore, we follow a different approach by modifying the above mentioned ANOVA-type statistic (MATS). It is motivated from the modified Dempster statistic proposed in [39] for high-dimensional one-way MANOVA. This statistic is also invariant under the change of units of measurements. However, until now, it has only been developed for a homoscedastic one-way setting assuming non-singularity and a specific distributional structure that is motivated from multivariate normality.

It is the aim of the present paper to modify and extend the [39] approach to factorial MANOVA designs, incorporating general heteroscedastic models. In particular, we only postulate existence of the group-wise covariance matrices, which may even be singular. The small sample behavior of our test statistic is enhanced by applying bootstrap techniques as in [23]. Thereby, the novel MATS procedure enables us to relax the usual MANOVA assumptions in several ways: While incorporating general heteroscedastic designs and allowing for potentially singular covariance matrices we postulate their existence as solely distributional assumption, i.e., only finite second moments are required. Moreover, we gain a procedure that is more robust against deviations from symmetry and homoscedasticity than the usual WTS approaches.

So far, only few approaches have been investigated which do not assume normality or equal covariance matrices (or both). Examples in the nonparametric framework include the permutation based nonparametric combination method [32, 33] and the rank-based tests presented in [6] and [7] for Split Plot Designs and in [2] and [27] for MANOVA designs. However, these methods are either not applicable for general MANOVA models or based on null hypotheses formulated in terms of distribution functions. In contrast we wish to derive inference procedures (tests and confidence regions) for contrasts of mean vectors. Here, beside all previously mentioned procedures, only methods for specific designs have been developed, see [11] for two-sample problems, [42, 43] for robust but homoscedastic one-way MANOVA or [18] for a particular two-way MANOVA. To our knowledge, mean-based MANOVA procedures allowing for possibly singular covariance matrices have not been developed so far.

The paper is organized as follows: In Section 2 we describe the statistical model and hypotheses. Furthermore, we propose a new test statistic, which is applicable to singular covariance matrices and is invariant under scale transformations of the data. In Section 3, we present three different resampling approaches, which are used for the derivation of statistical tests as well as confidence regions and simultaneous confidence intervals for contrasts in Section 4. The different approaches are compared in a large simulation study (Section 5), where we analyze different factorial designs with a wide variety of distributions and covariance settings. The motivating data example is analyzed in Section 6 and we conclude with some final remarks and discussion in Section 7. All proofs are deferred to the supplementary material, where we also provide further simulation results and the analysis of an additional data example.

2 Statistical Model, Hypotheses and Test Statistics

Throughout the paper, we will use the following notation. We denote by $\boldsymbol{I}_{d}$ the $d$-dimensional unit matrix and by $\boldsymbol{J}_{d}$ the $d\times d$ matrix of 1’s, i.e., $\boldsymbol{J}_{d}=\boldsymbol{1}_{d}\boldsymbol{1}_{d}^{\top}$, where $\boldsymbol{1}_{d}=(1,\dots,1)^{\top}$ is the $d$-dimensional column vector of 1’s. Furthermore, let $\boldsymbol{P}_{d}=\boldsymbol{I}_{d}-d^{-1}\boldsymbol{J}_{d}$ denote the $d$-dimensional centering matrix. By $\oplus$ and $\otimes$ we denote the direct sum and the Kronecker product, respectively.

In order to cover different factorial designs of interest, we establish the general model

[TABLE]

for treatment group $i=1,\dots,a$ and individual $k=1,\dots,n_{i}$, on which we measure $d$-variate observations. Here $\boldsymbol{\mu}_{i}=(\mu_{i1},\dots,\mu_{id})^{\top}\in\mathbb{R}^{d}$ for $i=1,\dots,a$. A factorial structure can be incorporated by splitting up indices, see, e.g., [23]. For fixed $1\leq i\leq a$, the error terms $\boldsymbol{\epsilon}_{ik}$ are assumed to be independent and identically distributed $d$-dimensional random vectors, for which the following conditions hold:

• (1)

${\rm E}(\boldsymbol{\epsilon}_{i1})=0,~{}i=1,\dots,a,$

• (2)

$0<\sigma_{is}^{2}=\operatorname{Var}(X_{iks})<\infty,~{}i=1,\dots,a,~{}s=1,\dots,d,$

• (3)

$\operatorname{Cov}(\boldsymbol{\epsilon}_{i1})=\boldsymbol{V}_{i}\geq 0,~{}i=1,\dots,a.$

Thus, we only assume the existence of second moments. For convenience, we aggregate the individual vectors into $\boldsymbol{X}=(\boldsymbol{X}_{11}^{\top},\dots,\boldsymbol{X}_{an_{a}}^{\top})^{\top}$ as well as $\boldsymbol{\mu}=(\boldsymbol{\mu}_{1}^{\top},\dots,\boldsymbol{\mu}_{a}^{\top})^{\top}$. Denote by $N=\sum_{i=1}^{a}n_{i}$ the total sample size. In order to derive asymptotic results in this framework, we will throughout assume the usual sample size condition:

[TABLE]

as $N\to\infty$.

An estimator for $\boldsymbol{\mu}$ is given by the vector of pooled group means $\boldsymbol{\overline{X}}_{i\cdot}=n_{i}^{-1}\sum_{k=1}^{n_{i}}\boldsymbol{X}_{ik},i=1,\dots,a$, which we denote by $\boldsymbol{\overline{X}}_{\bullet}=(\boldsymbol{\overline{X}}_{1\cdot}^{\top},\dots,\boldsymbol{\overline{X}}_{a\cdot}^{\top})^{\top}$. The covariance matrix of $\sqrt{N}~{}\boldsymbol{\overline{X}}_{\bullet}$ is given by

[TABLE]

where the group-specific covariance matrices $\boldsymbol{V}_{i}$ are estimated by the empirical covariance matrices

[TABLE]

resulting in

[TABLE]

In this semi-parametric framework, hypotheses are formulated in terms of the mean vector as $H_{0}:\boldsymbol{H}\boldsymbol{\mu}={\bf 0}$, where $\boldsymbol{H}$ is a suitable contrast matrix, i.e., $\boldsymbol{H}{\bf 1}_{ad}={\bf 0}$. Note that we can use the unique projection matrix $\boldsymbol{T}=\boldsymbol{H}^{\top}(\boldsymbol{H}\boldsymbol{H}^{\top})^{+}\boldsymbol{H}$, where $(\boldsymbol{H}\boldsymbol{H}^{\top})^{+}$ denotes the Moore-Penrose inverse of $\boldsymbol{H}\boldsymbol{H}^{\top}$. It is $\boldsymbol{T}=\boldsymbol{T}^{2},\boldsymbol{T}=\boldsymbol{T}^{\top}$ and $\boldsymbol{T}\boldsymbol{\mu}={\bf 0}\Leftrightarrow\boldsymbol{H}\boldsymbol{\mu}={\bf 0}$, see, e.g., [8].

A commonly used test statistic for multivariate data is the Wald-type statistic (WTS)

[TABLE]

which requires the additional assumption (3’) $\boldsymbol{V}_{i}>0,~{}i=1,\dots,a$. It is easy to show that the WTS has under $H_{0}:\boldsymbol{T}\boldsymbol{\mu}={\bf 0}$, asymptotically, as $N\to\infty$, a $\chi^{2}_{f}$-distribution with $f=\operatorname{rank}(\boldsymbol{T})$ degrees of freedom if (1) – (3’) holds. However, large sample sizes are necessary to maintain a pre-assigned level $\alpha$ using quantiles of the limiting $\chi^{2}$-distribution. [23] proposed different resampling procedures in order to improve the small sample behavior of the WTS for multivariate data. Therein, a parametric bootstrap approach turned out to be the best in case that the underlying distributions are not too skewed and/or too heteroscedastic. In the latter cases all considered procedures were more or less liberal. Moreover, assuming only (3) instead of (3’) for the WTS would in general not lead to an asymptotic $\chi^{2}_{f}$-limit distribution. The reason for this are possible rank jumps between $\boldsymbol{T}\widehat{\boldsymbol{\Sigma}}_{N}\boldsymbol{T}$, $\boldsymbol{T}\boldsymbol{\Sigma}\boldsymbol{T}$ and $\boldsymbol{T}$. To accept this, suppose that $\operatorname{rank}(\boldsymbol{T}\boldsymbol{\Sigma}\boldsymbol{T})=2$, while $\operatorname{rank}(\boldsymbol{T})=4$ (this corresponds to Scenario S5 in the simulation studies below). If additionally $\lim_{N\to\infty}\operatorname{rank}(\boldsymbol{T}\widehat{\boldsymbol{\Sigma}}_{N}\boldsymbol{T})=\operatorname{rank}(\boldsymbol{T}\boldsymbol{\Sigma}\boldsymbol{T})=2$, we have that the WTS follows, asymptotically, a $\chi^{2}_{g}$-distribution under the null hypothesis, where $g=\operatorname{rank}(\boldsymbol{T}\boldsymbol{\Sigma}\boldsymbol{T})=2$. The Wald-type test, however, compares $T_{N}$ to the quantile of a $\chi^{2}_{4}$-distribution. Thus, for a chosen significance level of $\alpha=0.05$ this results in a true asymptotic $(N\to\infty)$ type-I error of

[TABLE]

i.e., a strictly conservative behavior of the test. Here $\chi^{2}_{f;1-\alpha}$ denotes the $(1-\alpha)$-quantile of the $\chi^{2}_{f}$-distribution. Similarly, for $\alpha=0.1$ and $\alpha=0.01$ we obtain asymptotically inflated type-I error rates of $0.02$ and $0.0013$ (both again conservative), respectively. Moreover, the situation is even more complicated since $\lim_{N\to\infty}\operatorname{rank}(\boldsymbol{T}\widehat{\boldsymbol{\Sigma}}_{N}\boldsymbol{T})=\operatorname{rank}(\boldsymbol{T}\boldsymbol{\Sigma}\boldsymbol{T})$ is neither verifiable in practice nor holds in general.

We tackle this problem in the current paper, where we not only relax the assumption (3’) on the unknown covariance matrices but also gain a procedure that is more robust to deviations from symmetry and homoscedasticity. To this end, we consider a different test statistic, namely a multivariate version of the ANOVA-type statistic (ATS) proposed by [5] for repeated measures designs, which we obtain by erasing the Moore-Penrose term from (2.1):

[TABLE]

In the special two-sample case where we wish to test the null hypothesis $H_{0}:\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{2}={\bf 0}$, this is equivalent to the test statistic proposed by [15].

The drawback of the ATS for multivariate data is that it is not invariant under scale transformations of the components, e.g., under change of units ($cm\mapsto m$ or $g\mapsto kg$) in one or more components. We demonstrate this problem in a real data analysis given in the supplementary material, where we exemplify that a simple unit change can completely alter the test decision of the ATS. Thus, we consider a slightly modified version of the ATS, which we denote as MATS:

[TABLE]

Here, $\boldsymbol{\widehat{D}}_{N}=\operatorname{diag}\left(N/n_{i}\cdot{\widehat{\sigma}}_{is}^{2}\right),i=1,\dots,a,~{}s=1,\dots,d$, where ${\widehat{\sigma}}_{is}^{2}$ is the empirical variance of component $s$ in group $i$. A related test statistic has been proposed by [39] in the context of high-dimensional ($d\to\infty$) data for a special non-singular one-way MANOVA design. Here, we investigate in the classical multivariate case (with fixed $d$) how its finite sample performance can be enhanced considerably. We start by analyzing its asymptotic limit behavior.

Theorem 2.1**.**

Under Conditions (1), (2) and (3) and under $H_{0}:\boldsymbol{T}\boldsymbol{\mu}={\bf 0}$, the test statistic $Q_{N}$ in (2.2) has asymptotically, as $N\to\infty$, the same distribution as a weighted sum of $\chi^{2}_{1}$ distributed random variables, where the weights $\lambda_{is}$ are the eigenvalues of $\boldsymbol{V}=\boldsymbol{T}(\boldsymbol{T}\boldsymbol{D}\boldsymbol{T})^{+}\boldsymbol{T}\boldsymbol{\Sigma}$ for $\boldsymbol{D}=\operatorname{diag}\left(\kappa_{i}^{-1}\sigma_{is}^{2}\right)$ and $\boldsymbol{\Sigma}=\operatorname{diag}(\kappa_{i}^{-1}\boldsymbol{V}_{i})$, i.e.,

[TABLE]

*with $Z_{is}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\chi^{2}_{1}$ and ”$\xrightarrow{\mathcal{D}}$” denoting convergence in distribution.

Thus, we obtain an asymptotic level $\alpha$ benchmark test $\varphi_{N}={\mathds{1}}\{Q_{N}>c_{1-\alpha}\}$ for $H_{0}:\boldsymbol{T}\boldsymbol{\mu}={\bf 0}$, where $c_{1-\alpha}$ is the $(1-\alpha)$-quantile of $Z$. However, the distribution of $Z$ depends on the unknown variances $\sigma^{2}_{is},i=1,\dots,a,~{}s=1,\dots,d$ so that $\varphi_{N}$ is infeasible for most practical situations. For this reason we consider different bootstrap approaches in order to approximate the unknown limiting distribution and to derive adequate and asymptotically correct inference procedures based on $Q_{N}$ in (2.2). This will be explained in detail in the next section. Apart from statistical test decisions discussed in Section 4.1, a central part of statistical analyses is the construction of confidence intervals, which allows for deeper insight into the variability and the magnitude of effects. This univariate concept can be generalized to multivariate endpoints by constructing multivariate confidence regions and simultaneous confidence intervals for contrasts $\boldsymbol{h}^{\top}\boldsymbol{\mu}$ for any contrast vector $\boldsymbol{h}\in\mathbb{R}^{ad}$ of interest. Details on the derivation of such confidence regions for $\boldsymbol{h}^{\top}\boldsymbol{\mu}$ are given in Section 4.2 below.

3 Bootstrap Procedures

The first bootstrap procedure we consider is a parametric bootstrap approach as proposed by [23] for the WTS. The second one is a wild bootstrap approach, which has already been successfully applied in the context of repeated measures or clustered data, see [9, 10] or [17]. The third procedure is a group-wise, nonparametric bootstrap approach. All of these bootstrap approaches are based on the test statistic $Q_{N}$ in (2.2). Note that the procedures derived in the following can also be used for multiple testing problems, either by applying the closed testing principle [30, 38] or in the context of simultaneous contrast tests [19, 21].

3.1 A Parametric Bootstrap Approach

This asymptotic model based bootstrap approach has successfully been used in univariate one- and two-way factorial designs ([24, 46]), and has recently been applied to Wald-type statistics for general MANOVA by [23] and [37], additionally assuming finite fourth moments. The approach is as follows: Given the data, we generate a parametric bootstrap sample as

[TABLE]

The idea behind this method is to obtain an accurate finite sample approximation by mimicking the covariance structure given in the observed data. This is achieved by calculating $Q_{N}^{*}$ from the bootstrap variables $\boldsymbol{X}_{i1}^{*},\dots,\boldsymbol{X}_{in_{i}}^{*}$, i.e.,

[TABLE]

We then obtain a parametric bootstrap test by comparing the original test statistic $Q_{N}$ with the conditional $(1-\alpha)$-quantile $c^{*}_{1-\alpha}$ of its bootstrap version $Q_{N}^{*}$.

Theorem 3.1**.**

The conditional distribution of $Q_{N}^{*}$ weakly approximates the null distribution of $Q_{N}$ in probability for any parameter $\boldsymbol{\mu}\in\mathbb{R}^{ad}$ and $\boldsymbol{\mu}_{0}$ with $\boldsymbol{T}\boldsymbol{\mu}_{0}={\bf 0}$, i.e.,

[TABLE]

where $\Pr_{\boldsymbol{\mu}}(Q_{N}\leq x)$ and $\Pr_{\boldsymbol{\mu}}(Q_{N}^{*}\leq x|\boldsymbol{X})$ denote the unconditional and conditional distribution of $Q_{N}$ and $Q_{N}^{*}$, respectively, if $\boldsymbol{\mu}$ is the true underlying mean vector.

3.2 A Wild Bootstrap Approach

Another resampling approach, which is based on multiplying the fixed data with random weights, is the so-called wild bootstrap procedure. To this end, let $W_{ik}$ denote i.i.d. random variables, independent of $\boldsymbol{X}$, with ${\rm E}(W_{ik})=0,\operatorname{Var}(W_{ik})=1$ and $\sup_{i,k}{\rm E}(W_{ik}^{4})<\infty$. We obtain a bootstrap sample as

[TABLE]

Note that there are different choices for the random weights $W_{ik}$, e.g., Rademacher distributed random variables [13] or weights satisfying different moment conditions, see, e.g., [3, 28, 29, 45]. The choice of the weights typically depends on the situation. In our simulation studies, we have investigated the performance of different weights such as Rademacher distributed as well as $\mathcal{N}(0,1)$ distributed weights (see [26] for this specific choice). The results were rather similar and we therefore only present the results of our simulation study for standard normally distributed weights in Section 5 below.

Based on the bootstrap variables $\boldsymbol{X}_{ik}^{\star}$, we can calculate $Q_{N}^{\star}$ in the same way as described for $Q_{N}^{*}$ in (3.1) above. A wild bootstrap test is finally obtained by comparing $Q_{N}$ to the conditional $(1-\alpha)$-quantile of its wild bootstrap version $Q_{N}^{\star}$.

Theorem 3.2**.**

The conditional distribution of $Q_{N}^{\star}$ weakly approximates the null distribution of $Q_{N}$ in probability for any parameter $\boldsymbol{\mu}\in\mathbb{R}^{ad}$ and $\boldsymbol{\mu}_{0}$ with $\boldsymbol{T}\boldsymbol{\mu}_{0}={\bf 0}$, i.e.,

[TABLE]

3.3 A nonparametric bootstrap approach

The third approach we consider is the nonparametric bootstrap. Here, for each group $i=1,\dots,a$, we randomly draw $n_{i}$ independent selections $\boldsymbol{X}^{\dagger}_{ik}$ with replacement from the $i$-th sample $\{\boldsymbol{X}_{i1},\dots,\boldsymbol{X}_{in_{i}}\}$. The bootstrap test statistic $Q_{N}^{\dagger}$ is then calculated in the same way as described above, i.e., by recalculating $Q_{N}$ with $\boldsymbol{X}^{\dagger}_{ik},i=1,\dots,a,\,k=1,\dots,n_{i}$. Finally, a nonparametric bootstrap test is obtained by comparing the original test statistic $Q_{N}$ to the empirical $(1-\alpha)$-quantile of $Q_{N}^{\dagger}$. The asymptotic validity of this method is guaranteed by

Theorem 3.3**.**

The conditional distribution of $Q_{N}^{\dagger}$ weakly approximates the null distribution of $Q_{N}$ in probability for any parameter $\boldsymbol{\mu}\in\mathbb{R}^{ad}$ and $\boldsymbol{\mu}_{0}$ with $\boldsymbol{T}\boldsymbol{\mu}_{0}={\bf 0}$, i.e.,

[TABLE]

4 Statistical Applications

We now want to base statistical inference on the modified test statistic in (2.2) using the bootstrap approaches described above. A thorough statistical analysis should ideally consist of two parts: First, statistical tests give insight into significant effects of the different factors as well as possible interactions. We therefore consider important properties of statistical tests based on the bootstrap approaches in Section 4.1. Second, it is helpful to construct confidence regions for the unknown parameters of interest in order to gain a more detailed insight into the nature of the estimates. The derivation of such confidence regions is discussed in Section 4.2.

4.1 Statistical Tests

In this section, we analyze the statistical properties of the bootstrap procedures described above. For ease of notation, we will only state the results for the parametric bootstrap procedure, i.e., consider the test statistic $Q_{N}^{*}$ based on $\boldsymbol{X}_{ik}^{*}$ throughout. Note, however, that the results are also valid for the wild and the nonparametric bootstrap procedure, i.e., the test statistics $Q_{N}^{\star}$ and $Q_{N}^{\dagger}$.

As mentioned above, a bootstrap test $\varphi^{*}={\mathds{1}}\{Q_{N}>c^{*}_{1-\alpha}\}$ is obtained by comparing the original test statistic $Q_{N}$ to the $(1-\alpha)$-quantile $c_{1-\alpha}^{*}$ of its bootstrap version. In particular, $p$-values are numerically computed as follows:

• (1)

Given the data $\boldsymbol{X}$, calculate the MATS $Q_{N}$ for the null hypothesis of interest.

• (2)

Bootstrap the data with either of the bootstrap approaches described above and calculate the corresponding test statistic $Q_{N}^{*,1}$.

• (3)

Repeat step (2) a large number of times, e.g., $B=10,000$ times, and obtain values $Q_{N}^{*,1},\dots,Q_{N}^{*,B}$.

• (4)

Calculate the $p$-value based on the empirical distribution of $Q_{N}^{*,1},\dots,Q_{N}^{*,B}$ as

[TABLE]

Theorems 3.1 – 3.3 imply that the corresponding tests asymptotically keep the pre-assigned level $\alpha$ under the null hypothesis and are consistent for any fixed alternative $\boldsymbol{T}\boldsymbol{\mu}\neq{\bf 0}$, i.e., $E_{\boldsymbol{\mu}}(\varphi^{*})\to\alpha\cdot{\mathds{1}}\{\boldsymbol{T}\boldsymbol{\mu}={\bf 0}\}+{\mathds{1}}\{\boldsymbol{T}\boldsymbol{\mu}\neq{\bf 0}\}$. Moreover, for local alternatives $H_{1}:\boldsymbol{T}\boldsymbol{\mu}=\sqrt{N}^{-1}\boldsymbol{T}\boldsymbol{\nu},\boldsymbol{\nu}\in\mathbb{R}^{ad}$, the bootstrap tests have the same asymptotic power as $\varphi_{N}={\mathds{1}}\{Q_{N}>c_{1-\alpha}\}$, where $c_{1-\alpha}$ is the $(1-\alpha)$-quantile of $Z$ given in Theorem 2.1. In particular, the asymptotic relative efficiency of the bootstrap tests compared to $\varphi_{N}$ is 1 in this situation.

4.2 Confidence regions and confidence intervals for contrasts

In order to conduct a thorough statistical analysis, interpretation of the results should not be based on $p$-values alone. In addition, it is helpful to construct confidence regions for the unknown parameter. The concept of a confidence region is the same as that of a confidence interval in the univariate setting: We want to construct a multivariate region, which is likely to contain the true, but unknown parameter of interest. The aim of this section is to derive multivariate confidence regions and simultaneous confidence intervals for contrasts $\boldsymbol{h}^{\top}\boldsymbol{\mu}$ for any contrast vector $\boldsymbol{h}$ of interest. Such contrasts include, e.g., the difference in means $\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{2}$ in two-sample problems, Dunnett’s many-to-one comparisons, Tukey’s all-pairwise comparisons, and many more, see, e.g., [21] for specific examples. In this section, we will base the derivation of confidence regions on the bootstrap approximations given in Section 3, i.e., we will use one of the bootstrap quantiles. Again, we only formulate the results for $c^{*}_{1-\alpha}$.

For the derivation of a confidence region, first note that the results from Section 4.1 imply that the null hypothesis $H_{0}:\boldsymbol{H}\boldsymbol{\mu}=\boldsymbol{H}\boldsymbol{\mu}_{0}$ for a vector of contrasts $\boldsymbol{H}\boldsymbol{\mu}_{0},\boldsymbol{H}=(\boldsymbol{h}_{1}|\dots|\boldsymbol{h}_{q})^{\top}\in\mathbb{R}^{q\times ad},\boldsymbol{\mu}_{0}\in\mathbb{R}^{ad},$ is rejected at asymptotic level $\alpha$, if $N(\boldsymbol{H}\boldsymbol{\overline{X}}_{\bullet}-\boldsymbol{H}\boldsymbol{\mu}_{0})^{\top}(\boldsymbol{H}\boldsymbol{\widehat{D}}_{N}\boldsymbol{H}^{\top})^{+}(\boldsymbol{H}\boldsymbol{\overline{X}}_{\bullet}-\boldsymbol{H}\boldsymbol{\mu}_{0})$ is larger than the bootstrap quantile $c^{*}_{1-\alpha}$. Thus, a confidence region for the vector of contrasts $\boldsymbol{H}\boldsymbol{\mu}$ is determined by the set of all $\boldsymbol{H}\boldsymbol{\mu}$ such that

[TABLE]

A confidence ellipsoid is now obtained based on the eigenvalues $\widehat{\lambda}_{s}$ and eigenvectors $\boldsymbol{\widehat{e}}_{s}$ of $\boldsymbol{H}\boldsymbol{\widehat{D}}_{N}\boldsymbol{H}^{\top}$. As in [22], the direction and lengths of its axes are determined by going $\sqrt{\widehat{\lambda}_{s}\cdot c_{1-\alpha}^{*}/N}$ units along the eigenvectors $\boldsymbol{\widehat{e}}_{s}$ of $\boldsymbol{H}\boldsymbol{\widehat{D}}_{N}\boldsymbol{H}^{\top}$. In other words, the axes of the ellipsoid are given by

[TABLE]

Note that this approach is similar to the construction of confidence intervals in the univariate case, where we exploit the one-to-one relationship between CIs and tests. While we can calculate (4.1) for arbitrary dimensions $d$, we cannot display the joint confidence region graphically for $d\geq 4$. In the two-sample case with $d=2$ endpoints, however, the ellipse can be plotted: Beginning at the center $\boldsymbol{H}\boldsymbol{\overline{X}}_{\bullet}$ the axes of the ellipsoid are given by $\pm\sqrt{\widehat{\lambda}_{s}\cdot c_{1-\alpha}^{*}/N}\cdot\boldsymbol{\widehat{e}}_{s},~{}s=1,2$. That is, the confidence ellipse extends $\sqrt{\widehat{\lambda}_{s}\cdot c_{1-\alpha}^{*}/N}$ units along the estimated eigenvector $\boldsymbol{\widehat{e}}_{s}$ for $s=1,2$. Therefore, we get a graphical representation of the relation between the group-mean differences $\mu_{11}-\mu_{21}$ and $\mu_{12}-\mu_{22}$ of the first and second component, see Section 10 and Figure 6 in the supplementary material for an example.

Concerning the derivation of multiple contrast tests and simultaneous confidence intervals for contrasts, we consider the family of hypotheses

[TABLE]

As shown in Sections 2 and 3 a test statistic for testing the null hypothesis $H_{0}:\boldsymbol{H}\boldsymbol{\mu}={\bf 0}$ is given by $Q_{N}$ in (2.2). Consequently, working with a single contrast $\boldsymbol{h}_{\ell}$ as contrast matrix leads to the test statistic

[TABLE]

for the null hypotheses $H_{0}^{\ell}:\boldsymbol{h}^{\top}_{\ell}\boldsymbol{\mu}={\bf 0},\ell=1,\dots,q$. Here, $\boldsymbol{h}_{\ell}=(h_{\ell,11},\dots,h_{\ell,ad})^{\top}$.

To obtain a single critical value with one of the bootstrap methods we may, e.g., consider the usual maximum or sum statistics. We exemplify the idea for the latter. Thus, let

[TABLE]

and denote by ${q}_{1-\alpha}^{*}$ the conditional $(1-\alpha)$-quantile of its corresponding bootstrap version $S_{N}^{*}$. From the proofs of Theorem 3.1 – 3.3 given in the supplement it follows for any of the three bootstrap methods described in Section 3 that $\widehat{\sigma}_{is}^{*}$ are consistent estimates for ${\sigma}_{is}$ ($i=1,\dots,a,\,s=1,\dots,d$) and that $\sqrt{N}\boldsymbol{H}\boldsymbol{\overline{X}}_{\bullet}^{*}$ asymptotically mimics the distribution of $\sqrt{N}\boldsymbol{H}(\boldsymbol{\overline{X}}_{\bullet}-\boldsymbol{\mu})$. Thus, the continuous mapping theorem implies $\Pr\left(S_{N}\leq{q}^{*}_{1-\alpha}\right)\to 1-\alpha$ as $N\to\infty$ and therefore

[TABLE]

This implies, that simultaneous $100(1-\alpha)\%$ confidence intervals for contrasts \boldsymbol{h}_{\ell}^{\top}\boldsymbol{\mu},{\color[rgb]{0,0,0}\ell=1,\dots,q,} are given by

[TABLE]

In the supplement we additionally explain that the bootstrap idea also works for the usual maximum statistic.

5 Simulations

The procedures described in Section 3 are valid for large sample sizes. In order to investigate their behavior for small samples, we have conducted various simulations. In the simulation studies, the behavior of the proposed approaches was compared to a parametric bootstrap approach for the WTS as in [23] since this turned out to perform better than other resampling versions of the WTS and Wilk’s $\Lambda$. For comparison, we also included the asymptotic $chi^{2}$ approximation of the WTS. All simulations were conducted using R Version 3.3.1 [35] each with nsim = 5,000 simulation and nboot = 5,000 bootstrap runs. We investigated a one- and a two-factorial design.

5.1 One-way layout

For the one-way layout, data was generated as in [23]. We considered $a=2$ treatment groups and $d\in\{4,8\}$ endpoints as well as the following covariance settings:

[TABLE]

Setting 1 represents a compound symmetry structure, while setting 2 is an autoregressive covariance structure. Both settings 1 and 2 represent homoscedastic scenarios while settings 3 and 4 display two scenarios with unequal covariance structures. Data was generated by

[TABLE]

where $\boldsymbol{V}_{i}^{1/2}$ denotes the square root of the matrix $\boldsymbol{V}_{i}$, i.e., $\boldsymbol{V}_{i}=\boldsymbol{V}_{i}^{1/2}\cdot\boldsymbol{V}_{i}^{1/2}$. The mean vectors $\boldsymbol{\mu}_{i}$ were set to ${\bf 0}$ in both groups. The i.i.d. random errors $\boldsymbol{\epsilon}_{ik}=(\epsilon_{ik1},\dots,\epsilon_{ikd})^{\top}$ with mean ${\rm E}(\boldsymbol{\epsilon}_{ik})={\bf 0}_{d}$ and $\operatorname{Cov}(\boldsymbol{\epsilon}_{ik})=\boldsymbol{I}_{d\times d}$ were generated by simulating independent standardized components

[TABLE]

for various distributions of $Y_{iks}$. In particular, we simulated normal, $\chi^{2}_{3}$, lognormal, $t_{3}$ and double-exponential distributed random variables. We investigated balanced as well as unbalanced designs with sample size vectors $\boldsymbol{n}^{(1)}=(10,10)^{\top}$, $\boldsymbol{n}^{(2)}=(20,20)^{\top}$, $\boldsymbol{n}^{(3)}=(10,20)^{\top}$ and $\boldsymbol{n}^{(4)}=(20,10)^{\top}$, respectively. A major criterion concerning the accuracy of the procedures is their behavior in situations where increasing variances (settings 3 and 4 above) are combined with increasing sample sizes ($\boldsymbol{n}^{(3)}$, positive pairing) or decreasing sample sizes ($\boldsymbol{n}^{(4)}$, negative pairing).

In this setting, we tested the null hypothesis $H_{0}^{\mu}:\{(\boldsymbol{P}_{a}\otimes\boldsymbol{I}_{d})\boldsymbol{\mu}\ ={\bf 0}\}=\{\boldsymbol{\mu}_{1}=\boldsymbol{\mu}_{2}\}$, i.e., no treatment effect. The resulting type-I error rates (nominal level $\alpha=5\%$) for $d=4$ and $d=8$ endpoints are displayed in Table 1 (normal distribution) and Table 2 ($\chi^{2}_{3}$-distribution), respectively. Further simulation results for lognormal, $t_{3}$ and double-exponential distributed errors for the parametric bootstrap of WTS and MATS can be found in Tables 7 – 9 in the supplementary material.

As already noticed by [23], the WTS with the $\chi^{2}$ approximation is far too liberal, reaching type-I error rates of more than 50% in some scenarios (e.g., for $d=8$ with negative pairing, i.e., covariance setting S3 and $\boldsymbol{n}=(20,10)^{\top}$). Even in the scenarios with only $d=4$ dimensions and $\boldsymbol{n}=(20,20)^{\top}$, the error rates are around 9% instead of 5%. The parametric bootstrap of the WTS greatly improves this behavior for all situations. However, it still shows a rather liberal behavior with type-I error rates of around 10% in some situations, e.g., $d=8$ dimensions with S3 or S4 and $\boldsymbol{n}=(20,10)^{\top}$ in Tables 1 and 2.

The wild bootstrap of the MATS shows a rather liberal behavior across all scenarios and can therefore not be recommended in practice. In contrast, both the parametric and the nonparametric bootstrap of the MATS show a very accurate type-I error rate control. The nonparametric bootstrap is often slightly more conservative than the parametric bootstrap and thus works better in situations with negative pairing, especially for the $\chi^{2}_{3}$-distribution, i.e., for S3 and S4 with $\boldsymbol{n}=(20,10)^{\top}$ and $d=4$ or $d=8$ dimensions in Table 2. In most other scenarios, however, the parametric bootstrap yields slightly better results. The improvement of the parametric bootstrap MATS over WTS (PBS) and nonparametric bootstrap MATS is most pronounced for large $d$, i.e., in situations where $d$ is close to $\min(n_{1},n_{2})$.

However, in situations with negative pairing and skewed distributions (see Table 2 as well as Table 7 in the supplementary material), the parametric bootstrap MATS shows a slightly liberal behavior. For $t_{3}$ and double-exponentially distributed errors and negative pairing, in contrast, the parametric bootstrap MATS is slightly conservative, see Tables 8 and 9 in the supplementary material, respectively.

Surprisingly, the resampling approaches based on the MATS improve with growing $d$ in most settings, i.e., when the number of endpoints is closer to the sample size. The WTS approach, in contrast, gets worse in these scenarios. This might be an interesting approach for future research in high-dimensional settings such as in [31].

As a result, we find that the MATS with the parametric bootstrap approximation is the best procedure in most scenarios. Especially, it is less conservative than the nonparametric bootstrap approximation and less liberal than the WTS equipped with the parametric bootstrap approach over all simulation settings. Only in situations with negative pairing and skewed distributions, the new procedure shows a slightly liberal behavior.

5.1.1 Singular Covariance Matrix

In order to analyze the behavior of the discussed methods in designs involving singular covariance matrices, we considered the one-way layout described above with $a=2$ groups and $d\in\{4,8\}$ observations involving the following covariance settings (displayed for $d=4$):

[TABLE]

Setting 6 is based on an $AR(0.6)$ covariance matrix (see setting 2 above), where the last row and column have been replaced by half the row/column before, respectively. Setting 7 is based on $\boldsymbol{\widetilde{V}}_{1}=\operatorname{diag}(\sqrt{2^{s}}),s=0,\dots,d-1$, where the last row and column have been replaced by half the row/column before. We have considered the same sample size vectors as above.

The results are displayed in Tables 3 and 4. The parametric bootstrap of the MATS again yields the best results in almost all scenarios. The wild bootstrap, in contrast, is again rather liberal. For the $\chi^{2}$ approximation of the WTS, the results are in concordance with the theoretical reflections mentioned in Section 2: Covariance setting S5 corresponds to the case, where the rank of $\boldsymbol{T}$ and $\boldsymbol{T}\boldsymbol{\Sigma}\boldsymbol{T}$ differs and as calculated above, the $\chi^{2}$-approximation becomes very conservative here. In setting S6 and S7, in contrast, there is no rank jump despite the singular covariance matrices and the $\chi^{2}$-approximation shows its usual liberal behavior. Since the rank of $\boldsymbol{T}\boldsymbol{\Sigma}\boldsymbol{T}$ is not known in practice, the WTS should not be used for data with possibly singular covariance matrices. It turns out, however, that the parametric bootstrap of the WTS is relatively robust against singular covariance matrices. Its behavior is comparable to the scenarios above with non-singular covariance matrices. It is, however, rather liberal for $\boldsymbol{n}=(20,10)^{\top}$, especially with the $\chi^{2}_{3}$-distribution, see Table 4. This behavior is improved by the parametric bootstrap MATS, e.g., for $d=8$ and S7, the WTS (PBS) leads to a type-I error of 9%, whereas the MATS (PBS) is at 5.1%. The nonparametric bootstrap, in contrast, sometimes leads to strictly conservative test decisions. This is especially apparent for $d=8$ and covariance setting S7 in Tables 3 and 4.

5.2 Two-way layout

We have investigated the behavior of the methods in a setting with two crossed factors $A$ and $B$, which is again adapted from [23]. In particular, we simulated a $2\times 2$ designs with covariance matrices similar to the one-way layout above. A detailed description of the simulation settings as well as the results for the main and interaction effects are deferred to the supplementary material. Here we only summarize our findings: Since the total sample size $N$ is larger in this scenario, the asymptotic results come into play and therefore all methods lead to more accurate results than in the one-way layout. Nevertheless we find a similar behavior as in the one-way layout: Again, the MATS and the WTS with the parametric bootstrap approach control the type-I error very accurately, whereas the nonparametric bootstrap approach leads to slightly more conservative results. Both the WTS with $\chi^{2}$ approximation and the wild bootstrap MATS can not be recommended due to their liberal behavior. In situations with negative pairing (covariance setting 10 and 11 with sample size vector $\boldsymbol{n}^{(3)}$), the parametric bootstrap MATS improves the slightly liberal behavior of the WTS, see e.g., Table 10 for the normal distribution, where the WTS (PBS) leads to a type-I error of 6.1%, while the MATS (PBS) is at 4.9%.

5.3 Power

We have investigated the empirical power of the proposed methods to detect a fixed alternative in the simulation scenarios above. Data was simulated as described in Section 5.1 but now with $\boldsymbol{\mu}_{1}={\bf 0}$ and $\boldsymbol{\mu}_{2}=(\delta,\dots,\delta)^{\top}$ for varying shifts $\delta\in\{0,0.5,1,1.5,2,3\}$. Due to the liberality of the classical Wald-type test and the wild bootstrapped MATS, we have only considered the WTS with parametric bootstrap as well as the parametric and nonparametric bootstrap of the MATS. The results for selected scenarios are displayed in Figures 1 – 2. The plots show that both resampling versions of the MATS have a higher power for detecting the fixed alternative than the WTS. The parametric bootstrap of the MATS has a slightly higher power than the nonparametric bootstrap, a behavior that is more pronounced for the $\chi^{2}$-distribution (Figure 1). Moreover, the power analysis shows a clear advantage of applying the parametric bootstrap approach to the MATS over its application to the WTS. For example, in the scenario with normally distributed data, $d=8$ dimensions, covariance setting S4 and $\boldsymbol{n}=(10,20)^{\top}$ observations (Figure 2), the parametric bootstrap MATS has twice as much power as its WTS version in case of $\delta=0.5$ (34.4% as compared to 16.7%). Similar differences can also be observed in some of the other settings.

6 Application: Analysis of the Data Example

As a data example, we consider 7 demographic factors of US citizens in 43 states. Our aim is to investigate whether these factors differ between the states. The full data set ‘county_facts.csv‘ is available from kaggle (https://www.kaggle.com/joelwilson/2012-2016-presidential-elections). In order to have sufficient sample sizes for the analysis and to avoid a high-dimensional setting, we exclude all states with less than 15 counties. In particular, we removed Connecticut, Delaware, Hawaii, Massachusetts, New Hampshire, Rhode Island and Vermont. We consider the following demographic factors: the population estimate for 2014 (PST045214), the percentage of female citizens in 2014 (SEX255214) as well as the percentage of white (RHI125214), black or African American (RHI225214), American Indian and Alaska native (RHI325214), Asian (RHI425214) and native Hawaiian and other pacific islanders (RHI525214) citizens in 2014. This results in a one-way layout with $a=43$ levels of the factor ‘state‘ and $d=7$ dimensions. The sample sizes and mean values for the different states can be found in Table 13. Figure 3 exemplarily displays boxplots for the percentage of white citizens across the different states.

[FIGURE:]

We now want to analyze, whether there is a significant difference in the multivariate means for the different states. The null hypothesis of interest thus is $H_{0}:\{(\boldsymbol{P}_{43}\otimes\boldsymbol{I}_{7})\boldsymbol{\mu}={\bf 0}\}$. Since the empirical covariance matrix is computationally singular in this example (reciprocal condition number 1.7e-16), we cannot apply the Wald-type test. Thus, we consider the parametric bootstrap approach of the MATS which yielded the best results in the simulation study. Computation of the MATS results in a value of $Q_{N}=393.927$ and the parametric bootstrap routine with 1,000 bootstrap runs gives a $p$-value of $p<0.0001$, implying that there is indeed a significant difference between the states with respect to the 7 demographic measurements.

A confidence region for this effect can be constructed as described in Section 4. The analysis of this example, including the calculation of the confidence region, can be conducted using the R package MANOVA.RM.

7 Conclusions and Discussion

We have investigated a test statistic for multivariate data (MATS) which is based on a modified Dempster statistic. Contrary to classical MANOVA models, we incorporate general heteroscedastic designs and allow for singular covariance matrices while postulating their existence as solely distributional assumption. Moreover, our proposed MATS statistic is invariant under linear transformations of the response variables.

In order to improve the small sample behavior of the test statistic, we have investigated different bootstrap approaches, namely a parametric bootstrap, a wild bootstrap and a nonparametric bootstrap procedure. We have rigorously proven that they lead to asymptotically exact and consistent tests and even analyzed their local power behavior.

In a large simulation study, the parametric bootstrap turned out to perform best in most scenarios, even with skewed data and heteroscedastic variances. Although the type-I error control is still not ideal in the latter case, the method performed advantageous over the parametric bootstrap of the WTS proposed in [23] and has the additional advantage of being applicable to situations with singular covariance matrices. In situations with skewed distributions, the parametric bootstrap of the MATS yielded more robust results than the WTS. The wild bootstrap approach, in contrast, turned out to be very liberal in all scenarios, while the nonparametric bootstrap was mostly slightly more conservative than the parametric bootstrap. Power simulations showed a clear advantage of the parametric bootstrap MATS compared to the WTS (PBS) as well as the nonparametric bootstrap. All in all, we therefore recommend the parametric bootstrap based on the MATS for practical applications in a multivariate setting.

Furthermore, we have constructed confidence regions and simultaneous confidence intervals for contrasts $\boldsymbol{h}^{\top}\boldsymbol{\mu}$ based on the bootstrap quantiles. These confidence regions provide an additional benefit for the analysis of multivariate data since they allow for more detailed insight into the nature of the estimates.

In order to facilitate application of the proposed methods, the parametric bootstrap test and the calculation of confidence regions are implemented in the R package MANOVA.RM.

Following the idea of [39] we plan to extend our concepts to the high-dimensional setting, i.e., where the sample size $N$ may be less than the dimension $d$. This approach looks promising, since we have seen in the simulation study that the MATS with the parametric bootstrap approach exhibited an improved type-I error control with increasing $d$. However, the extension to high-dimensional data requires different techniques and will be part of future research.

Acknowledgment

The authors would like to thank Dr. Jan Paul and Prof. Dr. Volker Rasche for providing the cardiology data example used in the supplement. This work was supported by the German Research Foundation projects DFG PA 2409/3-1 and PA 2409/4-1.