 # Robust test statistics for the two-way MANOVA based on the minimum covariance determinant estimator

Robust test statistics for the two-way MANOVA based on the minimum covariance determinant (MCD) estimator are proposed as alternatives to the classical Wilks' Lambda test statistics which are well known to be very sensitive to outliers as they are based on classical normal theory estimates of generalized variances. The classical Wilks' Lambda statistics are robustified by replacing the classical estimates by highly robust and efficient reweighted MCD estimates. Further, Monte Carlo simulations are used to evaluate the performance of the new test statistics under various designs by investigating their finite sample accuracy, power, and robustness against outliers. Finally, these robust test statistics are applied to a real data example.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Two-way multivariate analysis of variance (MANOVA) deals with testing the effects of the two grouping variables, usually called factors, on the measured observations as well as interaction effects between the factors. It is the direct multivariate analogue of two-way univariate ANOVA and is able to deal with possible correlations between the variables under consideration. The classical MANOVA is based on the decomposition of the total sum of squares and cross-products (SSP) matrix, and test statistics then particularly compare some scale measures between two matrices, such as the determinant, the trace, or the largest eigenvalue of the matrices.

We now consider the two-way fixed-effects layout. Let us suppose we have independent -dimensional observations generated by the model

 yijk=μ+αi+βj+γij+eijk , (1)

with , , , where is an overall effect, is the -th level effect of the first factor with levels, and is the -th level effect of the second factor with levels. We will call the first and second factor row and column factor, respectively. is the interaction effect between the -th row factor level and the -th column factor level, and is the error term which is assumed under classical assumptions to be independent and identically distributed as for all , , . Additionally, we have . We require that the number of observations in each factor combination group should be the same, so that the total SSP matrix can be suitably decomposed.

We are interested in testing the null hypotheses of all being equal to zero, all being equal to zero, and all being equal to zero. One of the most widely used tests in literature to decide upon these hypotheses is the likelihood ratio (LR) test. Here, the LR test statistic is distributed according to a Wilks’ Lambda distribution with parameters , , and (Wilks, 1932). We will give the details in the following and start with the test for interactions. The hypothesis we want to test is

 HAB:γ11=γ12=…=γrc=0 (2)

against the alternative that at least one . Further, we may also test for main effects. E.g., the hypothesis we want to test is that there are no row effects, i.e.,

 HA:α1=α2=…=αr=0 (3)

against the alternative that at least one . A similar hypothesis may be tested for column effects . More details may be found in, e.g., Mardia et al. (1979).

It is widely known that test statistics based on sample covariances such as Wilks’ Lambda are very sensitive to outliers. Therefore, inference based on such statistics can be severely distorted when the data is contaminated by outliers. A common approach to robustify statistical inference procedures is to replace the classical non-robust estimates by robust ones. The underlying idea of this plug-in principle is that robust estimators reliably estimate the parameters of the distribution of the majority of the data and that this majority follows the classical model. Hence, robust test statistics require robust estimators of scatter instead of sample covariance matrices. Many such robust estimators of scatter have been proposed in the literature. See, e.g., Hubert et al. (2008) for an overview.

In this article we use highly robust and efficient reweighted minimum covariance determinant (MCD) estimates to replace the classical non-robust covariance estimates of the Wilks’ Lambda test statistic and, in this way, extend the approach of Todorov and Filzmoser (2010) to two-way MANOVA designs.

Further approaches to robustify one-way MANOVA have already been proposed in the literature to overcome the non-robustness of the classical Wilks’ Lambda test statistic in the case of data contamination. Instead of using the original measurements Nath and Pavur (1985) suggested to apply the classical Wilks’ Lambda statistic to the ranks of the observations. Another approach proposed by Van Aelst and Willems (2011) uses multivariate S-estimators and the related MM-estimators. The null distribution of the corresponding likelihood ratio type test statistic is obtained by a fast and robust bootstrap (FRB) method. Whereas Wilcox (2012) suggested a one-way MANOVA testing procedure based on trimmed means and Winzorized covariance matrices. Xu and Cui (2008)

proposed Wilks’ Lambda testing procedures for one-way and two-way MANOVA designs based on estimators of the covariance matrix where the mean is replaced by the median, taken component-wise over the indices. Although the vector of marginal medians is easily computed it is not affine equivariant.

The rest of this article is organized as follows. In Section 2 we propose a robust version of the Wilks’ Lambda test statistic based on the reweighted MCD estimator. As the distribution of this robust Wilks’ Lambda test statistic differs from the classical one, it is necessary to find a good approximation. Section 3 summarizes the construction of an approximate distribution based on Monte Carlo simulation. In Section 4 the design of the simulation study is described and Section 5 investigates the finite sample robustness and power of the proposed tests and compare their performance with that of the classical and rank transformed Wilks’ Lambda tests. In Section 6 the proposed tests are applied to a real data example and Section 7 finishes with a discussion and directions for further research. Additional numerical results of the simulation study are given in A.

## 2 The robust Wilks’ Lambda statistics

It is well known that the Wilks’ Lambda statistic—as it is based on SSP matrices—is prone to outliers. In order to obtain a robust procedure with high breakdown point in the two-way MANOVA model we construct a robust version of the Wilks’ Lambda statistic by replacing the classical estimates of all SSP matrices by reweighted MCD estimates. The MCD estimator introduced by Rousseeuw (1985) looks for a subset of observations with lowest determinant of the sample covariance matrix. The size of the subset defines the so called trimming proportion and it is chosen between half and the full sample size depending on the desired robustness and expected number of outliers. The MCD location estimate is defined as the mean of that subset and the MCD scatter estimate is a multiple of its covariance matrix. The multiplication factor is selected so that it is consistent at the multivariate normal model and unbiased at small samples (cf. Pison et al., 2002). This estimator is not very efficient at normal models, especially if is selected so that maximal breakdown point is achieved, but in spite of its low efficiency it is the most widely used robust estimator in practice. To overcome the low efficiency of the MCD estimator, a reweighed version is used. This approach is along the same lines as proposed by Todorov and Filzmoser (2010) for the one-way MANOVA model.

We start by computing initial estimates of the factor combination group means , , , and the common covariance matrix, . The means are the robustly estimated centers of the factor combination groups based on the reweighted MCD estimator. To obtain we first center each observation by its corresponding factor combination group mean . Then, these centered observations are pooled to get a robust estimate of the common covariance matrix, , by again using the reweighted MCD estimator.

Using the obtained estimates and we can now calculate the initial robust distances (Rousseeuw and van Zomeren, 1991)

 Δ0ijk=√(yijk−ˆμ0ij⋅)⊤C−10(yijk−ˆμ0ij⋅) . (4)

With these initial robust distances we are able to define a weight for each observation , , , by setting the weight to 1 if the corresponding robust distance is less or equal to and to 0 otherwise, i.e.,

 wijk=⎧⎨⎩1if Δ0ijk≤√χ2p;0.9750otherwise . (5)

With these weights we can calculate the final reweighted estimates which are necessary for constructing the robust Wilks’ Lambda statistics :

 ˆμij⋅ = 1wij⋅n∑k=1wijkyijk ,ˆμi⋅⋅=1wi⋅⋅c∑j=1n∑k=1wijkyijk ,ˆμ⋅j⋅=1w⋅j⋅r∑i=1n∑k=1wijkyijk , ˆμ⋅⋅⋅ = 1w⋅⋅⋅r∑i=1c∑j=1n∑k=1wijkyijk , WR = r∑i=1c∑j=1n∑k=1wijk(yijk−ˆμij⋅)(yijk−ˆμij⋅)⊤ , ER = r∑i=1c∑j=1n∑k=1wijk(yijk−ˆμi⋅⋅−ˆμ⋅j⋅+ˆμ⋅⋅⋅)(yijk−ˆμi⋅⋅−ˆμ⋅j⋅+ˆμ⋅⋅⋅)⊤ ,and RR = r∑i=1wi⋅⋅(ˆμi⋅⋅−ˆμ⋅⋅⋅)(ˆμi⋅⋅−ˆμ⋅⋅⋅)⊤ ,

where

 wij⋅=n∑k=1wijk ,wi⋅⋅=c∑j=1n∑k=1wijk ,w⋅j⋅=r∑i=1n∑k=1wijk ,

and

 w⋅⋅⋅=r∑i=1c∑j=1n∑k=1wijk .

The robust test statistic for testing the equality of the interaction terms, cf. (2), is

 ΛABR=|WR||ER| , (6)

where is the determinant of . We reject the hypothesis of no interactions for small values of . Further, the robust test statistic to test the equality of the being zero irrespective of the and , cf. (3), is

 ΛAR=|WR||WR+RR| . (7)

Equality is again rejected for small values of . We note that if significant interactions exist then it does not make much sense to test for row and column effects.

Alternatively, we might decide to ignore interaction effects completely. In such a case, the two-way MANOVA design without interactions,

 yijk=μ+αi+βj++eijk , (8)

we consider the error matrix instead of and the robust test statistic for testing for row effects is then

 ΛAR=|ER||ER+RR| . (9)

In similar manner we may define robust tests for column effects . Setting all weights equal to 1, the robust test statistics coincide with the classic ones.

For computing the MCD and related estimators the FAST-MCD algorithm of Rousseeuw and Driessen (1999) will be used as implemented in the R package rrcov (cf. Todorov and Filzmoser, 2009).

## 3 Approximate distribution of Λ⋅R

The Wilks’ Lambda distribution was originally proposed by Wilks (1932). For some special cases, functions of statistics that have a Wilks’ Lambda distribution may be expressed in terms of the  distribution (cf. Anderson, 1958, pp. 195-196). For other cases, provided is large, one of the most popular approximation is Bartlett’s  approximation (Bartlett, 1938) given by

 −(ν1−12(p−ν2+1))lnΛ(p,ν1,ν2)≈χ2pν2 . (10)

Analogously to this  approximation of the classical statistic and as proposed by Todorov and Filzmoser (2010) for the one-way MANOVA we will use the following approximation for :

 L⋅R=−lnΛ⋅R≈δχ2q , (11)

and express the multiplication factor

and the degrees of freedom of the

distribution, , through the expectation and variance of , i.e.,

 E(L⋅R)=δqandVar(L⋅R)=2δ2q . (12)

Hence, we get

 δ=E(L⋅R)1qandq=2E(L⋅R)2Var% (L⋅R) . (13)

and are determined by simulation as follows.

For a given dimension , number of levels , , and sample size in each factor combination group, samples of size from the

-variate standard normal distribution will be generated, i.e.,

, , , . For each sample the robust Wilks’ Lambda statistic based on the weighted MCD will be calculated. After performing trials, the sample mean and variance of will be obtained as

 ave(L⋅R) = 1m′m′∑ℓ=1L⋅(ℓ)Rand var(L⋅R) = 1m′−1m′∑ℓ=1(L⋅(ℓ)R−ave(L⋅R))2 .

Substituting these values into Eq. (13) we obtain estimates for the constants and which in turn will be used in Eq. (11) to obtain the approximate distribution of the robust Wilks’ Lambda statistic .

The above procedure to approximate the null distribution of the robust test statistic depends on the dimension , number of levels , , and sample size in each factor combination group as well as on the considered MANOVA model and the tested hypothesis. The resulting two parameters of the approximation, and , can be stored and reused whenever a new MANOVA problem with exactly the same parameters occurs. Alternatively, we may directly perform a simulation for the problem at hand in order to compute -values of the robust tests.

Todorov and Filzmoser (2010) have shown the accuracy of this approximation.

## 4 Monte Carlo simulations

In this section a Monte Carlo study is conducted to assess the performance of the proposed test statistics—the achieved significance level and the power of the tests. Additionally, the behavior of the robust test statistics is evaluated in the presence of outliers.

We compare the results of the MCD-based tests to the classical Wilks’ Lambda tests and to an alternative approach proposed by Nath and Pavur (1985) which uses the ranks of the observations. Although this approach was only proposed for one-way MANOVA models we extend it to the two-way MANOVA cases. The Wilks’ Lambda statistics are calculated on the ranks of the original data and are referred to as rank transformed statistics, . The distributions of the test statistics are approximated by those of the normal testing theory.

We considered several dimensions , number of levels , and number of samples assuming an equal sample size in each factor combination group. The different designs are given in Table 1. These 18 designs were used to simulate data according to both considered models, the two-way MANOVA model with main effects only as well as the one with main effects and interactions.

Here, the significance level is set to 0.05.

### 4.1 Finite-sample accuracy

Under the null hypotheses of , , and being zero in both models we assume that the observations come from identical multivariate distributions, i.e., , with and . Since the considered statistics are affine equivariant, without loss of generality, we can assume to be a null vector, i.e., , and the covariance matrix to be . Thus, for each design listed in Table 1 we generate -variate vectors distributed as and calculate the classical Wilks’ Lambda test statistics, , the rank transformed ones, , and the robust versions, , based on MCD estimates. This is repeated times and the percentage of values of the test statistics above the appropriate critical value of the corresponding approximate distribution is taken as an estimate of the true significance level. The classical Wilks’ Lambda and the rank transformed Wilks’ Lambda are compared to the Bartlett approximation given by Eq. (10) while the MCD Wilks’ Lambda is compared to the approximate distribution given in Eq. (11) with parameters and estimated by Eq. (13).

### 4.2 Finite-sample power comparisons

In order to assess the power of the robust Wilks’ Lambda statistic we will generate data samples under an alternative hypothesis and will examine the frequency of incorrectly failing to reject the null hypothesis, i.e., the frequency of Type II errors. The same combinations of dimensions

, number of levels , , and sample sizes in each factor combination group as in the experiments for studying the finite-sample accuracy will be used. There are infinitely many possibilities for selecting an alternative but for the purpose of the study we will borrow an idea from experimental design: we will set the values of the levels according to the least favorable case, i.e., we will choose the values of the levels in a way for which the null hypothesis is hardest to reject resulting in the lowest possible power (cf., e.g., Rasch et al., 2012).

Moreover, we will distinguish between both considered two-way MANOVA models. For the two-way MANOVA with interactions the data samples are generated from a -dimensional normal distribution, where each factor combination group , , , has a different mean and all of them have the same covariance matrix ,

 yijk∼Np(μij,Ip) ,i=1,…,r,j=1,…,c,k=1,…,n, (14)

with

 μ11=(d/4,0,…,0)⊤ , μr1=(−d/4,0,…,0)⊤ , μ1c=(−d/4,0,…,0)⊤ , μrc=(d/4,0,…,0)⊤ , μij=(0,0,…,0)⊤ , for all other i's and j's.

The parameter takes the following values: , 0.2, 0.5, 0.7, 1.0, 1.5, 2.0. We note that according to the chosen data generating procedure the null hypotheses and are true for the simulated data.

Similar, for the two-way MANOVA without interactions the data samples are generated again from a -dimensional normal distribution, where each factor combination group , , , has a different mean and all of them have the same covariance matrix ,

 yijk∼Np(μij,Ip) ,i=1,…,r,j=1,…,c,k=1,…,n, (15)

with

 μ1j=(d/2,0,…,0)⊤ , μ2j=(−d/2,0,…,0)⊤ , μij=(0,0,…,0)⊤ , i=3,…,r,j=1,…,c .

The parameter again takes the following values: , 0.2, 0.5, 0.7, 1.0, 1.5, 2.0. We note that according to the chosen data generating procedure the null hypothesis is true for the simulated data.

Then we calculate the classical Wilks’ Lambda test statistics, , the rank transformed ones, , and the robust versions, , based on MCD estimates. This is repeated times and the rejection frequency where the statistic exceeds its appropriate critical value is the estimate of the power for the specific configuration.

### 4.3 Finite-sample robustness comparisons

Here, for each design in Table 1 data samples are generated under the null hypotheses for all factor combination groups , , , i.e., all , , , are distributed as . The factor combination group follows the contamination model

 (1−ε)Np(0,Ip)+εNp(μ∗,0.252Ip) , (16)

with and . The amount of contamination is set to 0.1 and the outlier distance takes values 2.0, 5.0, 10.0. By adding to each component of the outliers we guarantee a comparable shift for different dimensions (cf. Rocke and Woodruff, 1996).

We then calculate the classical Wilks’ Lambda test statistics, , the rank transformed ones, , and the robust versions, , based on MCD estimates. This is repeated times and the percentage of values of the test statistics above the appropriate critical value of the corresponding approximate distribution is taken as an estimate of the true significance level.

## 5 Simulation results

In this section selected results of the simulation study are presented. Here we only consider two-way MANOVA designs with and without interactions with , , , and for the classical Wilks’ Lambda test (solid line), its rank-transformed version (dotted line), and the MCD-based test (dash-dotted line). The corresponding results for other designs and dimensions are similar. Additional numerical results may be found in A. All figures related to the two-way MANOVA with interactions contain three plots that correspond to the three possible hypotheses, (left), (middle), and (right), that are able to be tested. Whereas all figures related to the two-way MANOVA without interactions contain two plots that correspond to the two possible hypotheses, (left) and (right), that may be tested.

### 5.1 Finite-sample accuracy and power comparisons

First, we consider the two-way MANOVA with interactions. The estimated power is the observed percentage of samples for which the calculated -values are below 0.05. Fig. 1 shows the estimated power of the two-way MANOVA with interactions. Its power curves are plotted as a function of and include the case . The horizontal dashed line indicates the 5% significance level. We note that according to Section 4 the null hypotheses and are true for the simulated data. Hence, in the plots on the left and in the middle of Fig. 1 it is clearly visible that these tests keep the significance level independent of the value of . Further, it can be seen in the plot on the right of Fig. 1 that the power of the robust test procedure is almost as high as that of the classical Wilks’ Lambda test and its rank-transformed version.

To give a more complete picture of how test statistics follow the approximate distribution under the null hypothesis in simulated samples, we will make use of the value plots proposed by Davidson and MacKinnon (1998). Fig. 2 shows plots of the empirical distribution functions of the -values of the two-way MANOVA with interactions. The most interesting part of the value plot is the region where the size ranges from zero to 0.2 since in practice a significance level above 20% is never used. Therefore we limited both axis to -values . We expect that the results in the value plot follow the line since the -values are distributed uniformly on if the distribution of the test statistic is correct. As it can be seen in Fig. 2 all test follow closely the line.

Furthermore, the power of different test statistics can be visually compared by computing size-power curves under fixed alternatives, as proposed by Davidson and MacKinnon (1998). As the construction of size-power plots does not require knowledge of the exact distribution of the test statistic we may use size-power plots to examine to what extent the test statistic can differentiate between the null hypothesis and the alternative. Fig. 3 shows size-power curves for three different tests setting . The size-power curve should lie above the line, the larger the distance between the curve and the line the better. Again, the most interesting part of the size-power curve is the region where the size ranges from zero to 0.2 since in practice a significance level above 20% is never used. We again note that here, according to Section 4, the null hypotheses and are true for the simulated data. Hence, in the plots on the left and in the middle of Fig. 3 it can be seen that all curves follow the line. Further, it is clearly visible in the plot on the right of Fig. 3 that all curves are far above the line. The classical Wilks’ Lambda, being the likelihood ratio statistic under normality, obviously has the highest power, closely followed by its rank-transformed version. However, the robust MCD-based statistic is almost equally powerful.

Now, we consider the two-way MANOVA without interactions. Fig. 4 shows the estimated power of the two-way MANOVA without interactions. Its power curves are plotted as a function of and include the case . The horizontal dashed line again indicates the 5% significance level. We note that according to Section 4 the null hypothesis is true for the simulated data. Hence, in the plot on the right of Fig. 4 it is clearly visible that these tests keep the significance level independent of the value of . Further, it can be seen in the plot on the left of Fig. 4 that the power of the robust test procedure is almost as high as that of the classical Wilks’ Lambda test and its rank-transformed version.

Fig. 5 shows value plots of the two-way MANOVA without interactions. As it can be seen in both plots of Fig. 5 all test follow closely the line.

Fig. 6 shows size-power curves for the three different tests setting . We again note that here, according to Section 4, the null hypothesis is true for the simulated data. Hence, in the plot on the right of Fig. 6 it can be seen that all curves follow the line. Further, it is clearly visible in the plot on the left of Fig. 6 that all curves are far above the line. Again, the classical Wilks’ Lambda, being the likelihood ratio statistic under normality, obviously has the highest power, closely followed by its rank-transformed version. However, the robust MCD-based statistic is almost equally powerful here too.

### 5.2 Finite-sample robustness comparisons

First, we again consider the two-way MANOVA with interactions. Fig. 7

shows observed Type I error rates of the two-way MANOVA with interactions in the presence of outliers. The rates are plotted as a function of the outlier distance

, including the non-contaminated case of . The horizontal dashed line indicates the 5% significance level. It is clearly visible in all three plots of Fig. 7 that the MCD-based test keeps the significance level independent of the magnitude of the outlier distance . Thus, for the MCD-based test, the Type I error rates turn out to be quite robust and close to the nominal value. The rank-transformed test and even more the classical Wilks’ Lambda test are seen to be prone to outliers. Both yield very erroneous Type I error rates.

Further, Fig. 8 shows value plots for the three different tests in the presence of outliers setting the outlier distance to 5.0. As can be seen in all plots of Fig. 8 the MCD-based test follows the line fairly well whereas the classical Wilks’ Lambda test and its rank-transformed version deviates significantly from it. Figure 8: P value plot of two-way MANOVA with interactions (outlier distance 5) Figure 9: Type I error of two-way MANOVA without interactions

Now, we consider the two-way MANOVA without interactions. Fig. 9 shows observed Type I error rates of the two-way MANOVA without interactions in the presence of outliers. The rates are plotted as a function of the outlier distance , including the non-contaminated case of . The horizontal dashed line again indicates the 5% significance level. It is clearly visible in both plots of Fig. 9 that the MCD-based test keeps the significance level independent of the magnitude of the outlier distance . Thus again, for the MCD-based test, the Type I error rates turn out to be quite robust and close to the nominal value. The rank-transformed test and even more the classical Wilks’ Lambda test are again seen to be prone to outliers. Both yield very erroneous Type I error rates.

Further, Fig. 10 shows value plots for the three different tests in the presence of outliers setting the outlier distance to 5.0. As it can be seen in both plots of Fig. 10 the MCD-based test again follows the line fairly well whereas the classical Wilks’ Lambda test and its rank-transformed version deviates significantly from it. Figure 10: P value plot of two-way MANOVA without interactions (outlier distance 5)

## 6 Real data example

We will illustrate the application of the proposed robust test statistics using waste data collected in Salzburg, Austria (Lebersorger and Salhofer, 2013). In the years 2011, 2012, and 2013 waste bins of residual municipal solid waste were analyzed in two different districts. The waste bins were randomly chosen. The numbers of waste bins are given in Table 2. Due to regulatory reasons the data were anonymized.

In a residual waste analysis the waste of each bin is divided in up to 20 different fractions and each fraction’s portion (given in percentage of the total weight of the bin) is recorded. For our example we will only consider three main fractions, namely biogenic waste, recyclables, and residual waste, and aggregate the percentages of corresponding fractions. The first six observations of the raw data are given in the following:

  district year biogenic recyclables residual
1       XY 2011   0.2073      0.2493   0.5434
2       XY 2011   0.7065      0.1194   0.1741
3       XY 2011   0.1058      0.6923   0.2019
4       XY 2011   0.2537      0.2985   0.4478
5        A 2011   0.4793      0.1047   0.4160
6        A 2011   0.0966      0.1690   0.7345
...


So, our data matrix consists of rows and columns. However, we note that as each row sums up to 1 the observations are compositions being part of the 2-dimensional simplex (cf. Aitchison, 1986). Most methods from multivariate statistics developed for real valued data are misleading or inapplicable for compositional data (cf. van den Boogaart and Tolosana-Delgado, 2013). Hence, we use the isometric log-ratio (ilr) transformation which is an isometric linear mapping between the -dimensional simplex and to obtain a 2-dimensional data matrix for further analysis. The left panel of Fig. 11 shows the scatter plot matrix of the ilr-transformed data together with histograms of each variable. The middle and left panels display grouped boxplots for the different factor combination groups. In Fig. 12 the scatter plot matrices of each factor combination group are presented separately with classical (dashed lines) and robust (solid lines) 95% confidence ellipses in the upper triangle, classical (in parentheses) and robust correlations in the lower triangle and histograms on the diagonal. Figure 11: Scatter plot matrix of the isometric log-ratio (ilr) transformed waste data (left panel) and grouped boxplots of the same data (middle and right panel) Figure 12: Isometric log-ratio (ilr) transformed waste data: Scatter plot matrices of each factor combination group separately with classical (dashed lines) and robust (solid lines) 95% confidence ellipses in the upper triangle, classical (in parentheses) and robust correlations in the lower triangle and histograms on the diagonal

We now perform a two-way MANOVA with and without interactions using the ilr-transformed waste data. All possible hypotheses (on row effects, column effects, and, if applicable, interaction effects) were tested using the classical Wilks’ Lambda test statistic (cla), the rank transformed one (rnk), and the robust version based on MCD estimates (mcd). The -values of the corresponding statistics testing for main effects only as well as for main effects and interactions are given in Tables 3 and 4, respectively.

First, we consider MANOVA tests without interactions. For the tests based on the classical Wilks’ Lambda statistic both hypotheses testing for main effects cannot be rejected at a significance level of , whereas for the rank and MCD-based tests we can reject both hypotheses. The same is true for MANOVA tests with interactions. However, the hypothesis testing for interactions is rejected for classical and rank based test, whereas for the MCD-based test it cannot be rejected. These results coincide with the practitioners’ assumptions who expected significant main effects of the factors district and year but no interaction effect.

## 7 Conclusions

This paper considered robust test statistics for the two-way MANOVA model. Robust versions of the Wilks’ Lambda statistics testing main effects as well as interactions were introduced by replacing the classical estimates for mean and covariance with robust counterparts which extends the approach of Todorov and Filzmoser (2010). Here, the MCD estimator was chosen and approximate distributions of the robust test statistics were derived by simulations. The size and power of the new proposed tests were compared with the classical and rank transformed Wilks’ Lambda tests in Monte Carlo studies. Various simulations were performed considering different dimensions, number of levels, and sample sizes of factor combination groups. Although only a selection of the results is presented in the paper these are typical representative outcomes. Therefore it can be concluded that the significance level of the robust tests is reasonably precise in case of normally distributed errors as well as in the presence of outliers. In the latter case it turns out that the actual size of the robust tests are in general much closer to the nominal size than the classical and rank transformed Wilks’ Lambda tests. Furthermore, as indicated by the size-power curves, the robust tests do not lose much power compared to the classical Wilks’ Lambda test. Further research on extending the work of Van Aelst and Willems (2011) to two-way MANOVA designs and comparing the results to the robust tests introduced here is interesting and warranted.

All computations presented in Section 5 as well as the waste data example of Section 6 were performed using the statistical software environment R (R Core Team, 2018). The functions used to perform the two-way MANOVA tests can be obtained at http://short.boku.ac.at/km53zz. Moreover, the package compositions (van den Boogaart et al., 2014) was used to ilr-transform the waste data.

Here, only test statistics in the context of two-way MANOVA were introduced. However, by following the same principle of partitioning the total SSP matrix, these test statistics can easily be extended to many more complex and higher order MANOVA designs. This is also a topic for further research.

Moreover, the tests considered here focus on robustness against data contamination, but still assume that the different groups share a common covariance structure. Hence, although no improvement of the proposed tests compared to the classical Wilks’ Lambda tests with respect to robustness against heterogeneity of the covariance structure can be expected, this deserves further research. For tests, based on the Wilks’ Lambda statistic, the robustness against heterogeneity of the covariance has been extensively investigated. Rencher (1998) gives an overview and concludes that only severe heterogeneity seriously affects Wilks’ Lambda test statistics. An alternative test statistic, the Pillai’s trace statistic, is even more stable in the presence of heterogeneity of covariances. In future research robust versions of this test statistic may be studied and compared to the robust tests introduced here.

## Acknowledgments

We are grateful to Sandra Lebersorger and the state government of Salzburg for providing the waste data set as well as to Johannes Tintner and Karl Moder for helpful discussions and comments on an earlier draft of this paper.

## Appendix A Numerical results of the simulation study

In this section additional numerical results of the simulation study are presented. The performance of the MCD-based test (mcd) is compared to the classical Wilks’ Lambda test (cla) and its rank-transformed version (rnk). Tables related to the two-way MANOVA with interactions contain only results of testing for interaction effects whereas tables related to the two-way MANOVA without interactions contain only results of testing for row effects. The significance level was set to 0.05. The results for other significance levels were found to be similar.

### a.1 Finite-sample accuracy and power comparisons

First, we consider the two-way MANOVA with interactions. In Table 5 the results of the finite-sample accuracy and power comparison are given. In the column entitled the observed Type I error rates for a nominal level of 0.05 of testing are shown. It is clearly seen that all tests are capable to keep the significance level for all investigated designs and dimensions . The remaining columns of Table 5 give the estimated power of testing for different values of . Moreover, we note that the figures printed in bold in Table 5 and the plot on the right in Fig. 1 correspond to each other.

Now, we consider the two-way MANOVA without interactions. In Table 6 the results of the finite-sample accuracy and power comparison are given as before. In the column entitled the observed Type I error rates for a nominal level of 0.05 of testing are shown. It is clearly seen that all tests are capable to keep the significance level for all investigated designs and dimensions . The remaining columns of Table 6 give the estimated power of testing for different values of . Moreover, we note that the figures printed in bold in Table 6 and the plot on the left in Fig. 4 correspond to each other.

### a.2 Finite-sample robustness comparisons

First, we again consider the two-way MANOVA with interactions. In Table 7 the results of the finite-sample robustness comparison are given. We state the observed Type I error rates for a nominal level of 0.05 of testing in the presence of out outliers. Different outlier distances were considered. It is clearly seen that the robust MCD-based test is capable to keep the significance level for all investigated designs and dimensions whereas the classical Wilks’ Lambda test and its rank-transformed version fails to keep the significance level. Moreover, we note that the figures printed in bold in Table 7 and the plot on the right in Fig. 7 correspond to each other.

Now, we consider the two-way MANOVA without interactions. In Table 8 the results of the finite-sample robustness comparison are given. We state here the observed Type I error rates for a nominal level of 0.05 of testing in the presence of out outliers. Again, different outlier distances were considered. It is clearly seen that the robust MCD-based test is capable to keep the significance level for all investigated designs and dimensions whereas the classical Wilks’ Lambda test and its rank-transformed version fails to keep the significance level. Moreover, we note that the figures printed in bold in Table 8 and the plot on the left in Fig. 9 correspond to each other.

## References

• Aitchison (1986) Aitchison, J., 1986. The Statistical Analysis of Compositional Data. Chapman & Hall, London.
• Anderson (1958) Anderson, T., 1958. An Introduction to Multivariate Statistical Analysis. Wiley, New York.
• Bartlett (1938) Bartlett, M. S., 1938. Further aspects of the theory of multiple regression. Mathematical Proceedings of the Cambridge Philosophical Society 34 (1), 33–40.
• Davidson and MacKinnon (1998) Davidson, R., MacKinnon, J. G., 1998. Graphical methods for investigating the size and power of hypothesis tests. The Manchester School 66 (1), 1–26.
• Hubert et al. (2008) Hubert, M., Rousseeuw, P. J., van Aelst, S., 2008. High-breakdown robust multivariate methods. Statistical Science 23 (1), 92–119.
• Lebersorger and Salhofer (2013) Lebersorger, S., Salhofer, S., 2013. Evaluation der Auswirkungen des neuen Recyclinghofes in Gemeinde A.
• Mardia et al. (1979) Mardia, K., Kent, J., Bibby, J., 1979. Multivariate analysis. Academic Press.
• Nath and Pavur (1985) Nath, R., Pavur, R., 1985. A new statistic in the one-way multivariate analysis of variance. Computational Statistics & Data Analysis 2 (4), 297–315.
• Pison et al. (2002) Pison, G., Van Aelst, S., Willems, G., Apr 2002. Small sample corrections for lts and mcd. Metrika 55 (1), 111–123.
• R Core Team (2018) R Core Team, 2018. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
URL https://www.R-project.org/
• Rasch et al. (2012) Rasch, D., Spangl, B., Wang, M., 2012. Minimal experimental size in the three way anova cross classification model with approximate F-tests. Communications in Statistics - Simulation and Computation 41 (7), 1120–1130.
• Rencher (1998) Rencher, A., 1998. Multivariate Statistical Inference and Applications. Wiley, New York.
• Rocke and Woodruff (1996) Rocke, D. M., Woodruff, D. L., 1996. Identification of outliers in multivariate data. Journal of the American Statistical Association 91 (435), 1047–1061.
• Rousseeuw (1985) Rousseeuw, P. J., 1985. Multivariate estimation with high-breakdown point. In: Grossmann, W., Pflug, G., Vincze, I., Wertz, W. (Eds.), Mathematical Statistics and Applications. Reidel, Dortrecht, pp. 283–297.
• Rousseeuw and Driessen (1999) Rousseeuw, P. J., Driessen, K. V., 1999. A fast algorithm for the minimum covariance determinant estimator. Technometrics 41 (3), 212–223.
• Rousseeuw and van Zomeren (1991) Rousseeuw, P. J., van Zomeren, B. C., 1991. Robust distances: Simulations and cutoff values. In: Stahel, W., Weisberg, S. (Eds.), Directions in Robust Statistics and Diagnostics, Part II. Springer, New York, pp. 195–203.
• Todorov and Filzmoser (2009) Todorov, V., Filzmoser, P., 2009. An object-oriented framework for robust multivariate analysis. Journal of Statistical Software 32 (3), 1–47.
• Todorov and Filzmoser (2010) Todorov, V., Filzmoser, P., 2010. Robust statistic for the one-way manova. Computational Statistics & Data Analysis 54 (1), 37–48.
• Van Aelst and Willems (2011) Van Aelst, S., Willems, G., 2011. Robust and efficient one-way manova tests. Journal of the American Statistical Association 106 (494), 706–718.
• van den Boogaart et al. (2014) van den Boogaart, K. G., Tolosana, R., Bren, M., 2014. compositions: Compositional Data Analysis. R package version 1.40-1.
URL https://CRAN.R-project.org/package=compositions
• van den Boogaart and Tolosana-Delgado (2013) van den Boogaart, K. G., Tolosana-Delgado, R., 2013. Analyzing Compositional Data with R. Springer, Heidelberg.
• Wilcox (2012) Wilcox, R., 2012. Introduction to Robust Estimation and Hypothesis Testing. Academic Press, Boston.
• Wilks (1932) Wilks, S. S., 1932. Certain generalizations in the analysis of variance. Biometrika 24 (3/4), 471–494.
• Xu and Cui (2008) Xu, J., Cui, X., 2008. Robustified manova with applications in detecting differentially expressed genes from oligonucleotide arrays. Bioinformatics 24 (8), 1056–1062.