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. Furthermore, existing potential remedies (e.g., heteroskedasticity-robust approaches) become inappropriate in cases where the covariance matrices are singular or close to 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. Moreover, computational issues may also lead to singular covariance structures. 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 small-sample 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. Yet, for large sample sizes, Wald-type approaches are an attractive alternative option. Moreover, the choice between Wald-type and MANCATS-based tests might depend on the particular setting (e.g., the covariance structure, or the alternative under consideration).
