Multivariate analysis of covariance when standard assumptions are violated
Georg Zimmermann Markus Pauly Arne C. Bathke








Multivariate analysis of covariance when standard assumptions are violated
    Georg ZimmermannGeorg Zimmermann, Department of Mathematics, Paris Lodron University, Hellbrunner Strasse 34, A 5020 Salzburg, Austria (georg.zimmermann@pmu.ac.at)

Department of Mathematics, Paris Lodron University, Salzburg, Austria 

Department of Neurology, Paracelsus Medical University, Salzburg, Austria

Spinal Cord Injury and Tissue Regeneration Centre Salzburg,
 Paracelsus Medical University, Salzburg, Austria

	Markus Pauly 
 Department of Statistics, Ulm University, Germany 

	Institute for Mathematical Statistics and Applications in Industry, 
TU Dortmund, Dortmund, Germany

	Arne C. Bathke 

	Department of Mathematics, Paris Lodron University Salzburg, Austria 

	Department of Statistics, University of Kentucky, Lexington, KY, USA
    
====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




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.  


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




§ 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 <cit.>. 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 <cit.> and psychology <cit.>. 


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 <cit.>), were used. Analogously to their univariate counterparts, they rely on the assumption that the error vectors are i.i.d. multivariate normal, with expectation 0 and identical covariance matrix Σ 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 <cit.>, there are only few methods available that are potentially applicable in a more general setting. Recently, <cit.> 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 <cit.>, 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 <cit.> and later examined in combination with various bootstrap procedures <cit.> 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 <cit.> 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 <cit.>),  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 <cit.>. 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, <cit.>). 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 <ref>, 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 <ref>, 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 <ref>. 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 <ref>. Finally, Section <ref> 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.    
     



§ THE GENERAL MULTIVARIATE ANALYSIS OF COVARIANCE (MANCOVA) MODEL

 

In the sequel, let I_m denote the m-dimensional identity matrix, and let ⊕ and ⊗ denote the direct sum and the Kronecker product of matrices, respectively. 
Let Y_ij = (Y_ij1,Y_ij2,…,Y_ijp)' and z_ij = (z_ij^(1),z_ij^(2),…,z_ij^(c))' denote the p- and c-dimensional outcome and covariate vectors of subject j in group i, i ∈{1,2,…,a}, j ∈{1,2,…,n_i}. We assume that the covariates are fixed, whereas the outcomes are random variables, satisfying the model equation

    Y_ij = μ_i + (z_ij^'⊗)ν + ϵ_ij = μ_i + ∑_w=1^cz_ij^(w)ν^(w) + ϵ_ij.

Thereby, μ_i = (μ_i1,…,μ_ip)' denotes the vector of adjusted means in group i, and ν = (ν^(1)',…,ν^(c)')', where ν^(w) contains the p regression coefficients modelling the association between the w-th covariate and the components 1,…,p of the outcome, w ∈{1,2,…,c}. It should be noted that model ModelExplicit 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 ϵ_ij are independent, with E[ϵ_ij] = 0 and group-specific covariance matrix Σ_i:= Cov[ϵ_ij], i ∈{1,…,a}, j ∈{1,… n_i}. 


Model ModelExplicit can be expressed in a more compact form by using matrix notation: Let Y = (Y_11^',Y_12^',…,Y_an_a^')^', Z = (z_11,z_12,…,z_an_a)^', ϵ = (ϵ_11^',ϵ_12^',…,ϵ_an_a^')^', and μ = (μ_1^',…,μ_a^')^'. Moreover, let 1_m denote the m-dimensional vector containing all 1's. Then, model ModelExplicit is equivalent to 

    Y = M̃μ + Z̃ν + ϵ,

where M̃ = ⊕_i=1^a(1_n_i⊗), Z̃ = Z⊗ and Cov[ϵ] = ⊕_i=1^a(I_n_i⊗Σ_i). So, in particular, μ can be estimated by the ordinary least squares estimator μ̂ = (μ̂_1^',…,μ̂_a^')^', that is,

    μ̂_i = Y̅_i. - ∑_w=1^cz̅_i.^(w)ν̂^(w),

where the dot notation indicates averaging over all subjects in the respective group, and ν̂^(w) denotes the OLS estimator of ν^(w), w = 1,…,c. We are interested in testing hypotheses about the adjusted mean vectors μ_1,…,μ_a of the form H_0: Hμ = 0, where H denotes a contrast matrix (i.e.,H1_ap = 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 T := H'(HH')^+H. Thereby, (HH')^+ denotes the Moore-Penrose inverse of HH'. Note that T is uniquely defined, and Tμ = 0⇔Hμ = 0. For example, the hypothesis H_0: μ_1 = … = μ_a is equivalent to H_0: Tμ = 0, with T = P_a ⊗. Thereby, P_a = I_a - 1/aJ_a, where 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 <cit.>. 


Now, for testing the hypothesis H_0:Tμ = 0, one may consider the Wald-type statistic (WTS)

    W(T) = μ̂^'T^'(TΣ̂T)^+Tμ̂,

where Σ̂ denotes the upper-left block of the 2× 2 block diagonal matrix (X̃^'X̃)^-1X̃^'ŜX̃(X̃^'X̃)^-1. Thereby, X̃ = (M̃,Z̃), and Ŝ:= ⊕_i=1^a⊕_j=1^n_iΣ̂_ij (i.e., Ŝ is a block-diagonal matrix with matrices Σ̂_11,Σ̂_12,…,Σ̂_an_a on the diagonal), and Σ̂_ij = u_iju_ij^', u_ij = Y_ij - μ̂_i - (z_ij^'⊗)ν̂, i ∈{1,2,…,a}, j ∈{1,2,…,n_i}. Observe that this covariance matrix estimator basically represents the multivariate generalization of the approach proposed in <cit.> 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 <cit.>. 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 <ref>, we propose a multivariate generalization of the so-called HC4 version of Σ̂_ij, by multiplying the latter with 1/(1-p_ij,ij)^δ_ij, δ_ij:= min(4, p_ij,ij / (N^-1∑_r=1^a∑_s=1^n_rp_rs,rs)), where p_rs,rs denotes the diagonal element of the hat matrix X(X'X)^-1X' corresponding to individual s in group r, r ∈{1,2,…,a}, s ∈{1,2,…,n_r}. Thereby, X = (M,Z) = (⊕_i=1^a1_n_i, Z). This adjustment can be interpreted as a natural generalization of the univariate HC4 estimator, which was proposed by Cribari-Neto <cit.> because for p = 1, the covariance matrix estimator Σ̂_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 W(T) under H_0, the following assumptions are required, where all convergences are understood for N →∞:


  * The errors ϵ_ij are independent, with E[ϵ_ij] = 0, Cov[ϵ_ij] = Σ_i, and E(ϵ_ijk^4) ≤ C_1 < ∞, uniformly for all i ∈{1,2,…,a},j∈{1,2,…,n_i},k∈{1,2,…,p}.  

  * Σ_i > 0 for all i ∈{1,2,…,a}.

  * For all i ∈{1,2,…,a}: N/n_i→κ_i > 0. 

  * The columns of Z are linearly independent of each other and of the columns of ⊕_i=1^a1_n_i.

  * For all i ∈{1,2,…,a}, w ∈{1,…,c}: 1/n_i∑_j=1^n_iz_ij^(w)⟶ c_i^(w)∈ℝ.

  * For all i ∈{1,2,…,a}: 1/n_i∑_j=1^n_iz_ijz_ij^'⟶Ξ_i ∈ℝ^c× c.  


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



Assuming that model ModelMnot as well as WTSA1 – WTSA6 hold, the test statistic W(T) defined in WTS has an asymptotic (N→∞) χ_f^2 distribution under H_0: Tμ = 0, with f = rank(T). 


The assumptions WTSA1 and WTSA3 represent standard requirements in an asymptotic framework. Likewise, a violation of WTSA4 would mean that collinearity was present; in such a case, the validity of the results would be questionable anyway. Conditions WTSA5 and WTSA6 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, WTSA2 might be violated in a considerable number of practically relevant settings, as outlined in Section <ref>. 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) <cit.>, we consider a modified ANCOVA ANOVA-type statistic (MANCATS), which is defined as 

    A(T):=  μ̂^'T^'(TD̂T)^+Tμ̂.

Thereby,

    D̂ = ⊕_i=1^a1/n_i⊕_k=1^pσ̂_ik^2 = diag(n_1^-1σ̂_11^2,…,n_1^-1σ̂_1p^2,…,n_a^-1σ̂_ap^2),

where 

    σ̂_ik^2 = 1/n_i - c - 1∑_j=1^n_iu_ijk^2,

and u_ijk = Y_ijk - μ̂_ik - ∑_w=1^cz_ij^(w)ν̂_k^(w) is the residual of outcome k ∈{1,2,…,p} in group i ∈{1,2,…,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 - Y̅_i.k. Consequently, σ̂_ik^2 would be equal to the empirical variance estimator, which was used in <cit.>. So, Dest is a natural generalization of the MATS approach to the MANCOVA setting. As already mentioned before, when using A(T) instead of W(T), assumption WTSA2 can be replaced by the weaker requirement 


  (M2)^∗ Σ_i ≥ 0 and σ_ik^2 > 0, for all i ∈{1,2,…,a} and k ∈{1,2,…,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 D̂ in the MANCATS, the asymptotic chi-squared limit distribution from Theorem <ref> will not hold any more. Nevertheless, there is still at least a formal result regarding the asymptotic distribution of A(T).  



Under conditions WTSA1, MANCATSA1, WTSA3–WTSA6, and under H_0: Tμ = 0, the MANCATS test statistic A(T) defined in Mancats has, asymptotically (N →∞), the same distribution as the weighted sum 

    U = ∑_i=1^a∑_k=1^pλ_ikU_ik,

where U_iki.i.d.∼χ_1^2, and the weights λ_ik are the eigenvalues of T(TDT)^+TΛ_11. Thereby, D = ⊕_i=1^a⊕_k=1^pκ_iσ_ik^2, and Λ_11 denotes the upper-left block of the 2× 2 block matrix Λ = (Ξ^-1⊗)Ψ(Ξ^-1⊗) = lim_N →∞[(N^-1X̃^'X̃)^-1(N^-1X̃^'Cov(ϵ)X̃)(N^-1X̃^'X̃)^-1]. 


However, the weights λ_ik in Theorem <ref> 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.



§ 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 <cit.>, Liu <cit.> and Mammen <cit.>. 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 <cit.>, in regression models with clustered errors <cit.> and in multivariate factorial designs <cit.>, in order to improve the small-sample performance of the respective tests. In the present work, we use the following procedure: We generate N = ∑_i=1^an_i i.i.d. random variables T_ij independently from the data, with E(T_11) = 0, Var(T_11) = 1, and sup_i,jE(T_ij^4) < ∞, 1≤ i ≤ a, 1≤ j ≤ n_i. Then, the wild bootstrap observations are defined as 

    Y_ij^∗:= u_ijT_ij/√(1-p_ij,ij), 1≤ i ≤ a, 1≤ j ≤ n_i,

where p_ij,ij:=x_ij^'(X'X)^-1x_ij denotes the diagonal element of the hat matrix corresponding to individual j in group i, and u_ij is the p-dimensional residual vector of individual j in group i, i ∈{1,…,a}, j ∈{1,…,n_i}. The scaling factor (1-p_ij,ij)^-1/2 has been introduced by <cit.> 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 <ref> 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 β̂^∗:=(X̃'X̃)^-1X̃'Y^∗, the bootstrap covariance matrix estimator D̂^∗:= diag(n_1^-1σ̂_11^∗ 2,…,n_a^-1σ̂_ap^∗ 2) and, by plugging in the bootstrap versions instead of the original estimators in Mancats, the bootstrap analogon A^∗(T) of the MANCATS test statistic A(T) are calculated. The following theorem states that using the conditional (1-α)-quantile of the empirical distribution of the wild bootstrap test statistic A^∗(T) as the critical value indeed yields an asymptotically valid test.



For any parameter vector β:=(μ^',ν^')^'∈ℝ^(a+c)p and any β_0 = (μ_0^',ν_0^')' ∈ℝ^(a+c)p with Tμ_0 = 0, we have that given the data, as N →∞,

    sup_x∈ℝ|P_β(A^∗(T)≤ x | Y) - P_β_0(A(T)≤ x)|  0

holds, provided that WTSA1, MANCATSA1 and WTSA3–WTSA6 are fulfilled. Thereby,  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 <cit.>. Given the observed data, we impose a group-wise i.i.d. structure in the bootstrap world by drawing

    Y_ij^⋆i.i.d.∼𝒩(0,Σ̂_i), 1≤ i ≤ a, 1≤ j ≤ n_i,

where 
    Σ̂_i:=1/n_i-c-1∑_j=1^n_iu_iju_ij^', 1≤ i ≤ a.


Once the bootstrap observations have been generated, the parametric bootstrap version A^⋆(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:



For any parameter vector β:=(μ^',ν^')^'∈ℝ^(a+c)p and any β_0 = (μ_0^',ν_0^')' ∈ℝ^(a+c)p with Tμ_0 = 0, we have that given the data, as N →∞,

    sup_x∈ℝ|P_β(A^⋆(T)≤ x | Y) - P_μ_0(A(T)≤ x)|  0

holds, provided that  WTSA1, MANCATSA1 and WTSA3–WTSA6 are fulfilled.
  



§ SIMULATIONS



In order to investigate the finite-sample performance of the three methods, which have been discussed in Sections <ref> and <ref>, we have conducted simulations for a broad range of settings, using  version 3.5.1 <cit.>. 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,…,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 <cit.>. For each scenario, the data generation process was repeated n_sim = 10 000 times, and within each simulation run, n_boot = 5 000 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 ϵ_ij, at first we drew a random sample of Np independent observations ξ_111,ξ_112,…,ξ_an_ap from one out of several distributions (normal, χ_5^2, lognormal, double exponential). Then, we transformed them to standardized random variables ξ̃_ijk:= (ξ_ijk - E[ξ_ijk]) / (Var[ξ_ijk])^1/2 and subsequently calculated the products 

    ϵ_ij = Σ_i^1/2ξ̃_ij,
 
where Σ_i^1/2 denotes the matrix square root of the covariance matrix of group i, and ξ̃_ij=(ξ̃_ij1,…,ξ̃_ijp)^', i ∈{1,2,…,a}, j ∈{1,2,…,n_i}. Doing so, it is ensured that E[ξ̃_ij] = 0 and Cov[ξ̃_ij] = Σ_i. For each of the aforementioned error distributions, we considered three different choices of Σ_1 and Σ_2, namely Σ_1 =Σ_2 = I_2 + 0.5(J_2-I_2) (I), Σ_i = i·I_2 + 0.5(J_2-I_2) (II), and Σ_1 = Σ_2 = diag(1,0.25) + 0.5(J_2-I_2) (III). Observe that the latter covariance matrix is singular. Finally, the observations were obtained as Y_ij = μ_i + (z_ij^'⊗)ν + ϵ_ij, where ν = (-0.5,-1.0, 1.5, 1.0)'. 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 10 000 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 ν = (-0.5,-1.0, 1.5, 3.0)', 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) ∈{(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)^-δ_ij/2. Thereby, δ_ij:= min(4, p_ij,ij / (N^-1∑_r=1^a∑_s=1^n_rp_rs,rs)), where p_rs,rs denotes the diagonal element of the hat matrix corresponding to individual s in group r, r ∈{1,2,…,a}, s ∈{1,2,…,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 μ_i = 0, i∈{1,2}. The results are displayed in Table <ref>.  

 

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 <cit.>. 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 μ_1 = (0,0)' and μ_2 = (δ,0)', with δ∈{0.5,1.0,…,3.0}. For the singular setting III, we set μ = (δ,δ)', 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 <ref>. 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 <ref>). 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., δ∈{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 <ref>). 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 n_1 = (30,30,30,30), n_2 = (15,15,15,15), n_3 = (5,10,20,25), and n_4 = (25,20,10,5). The results are reported in Table <ref>. 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: Σ_i = I_4 + 0.5(J_4-I_4), II: Σ_i = i· I_4+ 0.5(J_4-I_4), and III: Σ_i = (s_1,s_2,s_3,s_4), where s_1 = s_2 = s_3 = (1,1,1,1/2)', s_4 = (1/2,1/2,1/2,1)', i ∈{1,2,3,4}. The vectors modelling the associations between a particular covariate and the components of the outcome were set to ν^(1) = (-0.5,-1.0,-2.0,-0.02)' and ν^(2) = (1.5,1.0,1.0,0.2)', respectively. For the singular setting III, we set ν^(1) = (-0.5,-1.0,-1.5,-0.02)' and ν^(2) = (1.5,3.0,4.5,0.2)'. Apart from that, all other specifications were the same as before. From Table <ref>, 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 <ref>, the results were similar. The parametric bootstrap MANCATS appeared to perform even slightly better, whereas the asymptotic Wald-type test was somewhat more conservative. 

 




 











§ 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 <cit.>. 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 

    Σ̂_h^(r) = (
    145.88     113.18 
    
    113.18     1073.21
    )  and  Σ̂_ℓ^(r) = (
    99.07     60.85 
    
    60.85     458.41
    ),

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 <cit.>. 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

    Σ̂_h^(s) = (
    145.88     113.18     259.06 
    
    113.18     1073.21     1186.39
    
    259.06     1186.39     1445.45
    )  and  Σ̂_ℓ^(s) = (
    99.07     60.85     159.92 
    
    60.85     458.41     519.26
    
    159.92     519.26     679.19
    ),

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 <ref>, we examined the empirical type I error rates, adopting the same specificiations for the data generation process as in the previous section, but with Σ_1 and Σ_2 replaced by Σ̂_h^(r) and Σ̂_l^(r), or Σ̂_h^(s) and Σ̂_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 <ref>, 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 <ref>. 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 <ref>, 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 <cit.>. 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. 

 



 



§ 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 <cit.> for uni- and multivariate linear models have been proposed in the literature (see, for example, <cit.> for an overview). We decided to use a multivariate generalization of the HC4 estimator <cit.>, since its performance turned out to be superior to other variants (HC0, HC2) in the univariate case <cit.>. 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 <cit.>. The reason is that in the context of heteroskedastic regression models, the nonparametric bootstrap has suboptimal properties with respect to variance estimation <cit.>. Moreover, in multivariate factorial designs, the nonparametic bootstrap performed somewhat worse than the respective competitors <cit.>.


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 <cit.>. 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 <ref>, efforts towards developing methods which allow for MANCOVA-like analyses of multivariate ordinal outcomes should be made in future research. 


To this end, applying probabilistic index models (PIM) might be a promising approach <cit.>. Apart from that, the authors conjecture that the ideas discussed in the present manuscript might be useful for ANCOVA models with clustered errors, too. Maybe, similarly to the MANCATS, it would suffice to apply a cluster-wise wild bootstrap procedure, even without subsequently using a cluster-robust covariance matrix estimator <cit.>. This would substantially facilitate calculations and most likely reduce the computational burden. Therefore, this topic might also be a promising goal of future research.           










  
 



Supplementary Material




Section <ref> 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 <ref>. 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 <ref>.  




§ PROOFS



 


 §.§ Notations and assumptions


Let I_m denote the m-dimensional identity matrix, and let ⊕ and ⊗ denote the direct sum and the Kronecker product of matrices, respectively.
In order to simplify notations, we consider Y = (⊗)  + ϵ, where E(ϵ)= 0, S:= Cov(ϵ) = ⊕_i=1^a(I_n_i⊗Σ_i), = (⊕_i=1^a1_n_i, Z), = (μ',ν')'. Moreover, let = (μ̂',ν̂')' denote the ordinary least squares estimator of , = (('/N)^-1⊗)(N^-1(^'⊗)S(⊗))(('/N)^-1⊗), = (('/N)^-1⊗)(N^-1(^'⊗)Ŝ(⊗))(('/N)^-1⊗), Ŝ = ⊕_i=1^a⊕_j=1^n_iΣ̂_ij, 
Σ̂_ij = u_iju_ij^', N = ∑_i=1^an_i. Thereby, u_ij denotes the p-dimensional residual vector corresponding to subject j in group i, 1≤ i ≤ a, 1≤ j ≤ n_i. 

Observe that assumption (M5) implies that the elements of  are uniformly bounded (at least for sufficiently large N). Moreover, (1/N') converges to, say, Ξ, due to (M3), (M5), and (M6). According to condition (M4),  is of full column rank, which in turn implies that 1/N' is positive definite. Consequently, (1/N')^-1 converges to Ξ^-1. Analogously, we get that 1/N('⊗)S(⊗) converges to a matrix Ψ, 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| = |x_ij'(X'X)^-1x_ij| ≤ C/N. 






 §.§ Proof of Theorem 1


The proof comprises two main steps:





  * √(N)(β̂-β) 𝒩(0,Λ), where Λ:=(Ξ^-1⊗)Ψ(Ξ^-1⊗), Ψ:= lim_N→∞(N^-1('⊗)S(⊗)), Ξ^-1:= lim_N →∞(N^-1')^-1. Note that both limits are well-defined (cf. Section <ref>).

  * - 0. 



Proof of <ref>:

Let = (x_11,x_12,…,x_an_a)^'. Doing some algebra yields 

    √(N)( - ) = √(N)((^')^-1^'⊗)ϵ = √(N)((^')^-1⊗)∑_i=1^a∑_j=1^n_i(x_ij⊗)ϵ_ij.

Hence, we define the triangular array η_ij:= √(N)((^')^-1⊗)(x_ij⊗)ϵ_ij, 1≤ i ≤ a, 1≤ j ≤ n_i, and apply the multivariate Lindeberg-Feller theorem.
Obviously, E(η_ij) = 0, and (η_ij)_i,j are independent, due to (M1). Secondly, we get

    ∑_i=1^a∑_j=1^n_iCov(η_ij)    = ∑_i=1^a∑_j=1^n_i N ((^')^-1x_ij⊗)Cov(ϵ_ij)(x_ij^'(^')^-1⊗)
       = N∑_i=1^a∑_j=1^n_i ((^')^-1x_ij⊗)Σ_i(x_ij^'(^')^-1⊗)
       = (^'/N)^-1⊗)(N^-1('⊗)S(⊗))(^'/N)^-1⊗)
       →Λ.

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

    ∑_i=1^a∑_j=1^n_iE(||η_ij||^2 1{||η_ij||^2 > δ}) → 0 ∀δ > 0.


Firstly, 

    ||η_ij||^2    = η_ij^'η_ij = N ϵ_ij^'(x_ij^'(^')^-1⊗) ((^')^-1x_ij⊗)ϵ_ij
       = 1/Nϵ_ij^'(x_ij^'(^'/N)^-1⊗) ((^'/N)^-1x_ij⊗)ϵ_ij.


Let η̃_ij:= ((^'/N)^-1⊗)(x_ij⊗)ϵ_ij (consequently, ||η_ij||^2 = N^-1 ||η̃_ij||^2). 







Using Lancaster's theorem, we get

    E(||η̃_ij||^2)    =tr((x_ij^'⊗) [1/N('⊗)]^-1[1/N('⊗)]^-1(x_ij^'⊗) Σ_i)
       = tr((x_ij⊗)Σ_i(x_ij^'⊗) [1/N('⊗)]^-1[1/N('⊗)]^-1).


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, 

    sup_i,j,NE(||η̃_ij||^2) ≤ C  for some  C>0.


Next, let A_ij:=(x_ij^'⊗)(('/N)^-1⊗)(('/N)^-1⊗)(x_ij⊗). Then, ||η̃_ij||^2 simplifies to ϵ_ij^'A_ijϵ_ij, which can be explicitely written as

    ||η̃_ij||^2 = ϵ_ij^'A_ijϵ_ij = ∑_k=1^p∑_ℓ = 1^pa_ij,kℓϵ_ijkϵ_ijℓ,

where a_ij,kℓ denotes the element in the k-th row and the ℓ-th column of A_ij. Consequently, the expectation E[||η̃_ij||^4] is equal to a sum of p^4 expectations of products ϵ_ijk_1ϵ_ijℓ_1ϵ_ijk_2ϵ_ijℓ_2, scaled with elements of the matrix 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[ϵ_ijk_1^4]E[ϵ_ijℓ_1^4]E[ϵ_ijk_2^4]E[ϵ_ijℓ_2^4], which is uniformly bounded, due to assumption (M1). So, all in all, the fourth moments E[||η̃_ij||^4] can be uniformly bounded from above by, say, C̃. Consequently,

    ∑_i=1^a∑_j=1^n_iE[||η_ij||^2 1{||η_ij||^2 > δ}]    = 
    1/N∑_i=1^a∑_j=1^n_iE[||η̃_ij||^2 1{||η̃_ij||^2 > δ N}] 
       ≤1/N∑_i=1^a∑_j=1^n_i√(E[||η̃_ij||^4])√(E[1{||η̃_ij||^2 > δ N}])
       ≤C̃/N∑_i=1^a∑_j=1^n_i√(E[1{||η̃_ij||^2 > δ N}])
       ≤C̃(C/δ N)^1/2,
 
where the last step was due to 

    E[1{||η̃_ij||^2 > δ N}]    = P(||η̃_ij||^2 > δ N)≤E[||η̃_ij||^2]/δ N≤C/δ N,

according to LindeExp. This completes the proof of the Lindeberg condition LindebergCond and, thus, the proof of <ref>.



 



Proof of <ref>:


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 ℓ and show that  

    1/N∑_ℓ=1^Nx_ℓx_ℓ^'⊗ (Σ̂_̂ℓ̂-Σ_ℓ)  0

holds element-wise, where Σ_ℓ and Σ̂_ℓ denote the covariance matrix and its estimator for subject ℓ, respectively, 1 ≤ℓ≤ N. Note that the MANCOVA setting (i.e., identical covariance matrices Σ_ℓ for all subjects within a particular group) is just a special case thereof. Due to the uniform boundedness of x_ℓ, it suffices to show that  

    1/N∑_ℓ=1^N(Σ̂_̂ℓ̂-Σ_ℓ)  0

holds element-wise. 
First, recall that u_ℓ = Y_ℓ - (x_ℓ^'⊗). Consequently, we have

    u_ℓu_ℓ^'   = (Y_ℓ - (x_ℓ^'⊗))(Y_ℓ - (x_ℓ^'⊗))^'
       = (Y_ℓ - (x_ℓ^'(X'X)^-1X'⊗)Y)(Y_ℓ - (x_ℓ^'(X'X)^-1X'⊗)Y)^'
       = (ϵ_ℓ - (x_ℓ^'(X'X)^-1X'⊗)ϵ)(ϵ_ℓ - (x_ℓ^'(X'X)^-1X'⊗)ϵ)^'
       = ϵ_ℓϵ_ℓ^' - ϵ_ℓϵ^'(X(X'X)^-1x_ℓ⊗) - (x_ℓ^'(X'X)^-1X'⊗)ϵϵ_ℓ^'
       + (x_ℓ^'(X'X)^-1X'⊗)ϵϵ^'(X(X'X)^-1x_ℓ⊗).
 

Therefore, we get (recalling that p_ℓ m = x_ℓ^'(^')^-1x_m)

    E[u_ℓu_ℓ^']    = E[ϵ_ℓϵ_ℓ^'] - ∑_m=1^NE[ϵ_ℓϵ_m^'](p_ℓ m⊗) - ∑_m=1^N(p_ℓ m⊗) E[ϵ_mϵ_ℓ^']
       + ∑_m=1^N∑_q=1^N(p_ℓ m⊗)E[ϵ_mϵ_q^'](p_ℓ q⊗),

where we have used that we can think of X as a column vector containing x_i^'⊗ in the i-th row, and ϵ = (ϵ_1^',…,ϵ_N^')^'. Now, since ϵ_i and ϵ_j are independent for i≠ j, and E[ϵ_i] = 0, the expressions from above can be further simplified to

    E[u_ℓu_ℓ^']    = E[ϵ_ℓϵ_ℓ^'] - E[ϵ_ℓϵ_ℓ^'](p_ℓℓ⊗) - (p_ℓℓ⊗) E[ϵ_ℓϵ_ℓ^']+ ∑_m=1^N(p_ℓ m⊗)E[ϵ_mϵ_m^'](p_ℓ m⊗)
       = E[ϵ_ℓϵ_ℓ^'] - E[ϵ_ℓϵ_ℓ^'](p_ℓℓ⊗) - (p_ℓℓ⊗) E[ϵ_ℓϵ_ℓ^']+ ∑_m=1^N(p_ℓ m^2⊗ E[ϵ_mϵ_m^'])
       = E[ϵ_ℓϵ_ℓ^'] - p_ℓℓE[ϵ_ℓϵ_ℓ^'] - p_ℓℓ E[ϵ_ℓϵ_ℓ^']+ ∑_m=1^Np_ℓ m^2 E[ϵ_mϵ_m^']
       = E[ϵ_ℓϵ_ℓ^'] - p_ℓℓE[ϵ_ℓϵ_ℓ^'] + ∑_m=1^Np_ℓ m^2 (E[ϵ_mϵ_m^']-E[ϵ_ℓϵ_ℓ^']),

where the latter step was due to p_ℓℓ = ∑_m=1^Np_ℓ m^2. Consequently, the k-th diagonal element of the difference E[u_ℓu_ℓ^'] -  Σ_ℓ is

    (E[u_ℓu_ℓ^'] -  Σ_ℓ)_k,k   =
    σ_ℓ k^2 - p_ℓℓσ_ℓ k^2 + ∑_m=1^Np_ℓ m^2(σ_mk^2-σ_ℓ k^2) - σ_ℓ k^2  = O(N^-1),
  
uniformly in ℓ, due to the uniform boundedness assumption on the variances, and because the uniform boundedness of the covariates and the elements of (N^-1')^-1 implies p_ℓ m = O(N^-1). 

Next, we consider Var[(u_ℓu_ℓ^' - Σ_ℓ)_k,k] = Var[(u_ℓu_ℓ^')_k,k] = Var[u_ℓ k^2]. Since u_ℓ = ϵ_ℓ - ∑_m=1^Np_ℓ mϵ_m, we have

    Var[u_ℓ k^2] = Var[ϵ_ℓ k^2 - 2∑_m=1^Np_ℓ mϵ_ℓ kϵ_mk + ∑_r,sp_ℓ rp_ℓ sϵ_rkϵ_sk].

Using Bienayme's equality yields that VarResid can be further simplified to 

    Var[u_ℓ k^2]    = Var[ϵ_ℓ k^2] 
       + 4 Var[∑_m=1^Np_ℓ mϵ_ℓ kϵ_mk]
       +Var[∑_r,sp_ℓ rp_ℓ sϵ_rkϵ_sk]
       + 2{ Cov[ϵ_ℓ k^2,-2∑_m=1^Np_ℓ mϵ_ℓ kϵ_mk]
       + Cov[ϵ_ℓ k^2,∑_r,sp_ℓ rp_ℓ sϵ_rkϵ_sk]
       + Cov[-2∑_m=1^Np_ℓ mϵ_ℓ kϵ_mk,∑_r,sp_ℓ rp_ℓ sϵ_rkϵ_sk]}

By using some algebra and the uniform boundedness conditions on E[ϵ_ℓ k^4] and the covariates, one can show that VarResid1 - VarResid6 are at least uniformly bounded for all ℓ and k (if not converging to 0). Analogously, we can show that for any ℓ_1≠ℓ_2, we have

    Cov[u_ℓ_1 k^2,u_ℓ_2 k^2] = O(N^-1),

uniformly in ℓ_1, ℓ_2 and k. So, summing up, we get that

    Var[1/N∑_ℓ = 1^Nu_ℓ k^2] = 1/N^2{∑_ℓ = 1^NVar[u_ℓ k^2] + ∑_ℓ_1 ≠ℓ_2Cov[u_ℓ_1 k^2,u_ℓ_2 k^2]}N →∞⟶ 0.

Combining convergenceE and convergenceVar and using Chebyshev's inequality, we get that

    1/N∑_ℓ = 1^N(u_ℓ k^2 - σ_ℓ k^2)  0,

for all k ∈{1,2,…,p}. 

Finally, we prove that an analogous result holds for the off-diagonal elements. Define σ_ℓ (jk):= Cov(ϵ_ℓ j,ϵ_ℓ k) for ℓ∈{1,2,…,N}, j, k ∈{1,…, p}, with j≠ k. Firstly, from errorRepresentation, one finds that

    (E[u_ℓu_ℓ^'] -  Σ_ℓ)_j,k   =
    σ_ℓ (jk) - p_ℓℓσ_ℓ (jk) + ∑_m=1^Np_ℓ m^2(σ_m(jk)-σ_ℓ (jk)) - σ_ℓ (jk)
       ≤ |p_ℓℓσ_ℓ (jk)| +∑_m=1^Np_ℓ m^2|σ_m(jk)-σ_ℓ (jk)| 
       = O(N^-1),

uniformly in ℓ,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

    Var[(u_ℓu_ℓ^'-Σ_ℓ)_j,k]   = Var[(u_ℓu_ℓ^')_j,k] = Var[u_ℓ ju_ℓ k] ≤ E[u_ℓ j^2u_ℓ k^2]
       ≤√(E[u_ℓ j^4])√(E[u_ℓ k^4])
       = √(Var[u_ℓ j^2] + E^2[u_ℓ j^2])√(Var[u_ℓ k^2] + E^2[u_ℓ k^2]).

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

Moreover,

    Cov[u_ℓ^(j)u_ℓ^(k),u_m^(j)u_m^(k)]    = E[u_ℓ^(j)u_ℓ^(k)u_m^(j)u_m^(k)] - E[u_ℓ^(j)u_ℓ^(k)]E[u_m^(j)u_m^(k)]
       ≤ |E[u_ℓ^(j)u_ℓ^(k)u_m^(j)u_m^(k)]| + O(N^-1)
       ≤√(E[u_ℓ^(j)2u_ℓ k^2])√(E[u_m^(j)2u_m^(k)2]) + O(N^-1)

(das O(N^-1) kommt daher, dass ja der zweite Teil, der subtrahiert wird, genau das Produkt der oben betrachteten Erwartungswerte des (j,k)-ten Elements ist!) Problem, das noch besteht: die oben verwendete Abschätzung für die Varianz ist wohl zu schwach!! → direkt ausrechnen (siehe Notizblaetter 1 und 2) → Cov[u_ℓ^(j)u_ℓ^(k),u_m^(j)u_m^(k)] ist mindestens von der Ordnung o(N^-1), uniformly in ℓ,m, passt also!



Moreover, one can show by tedious calculations that 

    Cov[u_ℓ_1 ju_ℓ_1 k,u_ℓ_2 ju_ℓ_2 k] = O(N^-1),

uniformly in ℓ_1,ℓ_2. 
Therefore, summing up, we get

    Var[1/N∑_ℓ = 1^Nu_ℓ ju_ℓ k] = 1/N^2{∑_ℓ = 1^NVar[u_ℓ ju_ℓ k] + ∑_ℓ_1 ≠ℓ_2Cov[u_ℓ_1 ju_ℓ_1 k,u_ℓ_2 ju_ℓ_2 k]}N →∞⟶ 0.


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 WhiteConsistency. As already mentioned, WTSConsistency 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  is weakly consistent for , because the deterministic matrix (N^-1')^-1 converges element-wise and has dimensions independent of N. This completes the proof of <ref>. 


Finally, the statements <ref> and <ref> are combined as follows: According to the assumptions from Section 2 of the manuscript, we have → (Ξ^-1⊗)Ψ(Ξ^-1⊗). Thus, <ref> implies that → (Ξ^-1⊗)Ψ(Ξ^-1⊗). Since  is of full rank for sufficiently large N, there is no rank jump in the convergence T̃T̃T̃(Ξ^-1⊗)Ψ(Ξ^-1⊗)T̃, which implies that the convergence result also holds for the Moore-Penrose inverse. Hence, applying Slutzky's theorem to <ref> yields the final result, namely that the quadratic form

    Q_N(T̃):= N ( - )'T̃^'(T̃T̃)^+T̃( - )

has, asymptotically, a central χ_f^2 distribution, where f = rank(T̃). Eventually, the statement of Theorem 1 follows by setting T̃ = diag(T,0) and applying H_0: T̃ = Tμ = 0. 











 §.§ Proof of Theorem 2




Again, the proof comprises two main steps:


  * √(N)(β̂-β) 𝒩(0,Λ), where Λ:=(Ξ^-1⊗)Ψ(Ξ^-1⊗), Ψ:= lim_N→∞(N^-1('⊗)S(⊗)), Ξ^-1:= lim_N →∞(N^-1')^-1.

  * (TND̂T)^+ (TDT)^+. 


Statement <ref> has been proven in the previous section; the proof of <ref> works analogously to the derivations in the second part of that section, because

    1/n_i∑_j=1^n_iu_iju_ij^' - Σ_i = 1/n_i∑_j=1^n_i(u_iju_ij^' - Σ_i) = 1/n_i∑_j=1^n_i(Σ̂_ij - Σ_i).

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

    N(T(μ̂-μ))^' (TDT)^+T(μ̂-μ)
       = N (T̃((μ̂-μ)^',(ν̂-ν)^')^')^'diag((TDT)^+,0)T̃((μ̂-μ)^',(ν̂-ν)^')^'
       = N(T̃(β̂-β))^'(T̃D̃T̃)^+T̃(β̂-β).



Now, according to the representation theorem and the continuous mapping theorem, statement <ref> implies that the quadratic form

    N(T(μ̂-μ))^' (TDT)^+T(μ̂-μ) = N(T̃(β̂-β))^'(T̃D̃T̃)^+T̃(β̂-β)

has, asymptotically, the same distribution as a weighted sum ∑_i=1^a+c∑_k=1^pλ_ikU_ik, where U_iki.i.d.∼χ_1^2, and the weights are the eigenvalues of the matrix of the quadratic form times the asymptotic covariance matrix of , that is, T̃(T̃D̃T̃)^+T̃·Λ. However, since T̃ = diag(T,0) and (T̃D̃T̃)^+ = diag((TDT)^+,0), that product is equal to diag(T(TDT)^+TΛ_11, 0), where Λ_11 denotes the upper-left part of the matrix Λ. Thus, the weighted some from above eventually reduces to

    U:= ∑_i=1^a∑_k=1^pλ_ikU_ik,

where the λ_ik are the eigenvalues of T(TDT)^+TΛ_11.  Note that per definition, Λ_11 is just the asymptotic covariance matrix of the vector of adjusted means. Now, using <ref> and Slutzky's theorem yields the desired result, namely that under H_0:Tμ = 0, 

    (Tμ̂)^' (TD̂T)^+Tμ̂ = N(Tμ̂)^' (TND̂T)^+Tμ̂ U.
 





Let T̃:=diag(T,0) (hypothesis matrix of general linear model representation of MANCOVA, in projector form, i.e., T = H'(HH')^-H⊗), X = (⊕_i=1^a1_n_i,Z), D̃:=diag(D,I_cp). Furthermore, let :=(μ^',ν^')^', and  be split up accordingly.
Then, we have

    N(T(μ̂-μ̂))^' (TDT)^+T(μ̂-μ̂)
       = N (T̃((μ̂-μ)^',(ν̂-ν)^')^')^'diag((TDT)^+,0)T̃((μ̂-μ)^',(ν̂-ν)^')^'
       = N(T̃(β̂-β))^'(T̃D̃T̃)^+T̃(β̂-β).

Consequently, it is sensible to consider the general case first and derive the MANCOVA as a special case, then. 




Proof of <ref>: 


Firstly, we do some simplifications: 

    √(N)(-)   = √(N)((^')^-1^'⊗)ϵ
       = √(N)((^')^-1⊗)∑_i,j(x_ij⊗)ϵ_ij


Next, we apply the multivariate Lindeberg-Feller theorem, by examining 

    η_ij:= √(N)((^')^-1⊗)(x_ij⊗)ϵ_ij, 1≤ i ≤ a, 1≤ j ≤ n_i


Firstly, E(η_ij) = 0 is obviously true. Secondly,  

    ∑_i,jCov(η_ij)    = N∑_i,j[((^')^-1⊗)(x_ij⊗)Cov(ϵ_ij)(x_ij^'⊗)((^')^-1⊗)] 
       = ((^'/N)^-1⊗)(1/N∑_i,j(x_ij⊗)Σ_i(x_ij^'⊗))((^'/N)^-1⊗),

which converges to (Ξ⊗)Δ(Ξ⊗) as N →∞, due to assumptions mancatsass2, mancatsass4 and mancatsass5. 


Thirdly, we shall examine the Lindeberg condition, i.e., 

    lim_N →∞∑_i,jE(||η_ij||^2 1{||η_ij||^2 > δ}) = 0 ∀δ >0.
 

Firstly, 

    ||η_ij||^2    = η_ij^'η_ij = N ϵ_ij^'(x_ij^'(^')^-1⊗) ((^')^-1⊗)(x_ij⊗)ϵ_ij


Let η_ij:= ϵ_ij^'(x_ij^'(^'/N)^-1⊗) ((^'/N)^-1⊗)(x_ij⊗)ϵ_ij (consequently, ||η_ij||^2 = N^-1 ||η_ij||^2). Then,

    ∑_i,jE(||η_ij||^2 1{||η_ij||^2 > δ}) = 1/N∑_i,jE(||η_ij||^2 1{||η_ij||^2 > Nδ})

Therefore, if we can show that E(||η_ij||^2 1{||η_ij||^2 > Nδ}) converges to 0, uniformly in i and j, we are done. Using Lancaster's theorem, we get

    E(||η_ij||^2)    = tr((x_ij^'⊗) ((^'/N)^-1⊗)((^'/N)^-1⊗)(x_ij⊗) Σ_i)
       = tr((x_ij⊗) Σ_i(x_ij^'⊗) ((^'/N)^-1⊗)((^'/N)^-1⊗)).

So, we have to calculate the trace of a product of matrices, which are uniformly bounded and have dimensions independent of N. Therefore, the trace is also uniformly bounded by, say, a constant C, for all i,j.  
Hence, by Markov's inequality, for any ξ>0, there is a positive constant M such that

    sup_i,j,N P(||η_ij||^2 > M) ≤sup_i,j,NE(||η_ij||^2)/M≤C/M < ξ,

that is, (||η_ij||^2) is uniformly tight. Analogously, we get that 1{||η_ij||^2>δ N} = O_P(N^-1). Consequently, 

    ||η_ij||^21{||η_ij||^2>δ N} 0.

This in turn implies that

    E(||η_ij||^21{||η_ij||^2>δ N})⟶ 0,

because 

    sup_i,jP(||η_ij||^21{||η_ij||^2>δ N} > M) ≤sup_i,jP(||η_ij||^2 > M) ≤sup_i,jE(||η_ij||^2)/M,

and by analogous arguments as above, we can show that the latter term converges to 0 as M→∞, which means that the sequence (||η_ij||^21{||η_ij||^2>δ N}) is uniformly integrable. 








 §.§ 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, 




  * √(N)β̂^∗𝒩(0,Λ) in probability, for any general linear model with coefficient vector β,where Λ:=(Ξ^-1⊗)Ψ(Ξ^-1⊗).  

  * (TND̂^∗T)^+ (TDT)^+ in probability. 




Proof of <ref>:



Firstly,

    √(N)β̂^∗   = √(N)((X'X)^-1'⊗)Y^∗ = √(N)((X'X)^-1⊗)∑_i=1^a∑_j=1^n_i(x_ij⊗) u_ijT_ij.


So, let ζ_ij:=√(N)((X'X)^-1⊗)(x_ij⊗) u_ijT_ij. Again, we use the Lindeberg-Feller theorem: E(ζ_ij|Y) = 0 is obviously fulfilled, and conditionally on the data, the ζ_ij's are independent, due to construction of the T_ij's. Next, 

    ∑_i,jCov(ζ_ij|Y)    = ∑_i,jN((X'X)^-1x_ij⊗)Cov(u_ijT_ij|Y) (x_ij^'(X'X)^-1⊗)
       = ∑_i,jN((X'X)^-1x_ij⊗)u_iju_ij^'Var(T_ij) (x_ij^'(X'X)^-1⊗)
       = N((X'X)^-1⊗)∑_i,j(x_ij⊗)u_iju_ij^' (x_ij^'⊗)((X'X)^-1⊗)
       = ((X'X/N)^-1⊗)1/N∑_i,j(x_ij⊗)u_iju_ij^' (x_ij^'⊗)((X'X/N)^-1⊗)
        (Ξ^-1⊗)Ψ(Ξ^-1⊗),


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

    1/N∑_i,j(x_ij⊗)u_iju_ij^' (x_ij^'⊗) - 1/N∑_i,j(x_ij⊗)Σ_i(x_ij^'⊗)  0.



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

    ||ζ_ij||^2 = ζ_ij^'ζ_ij   =1/Nu_ij^'T_ij(x_ij^'(X'X/N)^-1(X'X/N)^-1x_ij⊗)u_ijT_ij,

Now, analogously to the proof of Theorem 1, we let A_ij:=(x_ij^'(X'X/N)^-1(X'X/N)^-1x_ij⊗), and we denote the (k,l)-th element of A_ij by a_ij,kℓ. Consequently, we get

    ∑_i=1^a∑_j=1^n_iE[||ζ_ij||^2 1{||ζ_ij||^2 > δ}|Y] 
       = 
    1/N∑_i=1^a∑_j=1^n_iE[(T_ij^2u_ij^'A_iju_ij) 1{T_ij^2u_ij^'A_iju_ij > δ N}|Y] 
       ≤1/N∑_i,j(u_ij^'A_iju_ij)√(E[T_ij^4|Y])√(E[1{T_ij^2u_ij^'A_iju_ij > δ N}|Y]).

Now, using

    E[1{T_ij^2u_ij^'A_iju_ij > δ N}|Y] = P(T_ij^4(u_ij^'A_iju_ij)^2 > δ^2 N^2|Y) ≤E[T_ij^4|Y](u_ij^'A_iju_ij)^2/δ^2 N^2,

yields 

    1/N∑_i,j(u_ij^'A_iju_ij)√(E[T_ij^4|Y])√(E[1{T_ij^2u_ij^'A_iju_ij > δ N}|Y])
       ≤1/δ N^2∑_i,j(u_ij^'A_iju_ij)^2E[T_ij^4|Y].

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[(u_ij^'A_iju_ij)^2] is equal to the sum 

    ∑_k_1=1^p∑_ℓ_1=1^p∑_k_2=1^p∑_ℓ_2=1^pa_ij,k_1ℓ_1a_ij,k_2ℓ_2E[u_ijk_1u_ijℓ_1u_ijk_2u_ijℓ_2].

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 WBProofLindeberg2 converges to 0. Hence, 

    E[∑_i=1^a∑_j=1^n_iE[||ζ_ij||^2 1{||ζ_ij||^2 > δ}|Y]] N →∞⟶ 0.

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




Proof of <ref>:



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


  * lim_n_i →∞E[(n_i-c-1)^-1∑_j=1^n_i(u_ijk^∗ 2 - u_ijk^2)|Y] ⟶ 0 in probability, and

  * lim_n_i →∞Var[(n_i-c-1)^-1∑_j=1^n_i(u_ijk^∗ 2 - u_ijk^2)|Y] ⟶ 0 in probability,


for all i ∈{1,2,…,a} and k ∈{1,2,…,p}. Thereby, u_ijk^∗ denotes the k-th component of the wild bootstrap residual vector u_ij^∗:= Y_ij^∗ - (x_ij^'⊗)^∗ corresponding to subject j in group i. The proof works works completely analogously as in the univariate case <cit.> and is, therefore, omitted. 








 §.§ Proof of Theorem 4

Recall that the bootstrap observations are generated according to

    Y_ij^⋆i.i.d.∼𝒩(0,Σ̂_i), 1≤ i ≤ a, 1≤ j ≤ n_i,

where 
    Σ̂_i:=1/n_i-c-1∑_j=1^n_iu_iju_ij^', 1≤ i ≤ a.


Why is this sensible? - see handwritten notes in my notebook!

Then, we proceed as usual: Take least squares estimates and the White covariance matrix of the bootstrap observations. Observe, however, that we do not have any covariates to consider now, due to construction of the bootstrap observations! Consequently, it turns out that the least squares estimator is just the vector of group means (no adjusted means!). Likewise, the bootstrap version of the White covariance matrix estimator is just the empirical variance, because

    Σ̂_i^∗ = n_i^-1∑_j=1^n_iu_iju_ij^' = n_i^-1∑_j=1^n_i(Y_ij-Y̅_i.)(Y_ij-Y̅_i.)^'.

Analogously, if we take the MANCATS instead of the WTS, we just take the estimated variances in the diagonal! 

What is interesting to see here: We do not have n_i-1 in the denominator; so, we could perhaps try to make some adjustment to appropriately take the df's into account when calculating the average White estimator (also, we could do this already in the PB data generation process). Check this in simulations!


Therefore, we have that given the data,

    √(N)β̂^⋆∼𝒩(0,Λ̂_N),

where

    Λ̂_N = Cov(√(N)^⋆|Y) = 1/N∑_i,j((X'X/N)^-1x_ij⊗)Σ̂_i (x_ij^'(X'X/N)^-1⊗).


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

    1/N∑_i,j(x_ij⊗)Σ̂_i (x_ij^'⊗) - 1/N∑_i,j(x_ij⊗)Σ_i(x_ij^'⊗)  0,

which implies

    Λ̂_N  (Ξ^-1⊗)Ψ(Ξ^-1⊗).

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

    N μ̂^⋆'T(TDT)^+Tμ̂^⋆ = N ^⋆'T̃(T̃D̃T̃)^+T̃^⋆

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

What is left to show is that given the data, 

    N(D̂^⋆ - D̂) 0  in probability


Again, it suffices to show that given the data,


  * lim_n_i →∞E[(n_i-c-1)^-1∑_j=1^n_i(u_ijk^⋆ 2 - u_ijk^2)|Y] ⟶ 0 in probability, and 

  * lim_n_i →∞Var[(n_i-c-1)^-1∑_j=1^n_i(u_ijk^⋆ 2 - u_ijk^2)|Y] ⟶ 0 in probability,


for all i ∈{1,2,…,a} and k ∈{1,2,…,p}. Thereby, u_ijk^⋆ denotes the k-th component of the parametric bootstrap residual vector u_ij^⋆ = (u_ijk^⋆)_k=1^p:= Y_ij^⋆ - (x_ij^'⊗)^⋆ corresponding to subject j in group i. Firstly, observe that

    u_ijk^⋆ = Y_ijk^⋆ - ∑_i_1=1^a∑_j_1=1^n_i_1 p_ij,i_1j_1Y_i_1j_1k^⋆,


where p_ij,i_1j_1:= x_ij^'(')^-1x_i_1j_1. Consequently, we have

    u_ijk^⋆ 2 = Y_ijk^⋆ 2 - 2∑_i_1=1^a∑_j_1=1^n_i_1 p_ij,i_1j_1Y_ijk^⋆Y_i_1j_1k^⋆ + ∑_i_2,j_2,i_3,j_3 p_ij,i_2j_2p_ij,i_3j_3Y_i_2j_2k^⋆Y_i_3j_3k^⋆.


Since E[Y_ijk^⋆ 2|Y] = σ̂_ik^2, and the parametric bootstrap vectors are independent, we get

    0   ≤1/n_i-c-1∑_j=1^n_i|E[u_ijk^⋆ 2|Y]-σ̂_ik^2| 
       ≤1/n_i-c-1∑_j=1^n_i2|p_ij,ij|σ̂_ik^2 + 1/n_i-c-1∑_j=1^n_i∑_i_1=1^a∑_j_1=1^n_i_1p_ij,i_1j_1^2σ̂_i_1k^2.


Due to the fact that the elements of the hat matrix are bounded by C/N, uniformly in i,j, and σ̂_ik^2 σ_ik^2, the expression above converges to 0 in probability. So, the proof of Statement PBConsistency1 is complete. 

Next, we rewrite the conditional variance as follows:

    Var[1/n_i∑_j=1^n_i(u_ijk^⋆ 2-u_ijk^2)|Y]    = 1/n_i^2∑_j=1^n_iVar[u_ijk^⋆ 2|Y] + 1/n_i^2∑_j≠ lCov[u_ijk^⋆ 2, u_ilk^⋆ 2|Y]

Now, let

    A:=  1/n_i^2∑_j=1^n_iVar[u_ijk^⋆ 2|Y], B:=  1/n_i^2∑_j≠ℓCov[u_ijk^⋆ 2, u_iℓ k^⋆ 2|Y].

If both A and B converge to 0 in probability, we are done. Firstly, using BootResidRepresentation, Bienayme's equality and the fact that Y_ij and Y_ℓ m are independent for (i,j)≠ (l,m), we get

    A    = 1/n_i^2∑_j=1^n_iVar[u_ijk^⋆ 2|Y] 
       = 1/n_i^2∑_j=1^n_i[Var[Y_ijk^⋆ 2|Y] + 4∑_i_1,j_1p_ij,i_1j_1^2 Var[Y_ijk^⋆Y_i_1j_1k^⋆|Y].
       + ∑_i_2,j_2∑_i_3,j_3p_ij,i_2j_2^2p_ij,i_3j_3^2Var[Y_i_2j_2k^⋆Y_i_3j_3k^⋆|Y]
       + 2∑_i_2,j_2∑_i_3,j_3p_ij,i_2j_2p_ij,i_3j_3Cov[Y_i_2j_2k^⋆Y_i_3j_3k^⋆,Y_i_2j_2k^⋆Y_i_3j_3k^⋆|Y] 
       + 2 ( -2 p_ij,ij Cov[Y_ijk^⋆ 2,Y_ijk^⋆ 2|Y] + p_ij,ij^2  Cov[Y_ijk^⋆ 2,Y_ijk^⋆ 2|Y].
       ..- 4∑_i_1,j_1p_ij,i_1j_1^3 Cov[Y_ijk^⋆Y_i_1j_1k^⋆,Y_ijk^⋆Y_i_1j_1k^⋆|Y])].

Since given the data, Y_ijk^⋆ind.∼𝒩(0,σ̂_ik^2) for fixed k, we have Var[Y_ijk^⋆|Y] = σ̂_ik^2 and E[Y_ijk^⋆ 4|Y] = 3σ̂_ik^4, for i ∈{1,2,…,a}. Moreover, recall that the elements of the hat matrix are uniformly bounded by C/N. Therefore, each term of the sum converges to 0 in probability. For example,

    ∑_i_2,j_2∑_i_3,j_3p_ij,i_2j_2p_ij,i_3j_3Cov[Y_i_2j_2k^⋆Y_i_3j_3k^⋆,Y_i_2j_2k^⋆Y_i_3j_3k^⋆|Y]   = ∑_i_2,j_2∑_i_3,j_3p_ij,i_2j_2p_ij,i_3j_3Var[Y_i_2j_2k^⋆Y_i_3j_3k^⋆|Y]
       ≤3C^2/N^2∑_i_2,j_2∑_i_3,j_3σ̂_i_2k^2σ̂_i_3k^2
       = 3C^2/N^2∑_i_2=1^an_i_2σ̂_i_2k^2∑_i_3=1^an_i_3σ̂_i_3k^2
       = 3C^2 ∑_i_2=1^an_i_2/Nσ̂_i_2k^2∑_i_3=1^an_i_3/Nσ̂_i_3k^2,
 
and consequently,

    1/n_i^2∑_j=1^n_i∑_i_2,j_2∑_i_3,j_3p_ij,i_2j_2p_ij,i_3j_3Var[Y_i_2j_2k^⋆Y_i_3j_3k^⋆|Y]    ≤3C^2/n_i∑_i_2=1^an_i_2/Nσ̂_i_2k^2∑_i_3=1^an_i_3/Nσ̂_i_3k^2 0,

since n_i_2/Nσ̂_i_2k^2κ_i_2^-1σ_i_2k^2, n_i_3/Nσ̂_i_3k^2κ_i_3^-1σ_i_3k^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  0. So, summing up, the proof of Statement PBConsistency2 and, thus, of Theorem 4 is complete. 





§ GENERAL REMARKS



  * For all consistency proofs, things should go along the lines of the univariate case; so, tedious calculations, but nothing too special!

  * My σ_i estimation idea from the MANCOVA presentation slides (use group-specific instead of overall regression coefficients) is perhaps interesting, but a completely different story. So, postpone it (maybe, at some time point in the future, we could compare this estimator to the “averaged White estimator”, but at first, with respect to MSE, not with respect to hypothesis testing!).





§ IDEA OF PROOF CONCERNING THE DOMINATED CONVERGENCE PROBLEM ENCOUNTERED IN THE FIRST SECTION (WHITE WTS)


We have the situation that X_n≥ 0 (because X_n:=ϵ_i^'Γ_nϵ_i1{ϵ_i^'Γ_nϵ_i>δ N}. We have already shown that X_n 0. What is left to show is that E(X_n)⟶ 0. Firstly, the convergence of X_n in probability implies that the sequence is uniformly tight. Now, we consider an arbitrarily chosen subsequence X_n_j. Since this subsequence is uniformly tight, too, by Prohorov's theorem, there exists a sub-subsequence X_n_j_k such that X_n_j_k X. According to Portmanteau's lemma (see Van der Vaart; but no proof there, so cite another source!), we get 

    lim inf_N→∞E(X_n_j_k)≥ E(X) = 0,

where the latter is due to X = 0 a.s. (since weak limits are a.s. unique). Next, Fatou's lemma yields

    0 = E(lim sup_N→∞X_n_j_k)≥lim sup_N→∞E(X_n_j_k).

In combination with Portmanteau, we thus get 

    lim sup_N→∞E(X_n_j_k) = lim inf_N→∞E(X_n_j_k) = 0.

To sum things up, we have shown that for an arbitrarily chosen subsequence E(X_n_j) of E(X_n), one can find a further subsequence E(X_n_j_k), such that E(X_n_j_k) ⟶ 0. So, applying the subsequence principle yields the desired result. 

Warum nicht einfach von Anfang an Fatou verwenden? So geht das eigentlich ganz glatt durch! 









§ FURTHER SIMULATION RESULTS: TYPE I ERROR RATES






 



 


 






§ FURTHER SIMULATION RESULTS: ACHIEVED POWER



   


 
 













