A Quantile Shift Approach to Main Effects and Interactions in a 2-by-2 Design

When comparing two distributions, certainly the most common approach is to focus on a single measure of location, typically the mean or median. An alternative approach is to compare multiple quantiles with the goal of getting a more detailed understanding of where groups differ and by how much (Rousselet et al., 2017). Several methods have been derived for dealing with this issue (e.g., Doksum & Sievers, 1976; Goldman & Kaplan, 2018; Lombard, 2005; Wilcox, 1995; Wilcox et al., 2014). Extant results suggest how to generalize these methods to a between-by-between, 2-by-2 design. One goal here is to report results on two methods for controlling the family wise error rate (FWER), meaning the probability of one or more Type I errors.

Here, two distinct approaches are considered. The first defines interactions and main effects with means replaced by a collection of quantiles with an emphasis on the deciles. For example, if ${\theta }_{\mathrm{jk}}$ denotes the .2 quantile corresponding to level j of the first factor and level k of second factor, one goal is to test

which mimics the usual notion of an interaction in an obvious way. Here, however, the goal is to use multiple quantiles and to assess how well the FWER is controlled. Of course, a related issue is computing a reasonably accurate confidence interval for ${\theta }_{11}-{\theta }_{12}-{\theta }_{21}+{\theta }_{22}$. As is evident, main effects can be addressed as well. For example, one can test

for a collection of quantiles.

The second approach, when dealing with an interaction, has a certain connection to a rank-based method proposed by Patel and Hoel (1973) that in turn has a connection with the classic Wilcoxon–Mann–Whitney test. To describe the Patel–Hoel approach, let $X_{\mathrm{jk}}$ denote four independent random variables where X₁₁ and X₁₂ correspond to the first level of the first factor in a 2-by-2 design, while X₂₁ and X₂₂ correspond to the second level of the first factor. Let $p_1=P(X_{11}<X_{12})$ and $p_2=P(X_{21}<X_{22})$. The hypothesis of no interaction is

And there is the issue of computing a $1-\alpha$ confidence interval for $p_1-p_2$. Wilcox (2022, Section 7.9.2) describes a method for making inferences about this measure of effect size that performs well in simulations. For some related rank-based methods, see Gao and Alvo (2005), as well as De Neve and Thas (2017).

Let $X_{\mathrm{ijk}}$ denote a random sample from the jth level of the first factor and the kth level of the second factor ( $i=1,\dots ,n_{jk}$; j = 1, 2; k=1, 2). For convenience momentarily focus on p₁ and let $D_{ih}=X_{i11}-X_{h12}$ ( $i=1,\dots ,n_{11}$; $h=1,\dots ,n_{12}$). An estimate of p₁ is simply

where the indicator function $I(D_{ih})=1$ if $D_{ih}<0$; otherwise $I(D_{ih})=0$. The estimator ${\widehat{p}}_1$ is the estimator used by the classic Wilcoxon–Mann–Whitney test.

Let θ₁ denote the median of D and notice that

is the same as

The parameter θ₁ is defined based on Level 1 of the first factor. Let θ₂ denote the analog of θ₁ when dealing with Level 2 of the first factor. Then an analog of the Patel–Hoel interaction is ${\theta }_1-{\theta }_2$. Here, however, the goal is to consider the broader issue of comparing the deciles of these two distributions. More formally, let q₁ and q₂ denote the qth quantile of the distribution D for Level 1 of the first factor and Level 2 of first factor, respectively. The goal is to test

and to compute a $1-\alpha$ confidence interval for $q_1-q_2$ for $q=.1,.2,\dots ,.9$. A second goal is to control the FWER in a reasonably accurate manner.

Notice that when testing Equation (1), there is no distinction between ${\theta }_{11}-{\theta }_{12}={\theta }_{21}-{\theta }_{22}$ and ${\theta }_{11}-{\theta }_{22}={\theta }_{12}-{\theta }_{22}$. That is, interchanging the rows and columns does not alter the estimated effect size. However, when dealing with Equation (7), interchanging the rows and columns can yield different results.

To illustrate this point, consider, for example, a situation where for the first level of the first factor, both levels of the second factor have standard normal normal distributions, while for the second level of the first factor, the two levels of the second factor have lognormal distributions. Data were generated as just described based on sample sizes of 50 per group and the 50th quantile was estimated. The estimate of $q_1-q_2$ was .027. But interchanging the rows and columns, now the estimate was $-.088$. The code to reproduce this example can be found in the R notebook apd_ex.Rmd, which is available as part of the companion reproducibility package for this article (Wilcox & Rousselet, 2023) and on GitHub (Rousselet, 2023).

The Proposed Methods

First, consider testing Equation (1). The first issue is choosing a reasonable quantile estimator from among the many estimators that might be used. The focus here is on the estimator derived by Harrell and Davis (1982) for two fundamental reasons. The first has to do with tied values. When comparing two independent groups using the usual sample median, tied values can be accommodated using a slight generalization of a standard percentile bootstrap method. However, when comparing other quantiles using an estimator based on only one or two order statistics, such as those summarized by Hyndman and Fan (1996), this approach no longer performs in an adequate manner (e.g., Wilcox, 2022). Simulation results reported in the Simulation Results section indicate that this remains the case for the situation at hand.

In contrast to the estimators considered by Hyndman and Fan (1996), the Harrell–Davis estimator uses a weighted average of all the order statistics. Moreover, when using the Harrell–Davis estimator in conjunction with a percentile bootstrap method, all indications are that this approach does perform well when comparing two independent groups and there are tied values. An issue here is whether this continues to be the case when dealing with a two-way design. A second reason for using the Harrell–Davis estimator is that is has better efficiency under normality compared to using a weighted average of only two order statistics.

For completeness, it is noted that comparing quantiles can be accomplished using the quantile regression estimator derived by Koenker and Bassett (1978). However, this is tantamount to using one or two order statistics. When comparing the medians of two independent groups, for example, and the sample sizes are odd, in effect the usual sample median is being used. Evidently, there is no generalization of the Koenker–Bassett method that captures the spirit of the Harrell–Davis estimator.

There are several quantile estimators, in addition to the Harrell–Davis estimator, that use all of the order statistics (e.g., Liu et al., 2022; Navruz & Özdemir, 2020). A possible criticism is that these estimators, including the Harrell–Davis estimator, have a breakdown point of only 1/n. That is, the minimum number of order statistics that must be altered to make the estimate arbitrarily large is one. In practice this issue might not be a serious concern because the extreme order statistics get a relatively small weight. In a situation where the breakdown point is an issue, one possibility is to use the trimmed Harrell–Davis estimator derived by Akinshin (2022).

Let U be a random variable having a beta distribution with parameters $a=(n+1)q$ and $b=(n+1)(1-q)$. Let

The Harrell–Davis estimate of the qth quantile is

where $X_{(1)}\le \cdots \le X_{(n)}$ are the values written in ascending order. The beta weights used to calculate the deciles for a sample size of n = 50 are illustrated in Figure 1.

As an alternative to the Harrell–Davis estimator, we considered the default quantile estimator in R, called using the command quantile(x, probs = seq(0.1,0.9,0.1), Type = 7), to compute the deciles. This estimator, in addition to being widely used, relies on two order statistics (Hyndman & Fan, 1996), and could thus be more robust to outliers than the Harrell–Davis estimator in some situations.

To test Equation (1), we considered a percentile bootstrap method combined with the Harrell–Davis estimator and the quantile(Type = 7) estimator. The percentile bootstrap method begins by generating bootstrap samples by sampling with replacement $n_{\mathrm{jk}}$ values from the data associated with the jth level of the first factor and the kth level of second factor yielding $X_{ijk}^{*}(i=1,\dots ,n_{jk})$. Based on these bootstrap samples, compute the qth quantile using the Harrell–Davis estimator, or the quantile(Type = 7) estimator, yielding ${\theta }_{jk}^{*}$ followed by

Repeat this process B times yielding ${\Psi }_1^{*},\dots ,{\Psi }_B^{*}$. Let A denote the number of ${\Psi }^{*}$ values that are less than zero and let D denote the number of ${\Psi }^{*}$ values that are equal to zero. Let

A p-value, when testing Equation (1), is $2\min \{P,1-P\}$. The term D/B is important when dealing with tied values (e.g., Wilcox, 2022). To compute a $1-\alpha$ confidence interval, put ${\Psi }_1^{*},\dots ,{\Psi }_B^{*}$ in ascending order yielding ${\Psi }_{(1)}^{*}\le \cdots \le {\Psi }_{(B)}^{*}$. Let $\ell =\alpha B/2$, rounded to the nearest integer. Let $u=n-\ell$. Then a $1-\alpha$ confidence interval for Ψ is

Here, B = 2000 is used.

Note that the same bootstrap samples were used for each of the quantiles being compared. An alternative approach is to use separate bootstrap samples for each test to be performed (Wilcox, 1995). Both approaches were considered and there is no indication that separate bootstrap samples offer a practical advantage in terms of controlling the Type I error probability. Results of the simulations comparing the two bootstrap approaches are available in the R notebooks sim_fp_b1b9.Rmd and sim_fp_apd_b1b9.Rmd (Wilcox & Rousselet, 2023). Using the same bootstrap samples for all of the tests performed considerably reduces execution time, which is why it is assumed henceforth.

When performing C tests, there is the issue of controlling the FWER, meaning the probability of one or more Type I errors. Two approaches are considered here. The first is Hochberg’s (1988) improvement on the Bonferroni method. Let $p_1,\dots ,p_C$ be the p-values associated with the C tests. Put these p-values in descending order, and label the results $p_{[1]}\ge p_{[2]}\ge \cdots \ge p_{[C]}$. Set k = 0 and proceed as follows:

The second method for controlling the FWER was Benjamini and Hochberg (1995), which is aimed at controlling the false discovery rate (FDR). That is, the goal is to control the expected proportion of false positives among the rejected hypotheses. It is known that there are situations where the Benjamini–Hochberg method ensures that the false discovery rate is less than or equal to α, but it does not necessarily control the FWER (Hommel, 1988). However, the simulations described in the Simulation Results section, below, indicated that when using Hochberg's method, the FWER is always below the nominal level. Consequently, there was interest in how well the Benjamini–Hochberg method performs. It is readily verified that the Benjamini–Hochberg method always has as much or more power than Hochberg's method. Consequently, provided the Benjamini–Hochberg method controls the FWER for the situation at hand, it has a practical advantage over Hochberg's method.

The Benjamini–Hochberg method is applied simply by replacing $p_{[k]}\le \alpha /k$ in Step 1 of Hochberg’s method with

As for testing Equation (7), again both the Harrell–Davis and the quantile(Type = 7) estimators were considered, in conjunction with a percentile bootstrap method. Bootstrap samples are generated as before yielding estimates of q₁ and q₂, which are labelled $q_1^{*}$ and $q_2^{*}$. This process is repeated B times and the results are used to compute a $1-\alpha$ confidence interval for $q_1-q_2$, as well as a p-value, in essentially the same manner as done when testing Equation (1).

Simulation Results

Simulations were used to check the FWER when testing Equation (1) and Equation (7), as well as how the power of these methods compares to the classic ANOVA F test and a method for comparing 20% trimmed mean, which is described in Wilcox (2022); Section 7.4.3. Data were generated from four continuous distributions as well as three discrete distributions. The discrete distributions were a Poisson distribution having mean 9, and two beta-binomial distributions, one with parameter r = 1, the other with r = 9, and the other parameters set to s = 9, and 10 bins. The three discrete distributions were included as a partial check on the impact of tied values. The four continuous distributions were standard normal, mixed normal, lognormal and mixed lognormal. The distribution of the mixed normal is

where Φ(x) is the standard normal distribution. The mixed normal is a symmetric distribution with heavy tails, roughly meaning that outliers are likely to occur. The mixed lognormal distribution is given by Equation (13), but with Φ replaced by the lognormal distribution.

Based on over 1,500 estimates of skewness and kurtosis reported in various journal articles, Cain et al. (2017) report that 99% of the estimates were less than the skewness and kurtosis of a lognormal distribution. This suggests that if a method performs reasonably when dealing with a lognormal distribution, it is highly likely to perform reasonably well in practice. However, a possible concern is that when dealing with heavy-tailed distributions, the standard error of the usual kurtosis estimator can be quite high even when the sample size is fairly large. Moreover, the usual estimate of kurtosis can grossly under estimate the true value.

Consider, for example, a lognormal distribution, which has kurtosis 113.9. Based on a sample of 100,000, the kurtosis of the lognormal distribution was estimated and this process was repeated 1000 times. It was found that 79% of the estimates were less than the true value. The median estimate was 82. This process was repeated using the mixed lognormal distribution which is skewed and very heavy-tailed. The kurtosis of this distribution was estimated based on a sample size of one million. This process was repeated 1000 times yielding estimates ranging between 242 and 16400. The median estimate was 429. Consequently, the mixed lognormal distribution was used here as an additional check on how well the methods perform. The code for these simulations is available in the notebook kurtosis_estimation.Rmd (Wilcox & Rousselet, 2023).

Simulations were performed for sample sizes 20, 30, 40, ..., 100 per group, using 10,000 iterations. This was done using both types of quantile estimators and for the two main effects and the interaction. A p-value was computed for each of the nine quantiles to be compared. In terms of controlling the FWER, both Hochberg's method and the Benjamini–Hochberg method were considered. And as previously indicated, this was done for seven distributions. So in terms of Type I errors there is a total of 3402 results (7 distributions * 2 quantile estimators * 3 contrasts * 9 deciles * 9 sample sizes).

Finally, simulations were performed dealing with power. Various situations were considered, including shifting the four continuous distributions used to estimate Type I errors, varying the skewness of g-and-h distributions (Hoaglin, 1985), and the skewness of Poisson and beta-binomial distributions. In all situations, the distributions were adjusted to provide high power at the maximal sample size, to better differentiate the various methods. Complete details, including the code that was used, are reported in files available on figshare (Wilcox & Rousselet, 2023) and on GitHub (Rousselet, 2023). The main R packages used to perform the simulations and illustrate the results were Rcpp (Eddelbuettel & Francois, 2011), ggplot2 (Wickham, 2016) and cowplot (Wilke, 2017).

First we consider the results for testing Equation (1), then for testing Equation (7).

An Illustration

Both methods are illustrated using data dealing with perceived health (PH) among older adults (Clark et al., 2011). The first factor consists of two educational groups: those who did not complete high school and those who have some college or technical training. The other two groups are based on a measure of depressive symptoms (CESD). One group corresponds to participants with a CESD score greater than 15, which is often taken to indicate mild depression or worse. The other level consists of participants with CESD scores less than or equal to 15. The four groups are defined like this:

Perceived health results are illustrated for the four groups in Figure 11A. The figure was generated using the notebook examples.Rmd.

A 2x2 ANOVA on means returns these p-values: main effect of A (education) = 0.001; main effect of B (depression) < 0.0001; interaction = 0.09. Should we conclude that we have failed to obtain sufficient evidence about the presence of an interaction? This conclusion would be appropriate if the populations were symmetric and differed only in central tendency. However, the plot of marginals suggests differences in skewness and spread (Figure 11A). In keeping with this observation, considering the deciles reveals a more complex picture, as illustrated in Panels B-F of Figure 11, with patterns of non-uniform group differences. Panels B and C illustrate the shift functions comparing B1 and B2 at each level of A: the decile differences between two groups are plotted as a function of the deciles in one group. Panels D and E illustrate the main effects. The values along the x axis (A1 and B1) correspond to the deciles of observations pooled across groups (for instance A1 = (A1B1, A1B2). Computing the average of the deciles leads to very similar graphs. Finally, Panel F illustrates the interaction, which is highly non-linear, growing from the first decile to the median, and then decreasing and reversing sign.

The output of the decinter R function that tests Equation (1) is shown in Table 1. As can be seen, the unadjusted p-values suggest that there is an interaction for the .2-.6 quantiles. Shown are the adjusted p-values based on Hochberg’s method (using the R function p.adjust) in order to underscore the practical advantage of the Benjamini–Hochberg method. As indicated, no significant difference is found at the .05 level using Hochberg’s method. But using the Benjamini–Hochberg correction when testing Equation (1), the adjusted p-values, when dealing with the .3, .4 and .5 quantiles are all equal to .030.

When testing Equation (7), plots of the difference scores help provide perspective. Figure 12A illustrates this point. For each level of education (A1 = lower level; A2 = higher level), every participant with a lower CESD score was compared to every participant with a higher CESD score. As previously noted, when two distributions are identical, the distribution of D is symmetric about zero. Figure 12A suggests that for both groups the distribution of the difference scores is shifted to the right, but with a stronger shift for Group 2 (completed high school – lower panel in Figure 12A). A positive difference indicates higher perceived health in not depressed participants relative to depressed participants. For the second group, testing the hypothesis that the median of the difference scores is zero, the estimate is 10.5, with a [4.2, 14.7] 95% confidence interval, p = 0. In the first group the median is 4.21 [0, 8.4], p = 0.0765.

Table 2 summarizes the results when testing Equation (7). Here, the .1, .25, .5, .75 and .9 quantiles are used, which is the default for the R function iband that was used. Now the unadjusted p-values indicate an interaction in the upper tails of the two distributions. For example, the estimates of the 0.75 quantile indicate that for the first level of the first factor (did not complete high school), when comparing the not depressed group to the depressed group, there is 25% chance of getting a difference between perceived health values greater than 14.75, while for the second group there is 25% chance of getting a difference greater than 22. But using the Hochberg adjustment shown here, no significant difference is found at the .05 level and this remains the case when using the Benjamini–Hochberg method.

Concluding Remarks

The two new methods presented in this article help gain a deeper understanding of where and by how much groups differ in a 2x2 factorial design. With a sample size of at least 30 per group, all indications are that the methods perform well even when dealing with distributions that have a relatively high amount of skewness and kurtosis. These methods are implemented in the functions decinter to test Equation (1) and iband to test Equation (7). The R functions decinter and iband can be found in the file Rallfun-v41.txt, downloadable from Wilcox (2024). These functions, the R code for the simulation and the figures are available in the reproducibility package for the article (Wilcox & Rousselet, 2023) and on GitHub (Rousselet, 2023). The reproducibility package also contains Rcpp code that is much faster to execute than the base R version (see notebook examples.Rmd). By default, these functions use 2000 bootstrap samples and correct for multiple comparison using the FDR correction from Benjamini & Hochberg (1995), which outperformed the FWER correction of Hochberg (1988) in our simulations. Even with the FDR correction, the two methods presented here remain conservative, so it would be worthwhile to explore other correction strategies. We have already considered several alternative methods, but they do not improve matters (Benjamini et al., 2006; Benjamini & Yekutieli, 2001; Blanchard & Roquain, 2008). Preliminary investigations of yet other approaches to control the FWER, such as using a maximum statistic distribution (Nichols & Holmes, 2002), have not revealed any method that could significantly improve power in all situations, making recommendations difficult. While this issue requires further work, for the moment we recommend to use the FDR correction from Benjamini and Hochberg (1995) by default.

More generally, the obvious concern with comparing multiple quantiles is that power might be negatively impacted due to controlling the FWER. But power might also be negatively impacted when focusing on a single quantile simply because other quantiles are being ignored. The example in the previous section demonstrated the risk of drawing conclusions from a single measure of central tendency or quantile. It is also worth keeping in mind that there is no free lunch in inference: methods that can reveal more complex patterns in the data, such as those proposed here, necessarily require larger sample sizes to reveal where and by how much distributions differ. For applications to the deciles considered here, given the increased uncertainty associated with the estimation of the 1st and 9th deciles, we recommend sample sizes of at least 30, echoing earlier recommendations (Wilcox et al., 2014).

Several strategies are worth considering to boost power, starting with testing more specific hypotheses involving only a subset of quantiles. For instance, one could imagine a cross-validation approach in which a large dataset is split between a discovery set and a testing set. Instead of testing each quantile individually, one could also look for specific patterns across quantiles, such as stochastic dominance, which is characterised by all quantile differences having the same sign, or differences in spread, which are characterised by monotonic trends across quantile differences (Rousselet et al. 2017). Other obvious strategies are to consider different quantile estimators and bootstrap methods. Preliminary investigations suggest that using the quantile(Type = 8) estimator recommended by Hyndman and Fan (1996) improves Type I error rates and power relative to the quantile(Type = 7) estimator for instance, but is still outperformed by the Harrell–Davis estimator in all situations. As for parametric methods, Goldman and Kaplan (2018) have proposed a fast and powerful extension of the Kolmogorov–Smirnov test to compare all the quantiles of two independent groups (Kolmogorov, 1992; Stephens, 1992). However, their approach assumes no tied values, and it is unclear how it could be generalised to deal with interactions.

A bootstrap approach provides enough flexibility to deal with a variety of experimental designs. As such, this work could be extended to mixed designs with within- and between-subject factors. A covariate could be handled using different strategies, including non-parametric methods with smoothers (Wilcox, 1997; Wilcox, 2022, Chapter 12).

A referee inquired about how applied researchers might report results based on the methods studied here. We suggest reporting results as done in Table 1 and Table 2.