One-Way and Two-Way ANOVA: Inferences About a Robust, Heteroscedastic Measure of Effect Size

Consider the goal of comparing independent groups. A broad issue is how many methods are needed to get a good understanding of how, and by how much, groups differ. Certainly, the best-known approach is to use some measure location such as the mean, median, 20% trimmed mean or even an M-estimator. Classic methods based on means provide some information regarding how groups differ, but it is fairly evident that this approach can miss details that are clinically important. For example, the difference between means might be relatively small when the difference between the medians is large. Of course, the reverse can happen where the difference between the medians is small but the difference between the means is large.

A limitation of methods that test hypotheses based on measures of location only is that they ignore the variation within the groups. The goal in this paper is to suggest methods that are based on a robust measure of effect size that is based in part on a robust measure of dispersion. As is illustrated later in this paper, taking into account the variation within groups can increase power substantially in some situations.

For $J$ independent groups, let ${\tau }_j$ ( $j=1$, $\dots ,J$) denote some measure of dispersion associated with the $j$th group. The better-known measures of effect size assume homoscedasticity, meaning that ${\tau }_1=\dots ={\tau }_J$ is assumed. The measures of effect size used here avoid this assumption. That is, they allow heteroscedasticity. The idea is that comparing groups based on a heteroscedastic measure of effect size provides an alternative perspective that helps provide a deeper understanding of how groups differ.

Two situations are considered. The first is where the goal is to perform all pairwise comparisons among $J$ independent groups. More formally, let ${\eta }_{jk}$ denote a measure of effect size when comparing groups $j$ and $k$. The goal is to test

and to compute a confidence interval for ${\eta }_{jk}$ for all $j<k$. Typically, as is the case here, ${\eta }_{jk}$ is based in part on the difference between two measures of location, say ${\theta }_j-{\theta }_k$. Consequently, a decision can be made about whether ${\eta }_{jk}$ is greater than or less than zero by simply testing the hypothesis $H_0:{\theta }_j={\theta }_k$. But here, ${\eta }_{jk}$ is also based on a measure of dispersion. Consequently, power when testing $H_0:{\theta }_j={\theta }_k$ can differ from the power of testing (1) simply because the method for testing (1) is sensitive to different features of the data. And even when these two approaches have similar power, the problem of computing a confidence interval for ${\eta }_{jk}$ remains.

The second goal deals with comparing the levels of the first (or second) factor of a $J$-by- $K$ ANOVA design. To elaborate a bit, momentarily focus on a 2-by- $K$ design where the goal is to compare the first and second levels of the first factor. To be a bit more precise, consider $2K$ independent random variables. Let ${\xi }_1$ be some global measure of effect size associated with the first $K$ random variables and let ${\xi }_2$ denote the measure of effect size for the remaining $K$ variables. As will be made evident, the choice for $\xi$ used here is related to the choice for $\eta$ but it differs in a fundamental way. The goal is to test

And there is the related goal of computing a confidence interval for ${\xi }_1-{\xi }_2$. The proposed method is readily extended to a $J$-by- $K$ design where the goal is to make inferences based on all pairwise levels of the first or second factor. The results here extend the results in Wilcox (2022b), which was limited to two-by-two design.

When testing (1), a seemingly natural guess is that a percentile bootstrap method will perform well in terms of controlling the Type I error probability. Simulations related to this issue are reported here. As for testing (2), a percentile bootstrap method was found to perform poorly; the actual level was estimated to be well below the nominal level. An alternative approach is proposed here and studied via simulations.

The paper is organized as follows. The next section describes the robust heteroscedastic measures of effect size that will be used. Both are based on a simple extension of measures of effect size proposed by Kulinskaya et al. (2008) and Kulinskaya and Staudte (2006). This is followed by proposed methods for testing (1) and (2). Then some simulations results are reported followed by two illustrations.

Robust, Heteroscedastic Measures of Effect Size

Momentarily consider two independent groups and let ${\overline{X}}_j$ denote the sample mean for the $j$th group. Let $n_j$ denote the corresponding sample size. Let ${\mu }_1$ and ${\mu }_2$ denote the population means and let ${\sigma }_1$ and ${\sigma }_2$ denote the population standard deviations. Certainly, one of the better-known measures of effect size is

where it is assumed that $\sigma ={\sigma }_1={\sigma }_2$. That is, homoscedasticity is assumed. Cohen’s d, (e.g., Cohen, 1988) provides an estimate of $\Delta$ that is biased. Hedge’s g (Hedges & Olkin, 1985) provides an unbiased estimator for $\Delta$. Glass et al. (1981) avoid the homoscedasticity assumption by specifying one of the groups as the control group and then use only the estimate of the variance for the control group to measure effect size. Kulinskaya et al. (2008) avoid the homoscedasticity assumption in the following manner. Note that the standard error of ${\overline{X}}_1-{\overline{X}}_2$ can be written as ${\varsigma }^2/N$ where

$N=n_1+n_2$ and $q=n_1/N$. As a result, they use

as a measure of effect size. Kulinskaya and Staudte (2006, p. 101) note that a natural generalization of $\Delta$ to the heteroscedastic case does not appear to be possible without taking into account the relative sample sizes. Under normality and when the population variances are equal, $\Delta =2{\Delta }_{KMS}$. For $J\ge 2$ independent groups, Kulinskaya and Staudte (2006) generalized this approach by estimating effect size with

where $\tilde{X}=\sum q_j{\overline{X}}_j/\sum q_j$ and $s_j^2$ is the sample variance associated with the $j$th group. Their results include a method for computing an interval estimate of the population analog of (4) assuming normality.

A parameter is said to be non-robust if a small change in a distribution has a large impact on its value. Formal mathematical methods for characterizing robust parameters are summarized by Hampel et al. (1986), Staudte and Sheather (1990) as well as Huber and Ronchetti (2009). The point here is that the population mean and variance are not robust. In practical terms, when using the measures of effect size just reviewed, even a small departure from a normal distribution can mask a large effect size among the bulk of the participants. Algina et al. (2005) illustrate this point when using $\Delta$.

Here, robust versions of (3) and (4) are used, which are based on a trimmed mean and a Winsorized variance. This mimics the basic approach used by Algina et al. (2005) to derive a robust version of $\Delta$. For notational convenience, momentarily focus on a single random sample $X_1,\dots ,X_n$ and let $X_{\left(1\right)}\le \dots \le X_{\left(n\right)}$ denote the values written in ascending order. Let $\gamma$ denote some specified constant, $0\le \gamma <0.5,$ and let $g=\left[\gamma n\right]$, where $\left[\gamma n\right]$ is the value of $\gamma n$ rounded down to the nearest integer. The sample $\gamma$ trimmed mean is computed by removing the g largest and g smallest observations and averaging the values that remain. More formally, the sample trimmed mean is

The choice $\gamma =0.2$ has been studied extensively (e.g., Wilcox, 2022a) and is used here. It has good efficiency under normality and its standard error can be substantially smaller than the standard error of the mean when dealing with heavy-tailed distributions where outliers are likely to occur.

Winsorizing $X_1,\dots ,X_n$ means that the g smallest values are set equal to $X_{\left(g+1\right)}$ and the g largest are set equal to $X_{\left(n-g\right)}$. The γ Winsorized mean, ${\overline{X}}_w$, is the mean of the Winsorized values and the Winsorized sample variance, $s_w^2$, is the usual sample variance based on the Winsorized data. When sampling from a normal distribution and when $\gamma =0.2$, $s_{wN}^2=s_w^2/0.642$ estimates the population variance. Let ${\overline{X}}_{tj}$ denote the trimmed mean for the $j$th group. Here, the measure of effect size that will be used is estimated with

where $\tilde{X}=\sum q_j{\overline{X}}_{tj}/\sum q_j$ and $s_{jwN}^2$ is the Winsorized variance of the $j$th group rescaled to estimate the variance when dealing with a normal distribution. This will be called the KMS measure of effect size henceforth. The robust version of (3) is estimated with

where

Consider the case where $J-1$ of the groups have a common population trimmed mean and the other group has a population trimmed mean that is larger than the other $J-1$ trimmed means by some specified amount. If the term $J/2$ is excluded from (5), the resulting measure of effect size decreases as $J$ increases. By including this term, $\widehat{\xi }$ remains similar to what would be obtained when $J=2$.

Simulation Results

Simulations were used to assess the small sample properties of methods M1 and M2. Data were generated from four types of distributions: normal, symmetric and heavy-tailed symmetric and relatively light-tailed, and asymmetric and relatively heavy-tailed, roughly meaning that outliers tend to be common. More specifically, data were generated from four g-and-h distributions. The four distributions used here are the standard normal ( $g=h=0$), a symmetric heavy-tailed distribution ( $h=0.2$, $g=0$), a skewed distribution with relatively light tails ( $g=1$, $h=0$), and a skewed distribution with heavy tails ( $g=1$, h = 0.2). Table 1 summarizes the skewness ( ${\kappa }_1$) and kurtosis ( ${\kappa }_2$) of these distributions. All g-and-h distributions have a median equal to zero. For $g=0$ and $h=0.2$ the variance is 2.15. For g $=1$ and h $=0$, the mean and variance are 0.648 and 4.67, respectively. For g $=1$ and h $=0.1$, the mean and variance are 0.97 and 30.6. Figure 1 shows plots of these four distributions.

The range of distributions used here is motivated by a review of several studies aimed at charactering the extent distributions differ from normality (Wilcox, 2022a, section 4.2). It is noted that the reported skewness and kurtosis for $g=h=0.2$ are estimates of the actual skewness and kurtosis based on a sample size of one million. Repeating this process yielded even larger estimates. This particular distribution might appear to be extreme. The point here is that if a method performs reasonably well when dealing with this distribution, this provides some assurance that it will perform reasonably well for distributions that are likely to be encountered.

For method M1, simulations were run for $J=2$ and 4 groups. For $J=4,$ three sample size configurations were used: (20, 20, 20, 20), (40, 40, 40, 40) and (20, 20, 40, 40). For $J=2$ the sample sizes were (20, 20), (40, 40) and (20, 40). For convenience these three sample size configurations are labeled N1, N2 and N3, respectively. The results for N1 and N2 are reported in Table 2 and are based on 1000 replications. The estimates for N3 were very similar and so for brevity are not reported. Note that the largest estimate occurs for N2, $J=4$, $g$= 0 and $h=0.2$. Increasing the common sample size to 60, the estimate is 0.040. Some additional simulations were run where the first group has a standard deviation four times larger than the other groups. Again, control over the Type I probability was similar to the results in Table 2. Bradley (1978) has suggested that as a general guide, when testing at the 0.05 level, the actual level should be between 0.025 and 0.075. All indications are that M1 satisfies this criterion. Replacing the KMS method with the robust version of $\Delta$ derived by Algina et al., 2005, gives very similar results.

There is a well-established heteroscedastic method for performing all pairwise comparisons based on trimmed means (e.g., Wilcox, 2022a, section 7.4.1). For convenience, this method is labeled T1. In some cases, T1 will have about the same amount of power as M1, but the expectation is that the choice of method might make a practical difference simply because they are sensitive to different features of the data. Consider, for example, normal distributions having variance one and suppose the mean of the first group is 0.5 while the remaining groups have a mean of zero. For N2, the probability of one or more significant results was estimated to be 0.553 for method M1 and 0.543 for method T1: Hochberg’s method was also used to control FWE when using T1. Increasing the common sample size to 60, the estimates were 0.734 and 0.693, respectively. If instead the first random variable is divided by 4 and 0.8 is added to the first and second random variables, the estimates are 0.945 and 0.865.

As for method M2, simulations were run for $J=2$ and $K=4$. Four sample size configurations were used. The first three had a common sample size of 20, 50 and 100. Now N1, N2 and N3 are used to denote these three sample configurations, respectively. The fourth, N4, consisted of a common sample size of 20 for the first four groups (corresponding to the first level of the first factor) and a common sample size of 50 for the remaining four groups. The Type I error probability was estimated with 5000 replications. Table 3 shows the estimated probability of a Type I error, $\widehat{\alpha }$, when testing at the $\alpha$= 0.05 level.

As can be seen, the estimates satisfy Bradley’s criterion in all situations except N1 and when dealing with a symmetric distribution; the estimates are less than 0.025, the lowest estimate being 0.016.

Method T1 is readily extended to a two-way ANOVA design where the goal is to perform all pairwise comparisons of the levels of the first or second factor, which is called method T2 henceforth. The details are in Wilcox (2022a, section 7.2.1; here, the R function t2way in the R package WRS was used to apply method T2). As was the case for a one-way design, M2 and T2 are sensitive to different features of the data, so which method has more power depends on the situation. A few simulations are reported here to provide some indication of the extent the choice of method matters, where again $J=2$ and $K=4$ and the goal is to compare the two levels of the first factor. Power was estimated for situations where, excluding the first group (the group corresponding to the level 1 of both factors), the remaining groups have a standard normal distribution. The first group also had a normal distribution but a mean $\mu >0$. Values for the standard deviation for the first group were taken to be $\sigma =0.5$, 1 and 2. Table 4 reports the results. As can be seen, M2 generally has the highest power. But as indicated, this is not always the case. The only point is that the power of the two methods can differ substantially. Moreover, the method that has the most power depends on how the groups differ which is not known. But perhaps the more important point is that the two methods provide different perspectives on how the groups differ.

Concluding Remarks

Steegen et al. (2016) discuss and illustrate an extremely important issue: no single method reveals everything of interest when comparing groups. Multiple methods can be required to get a deep and nuanced understanding of data. The goal here is to suggest a method for comparing groups, via a robust, heteroscedastic measure of effect size, that helps achieve this goal. Here, the illustration involving a 2-by-2 ANOVA design demonstrates this point.

Perhaps it should be stressed that it is not being suggested that methods based on measures of location only should be abandoned. Surely, they provide useful information about how groups compare. Again, the only suggestion is that comparing groups based on a robust measure of effect size that includes some measure of variation can provide insights about how groups compare.

There are several other ways for dealing with an interaction, based on a heteroscedastic measure of effect size, beyond the approach used here (Wilcox, 2022a, section 7.4.17). One approach is to use a quantile shift measure of effect size and another is to use the notion of explanatory power. Inferences based on these methods are under investigation.