2.1 Permuting independent measured observations: the Pitman permutation test

Appropriate permutation strategies can be inferred from knowledge of the experimental design and what the null hypothesis would mean for counterfactual sample draws. Here, the null hypothesis of no treatment effect would mean that average prices in every session of the experiment were independent draws from a common average-price distribution. By the same logic, each of the observed price draws would be just as likely to have been assigned to the cash-rich treatment as the asset-rich treatment, so every permutation of the data between these treatments would have an equal ex ante probability of having been observed. With samples of 3 observations in the cash-rich treatment and 4 observations in the asset-rich treatment, there are $\left(\frac{7}{3}\right)=35$ equally probable ways that these data could have been assigned to the two treatments if the null hypothesis were true. The null distribution of the test statistic can be constructed by computing the value of the test statistic under each of these 35 permutations of the sample data, a subset of which are reproduced in Table 3.

Table 3

Computing the null distribution of the test statistic

Permutation index (i)	Cash-rich session-average price (x)	Asset-rich session-average price (y)	Test statistic $T_i={\overline{x}}_i-{\overline{y}}_i$
1	3.03 2.55 3.06	3.43 3.73 3.97 3.32	− 0.73
2	3.03 3.32 2.55	3.43 3.73 3.97 3.06	− 0.58
3	3.32 2.55 3.06	3.43 3.73 3.97 3.03	− 0.56
4	3.43 3.03 2.55	3.73 3.97 3.32 3.06	− 0.52
5	3.43 2.55 3.06	3.73 3.97 3.03 3.32	− 0.50
…	…	…	…
31	3.43 3.97 3.32	3.73 3.03 2.55 3.06	0.48
32	3.73 3.97 3.03	3.43 3.32 2.55 3.06	0.49
33	3.73 3.97 3.06	3.43 3.03 3.32 2.55	0.50
34	3.73 3.97 3.32	3.43 3.03 2.55 3.06	0.66
35*	3.43 3.73 3.97	3.03 3.32 2.55 3.06	0.72

*Observed samples

After permutation values of the test statistic are computed, hypothesis testing becomes straightforward. In conventional hypothesis testing, the p-value for a two-sided test represents the probability of drawing a value of the test statistic, T, at least as extreme as the observed value, $T_{\mathrm{obs}}$ , if the null hypothesis were true:

(2) $\text{two-sided}p\text{-value}=P_{H_0}\left(\left(T\right)\ge \left(T_{\mathrm{obs}}\right)\right)$

This probability is usually calculated by assuming that the test statistic follows a known distribution. Instead of assuming a distribution, permutation testing constructs the empirical null distribution of the test statistic from the sample data. Here, $T_{\mathrm{obs}}=0.72$ . There are only two ways that these data could be rearranged to result in a test statistic as extreme or more extreme than $T_{\mathrm{obs}}$ . These are the permutations in the first and last rows of Table 3. Since the null hypothesis implies that all permutations are equally probable, the p-value for a two-sided test of the null hypothesis is $2/35\approx 0.057.$ More generally, the two-sided p-value for this permutation test is the proportion of all test-statistic values greater than or equal to the observed value among all $\left(\frac{m}{n}\right)$ ways of permuting the data between the two samples:

(3) $\text{Pitman}\text{permutation}\text{test},\text{two-sided}p\text{-value}=\frac{{\sum }_{i=1}^{\left(\frac{m}{n}\right)}1\left(\left(T_i\right)\ge \left(T_{\mathrm{obs}}\right)\right)}{\left(\frac{m}{n}\right)}$

where $T_i$ is the value of the test statistic for the ith permutation and $1\left(·\right)$ is the indicator function with a value 1 if its argument is logically true and 0 otherwise.

One-sided versions of this permutation test are computed similarly; the only differences are that the numerator consists of signed values of the test statistic and the inequality identifies more extreme observations in the hypothesized direction. A one-sided test for the Caginalp et al. (Reference Caginalp, Porter and Smith1998) data could be defended because the experiment was designed to evaluate the excess cash hypothesis that high cash-to-asset ratios would tend to increase share prices. The p-value for this one-sided test is the probability of observing a signed value of the test statistic greater than or equal to the observed value under the null. Reviewing Table 3, only the observed permutation meets this criterion, so the one-sided p-value is $1/35\approx 0.029$ .

In our experience, students often find the methodology of permutation testing more intuitive than normal-theory testing. Some find it so appealing that they never look back. Even so, it is instructive to consider some of the similarities and differences between the permutation method and the familiar approach of Student’s two-sample t test.

First, unlike the t or normal distributions, the null distribution of a permutation test can be highly discrete. In the above one-sided example, the p-value of 0.029 is the strongest possible rejection of the null hypothesis for these sample sizes. Intuitively, no configuration of the data can be more extreme than the case where all 3 observations in the cash-rich sample are greater than all 4 observations in the asset-rich sample, so the p-value for this permutation test could never fall below 1/35 no matter how extreme the difference between the samples.

Second, while having a less discrete null distribution can allow Student’s t-test to reach lower p-values (0.025 for the two-sided t test versus 0.057 in the above two-sided permutation test), it does so at the cost of assuming a specific distribution for the test statistic. This distributional assumption is not innocuous. Inaccurate distributional assumptions can invalidate a parametric test. The permutation test imposes no distributional assumptions in the sense of requiring the sample data to come from a specific population distribution. While the permutation test is not equally powerful for all possible distributions, its validity is not based upon distributional assumptions.

Third, the similarity of this permutation test and Student’s t test is not superficial. Both tests are based on permutation-equivalent test statistics (functions of the difference of sample averages). As sample sizes become large, the efficiency of the permutation test converges with that of the t-test (Hoeffding, Reference Hoeffding1952; Miller, Reference Miller1997: p. 54). Thus, the two-independent-sample permutation test can be seen as a distributionally robust version of the t test.

Fourth, while the intuition of permutation testing fits most naturally to hypothesis testing, the methodology can also be used to generate a form of confidence interval (Manly, Reference Manly2007: pp. 18–20; Miller, Reference Miller1997: pp. 24, 28). To generate this interval, let $y\left(\theta \right)=y+\theta$ be the observed values of sample $y$ , shifted up or down by the value $\theta$ . By inspection or numeric search, find the smallest and largest values of $\theta$ such that the permutation test would fail to reject the null of no difference between $x$ and $y\left(\theta \right)$ at level $\alpha$ in a two-sided test. These end points constitute a $100\left(1-\alpha \right)\%$ confidence interval: a range of shift values that could not be rejected at the $\alpha$ level in a two-sided test. Applied to the Caginalp et al. (Reference Caginalp, Porter and Smith1998) data, for example, this process yields a 90% confidence interval of (0.244, 1.182). This interval is wider than the (0.266, 1.174) confidence interval produced by the standard t test but, unlike the t test, does not depend on distributional assumptions.

If everything in this section seems intuitive, one might wonder why this type of permutation testing is not more common in the experimental literature. Indeed, the Pitman subcaption of this section is a nod to one of the earliest proponents of this form of permutation testing (Pitman Reference Pitman1937a, Reference Pitman1937b, Reference Pitman1938), yet few experimenters would recognize that name today.Footnote ⁴ The explanation is that, while the theory of permutation-based inference has been understood for more than 75 years, computing power has only recently made this type of case-by-case construction of null distributions a convenient option (Berry et al., Reference Berry, Johnston and Mielke2019: ch. 2). Recent clearance of computational barriers might raise some doubts about the ease or practicality using permutation tests in published work, but none of them withstand scrutiny.

For example, one might worry that, for large sample sizes, the computational burdens of permutation testing would still become prohibitive even with modern computers. This does happen but it is a problem of little practical importance. For one thing, numerical simulation methods carry the theory of permutation testing past the computational horizon. An approximate permutation test simply replaces exhaustive permutation with many random reshuffles of the observed data.Footnote ⁵ For another, normal theory and asymptotic testing become defensible as sample sizes grow, meaning the practical need for permutation testing fades just as its computational burdens begin to turn troublesome.

In a similar vein, one might worry that considerable scripting experience would be needed to perform tests that are not already standard routines in common statistical software. Some comfort with scripting is indeed helpful, especially in stratified permutation testing (discussed in Sect. 6) or when designing tests for complicated randomization routines (Heß, Reference Heß2017). But common statistical programs provide standard routines for permuting sample data; in some cases, available permutation testing libraries make scripting entirely unnecessary. We collect relevant software and routines for permutation testing in online appendixes.

Finally, a pragmatic researcher might hesitate to use an uncommon statistical test for fear that editors and reviewers would question it. We have encountered no evidence of this ourselves. The theory of permutation testing is old and well established. Methods papers have specifically recommended it for experimental work (Moir, Reference Moir1998). And examples of the use of permutation testing in published work, while not common, are available to cite as precedent if needed. Examples include Sherstyuk (Reference Sherstyuk1999), Abbink (Reference Abbink2004), Orzen (Reference Orzen2008), Anderson et al. (Reference Anderson, DiTraglia and Gerlach2011), Sieberg et al. (Reference Sieberg, Clark, Holt, Nordstrom and Reed2013), Nosenzo et al. (Reference Nosenzo, Quercia and Sefton2015), Rosokha and Younge (Reference Rosokha and Younge2020), Erkal et al. (Reference Erkal, Gangadharan and Koh2020), Kujansuua et al. (Reference Kujansuua and Schram2021), Stephenson and Brown (Reference Stephenson and Brown2021), and Schram et al. (Reference Schram, Zheng and Zhuravleva2022).

With no remaining obstacles to excuse its unfamiliarity, the Pitman permutation test deserves a more prominent role in the experimenter’s toolkit than it commands today.

2.2 Permuting independent ranked observations: the Mann–Whitney test

While the Pitman permutation test is unfamiliar, related tests devised by Wilcoxon (Reference Wilcoxon1945) and Mann and Whitney (Reference Mann and Whitney1947) are commonly employed in the experimental economics literature. The Wilcoxon and Mann–Whitney approaches describe different but equivalent tests and are sometimes referenced jointly as the Wilcoxon-Mann–Whitney test. For ease of exposition, we refer to both tests as the Mann–Whitney test, which helps to distinguish this permutation strategy from the Wilcoxon Signed Rank test (discussed in relation to matched-pair samples in Sect. 3.2). Just as lack of familiarity with the Pitman permutation test owes to historic and now outdated computational difficulties, the popularity of the Mann–Whitney test owes mainly to inertia from computational shortcuts of little modern importance.

The Mann–Whitney procedure presented in most non-parametric books involves ranking all sample data (both samples combined) and replacing each observation with its ordinal rank in the combined sample. Both the Wilcoxon and Mann–Whitney versions of the test then compute special test statistics with computationally convenient null distributions (Gibbons & Chakraborti, Reference Gibbons and Chakraborti2003; Miller, Reference Miller1997; Siegel, Reference Siegel1956). The specific definition of the test statistic is not important for present purposes. The thing to note is that the approach is equivalent to running a permutation test on the rank-transformed sample data.

Specifically, the Mann–Whitney test is a two-independent-sample permutation test (covered in Sect. 2.1) in which the data being permuted are not the measured observations but the ordinal ranks of those observations in the combined sample (Siegel, Reference Siegel1956: p. 155). Instead of permuting observed average prices in the Caginalp et al. (Reference Caginalp, Porter and Smith1998) experiment, the Mann–Whitney test would permute the ranked values of each session-average price in the combined sample. The cash-rich average prices of 3.43, 3.73, and 3.97 were the three highest prices observed in this experiment, with ranks 5, 6, and 7 respectively; the asset-rich average prices of 3.03, 3.32, 2.55, and 3.06 were the four lowest prices, with ranks of 2, 4, 1, and 3 respectively. The Mann–Whitney test is equivalent to running a Pitman permutation test on these rank-transformed observations.

It happens to make no difference, in this example, whether one works in ranks or level data. Since the ranked versions of the samples, x = {5, 6, 7} and y = {2, 4, 1, 3}, are more extreme than any other possible permutation, the one-sided p-value is still $1/35\approx 0.029$ and the two-sided p-value is still $2/35\approx 0.057$ . This is a special case, however, and permutation tests based on ranks are not generally the same as those based on observed values.

To illustrate the potential difference between these tests, consider an experiment reported by Bohr et al. (Reference Bohr, Holt and Schubert2019), involving asset market performance with saving decisions over a simulated lifetime. Subjects, in this experiment, were permitted to buy and sell asset “shares” that paid dividends each period. The dividends and interest paid on cash induced a flat fundamental share value (present value of future dividends) of $20 per share. In six sessions of the experiment (the “private-savings” treatment), subjects traded assets while also deciding how much of their incomes to save for low-income “retirement” years. In another 6 sessions (the “government-savings” treatment), a fixed portion of each subject’s income was withheld by the “government” for retirement years. As a result of the difference in savings policy, subjects in the private-savings treatment carried cash amounts that were about twice as large as those in the government-savings treatment. Peak price data are provided for both treatments in Table 4, below.

Table 4

Peak price data and ranks for asset shares

Treatment	Session peak prices	Mean
Private savings	42	36	53	61.5	38.5	70	50.2
Government savings	42.5	21.25	30	26	43	38	33.5

Bohr et al. (Reference Bohr, Holt and Schubert2019): selected data from Table 1

While major price bubbles (above the $20 fundamental value) were observed in most sessions of the experiment, an interesting research question is whether peak asset prices were greater under the cash-rich private-savings treatment than under the government-savings treatment.Footnote ⁶ The data are qualitatively consistent with this hypothesis but there is overlap between the samples. When permuting the measured observations under the Pitman test, there are $\left(\frac{12}{6}\right)=924$ possible permutations of the 12 observations across treatments, 36 of which yield a treatment difference at least as extreme as the observed value, for a two-sided p-value of about 0.039. In contrast, when permuting rank-transformed observations under the Mann–Whitney test, there are 86 permutations in which the treatment difference (in average ranks) is at least as extreme as the observed value, yielding a two-sided p-value of about 0.093. The rank-based Mann–Whitney test barely supports rejection of the null hypothesis at the 10% level, while the Pitman test’s greater use of the sample information allows for rejection at the 5% level.Footnote ⁷

Since the Pitman permutation test is sensitive to the magnitudes of observation differences, not just their rank comparisons, one might think that the p-values produced by Pitman tests will always be lower than those of equivalent Mann–Whitney tests. This intuition is incorrect. To illustrate, consider an experiment reported by Holt and Smith (Reference Holt and Smith2022) in which subjects competed for a fixed monetary prize by expending lab dollars on “effort” that has a specified cost per unit. One treatment contrast compared average total expenditures on “rent-seeking” activity with low effort cost (77, 83, 130, 81) against expenditures with high effort cost (132, 111, 93, 87).Footnote ⁸ Of the $\left(\frac{8}{4}\right)=70$ possible ways of permutating the sample observations, 28 yield a treatment difference more extreme than what was observed, making the p-value of a two-sided Pitman test $28/70=0.4$ . After rank-transforming the sample observations, however, only 14 of 70 permutations yield an average difference more extreme than the observed value, meaning the p-value of a two-sided Mann–Whitney test is $14/70=0.2$ .

Another way that the Pitman and Mann–Whitney tests often differ is in confidence intervals. Recall that the Pitman test produced a 90% confidence interval of (0.244, 1.182) when assessing the difference between cash-rich and asset-rich average prices in the Caginalp et al. (Reference Caginalp, Porter and Smith1998) experiment. Confidence intervals can also be generated for rank-transformed tests (Bauer, Reference Bauer1972). Computing the end points of an interval by searching for the largest and smallest shift-terms where a two-sided Mann–Whitney test would not reject the null hypothesis at the 90% level yields a 90% confidence interval of $\left(0.111,1.420\right)$ .

If the only difference between the Mann–Whitney and Pitman permutation tests is that the Mann–Whitney tests drops information when sample data are measured as interval values, then why is the Mann–Whitney test so popular? It once enjoyed the important advantage of having a null distribution that could be pre-computed and printed in critical-value tables in statistics textbooks. But modern computing power makes this all but irrelevant. The rank-based test is still appropriate when the sample data are ordinal as measured.Footnote ⁹ And, because the rank conversion suppresses distortions caused by outliers, the Mann–Whitney test may also be preferable to the Pitman test when experimenters are concerned that outliers might be introduced by measurement errors, procedural issues, or other artificial influences.

In most cases involving interval measurements, however, the Pitman test presents the more compelling option. The Pitman test uses more of the sample information and is more sensitive to differences between measured observations. In a Monte Carlo comparison of the Pitman test against Student’s t-test and the Mann–Whitney test, Moir (Reference Moir1998) finds the Pitman test to match or outperform the more familiar tests in most applications.Footnote ¹⁰ Where data are interval as measured (prices, auction revenues, market efficiencies), and where policy interest concerns the magnitudes of treatment effects, it is hard to justify passing over the Pitman test in favor of an alternative that is less sensitive to these aspects of the data.

3.1 Permuting pairs of measured observations: the Fisher permutation test

As before, the appropriate permutation strategy for constructing the null distribution of a test statistic can be inferred from knowledge of the experimental design and what the null hypothesis would mean for counterfactual sample draws. Here, the null hypothesis of no treatment effect implies that the vector-difference of matched pairs should have zero mean, suggesting the use of the average difference as an appropriate test statistic:Footnote ¹⁴

(4) $T=\overline{d}$

As noted above, the null hypothesis implies that the difference distribution should be symmetric around zero. This motivates a simple permutation strategy for constructing the null distribution of the test statistic.

Since every difference is equally likely to be either positively or negatively signed under the null hypothesis, the null distribution of the test statistic can be computed by recalculating the value of the test statistic under all possible permutations of the signs of the sample differences. Equivalently, and perhaps more intuitively, since $F_x$ and $F_y$ are equal under the null hypothesis, the observed values in every matched pair are equally likely to have been assigned to opposite treatments under the null—which corresponds to simply swapping the sign of their difference. This strategy of permuting signs is often attributed to Fisher (Reference Fisher1935). We refer to it as the Fisher permutation test.

For a sample of $n$ matched pairs, there are $2^n$ ways that the signs of all matched pairs could be permuted under the null. Taking the p-value to be the probability of seeing a value of the test statistic, T, as or more extreme than the observed value, $T_{\mathrm{obs}}$ , a permutation p-value is as follows:

(5) $\text{Fisher}\text{permutation}\text{test},\text{two-sided}p\text{-value}=\frac{{\sum }_{i=1}^{2^n}1\left(\left(T_i\right)\ge \left(T_{\mathrm{obs}}\right)\right)}{2^n}$

where $1\left(·\right)$ is again the indicator function.

The Davis and Holt experiment provides a concrete context for illustration. Recall that the first question is whether increased market power confers pricing power in this setting. For the 6 matched pairs in sessions S1–S6 of Table 5, the average price is higher for the market power treatment in every case. That is, the sample-difference vector has only positive signs. Since there are $2^6=64$ possible permutations of the signs of these matched pairs under the null, only two of which result in a test statistic as or more extreme than the observed test statistic, the two-sided p-value is $2/64\approx 0.031$ . Since economic theory predicts a positive sign for this treatment effect, it is arguably more appropriate to use a one-sided test in this setting. The one-sided p-value corresponding to a test of the alternative hypothesis that prices are higher with greater market power seeks only those permutations with test statistic values greater than or equal to the observed value. Since no alternative permutation of the signs of these data provides an average value larger than what was observed, the one-sided p-value is about $0.016$ . This example is, however, misleading in its simplicity. Since all differences were of the same sign, the magnitudes of differences did not need to be considered in calculating the $1/64$ probability for each tail. Next, we consider a more realistic example with some overlap.

The second question is whether a pure numbers effect gives smaller numbers of competitors greater pricing power, even when holding technological market power constant. For the 6 matched pairs in sessions S7–S12 of Table 5, the average price is greater under the 3-seller treatment in all but one case. Focusing on the one-sided test, 2 possible permutations yield an average difference greater than or equal to the observed value. The p-value of a one-sided test is thus $2/64\approx 0.031$ ; the corresponding p-value for a two-sided test would be $4/64\approx 0.062$ .

Note that the previous conclusion would not have been reached if testing had not accounted for matched pairs in this design. If the 3-seller power and 5-seller power samples on the right side of Table 5 had been treated as independent samples, a Pitman permutation test would give a one-sided p-value of about 0.273, which would not support rejection of the null hypothesis of no numbers effect. The intuitive explanation for this is that groups of subjects can differ in competitiveness. The most collusive outcome observed in any 3-seller market was produced by the same group that produced the most collusive outcome in any 5-seller market. Using each group as its own control helps to mitigate the effects of subject heterogeneity and thus helps to reveal treatment effects that might otherwise be difficult to distinguish from noise in the data.

Like the Pitman test, the Fisher permutation test has appeared in published articles. Examples include Baillon, Schlesinger, and van de Kuilen (Reference Baillon, Schlesinger and van de Kuilen2018) and Bao et al. (Reference Bao, Kalaycı, Leibbrandt and Oyarzun2020).

3.2 Permuting ranked observations: the Wilcoxon signed-rank test

Just as the ability to perform permutation tests on measured data is often overlooked in experimental analysis of independent samples, so is it overlooked in the matched-pairs context. By far, the most common small-sample test in the matched-pairs setting is the Wilcoxon (Reference Wilcoxon1945) signed-rank test. The name of this test reflects the peculiar transformation that it applies to the measured data. All matched-pair differences are ranked from smallest to largest in absolute value, and these ranks are then reassigned the signs of the original difference data. To illustrate, the vector of sample treatment differences $d=\left(-6,4,0,-3\right)$ becomes $S{R}_d=\left(-4,3,1,-2\right)$ .

As was the case for the Mann–Whitney test, the original motivation for using the signed-rank transformation was primarily to enable reliance on a test statistic for which a precompiled null distribution could be provided in printed form. The specifics of the relevant test statistic and its distribution are interesting but are not belabored here. Accessible presentations are provided by Wilcoxon (Reference Wilcoxon1945) and most introductory texts on nonparametric statistics (Gibbons & Chakraborti, Reference Gibbons and Chakraborti2003; Miller, Reference Miller1997; Siegel, Reference Siegel1956). It is sufficient to note that the Wilcoxon signed-rank test is a Fisher permutation test conducted on the signed ranks of the sample difference vector.

How do the Fisher permutation test and the Wilcoxon signed-rank test compare? For the Davis and Holt experiment, both yield the same p-values. Fisher’s permutation test has strong intuitive appeal as the permutation analog of Student’s one-sample t-test (Miller, Reference Miller1997). Monte Carlo evidence demonstrates the superior power of the Fisher permutation test over the Wilcoxon signed rank test for a variety of sample sizes and distributions (Kempthorne & Doerfler, Reference Kempthorne and Doerfler1969). This is intuitive, as a Fisher permutation test run on interval data uses more sample information than does a Wilcoxon signed-rank test. In our opinion, the Fisher permutation test should be the experimenter’s default choice when differences between observed measurements reflect important outcomes—with the Wilcoxon test reserved for situations in which data are ordinal as measured or when outliers are a concern.

To illustrate the last point, suppose that the first session of the 5 seller no-market-power treatment had yielded an observed price of 700, rather than 329. Since the Fisher permutation test is based on measured values instead of ranks, the magnitude of this observation affects the test statistic quite a bit. When considering the one-sided alternative, there are now 26 possible sign permutations that yield average differences greater than or equal to the observed difference, for a p-value of $26/64\approx 0.406$ . The magnitude of the high-price observation has less effect when converted to a signed rank in the Wilcoxon test. When considering the one-sided alternative, there are only 14 possible sign permutations that would yield an average difference of signed ranks greater than or equal to the observed value, for a p-value of $14/64\approx 0.219$ . The Fisher permutation test’s sensitivity to extreme observations is a strength of the test when these observations reflect underlying economic conditions but is a weakness when these observations are artificial outliers.

4.1 Permuting measured observations: the F test

The test of multiple unordered treatment effects is a simple generalization of the test of locational difference between two independent samples. In the two-sample setting, samples $x$ and $y$ are compared to see whether their means are statistically distinguishable. In the more general setting, $k>2$ samples, $x_1,\dots ,x_k,$ of size $n_1,\dots ,n_k$ are compared simultaneously to see whether any of their means are statistically distinguishable. If experiment sessions are run at k different universities, for example, a researcher might want to determine whether any subject pool differences are revealed. In statistics textbooks, this procedure falls under the heading of one-factor analysis of variance (ANOVA) testing, and statistical inference would most likely be based on the F statistic and its comparison to the F distribution:

(6) $F=\frac{{\left(k-1\right)}^{-1}{\sum }_{j=1}^k{n}_j{\left({\overline{x}}_j-\overline{x}\right)}^2}{{\left(N-k\right)}^{-1}{\sum }_{j=1}^k{\sum }_{i=1}^{n_j}{\left(x_{i,j}-{\overline{x}}_j\right)}^2}$

where $x_{i,j}$ denotes the ith observation in sample $x_j$ , ${\overline{x}}_j$ and $n_j$ denote the mean and number of observations in sample $x_j$ , $\overline{x}$ denotes the grand mean when all samples are pooled together, and N is the total number of observations. This test is known to be robust against non-normality (Miller, Reference Miller1997: 80–81), but experimental data may push the boundaries of what could plausibly justify the assumption of a specific, parametric null distribution.

Applications of permutation testing in one-way ANOVA go back nearly as far as the use of permutation testing in the two-sample setting (Pitman, Reference Pitman1938). The reason is that the two-sample permutation process generalizes in an intuitive way to the higher order setting. To illustrate, consider an experiment that Holt and Sprott (Reference Holt and Sprott2022) use to evaluate different methods for auctioning off-shore wind energy leases. One treatment implements a clock auction used by the US Department of Interior. The initial clock price is set low to generate excess demand for leases, and the clock price is raised in a sequence of rounds until there is no excess demand, which determines the uniform sale price. A second treatment models the multi-round sealed-bid procedure used in the UK. In that auction, the high bidder in the first round receives any leases they bid for at their own bid price, and these leases are removed from the auction. The high bidder in the second round then gets any leases they bid for at their bid price, and so on until all leases are sold.Footnote ¹⁵ Finally, a third treatment implements a single-round sealed-bid auction. After the bids are submitted, the auction is cleared by going down the bid list from high to low and accepting bids until no leases remain unsold. Winning bidders pay their own bid amounts.

Table 6 shows average earnings over each of 4 sessions of each treatment. Since the costs of wind-energy generation are difficult to estimate precisely in advance, the experiment was structured so that each bidder received independent noisy signals of the underlying lease value. This common-value setting entails a risk that the winning bidder might have bid high because they overestimated the lease value. This well-known “winner’s curse” effect is clearly indicated by the negative average earnings for the UK auction treatment in the top row.

Table 6

Wind energy auction earnings

Treatment	Session-average earnings				Average
UK (right to choose auctions)	− 0.61	− 3.64	0.82	1.39	− 0.51
US (clock auctions)	3.84	1.04	3.17	3.97	3.00
Single-round pay-as-bid auction	0.22	− 0.24	0.42	1.13	0.38

Holt & Sprott (Reference Holt and Sprott2022): data on file with authors

To illustrate how a permutation test could be applied to these data, consider a test of the null hypothesis of no difference between any of the auction formats, against the alternative of a locational difference between at least two of the mechanisms. The observed data yield an F statistic of 5.54. All that remains is to construct the null distribution and to see how the observed value of the test statistic compares with it.

As before, the null distribution can be derived from knowledge of the experimental design and the implications of the null hypothesis. If the null of no treatment effect were true, then all observations in the experiment would be equally likely to have been assigned to any of the treatments. If there are $k>2$ treatments with $n=\sum n_k$ total observations, then there are $\left(\frac{n}{n_1}\right)$ ways that observations could have been assigned to the first treatment. For each of these, there are $\left(\frac{n-n_1}{n_2}\right)$ ways that observations could have been assigned to the second treatment. For each of these, there are $\left(\frac{n-n_1-n_2}{n_3}\right)$ ways that observations could have been assigned to the third treatment. And so on, for a total of $n!/\left(n_1!\times \dots \times n_k!\right)$ equally likely permutations of the sample data under the null. The p-value for a permutation-based test of the null hypothesis (equality of all samples) comes from computing the F statistic for every permutation of the data and counting the proportion of these F statistics that are greater than or equal to the observed value:

(7) $\text{permutation}F\text{test},p\text{-value}=\frac{{\sum }_{i=1}^{N}1\left(F_i\ge F_{\mathrm{obs}}\right)}{N}$

where $N=n!/\left(n_1!\times \dots \times n_k!\right)$ . For the wind auction experiment, three groups of four observations yields $N=12!/\left(4!\times 4!\times 4!\right)=34,650$ possible permutations of the observed data. As noted above, the observed value of the test statistic is 5.54. Of the 34,650 possible permutations of the sample data, 642 yield F statistics greater than or equal to 5.54, so the p-value for a permutation F test of the null hypothesis is about 0.019. This p-value is smaller than the p-value of $4/70\approx 0.057$ that comes from pairwise comparison of even the highly contrasting US and UK treatments using the Pitman test discussed in Sect. 2.1. Of course, specific pairwise comparisons would be relevant for the actual choice of an auction format.

4.2 Permuting ranked observations: the Kruskal–Wallis test

The above description of the permutation F test differs from the familiar Kruskal–Wallis test—the current default choice of most experimenters when seeking a nonparametric test of locational difference among k samples (Kruskal & Wallis, Reference Kruskal and Wallis1952). Like the Wilcoxon and Mann–Whitney tests, the Kruskal–Wallis test statistic is a function of rank-transformed observations:

(8) $H=\left(N-1\right)\frac{{\sum }_{j=1}^k{n}_j{\left({\overline{r}}_j-\overline{r}\right)}^2}{{\sum }_{j=1}^k{\sum }_{i=1}^{n_j}{\left(r_{i,j}-{\overline{r}}_j\right)}^2}$

where $r_{i,j}$ denotes the rank (among all observed values) of the ith observation in sample $x_j$ , ${\overline{r}}_j$ and $n_j$ denote the mean rank and number of observations in sample $x_j$ , $\overline{r}$ denotes the average of all ranks, and N is the total number of observations.

How does the Kruskal–Wallis Test compare to the permutation F test? Consistent with our discussion in the two-sample context, rank-based tests make sense when the sample data are ordinal as measured, or when the experimenter has reason to be concerned about serious outliers not connected to fundamentals of the subject being investigated. There is, however, an important difference between the Kruskal–Wallis test and the previously discussed two-sample rank-based tests. Unlike its two-sample analogs, the Kruskal–Wallis test does not have an easily computed null distribution. Implementations in common statistics software substitute approximate null distributions, but these are known to be poor approximations for small sample sizes (Meyer & Seaman, Reference Meyer and Seaman2013). The exact null distribution can be constructed by permutation, but then the Kruskal–Wallis test lacks any computational advantage over the permutation F test. For these reasons, rote reliance on the Kruskal–Wallis test is not advisable.

To illustrate how these considerations play out in an actual example, return to the Wind Auction experiment and the null hypothesis of no difference in treatment effect between any of the auction mechanisms. Recall that the permutation F test rejected the null hypothesis with a p-value of 0.019. An exact p-value for the Kruskal–Wallis test can be computed by following the same process, only replacing the F statistic with the H statistic. For the Wind Auction data, the observed value of the test statistic is $H=5.69$ . Of the 34,650 possible permutations of the sample data, 1,686 yield H statistics greater than or equal to 5.69, making the exact p-value of the rank-based Kruskal–Wallis test $1,686/34,650\approx 0.049$ .

5.1 Permuting ranked observations: the Jonckheere–Terpstra test

To ground discussion of ordered alternative hypotheses, consider the classic Smith (Reference Smith1964) experiment comparing three variations of the double-auction trading institution. Based on observations of a prior pilot experiment, Smith conjectured that prices would be lower when sellers competed against and undercut each other and that prices would be higher when buyers bid against each other to compete for purchases. To test this, Smith ran two sessions for each of three treatments—each session lasting 5 rounds—with infra-marginal buyer values and seller costs creating identical supply–demand arrays symmetrically configured around a unique equilibrium price in every round. In a sellers-offer treatment, sellers could make price offers during a trading period, whereas buyers could only observe offers and decide whether to accept them. The situation was reversed in a buyers-bid treatment: only buyers could make bids to purchase, and sellers were constrained to observe and decide whether to accept bids. Finally, the double-auction treatment is symmetric: buyers could make bids, sellers could make offers, and either side could accept a bid or offer. Sessions in groups A and B were run with 20 and 24 traders, respectively. Here, we rely on the theoretical invariance of equilibrium outcomes to group size to pool the two session-average prices in each column of Table 7 (yielding 2 observations per treatment), but we will return to this group-size nuisance variable in Sect. 6.2.Footnote ¹⁶ Table 7 shows average prices for the final two trading periods in each session. The null hypothesis is that the trading rules have no effect on price levels. The alternative hypothesis is that prices will be lowest in the sellers-offer treatment and highest in the buyers-bid treatment.

Table 7

Session-average contract price by trading condition

Group	Sellers offer	Double auction	Buyers bid
Group A (20 traders)	208	213	217
Group B (24 traders)	195	209	213

Smith (Reference Smith1964): selected data from Table 3

The most commonly used nonparametric test when testing for this type of ordered treatment effects is the Jonckheere–Terpstra test (Jonckheere, Reference Jonckheere1954; Terpstra, Reference Terpstra1952. When treatment categories are ordered so that the predicted effect increases from left to right (as in Table 7), the test statistic J is the sum of all “binary wins” in the predicted direction. In other words, J is the total number of larger observations in columns to the right of each observation:

(9) $J=\sum _{s=1}^{k-1}\sum _{t=s+1}^k\sum _{i=1}^{n_s}\sum _{j=1}^{n_t}1\left(x_{i,s}<x_{j,t}\right)$

where $x_{i,s}$ is the ith observation in the sample data from ordered treatment s, $x_{j,t}$ is the jth observation in the sample data from ordered treatment $t>s$ , $1\left(·\right)$ is the indicator function, and $n_s$ , and $n_t$ are respective sample sizes. Notice that the first summation, indexed by s, is over all columns except the final column, k. The second summation, indexed by t, is over all columns to the right of s, up to and including column k. The final two summations are over all pairs of observations in columns s and t, which are used to obtain a “less-than” count via the indicator function. Here, the average contract price of 208 for the Group B session in the Seller Offers treatment is smaller than all 4 numbers in columns to its right, so the first term in the sum for J would be a 4. Similarly, the contract price of 195 is smaller than all 4 numbers to its right, so the second term in the sum for J is also 4. In the middle column, the contract price of 209 is smaller than 2 numbers to its right, so the third term in the sum is 2, but the contract price of 213 in the middle column is only smaller than 1 number to its right, so the final term in the sum is 1. The sum of these 4 terms yields the test statistic, $J_{\mathrm{obs}}=4+4+2+1=11$ .

For large sample sizes, approximate null distributions for J are available. For small samples like this, the researcher is left to search for a precomputed null distribution or to generate one via permutation. By the logic used in Sect. 4.1, the null hypothesis of no treatment effect implies that there are $n!/\left(n_1!\times \dots \times n_k!\right)$ equally likely permutations of the observed sample data. The permutation p-value for this test is thus computed the same way as the F test described in Eq. (7) but with $J_i$ and $J_{\mathrm{obs}}$ substituted in place of $F_i$ and $F_{\mathrm{obs}}$ . Here, the observed value of the test statistic is $J_{\mathrm{obs}}=11$ and, of the $N=6!/{\left(2!\right)}^3=90$ possible permutations of the data, only 2 yield a value of J as or more extreme than the observed value of 11, so the p-value for the Jonckheere–Terpstra test is $2/90\approx 0.022$ .

The Jonckheere–Terpstra test is an attractive option when sample data are measured as ordinal values. It is also a test that appears with some regularity in the literature. Smith (Reference Smith1964) used a Jonckheere–Terpstra test in his analysis of these data.Footnote ¹⁷ Other recent examples include Gürerk and Selten (Reference Gürerk and Selten2012) and Conrads et al. (Reference Conrads, Irlenbusch, Reggiani, Rilke and Sliwka2016). When sample data are interval as measured, however, the Jonckheere–Terpstra test operates like other rank-based tests in discarding sample information. The obvious question is whether a different test could be constructed that does not involve this loss of information.

5.2 Permuting measured observations: a directional difference test

How might a magnitude-sensitive version of the Jonckheere–Terpstra test be constructed? A simple approach would be to replace the “binary win” count with a sum of differences.Footnote ¹⁸ With the treatment vectors of observations still listed left to right in increasing order of predicted effect, let the test statistic D be the sum of all differences between each observation and all observations in columns to the left in the ordered array:

(10) $D=\sum _{s=1}^{k-1}\sum _{t=s+1}^k\sum _{i=1}^{n_s}\sum _{j=1}^{n_t}\left(x_{j,t}-x_{i,s}\right)$

where all terms are defined as above. For the Smith (Reference Smith1964) data, the observed test statistic is $D=108$ , which is equaled or exceeded by only 2 test statistics in the 90 possible permutations of the data, for a p-value of $2/90\approx 0.022$ . Here, the Jonckheere–Terpstra test (based on binary win comparisons) and the Direction Difference test [based on calculated directional differences in (10)] yield the same p-values. This is a special case, however, and these tests will not generally yield the same results when applied to the same data. Rather than belabor this difference, which we have already seen in other applications, we note two points.

First, as in other comparisons, the availability of a test that operates on measured observations demands an explanation for preferring a rank-based test when sample data are interval as measured. Unless a precomputed null distribution is available, the computational cost of calculating exact p-values for the Jonckheere–Terpstra test is the same as for the Directional Difference test. We conjecture that the Jonckheere–Terpstra test could have superior properties when the samples contain outliers, though we are aware of no Monte Carlo analysis or proof to this effect. Absent special considerations, the Directional Difference test provides a more natural and intuitive basis for statistical inference.

Second, there is potentially great research value in designing experiments with the type of treatment variation illustrated in this example. Suppose that, instead of collecting two data points from each of three treatments, the researchers had followed the more conservative approach of collecting three data points for each of two treatments. The strongest possible rejection by a permutation test in the two-treatment context would be at a p-value of $1/20=0.05$ for a one-tailed test, assuming no reversals between the treatments. In Table 7, with the same number of data points spread across three treatments, the p-value is smaller even in the presence of a tied “reversal” for the 213 price in each of the double-auction and buyers-bid treatments.

Several generalizations of the Jonckheere–Terpstra test have been used for physical systems in which too much of a treatment may have a negative effect. The “umbrella test” is for the case where the alternative hypothesis involves a “hill-shaped” data pattern as treatment intensity is increased (Mack & Wolfe, Reference Mack and Wolfe1981). This test essentially combines two directional tests: one for data to the left of the mode, and one for data to the right. Analogous test statistics could be devised for the directional difference test. This illustrates an intriguing property of permutation testing. As Pearson (Reference Pearson1937) observed at the dawn of this methodology, the permutation method’s decoupling of computation of the null distribution from choice of test statistic frees researchers to select whatever test statistic is most sensitive to a hypothesized relationship.

6.1 Stratified permutation tests for k = 2 treatments with nuisance variables

A common problem in experimental data analysis is dealing with procedural differences in data groupings that are unrelated to the difference of interest. These procedural groupings are essentially nuisance variables: secondary treatments that ideally would be held constant when evaluating the effects of the primary treatments. For example, suppose an experiment involves two treatments, each applied to subjects from two different pools. If the experimenter has reason to believe that these subject pools are interchangeable, then observations can be pooled by treatment and a simple two-sample test can be employed. If subject pools cannot be assumed to be interchangeable, however, then the need to perform separate tests for each subject pool can present both narrative and statistical challenges (a point we illustrate below).

Stratified permutation testing may be a more attractive option. To explain what we mean, let the primary treatments be indexed by j, and let the secondary groupings be indexed by g. Thus, $x_{\mathrm{ijg}}$ denotes experimental observation i taken when treatment j is applied to subjects from group g. The idea behind stratified permutation testing is to construct the null sampling distribution by permuting the primary treatment labels (j) within groups but not between groups. This procedure captures the null hypothesis—that observed measurements are equally likely to be seen under any treatment—without imposing the additional assumption that observed measurements are equally likely to be seen under any secondary grouping.

As a concrete example, consider the two-treatment asset market experiment reported by Holt et al. (Reference Holt, Porzio and Song2017). One treatment, applied to 14 sessions, involved a 25-period trading sequence; the other treatment, applied to 10 sessions, involved a 15-period trading sequence. These markets were blocked on gender: half of the sessions in each treatment were female-only and half were male-only. In all treatments, the fundamental (present) value of asset shares was constant, at $28, for all trading periods. Price bubbles, with peaks well above $28, were observed in all sessions. Table 8 shows peak asset prices for female-only sessions (top row) and male-only sessions (bottom row). Shorter sessions afforded subjects less opportunity to accumulate large cash balances, resulting in smaller cash-asset-value ratios. This difference in cash-asset values motivated a question whether the peak asset prices were also lower in the shorter markets. For purposes of testing the null hypothesis of no trading-length effect, trading-length is the control variable of primary interest and gender groupings are a nuisance variable.Footnote ¹⁹

Table 8

Peak prices by market with gender sorting

Market pool	25-period markets	15-period markets	Mean
Female only	87 95 61 177.5 75.5 37 152	66 36 58 64 42	79.3
Male only	55 48 68 85 65 56.5 50	50 70 45 43 53	57.4
Mean	79.5	52.7	68.3

Holt et al. (Reference Holt, Porzio and Song2017): selected data from Tables 1 and 5

In order to evaluate the effect of the number market periods, while controlling for gender, a stratified permutation test permutes peak price observations across session-length treatments (the columns of Table 8), but not across gender-groups (the rows of Table 8). A one-sided permutation test of the null of no treatment effect against the alternative of higher price peaks in longer markets can be based on the observed difference in the bottom row of the table, $T_{\mathrm{obs}}=79.5-52.7=26.8.$ The null distribution of this test statistic is constructed by computing the test statistic over each of the constrained set of permutations in which observations are moved between market-length treatments but not across gender labels. There are $\left(\frac{12}{5}\right)=792$ ways that treatment labels could be assigned to the numbers in the top row. For each of these top row permutations, there are another 792 ways that treatment labels could be reassigned in the bottom row. Thus, there are 627,264 total permutations to consider. Of these, 6259 involve a treatment effect greater than or equal to the observed value, yielding a one-sided p-value of about 0.01. The exact same approach could be used with the Mann–Whitney test, in which case 11,577 permutations are greater than or equal to the observed test statistic, for a p-value of about 0.018.

How does this approach compare to conducting two separate tests, one for each gender group? Conducting separate Pitman tests for each gender group yields one-sided p-values of 0.12 for the male-only group, and 0.028 for the female-only group. This example illustrates the previously discussed problems with the multiple-comparison approach. Can the experimenter conclude, in this multiple-comparison exercise, that the null hypothesis of no treatment effect is rejected? Moreover, if the experimenter really was seeking to test every combination of treatment effect and gender grouping simultaneously, then the testing procedure should be constructed in a way that controls the family-wise error rate of these tests (the probability of falsely rejecting at least one null hypothesis among the two tests). While a detailed treatment of multiple-comparison adjustments is beyond the scope of this discussion, it is helpful to note that Bonferroni-adjusted p-values are 0.24 for the male-only group and 0.056 for the female-only group, substantially greater than the 0.01 value of the stratified test.Footnote ²⁰

Stratified permutation testing provides expositional simplicity and statistical power when nuisance variables are discrete. It does so by tailoring the permutation strategy to the underlying randomization of the experiment’s design. It bears emphasis that the null hypothesis for the stratified permutation test—that observations are drawn from the same distribution within strata—does not restrict observations to share a common distribution across strata. Distributions should be the same within strata, though, apart from potential differences in location. Treatment effects should also be the same across strata; otherwise, separate tests would be appropriate.Footnote ²¹

Before moving on, note that this stratification process generalizes easily to situations with multiple nuisance variables. A simple illustration is provided by Comeig et al. (Reference Comeig, Holt and Jaramillo2022), who report an experiment designed to study risk appetite as a function of framing (“downside risk” vs “upside risk”), payoff scale, and subject gender. The experiment involved 256 subjects, half male and half female, each tasked with making a choice between a risky lottery and a safe lottery.Footnote ²² Half of the subjects were exposed to treatments in which the risky option was presented as downside risk (a small probability of a low payoff); the other half were presented the risky option as upside risk (a small probability of a large payoff). If small probabilities are “over-weighted,” as prospect theory predicts, then subjects would tend to shy away from the downside risk of a low payoff and be attracted to the upside risk of a low probability of a high payoff. In every pairwise choice, both experimental lotteries were constructed so that the expected payoff from the risky lottery exceeded the payoff from the safe lottery by the same fixed amount. The standard deviations of the two lotteries were also the same. Finally, these treatments were blocked on payoff scale, with half of subjects presented payoffs five times larger than the other half. Of the 32 male and 32 female subjects exposed to each combination of risk profile and payoff scale, the number of subjects choosing the risky option is presented in Table 9, below.

Table 9

Risky option choice proportions by treatment

Gender	Payoff scale	Downside risk (%)	Upside risk (%)
Male	1×	25/32 = 78	30/32 = 94
Male	5×	17/32 = 53	28/32 = 88
Female	1×	19/32 = 59	29/32 = 91
Female	5×	4/32 = 13	27/32 = 84

Comeig et al. (Reference Comeig, Holt and Jaramillo2022): Table 2

While these data may be used to explore various hypotheses, perhaps the most interesting prediction is that subjects should be more willing to take upside risks than downside risks, even when the expected payoff and standard deviation of the safe and risky options are the same. In testing this hypothesis, both gender and payoff scale are nuisance variables.Footnote ²³

A stratified permutation test of the effect of risk type on lottery choice would involve permuting risk-type labels across each of the 256 lottery choices, subject to the constraint that labels are not reassigned across any of the strata in the different rows of the table. With two crossed nuisance variables, this equates to preserving four separate strata during the permutation process: male/1×, male/5×, female/1×, and female/5×. Otherwise, the procedure is the same as with one nuisance parameter.

Given the large number of possible permutations for a sample of this size, an exact permutation test would be computationally costly. An approximate permutation test can be conducted by randomly sampling many possible permutations, subject to stratification constraints, and comparing the test statistic values of these random draws against the observed value. Here, pooling across gender and payoff scale, 49 more subjects selected the risky lottery in the case of upside risk than in the case of downside risk. Over 999,999 random permutations of risk-type labels within the respective strata, none yielded a pooled difference between the upside and downside risk treatments as extreme as the observed difference, implying a p-value of less than 0.001. Similar tests, omitted here, could be used to evaluate the effects of payoff scale or subject gender.

One of the most attractive properties of this stratified permutation approach is the intuitive nature of randomizing across only the dimension of the data at focus. Indeed, the stratification tests we consider in this section are really just generalizations of the matched-pairs permutation strategy. In the matched-pairs context, each pair of observations is treated as its own stratum. In the more general stratified permutation testing context, multiple observations may fall within each stratum, and different strata may have different numbers of observations. The intuition behind stratified permutation also generalizes to other influences that might be accounted for in the permutation process (Ehlert et al., Reference Ehlert, Kindschi, Algesheimer and Rauhut2020; Heß, Reference Heß2017). As a general strategy for conducting statistical inference in the presence of nuisance variables, stratified permutation testing strikes an attractive balance of analytical flexibility and ease of presentation.

6.2 Stratified permutation tests for k > 2 treatments with nuisance variables

The idea of stratified permutation testing generalizes in an intuitive way to higher order comparisons with $k>2$ treatments. As a concrete example, return to the Smith (Reference Smith1964) experiment, selected data from which is reproduced in Table 7. When using the Jonckheere–Terpstra test to look for an ordered treatment effect in Sect. 5.1, we relied on economic theory to justify pooling results across sessions with different group sizes: 20 traders in Group A sessions versus 24 traders in Group B. Instead of combining the data for the two rows of each treatment column as was done in the prior section, an alternative would be to treat the presence of extra-marginal units as a nuisance variable.

In other words, the rows of Table 7 can be seen as strata with a one average price observation for each of three treatments. In general, a stratified Jonckheere–Terpstra test would involve permuting observations across treatments but not across strata. Here, since there are no ties or reversals for treatment comparisons within a row of the table, the p-value can be computed analytically. There are $3!$ ways to permute the observations in each row, so there are ${\left(3!\right)}^2=36$ possible permutations of the sample data that do not move observations across the separate strata. Since the most extreme outcome in the conjectured direction is what is observed, the p-value of a stratified Jonckheere–Terpstra test is $1/36\approx 0.028$ , only a little larger than the 0.022 p-value obtained when ignoring the nuisance variable. In contrast, running separate tests for each row would limit the researcher to minimum p-values of $1/3\approx 0.333$ for each test.

The similarity of the previous example to the analysis of matched-pairs data is no coincidence. The matched-pair design (in which each subject is exposed to two treatment conditions) is a special case of a more general design in which a subject, or group of subjects, is exposed to $k>2$ treatment conditions. In the language of classical statistics, the resulting data invite two-way analysis of variance: the experimenter may be primarily interested in the k treatment effects, but analysis also needs to account for dependence relationships among the multiple measurements taken from each subject or group.

A clean example of this k-treatment data structure is provided by Ma et al. (Reference Ma, Noussair and Renneboog2022), who report a laboratory experiment designed to measure valuation of paintings with 6 different color configurations. In a laboratory setting, the authors showed 465 unique subjects sets of 6 constructed Rothko-like paintings, each with a different primary color (blue, red, green, yellow, etc.). Subjects were given an opportunity to purchase each painting under a bidding process that incentivized private-value revelation.Footnote ²⁴ Results indicated substantial differences in average valuation by color, with bids for red and blue paintings exceeding the average bid by about 17–19 percent.

With each subject viewing and bidding on multiple paintings, subject heterogeneity could be pronounced in this design. A subject who particularly liked art, or who happened to need a painting for decorative purposes, might bid higher on all 6 colors than another subject. If this experiment had involved only a few subjects, then a stratified permutation test for differences in valuation by color could be based upon the two-way ANOVA F statistic. For the null hypothesis of no treatment effect between any of the 6 colors, the observed F statistic could be compared to a null distribution generated by permuting color labels within each subject’s bids (controlling heterogeneity by constraining permutations never to cross subject-strata in the sample data). With a sample consisting of the price bids of 465 subjects, however, standard parametric tests are a better option.

Interestingly, a rank-based analogue of the permutation ANOVA test just described was long ago proposed by economist Milton Friedman (Reference Friedman1937, Reference Friedman1939, Reference Friedman1940), who later won a Nobel Prize in Economics for other works. Friedman argued that normality is “likely to be the exception rather than the rule” when studying socio-economic data and devised a rank-based test because computations using ranks were “less arduous… requiring but a fraction of the time” when working with “large scale collections of social and economic data” in the years leading up to World War II. (Friedman, Reference Friedman1937: p. 675). Times have changed and, like other rank-based tests, the advantages of the Friedman test are now limited. Importantly, the difficulty of computing exact p-values under the Friedman test is the same as the difficulty of computing exact p-values for analogous ANOVA models based on measured observations.

7.1 Tests of correlation coefficients

Correlation studies arise with some frequency in experimental settings. The data in a correlation study consist of $n$ pairs of observations $\left(\left(x_1,y_1\right),\dots ,\left(x_n,y_n\right)\right)$ . As explained in every basic course on applied statistics, the Pearson correlation coefficient, $r$ , is a function of the products of deviations of the sample observations from their respective sample means:

(11) $r=\frac{{\sum }_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$

The null hypothesis of no correlation corresponds to a situation in which there is no monotonic relationship between the $x$ and $y$ values.

In conventional hypothesis testing, the observed value of the correlation coefficient, $r_{\mathrm{obs}}$ , would be compared to the null distribution of this statistic—something that could be inferred from assumptions about the distributions of the $x$ and $y$ values under the null hypothesis, or approximated via reference to a limit theorem for large enough sample sizes. As in the previously discussed locational tests, however, these approaches are difficult to defend for the small samples typical of many experiments.

Fortunately, permutation tests can be constructed around correlation coefficients. If there were no association between the variables, as the null hypothesis insists, then observed covariation would reflect exchangeable error in the experimental design. Because each observed value of the measured outcome would be equally likely to have been observed under every value of the treatment variable, the null distribution of the correlation test-statistic can be constructed by computing the value of $r$ under each of the $n!$ ways that one of the two variables could be reordered, holding the order of the other variable constant. For a two-sided test, the permutation p-value would be computed as follows:

(12) $\text{permutation}\text{test}\text{of}\text{correlation},\text{two-sided}p\text{-value}=\frac{{\sum }_{i=1}^{n!}1\left(\left(r_i\right)\ge \left(r_{\mathrm{obs}}\right)\right)}{n!}$

To illustrate, consider the traveler’s dilemma experiment reported by Capra et al. (Reference Capra, Gomez, Goeree and Holt1999). Subjects were randomly paired and tasked with making simultaneous “claims,” subject to the following rules. Each subject’s claim had to fall within the range from 80 to 200. After both claims were made, each subject would earn the smaller of the two claims, minus a penalty of R if the subject’s claim was the larger claim and plus a reward of R if the subject’s claim was the smaller claim. No penalty or reward was applied if the two claims were equal. Since each subject has a unilateral incentive to undercut anything greater than the minimum in this game, the unique Nash equilibrium is for both players to make the minimum claim regardless of the size of the incentive parameter. This prediction is counterintuitive. Common sense suggests that claims should fall with increases in the size of the penalty (for being high) and reward (for being low).

The experiment was conducted with session groups of 9–12 subjects, randomly matched into pairs to play 10 rounds of the game with one incentive parameter in Part A, followed by additional rounds with a different incentive parameter in Part B. There was one session for each pair of A and B treatments, as shown in Table 10. Here, we focus on the Part A data; we will turn to Part B in Sect. 7.2.

Table 10

Session-average claims for a traveler’s dilemma game

Treatments and results	S1	S2	S3	S4	S5	S6
Part A
Incentive term (R)	80	10	50	20	25	5
Average claim	82	186	92	116	146	196
Part B
Incentive term (R)	10	80	20	50	5	25
Average claim	163	99	86	82	171	170

Capra et al. (Reference Capra, Gomez, Goeree and Holt1999): average claims over final 5 rounds of treatment

The most salient feature of the Part A data in the top part of the table is the apparent sensitivity of claims to the size of the incentive parameter. Indeed, the Part A data appear more consistent with the intuitive prediction that average claims would vary inversely with the incentive parameter than with the Nash prediction of no effect. In fact, the observed value of the correlation between penalties and average claims in Part A is $r_{\mathrm{obs}}=-0.873$ . To test whether that observation is statistically significant, the null distribution of the correlation coefficient can be constructed by computing the correlation coefficient under every possible permutation of the observed claims data. Of the $6!=720$ possible permutations to be considered in this way, only 3 result in a correlation coefficient equal to or smaller than the observed value. Thus, a one-sided p-value for a permutation test of negative correlation against the null is $3/720\approx 0.004,$ a significant rejection of the Nash hypothesis, despite the limited sample size.Footnote ²⁵

Before moving on, note that nothing about the permutation strategy just described depends on use of the Pearson correlation coefficient as the test statistic. The Pearson statistic made sense for the theoretically continuous data under investigation, but other experimental contexts could motivate the use of other statistics. In an “bomb” risk aversion experiment (Crosetto & Filippin, Reference Crosetto and Filippin2013), for example, a subject is shown 12 boxes and allowed to check as many boxes as desired, understanding that each box checked earns the subject $1 unless and until a randomly hidden bomb is encountered—in which case nothing is earned. This setup can be used to elicit risk aversion, but the mapping from number-of-boxes-checked to some measure of risk aversion depends on the measure of risk aversion used (e.g., constant relative risk aversion) and is nonlinear in any event. In this application, it might make sense to rank subjects by the number of boxes they choose to check and by some other proxy for risk aversion, such as amount saved for retirement. In this ordinal context, a rank-based measure of association—such as Kendall’s $\tau$ or Spearman’s $\rho$ statistic—would be a better correlation concept. The same approach described above could be used to arrive at a permutation p-value in this setting. The only difference would be the measure of association used as the test statistic.

7.2 Tests of linear regression models

Moving beyond correlation analysis of a bivariate pattern, permutation methods can also be used to conduct statistical inference for linear models with more independent variables. There are, however, significant limitations to permutation inference in this setting. Obvious permutation strategies are forthcoming only for a few special regression models. In most cases, the experimenter will face a choice of different permutation strategies, and effort may be required to identify the appropriate strategy for the application. Standard parametric tests do not exhibit these difficulties and, for moderate sample sizes, may be robust enough to deviations from parametric assumptions to support credible inference. For small sample sizes, however, permutation tests will still constitute a more reliable basis for inference.

Starting with one of the lucky special cases for permutation testing, consider the following bivariate data generating process:

(13) $y=\alpha +\delta z+ϵ$

For $ϵ$ a mean-zero error term unrelated to the value of the regressor $z$ . Suppose interest is in testing the null hypothesis $\delta =0.$ Under the null hypothesis, all variation in the elements of $y$ is attributable to random error, $y_0=\alpha +ϵ$ , so every element of the $y$ vector is equally likely to have been paired with every element of the $z$ vector. This observation suggests a simple permutation strategy: compare the observed t statistic for the least squares estimate of $\delta$ against the set of t statistics calculated under all possible permutations of the order of elements in the $y$ vector while holding fixed the order of elements in the $z$ vector.Footnote ²⁶ This should look familiar. It is the permutation strategy for correlation coefficients described in Sect. 7.1.

As a concrete illustration, return to the Capra et al. (Reference Capra, Gomez, Goeree and Holt1999) traveler’s dilemma experiment, with data reproduced in Table 10. First consider a simple linear regression of average claim ( $y$ ) on penalty term ( $z$ ) using data from Part A of the experiment. Fitting the linear model $y=\alpha +\delta z+ϵ$ via OLS yields the following parameter estimates and standard errors:

$\begin{array}{cc}\hat{\alpha }=182.917& \hat{\delta }=-1.4711\\ \left(16.814\right)& \left(0.411\right)\end{array}$

The t statistic for this estimate of $\delta$ is $t_{\mathrm{obs}}=-1.471/0.411=-3.581$ . Under the null hypothesis that $\delta =0$ (the size of the incentive parameter has no effect on average claims), all variation in the $y$ vector is attributable to exchangeable error. This means that the null distribution of the t statistic can be computed by permuting the order of the $y$ vector and recalculating the value of the t statistic at each permutation. With 6 observations, there are $6!=720$ ways to permutate the order of elements in the $y$ vector under the null hypothesis. Of these 720 permutations, 12 yield t statistic values at least as extreme as $t_{\mathrm{obs}}$ , yielding a two-sided p-value of $12/720\approx 0.017$ for the effect of the incentive parameter on average claims in Part A.

Next consider the same model applied to Part B. Fitting the same model via OLS when using the data from Part B of the experiment yields the following parameter estimates and standard errors:

$\begin{array}{cc}\hat{\alpha }=159.434& \hat{\delta }=-0.977\\ \left(24.345\right)& \left(0.595\right)\end{array}$

The t statistic for this estimate of $\delta$ is $t_{\mathrm{obs}}=-1.642$ and—by the same permutation process as above—a two-sided p-value for the effect of the incentive parameter on average claims in Part B is found to be $150/720\approx 0.208$ .

Visual examination of the Part A and Part B outcomes in each treatment column of Table 10 suggests the presence of order effects—that is, experiences in Part A seem like they could be introducing confounding variation in the Part B data. To explore this possibility, consider the more complicated linear model that allows for average claims in Part B of the experiment ( $y$ ) to depend on both the Part B incentive parameter ( $z$ ) and the nuisance influence captured by the average claim in Part A ( $x$ ). Fitting the multiple regression model $y=\alpha +\beta x+\delta z+ϵ$ via OLS is trivial. Conducting permutation inference on this model is not.

To see why multiple regression with nuisance variables presents a harder problem, consider the general case of a data generating process with the following form:

(14) $y=X\beta +Z\delta +ϵ$

For $ϵ$ an error term as before, $X$ a matrix of nuisance variables which may include a constant term, and $Z$ a matrix of regressors of interest. Under the null hypothesis that $\delta =0$ , each element of the response vector is now more than random error: $y_0=X\beta +ϵ$ . This connection makes the simple permutation strategy of reordering $y$ generally indefensible. Variation in the $y$ vector is partly attributable to the influence of the nuisance variables in X . One might initially think that this could be solved by either subtracting $X\beta$ from both sides of the equation or by permuting the rows of the $X$ matrix in lockstep with the rows of the $y$ vector but neither of these strategies is attractive. The first requires knowledge of $\beta$ , which is unavailable in most interesting cases; the second fails to preserve collinearity between $X$ and $Z$ .

In fact, while many permutation strategies have been suggested for the multiple regression context, there remains no generally accepted permutation approach for this problem. A full survey of this literature is beyond the scope of discussion, but helpful commentaries are provided by Kennedy (Reference Kennedy1995), Kennedy and Cade (Reference Kennedy and Cade1996), Manly (Reference Manly2007: ch. 8), Anderson and Robinson (Reference Anderson and Robinson2001), and Winkler et al. (Reference Winkler, Ridgway, Webster, Smith and Nichols2014). To illustrate one intuitive option, consider the following permutation strategy due to Freedman and Lane (Reference Freedman and Lane1983):

1. 1. Fit the full model by OLS, $y=X\beta +Z\delta +ϵ$ , and compute the observed F statistic for testing the null hypothesis that $\delta =0$ .

2. 2. Fit the reduced model by OLS, $y=X\beta +ϵ$ , and use this model to compute a reduced-model prediction vector $\hat{y}=X\hat{\beta }$ and a reduced-model residual vector $r=y-\hat{y}$ .Footnote ²⁷

3. 3. Permute the order of the reduced-model residual vector $r$ and add each permutation of the residual vector to the reduced-model prediction vector $\hat{y}$ to generate a new permutation of the $y$ vector. For the observed data, this reconstructs the observed $y$ vector. For all other permutations, it constructs a new $y_p$ vector in which only the variation not explained by the nuisance variables is being permuted.

4. 4. For each such permutation, fit the full model $y_p=X\beta +Z\delta +ϵ$ and compute the F statistic for testing the null hypothesis that $\delta =0$ . For testing individual parameters, t statistics could be used instead.

5. 5. Compare the observed value of the F statistic (or t statistic) from step 1 to the permutation distribution from step 4 to compute a p-value for this test.

This procedure can be made concrete by applying it to the multiple regression model described above for Part B of the Capra et al. (Reference Capra, Gomez, Goeree and Holt1999) data. Specifically, consider regressing average claim in Part B ( $y$ ) on both average claim in Part A ( $x$ ) and incentive term in Part B ( $z$ ). Fitting the model $y=\alpha +\beta x+\delta z+ϵ$ via OLS yields the following parameter estimates and standard errors:

$\begin{array}{ccc}\hat{\alpha }=93.300& \hat{\beta }=0.589& \hat{\delta }=-1.425\\ \left(39.033\right)& \left(0.305\right)& \left(0.514\right)\end{array}$

For purposes of testing the joint null hypothesis that $\beta =\delta =0$ , the simple approach used in univariate regression could be repeated using the F statistic as a test statistic. For testing the less restrictive null that $\delta =0,$ a more involved permutation scheme is needed.

Following the Freedman and Lane (Reference Freedman and Lane1983) procedure outlined above, we can start by noting the observed t statistic value of $t_{\mathrm{obs}}=-1.425/0.514=-2.773$ for the estimate of $\delta$ . Fitting the reduced model $y=\alpha +\beta x+ϵ$ via OLS then yields the two components needed to construct approximate permutations of the $y$ vector: (1) a vector of predicted values from this reduced regression and (2) a vector of residuals from the reduced regression. To construct the null distribution, we permute the order of elements in the residual vector, each time adding it to the predicted value vector to form a new permutation of the $y$ vector, and then fit that new $y$ vector to the full model, recording the t statistic associated with the estimated value of $\delta$ . Of the 720 permutations of the t statistic computed in this manner, 35 yield t statistics with equal or greater absolute value than $t_{\mathrm{obs}}$ , resulting in a two-sided test p-value of about $0.049$ . This result is consistent with the guess that order effects could be muddying relationships in the Part B data.