Multi-Sample Nonparametrics, Rank Correlation & Resampling: Kruskal-Wallis, Friedman, Kendall & Permutation

Consider a completely randomized one-way experimental design comparing $k \ge 3$ mutually independent treatment groups or populations. For each group $i \in \{1, 2, \dots, k\}$, we observe an independent random sample $X_{i1}, X_{i2}, \dots, X_{in_i}$ of size $n_i$ drawn from a continuous cumulative distribution function $F_i(x)$, with total combined sample size $N = \sum_{i=1}^k n_i$. Under the location-shift parameterization $F_i(x) = F(x - \theta_i)$, where $F(\cdot)$ is an arbitrary unspecified continuous distribution with median zero, we wish to test the omnibus null hypothesis of identical population distributions against the general shift alternative:

Pooling, Ranking, and Exact Moments of Group Rank Sums Under $H_0$

To construct a distribution-free test statistic that is invariant under any strictly increasing monotonic transformation of the measurement scale, William Kruskal and W. Allen Wallis (1952) pooled all $N = \sum_{i=1}^k n_i$ observations into a single combined sample and ordered them from smallest to largest, assigning ranks $1, 2, \dots, N$. Let $R_{ij} = \operatorname{rank}(X_{ij})$ denote the joint rank of the $j$-th observation in the $i$-th group within the combined sample of size $N$. Define the $i$-th Group Rank Sum $R_{i\cdot}$ and the $i$-th Group Mean Rank $\bar{R}_{i\cdot}$:

Under the null hypothesis $H_0$ and continuity of $F$, all $N!$ rank permutations of $\{1, 2, \dots, N\}$ are equally likely, and the $\frac{N!}{n_1! n_2! \cdots n_k!}$ partitions of the integer set $\{1, 2, \dots, N\}$ into $k$ subsets of sizes $n_1, \dots, n_k$ occur with uniform probability $\frac{n_1! n_2! \cdots n_k!}{N!}$. Thus, for each group $i$, the set of ranks $\{R_{i1}, \dots, R_{in_i}\}$ constitutes a Simple Random Sample Without Replacement (SRSWOR) of size $n_i$ drawn from the finite population $\{1, 2, \dots, N\}$. Using the exact first and second moments of a finite population with mean $\mu_R = \frac{N+1}{2}$ and finite-population variance $\sigma_R^2 = \frac{1}{N}\sum_{r=1}^N \left(r - \frac{N+1}{2}\right)^2 = \frac{N^2 - 1}{12}$, we obtain the exact finite-sample expectation, variance, and covariance of the group mean ranks $\bar{R}_{i\cdot}$ under $H_0$:

Derivation of the Kruskal–Wallis $H$-Statistic and Its Computational Shortcut

If the classical one-way ANOVA between-group sum of squares $\text{SSB}_R = \sum_{i=1}^k n_i \left(\bar{R}_{i\cdot} - \bar{R}_{\cdot\cdot}\right)^2 = \sum_{i=1}^k n_i \left(\bar{R}_{i\cdot} - \frac{N+1}{2}\right)^2$ is scaled by the total rank variance divisor $\frac{N(N+1)}{12}$ (noting that $\operatorname{Var}_{H_0}\left(\sqrt{n_i}(\bar{R}_{i\cdot} - \frac{N+1}{2})\right) = \frac{N(N+1)}{12}\left(1 - \frac{n_i}{N}\right)$), we obtain the Kruskal–Wallis $H$-Statistic:

Wait: Why does the coefficient $\frac{12}{N(N+1)}$ appear rather than the naive standardized sum $\sum_{i=1}^k \frac{(\bar{R}_{i\cdot} - \mathbb{E}[\bar{R}_{i\cdot}])^2}{\operatorname{Var}(\bar{R}_{i\cdot})}$? Because the vector of centered rank averages $\mathbf{Z} = \left(\bar{R}_{1\cdot} - \frac{N+1}{2}, \dots, \bar{R}_{k\cdot} - \frac{N+1}{2}\right)^\top$ satisfies the linear dependency $\sum_{i=1}^k n_i Z_i = 0$, its $k \times k$ covariance matrix $\boldsymbol{\Sigma}_Z$ has rank $k - 1$ with generalized inverse $\boldsymbol{\Sigma}_Z^- = \frac{12}{N(N+1)}\operatorname{diag}(n_1, n_2, \dots, n_k)$! Thus, $H = \mathbf{Z}^\top \boldsymbol{\Sigma}_Z^- \mathbf{Z}$ is the exact Mahalanobis quadratic form of the rank means. Expanding the square inside $H$ yields the computationally convenient hand-calculation shortcut formula:

Multiplying through by $\frac{12}{N(N+1)}$ gives:

By the multivariate Hájek–Wald–Wolfowitz Central Limit Theorem for linear rank statistics, as $\min_{1 \le i \le k} n_i \to \infty$ with $\frac{n_i}{N} \to \lambda_i \in (0, 1)$, the quadratic form $H$ converges in distribution under $H_0$ to a central Chi-square distribution with $k - 1$ degrees of freedom:

We reject $H_0$ at significance level $\alpha$ if $H \ge \chi^2_{\alpha, \, k - 1}$ (or if $H$ exceeds the exact permutation critical value for small $n_i \le 5$). Under Gaussian distributions, the Pitman Asymptotic Relative Efficiency (ARE) of the Kruskal–Wallis test relative to the parametric one-way ANOVA $F$-test is $\frac{3}{\pi} \approx 0.9549$ ($95.5\%$), while for Laplace (double exponential) errors its ARE is $1.50$ ($150\%$) and for heavy-tailed distributions it never falls below $0.864$!

Exact Variance Adjustment for Tied Observations ($H_{\text{adj}}$)

When ties occur among the $N$ observations and each tied group of size $t_g$ ($g = 1, \dots, G$, where $t_g$ is the number of observations tied at the $g$-th distinct value) is assigned the average midrank, the mean of the ranks remains $\frac{N+1}{2}$, but the finite-population variance of the $N$ midranks is strictly reduced from $\frac{N^2 - 1}{12}$ by $\frac{\sum_{g=1}^G (t_g^3 - t_g)}{12 N}$:

Dividing the unadjusted $H$ statistic by the tie-correction factor $C_{\text{ties}} = 1 - \frac{\sum_{g=1}^G (t_g^3 - t_g)}{N^3 - N} \in (0, 1]$ yields the Tie-Adjusted Kruskal–Wallis Statistic:

Because $C_{\text{ties}} < 1$ whenever at least one tie is present ($t_g \ge 2 \implies t_g^3 - t_g \ge 6$), $H_{\text{adj}} > H$, so failing to correct for ties makes the unadjusted $H$ test slightly conservative.

Dunn's (1964) Post-Hoc Multiple Comparisons Procedure

When the omnibus Kruskal–Wallis test rejects $H_0$, we identify which specific pairs of groups $(i, \ell)$ differ significantly using Olive Jean Dunn's (1964) joint-rank multiple comparison test. Under $H_0$, the difference between the $i$-th and $\ell$-th group mean ranks $\bar{R}_{i\cdot} - \bar{R}_{\ell\cdot}$ has exact variance:

Thus, Dunn's Standardized $Z$-Statistic for comparing Group $i$ and Group $\ell$ is:

To control the Family-Wise Error Rate (FWER) at level $\alpha$ across all $M = \binom{k}{2} = \frac{k(k-1)}{2}$ pairwise comparisons, Dunn applied the Bonferroni inequality (or Holm's step-down procedure), declaring Group $i$ and Group $\ell$ significantly different if $Z_{i\ell} \ge z_{\alpha / (2M)} = z_{\alpha / [k(k-1)]}$.

---

In many agricultural, industrial, clinical, and psychological experiments, observations are arranged in a Randomized Complete Block Design (RCBD) (or a repeated-measures design) in which $b$ homogeneous blocks (or subjects) are each evaluated across all $k$ treatments. Let $X_{ij}$ denote the response of block $i \in \{1, \dots, b\}$ under treatment $j \in \{1, \dots, k\}$:

where $\beta_i$ is the nuisance block effect and $\tau_j$ is the treatment effect of interest. We wish to test $H_0: \tau_1 = \tau_2 = \dots = \tau_k = 0$ against $H_1: \text{not all } \tau_j \text{ are equal}$, without assuming normality of the errors $\varepsilon_{ij}$.

Intra-Block Ranking and Derivation of Friedman's $Q_F$ Statistic (1937)

Because block effects $\beta_i$ vary arbitrarily across blocks, pooling all $b k$ observations together (as in Kruskal–Wallis) would confound treatment differences with inter-block heterogeneity. Instead, Milton Friedman (1937) eliminated the additive block effect $\beta_i$ completely by ranking observations separately within each block:

• Within block $i$, order $X_{i1}, X_{i2}, \dots, X_{ik}$ from smallest to largest and assign intra-block ranks $r_{i1}, r_{i2}, \dots, r_{ik} \in \{1, 2, \dots, k\}$ (assigning midranks if ties occur within block $i$).
• Compute the $j$-th Treatment Rank Sum $R_{\cdot j}$ and Treatment Mean Rank $\bar{R}_{\cdot j}$ across the $b$ blocks:

Under $H_0$, within each independent block $i$, the rank vector $(r_{i1}, \dots, r_{ik})$ is uniformly distributed over all $k!$ permutations of $\{1, \dots, k\}$, with $\mathbb{E}[r_{ij}] = \frac{k+1}{2}$, $\operatorname{Var}(r_{ij}) = \frac{k^2 - 1}{12}$, and $\operatorname{Cov}(r_{ij}, r_{i\ell}) = -\frac{k+1}{12}$ ($j \neq \ell$). Summing across the $b$ mutually independent blocks yields:

Scaling the between-treatment sum of squares $\sum_{j=1}^k b\left(\bar{R}_{\cdot j} - \frac{k+1}{2}\right)^2 = \frac{1}{b}\sum_{j=1}^k \left(R_{\cdot j} - \frac{b(k+1)}{2}\right)^2$ by $\frac{12}{k(k+1)}$ yields Friedman's Two-Way Analysis of Variance by Ranks Statistic ($Q_F$ or $\chi^2_r$):

If intra-block ties occur in block $i$ with tie group sizes $t_{ig}$ ($g = 1, \dots, G_i$), the tie-adjusted Friedman statistic is:

The Iman–Davenport (1980) $F$-Distribution Improvement

Because the $\chi^2_{k-1}$ approximation to $Q_F$ is noticeably conservative in small-to-moderate block sizes $b$ (since $Q_F$ is bounded above by $b(k - 1)$), Ronald Iman and James Davenport (1980) demonstrated that transforming $Q_F$ into the classical two-way ANOVA $F$-statistic computed directly on the intra-block ranks produces a dramatically more accurate finite-sample approximation:

Notice that as $Q_F \to b(k - 1)$ (complete agreement of rankings across all $b$ blocks), $F_{\text{ID}} \to \infty$, faithfully reflecting zero residual error variance!

Cochran's $Q$-Test for Binary ($0/1$) Matched Responses across $k \ge 3$ Treatments (1950)

When the response $X_{ij} \in \{0, 1\}$ is dichotomous ($1 = \text{Success}$, $0 = \text{Failure}$) across $b$ matched subjects or blocks and $k \ge 3$ treatments (extending McNemar's $k = 2$ paired test to $k \ge 3$ treatments), William Cochran (1950) derived the exact conditional randomization test statistic. Let $C_j = \sum_{i=1}^b X_{ij}$ denote the total number of successes under treatment $j$ (column total), and let $L_i = \sum_{j=1}^k X_{ij}$ denote the total number of successes for block $i$ (row total). Cochran's $Q$-Statistic is identical to the tie-adjusted Friedman statistic $Q_{F,\text{adj}}$ applied to binary $0/1$ scores:

Notice that blocks with $L_i = 0$ (all $0$s) or $L_i = k$ (all $1$s) contribute $k L_i - L_i^2 = 0$ to the denominator and equal increments to all $C_j$, so concordant rows automatically drop out of Cochran's $Q$, exactly generalizing McNemar's test!

When Friedman's omnibus test (or Cochran's $Q$-test) rejects the null hypothesis $H_0: \tau_1 = \dots = \tau_k$, post-hoc pairwise comparisons between Treatment $j$ and Treatment $\ell$ are conducted using the Nemenyi–Wilcoxon–Wilcox Multiple Comparison Procedure or Conover's (1999) $t$-Test on Intra-Block Ranks. Under $H_0$, the difference between two treatment rank sums $R_{\cdot j} - R_{\cdot \ell}$ has exact variance $\operatorname{Var}_{H_0}(R_{\cdot j} - R_{\cdot \ell}) = 2\operatorname{Var}_{H_0}(R_{\cdot j}) - 2\operatorname{Cov}_{H_0}(R_{\cdot j}, R_{\cdot \ell}) = \frac{2b(k^2 - 1)}{12} + \frac{2b(k + 1)}{12} = \frac{b k(k + 1)}{6}$, so the standard error of the mean-rank difference $\bar{R}_{\cdot j} - \bar{R}_{\cdot \ell}$ is $\sqrt{\frac{k(k + 1)}{6b}}$. Applying a Bonferroni correction across all $M = \binom{k}{2}$ pairs declares Treatments $j$ and $\ell$ significantly different at family-wise error rate $\alpha$ whenever:

where $q_{\alpha, k, \infty}$ is the upper-$\alpha$ critical value of the Studentized Range distribution for $k$ groups.

---

Suppose $b$ independent judges (or clinical raters, consumer panels, or machine learning ranking algorithms) each rank $k$ candidate objects (e.g., job applicants, medical images, or competing products) from $1$ to $k$. Rather than testing whether treatments differ as an experimental intervention, our goal is to measure the degree of inter-rater agreement (concordance) among the $b$ judges.

Let $R_{\cdot j} = \sum_{i=1}^b r_{ij}$ denote the sum of ranks assigned to object $j$ by the $b$ judges, with mean $\bar{R}_{\cdot} = \frac{b(k+1)}{2}$. Define the Sum of Squared Deviations of Object Rank Sums ($S$):

Let us evaluate the two extreme boundary regimes for $S$:

1. Complete Disagreement / Maximum Discordance ($S_{\min} = 0$): If the $b$ judges disagree in such a completely balanced fashion that every object $j = 1, \dots, k$ receives the exact same total rank sum $R_{\cdot j} = \frac{b(k+1)}{2}$, then $S = 0$.

2. Perfect Agreement / Complete Concordance ($S = S_{\max}$): If all $b$ judges assign identical rankings to the $k$ objects, then after relabeling the objects in increasing order of preference, object $1$ receives rank $1$ from every judge ($R_{\cdot 1} = b$), object $2$ receives rank $2$ from every judge ($R_{\cdot 2} = 2b$), $\dots$, and object $k$ receives rank $k$ from every judge ($R_{\cdot k} = b k$). Therefore, the maximum possible sum of squared deviations is:

Definition of Kendall's Coefficient of Concordance $W$ (1939)

Maurice Kendall and Bernard Babington Smith (1939) defined Kendall's Coefficient of Concordance ($W$) as the ratio of the observed sum of squares $S$ to its theoretical maximum $S_{\max}$:

When intra-judge ties occur, $S_{\max}$ is reduced by $\frac{b}{12}\sum_{i=1}^b \sum_{g=1}^{G_i}(t_{ig}^3 - t_{ig})$, yielding the tie-corrected coefficient $W_{\text{adj}} = \frac{12 S}{b^2 k(k^2 - 1) - b \sum_{i=1}^b \sum_{g=1}^{G_i}(t_{ig}^3 - t_{ig})}$.

Exact Algebraic Identities Linking $W$, Friedman's $Q_F$, and Average Pairwise Spearman Correlation $\bar{\rho}_s$

Comparing the numerator and denominator of $W = \frac{12 S}{b^2 k(k^2 - 1)}$ with Friedman's statistic $Q_F = \frac{12 S}{b k(k + 1)}$ reveals a remarkably clean exact algebraic identity:

Consequently, testing the null hypothesis of random, independent judge rankings $H_0: W = 0$ is mathematically identical to executing Friedman's $\chi^2_{k-1}$ test via $Q_F = b(k - 1)W \sim \chi^2_{k-1}$ (or the Iman–Davenport test $F_{\text{ID}} = \frac{(b - 1)W}{1 - W} \sim F_{k-1, (b-1)(k-1)}$)!

Furthermore, let $\rho_{s}^{(i, m)} = 1 - \frac{6\sum_{j=1}^k (r_{ij} - r_{mj})^2}{k(k^2 - 1)}$ denote Spearman's rank correlation coefficient between judge $i$ and judge $m$, and let $\bar{\rho}_s = \frac{2}{b(b - 1)}\sum_{1 \le i < m \le b} \rho_s^{(i, m)}$ denote the average of all $\binom{b}{2}$ pairwise Spearman rank correlations among the $b$ judges. Expanding $\sum_{j=1}^k \left(\sum_{i=1}^b (r_{ij} - \frac{k+1}{2})\right)^2$ into diagonal and off-diagonal terms establishes the exact relationship:

Notice why $W \in [0, 1]$ cannot be negative even when judges negatively correlate: because $W \ge 0$, the average correlation among $b$ rankings is bounded below by $\bar{\rho}_s \ge -\frac{1}{b - 1}$! Moreover, when $H_0: W = 0$ is rejected and the judges exhibit statistically significant concordance, Kendall proved that the optimal consensus ranking—defined as the permutation $\hat{\boldsymbol{\pi}} \in \mathcal{S}_k$ that minimizes the total sum of squared Spearman rank distances $\sum_{i=1}^b \sum_{j=1}^k (r_{ij} - \pi_j)^2$ across all $b$ judges—is achieved by ordering the $k$ objects in ascending order of their aggregate rank sums $R_{\cdot 1}, R_{\cdot 2}, \dots, R_{\cdot k}$.

---

Both the Kruskal–Wallis $H$-test and Friedman's $Q_F$-test are omnibus tests whose test statistics are symmetric functions of the group rank sums $(R_{1\cdot}, \dots, R_{k\cdot})$. If we permute the labels of the $k$ treatment groups arbitrarily, $H$ and $Q_F$ remain completely unchanged! However, in dose-response toxicology, educational intervention trials, and economic development studies, the $k$ treatments often correspond to increasing dosage levels $d_1 < d_2 < \dots < d_k$, and the investigator specifically hypothesizes a monotonic ordered trend:

Spreading omnibus test power across all $k!$ possible orderings wastes substantial statistical sensitivity when the alternative is restricted to the single monotonic cone $\theta_1 \le \dots \le \theta_k$.

1. The Jonckheere–Terpstra ($JT$) Test for Independent Samples (1952, 1954)

Independently proposed by T. J. Terpstra (1952) and A. R. Jonckheere (1954), the Jonckheere–Terpstra Test aggregates all $\binom{k}{2} = \frac{k(k-1)}{2}$ pairwise Mann–Whitney count statistics $U_{i\ell}$ comparing each lower-indexed group $i$ against each higher-indexed group $\ell$ ($1 \le i < \ell \le k$):

Here $U_{i\ell}$ counts the number of pairs $(X_{ia}, X_{\ell b})$ with $i < \ell$ for which the observation in the higher-indexed group $\ell$ exceeds the observation in the lower-indexed group $i$. Under the strictly increasing trend $\theta_1 < \theta_2 < \dots < \theta_k$, each $U_{i\ell}$ exceeds its null expectation $\frac{n_i n_\ell}{2}$, reinforcing constructively across all $\binom{k}{2}$ pairs! Under $H_0$ (with no ties and total sample size $N = \sum_{i=1}^k n_i$), the exact mean and variance of $J$ are:

As $N \to \infty$, the standardized statistic $Z_{\text{JT}} = \frac{J - \mathbb{E}_{H_0}[J]}{\sqrt{\operatorname{Var}_{H_0}(J)}} \xrightarrow{d} \mathcal{N}(0, 1)$, and we reject $H_0$ in favor of an increasing trend if $Z_{\text{JT}} \ge z_\alpha$. Notice that whereas the omnibus Kruskal–Wallis test converts deviations into a $k-1$ degree-of-freedom quadratic form whose non-centrality parameter is diluted by $k - 1$, the Jonckheere–Terpstra statistic concentrates all monotonic trend signal into a single $1$-degree-of-freedom directional normal score $Z_{\text{JT}}$ (or equivalently $Z_{\text{JT}}^2 \sim \chi^2_1$). Under linear shift alternatives $\theta_i = i \cdot \Delta / \sqrt{N}$ with Gaussian disturbances, the Pitman Asymptotic Relative Efficiency of the Jonckheere–Terpstra test relative to parametric linear regression on group indices is $\frac{3}{\pi} \approx 0.955$, and its power advantage over Kruskal–Wallis increases monotonically with the number of groups $k$.

2. Page's $L$-Test for Ordered Alternatives in Randomized Block Designs (1963)

In a randomized complete block design with $b$ blocks and $k$ treatments ordered a priori as $1, 2, \dots, k$ under the monotonic alternative $H_1: \tau_1 \le \tau_2 \le \dots \le \tau_k$, Ellis Batten Page (1963) defined the linear rank contrast statistic $L$ as the weighted sum of the treatment rank sums $R_{\cdot j} = \sum_{i=1}^b r_{ij}$ weighted by their hypothesized integer order $j \in \{1, 2, \dots, k\}$:

Why does $L$ focus power on the ordered alternative? Notice that within each block $i$, the linear combination $\sum_{j=1}^k j r_{ij}$ is linearly equivalent to Spearman's rank correlation between the observed intra-block ranks $(r_{i1}, \dots, r_{ik})$ and the hypothesized order $(1, 2, \dots, k)$! Using $\mathbb{E}_{H_0}[R_{\cdot j}] = \frac{b(k+1)}{2}$ and the covariance matrix of $(R_{\cdot 1}, \dots, R_{\cdot k})$, the exact null expectation and variance of $L$ under $H_0$ are:

For moderate or large $b$, the standardized statistic (with continuity correction $-0.5$) is asymptotically standard normal:

We reject $H_0$ in favor of the increasing monotonic trend $H_1: \tau_1 \le \dots \le \tau_k$ at level $\alpha$ if $Z_L \ge z_\alpha$. Notice that $Z_L^2 = \frac{[12 L - 3 b k(k+1)^2]^2}{b k^2(k+1)^2(k-1)} \le Q_F$ represents the exact single-degree-of-freedom linear trend component extracted from Friedman's $(k-1)$-degree-of-freedom statistic $Q_F$, while the remaining difference $Q_{\text{rem}} = Q_F - Z_L^2 \sim \chi^2_{k-2}$ tests for higher-order non-linear departures from monotonicity!

---

Let $(X_1, Y_1), (X_2, Y_2), \dots, (X_n, Y_n)$ be an independent random sample of size $n$ from a bivariate population with joint cumulative distribution function $F_{X,Y}(x, y)$ and marginal distributions $F_X(x)$ and $F_Y(y)$. While Pearson's product-moment correlation $r = \frac{\sum (X_i - \bar{X})(Y_i - \bar{Y})}{\sqrt{\sum (X_i - \bar{X})^2 \sum (Y_i - \bar{Y})^2}}$ measures strictly linear dependence and is non-robust to outliers or non-linear monotonic transformations (e.g., $Y = e^X$), nonparametric correlation coefficients measure monotonic association and are strictly invariant under any order-preserving transformations $X \mapsto g_1(X), Y \mapsto g_2(Y)$.

1. Spearman's Rank Correlation Coefficient $\rho_s$ (1904)

Let $R_i = \operatorname{rank}(X_i) \in \{1, \dots, n\}$ and $S_i = \operatorname{rank}(Y_i) \in \{1, \dots, n\}$ denote the marginal ranks of $X_i$ and $Y_i$ within their respective samples. Charles Spearman (1904) defined $\rho_s$ (or $r_s$) as Pearson's product-moment correlation evaluated directly on the rank pairs $(R_i, S_i)$:

When there are no ties in either $X$ or $Y$, both $\{R_1, \dots, R_n\}$ and $\{S_1, \dots, S_n\}$ are permutations of $\{1, 2, \dots, n\}$ with $\bar{R} = \bar{S} = \frac{n+1}{2}$ and $\sum_{i=1}^n (R_i - \bar{R})^2 = \sum_{i=1}^n (S_i - \bar{S})^2 = \frac{n(n^2 - 1)}{12}$. Using the algebraic identity for the squared rank differences $d_i = R_i - S_i$:

we solve for $\sum_{i=1}^n (R_i - \bar{R})(S_i - \bar{S}) = \frac{n(n^2 - 1)}{12} - \frac{1}{2}\sum_{i=1}^n d_i^2$ and divide by $\frac{n(n^2 - 1)}{12}$ to obtain Spearman's Shortcut Formula:

Under the null hypothesis of bivariate independence $H_0: F_{X,Y}(x, y) = F_X(x)F_Y(y)$, all $n!$ pairings of $\mathbf{S}$ with $\mathbf{R}$ are equally likely, yielding exact moments:

Thus, for large $n$, $Z_s = \rho_s \sqrt{n - 1} \xrightarrow{d} \mathcal{N}(0, 1)$, or even more accurately in finite samples ($n \ge 10$), Student's $t$-approximation $t_s = \rho_s \sqrt{\frac{n - 2}{1 - \rho_s^2}} \sim t_{n - 2}$. Under bivariate normality, the Pitman ARE of Spearman's test relative to Pearson's $r$ test is $\left(\frac{3}{\pi}\right)^2 = \frac{9}{\pi^2} \approx 0.9119$ ($91.2\%$).

2. Kendall's Rank Correlation $\tau$ (1938): Concordant/Discordant Pairs and U-Statistic Variance

Maurice Kendall (1938) formulated a direct probabilistic measure of ordinal concordance based on all $\binom{n}{2} = \frac{n(n-1)}{2}$ distinct pairs of observations $\{(X_i, Y_i), (X_j, Y_j)\}$ ($1 \le i < j \le n$):

• A pair $(i, j)$ is Concordant ($C$) if $(X_i - X_j)(Y_i - Y_j) > 0$ (both $X_i > X_j, Y_i > Y_j$ or both $X_i < X_j, Y_i < Y_j$; the two variables move in the same direction).
• A pair $(i, j)$ is Discordant ($D$) if $(X_i - X_j)(Y_i - Y_j) < 0$ ($X_i$ and $Y_j$ move in opposite directions).
• A pair $(i, j)$ is Tied if $X_i = X_j$ or $Y_i = Y_j$.

In the population, let $\pi_c = \mathbb{P}\left((X_1 - X_2)(Y_1 - Y_2) > 0\right)$ and $\pi_d = \mathbb{P}\left((X_1 - X_2)(Y_1 - Y_2) < 0\right)$. Kendall's Population Parameter $\tau = \pi_c - \pi_d \in [-1, +1]$ is the difference between the probability of concordance and the probability of discordance for two randomly drawn bivariate observations! Its sample estimator (Kendall's $\tau_a$ when no ties exist) is a symmetric Hoeffding $U$-Statistic of degree $m = 2$ with bounded kernel $h\left((X_i, Y_i), (X_j, Y_j)\right) = \operatorname{sgn}(X_i - X_j)\operatorname{sgn}(Y_i - Y_j) \in \{-1, 0, +1\}$:

where $N_c$ is the total number of concordant pairs and $N_d$ is the total number of discordant pairs ($N_c + N_d = \binom{n}{2}$ when no ties exist). Under $H_0$ (independence), by Hoeffding's $U$-statistic projection theorem, the exact finite-sample expectation and variance of $\hat{\tau}_a$ and of $S_K = N_c - N_d$ are:

Thus, the standardized asymptotic normal test statistic under $H_0$ is:

Under bivariate normality with correlation $\rho$, Greiner's relation states $\tau = \frac{2}{\pi}\arcsin(\rho)$, and the Pitman ARE of Kendall's $\tau$ test relative to Pearson's $r$ is identical to Spearman's test: $\frac{9}{\pi^2} \approx 0.9119$.

3. Adjusting for Ties: Kendall's $\tau_b$ vs. Goodman–Kruskal $\gamma$ (1954)

When ties occur ($T_X = \frac{1}{2}\sum_g t_{X,g}(t_{X,g} - 1)$ pairs tied on $X$, and $T_Y = \frac{1}{2}\sum_h t_{Y,h}(t_{Y,h} - 1)$ pairs tied on $Y$), two primary adjustments exist:

1. Kendall's $\tau_b$: Normalizes $N_c - N_d$ by the geometric mean of the number of untied $X$-pairs and untied $Y$-pairs:

2. Goodman and Kruskal's Gamma ($\gamma$, 1954): Completely conditions on the event that a pair is untied on both variables, measuring the conditional excess of concordant over discordant pairs among all strictly untied pairs:

In ordinal contingency tables with many ties, $|\hat{\gamma}| \ge |\hat{\tau}_b| \ge |\hat{\tau}_a|$, and $\hat{\gamma} = \pm 1$ under weak monotonicity (e.g., upper-triangular or staircase tables) even when $\hat{\tau}_b < 1$.

---

Whereas rank tests achieve distribution-free exactness by replacing raw observations $X_i$ with their integer ranks $R_i$ (thereby discarding exact metric distances between observations), Ronald Fisher (1935) and Edwin Pitman (1937) demonstrated that **exact distribution-free tests can be constructed directly from the original metric data** by conditioning on the observed order statistics! ### Exchangeability Under $H_0$ and Exact Two-Sample Fisher–Pitman Permutation Tests Let $\mathbf{Z} = (Z_1, Z_2, \dots, Z_N) = (X_1, \dots, X_m, Y_1, \dots, Y_n)$ denote the combined vector of $N = m + n$ observations from two independent samples $X_1, \dots, X_m \sim F_X$ and $Y_1, \dots, Y_n \sim F_Y$. Under the null hypothesis of identical distributions $H_0: F_X = F_Y = F$, the joint distribution of $(Z_1, \dots, Z_N)$ is **exchangeable**: for every permutation $\pi \in \mathcal{S}_N$ of $\{1, 2, \dots, N\}$, the permuted vector $(Z_{\pi(1)}, \dots, Z_{\pi(N)})$ has the exact same joint probability distribution as $(Z_1, \dots, Z_N)$. Let $\mathbf{Z}_{(\cdot)} = (Z_{(1)} \le Z_{(2)} \le \dots \le Z_{(N)})$ denote the vector of pooled order statistics. Under $H_0$, $\mathbf{Z}_{(\cdot)}$ is a **complete sufficient statistic** for the unknown common distribution $F$. Conditioning on the observed order statistics $\mathbf{Z}_{(\cdot)} = \mathbf{z}_{(\cdot)}$, all $\binom{N}{m} = \frac{N!}{m! n!}$ distinct assignments of $m$ values to the treatment group $\mathcal{X}^*$ and $n$ values to the control group $\mathcal{Y}^*$ have **exact conditional probability $1 / \binom{N}{m}$**, regardless of what the true underlying distribution $F$ is! Let $T(\mathbf{Z}) = \bar{X} - \bar{Y}$ (or any test statistic, such as the Studentized $t$-statistic or trimmed mean difference) be evaluated on the observed sample, yielding $t_{\text{obs}} = T(\mathbf{z})$. Note that because the grand sum $S_{\text{total}} = \sum_{i=1}^m X_i + \sum_{j=1}^n Y_j$ is strictly constant across all $\binom{N}{m}$ partitions, $\bar{X} - \bar{Y} = \frac{1}{m}\sum_{i=1}^m X_i - \frac{1}{n}\left(S_{\text{total}} - \sum_{i=1}^m X_i\right) = \frac{N}{m n}\sum_{i=1}^m X_i - \frac{S_{\text{total}}}{n}$ is a strictly increasing linear function of the first-sample sum $S_X = \sum_{i=1}^m X_i$! Thus, the **Exact One-Sided Fisher–Pitman Permutation $p$-Value** is simply the fraction of all $M = \binom{N}{m}$ partitions $\mathcal{P}_r$ ($r = 1, \dots, M$) whose test statistic $T(\mathbf{z}_r^*)$ is at least as extreme as $t_{\text{obs}}$: $$p_{\text{exact}} = \mathbb{P}_{H_0}\left( T(\mathbf{Z}^*) \ge t_{\text{obs}} \;\middle|\; \mathbf{Z}_{(\cdot)} = \mathbf{z}_{(\cdot)} \right) = \frac{1}{\binom{N}{m}} \sum_{r=1}^{\binom{N}{m}} \mathbf{1}_{\left\{ T(\mathbf{z}_r^*) \ge t_{\text{obs}} \right\}}$$ By Lehmann and Stein's (1949) theorem, any test that maintains an exact Type I error rate $\alpha$ across all continuous distributions $F$ **must** be a conditional permutation test given $\mathbf{Z}_{(\cdot)}$, and the Fisher–Pitman permutation test has a Pitman ARE of **$1.000$ ($100\%$)** relative to Student's $t$-test under Gaussian populations! In matched-pairs designs where $n$ paired differences $D_i = X_i - Y_i$ ($i = 1, \dots, n$) are observed under the null hypothesis that each $D_i$ is symmetrically distributed around zero ($f_D(d) = f_D(-d)$), the conditional randomization group consists of all $2^n$ **sign-flippings** $\mathbf{D}^*(\mathbf{s}) = (s_1 |D_1|, s_2 |D_2|, \dots, s_n |D_n|)$ with independent signs $s_i \in \{-1, +1\}$. Conditioning on the observed absolute differences $(|D_1|, \dots, |D_n|)$, each of the $2^n$ sign assignments occurs with exact probability $2^{-n}$, and the Fisher–Pitman paired randomization test computes the fraction of the $2^n$ signed sums $\sum_{i=1}^n s_i |D_i|$ that equal or exceed the observed sum $\sum_{i=1}^n D_i$. Moreover, by exploiting the **Neyman–Pearson Duality Between Hypothesis Tests and Confidence Regions**, we can invert the two-sample Fisher–Pitman permutation test of $H_0(\Delta_0): F_X(x - \Delta_0) = F_Y(x)$ by subtracting the hypothesized shift $\Delta_0$ from every treatment observation ($X_i(\Delta_0) = X_i - \Delta_0$) and collecting all values of $\Delta_0 \in \mathbb{R}$ that are not rejected at level $\alpha$. This test inversion produces an exact $100(1 - \alpha)\%$ distribution-free confidence interval $[\hat{\Delta}_L, \hat{\Delta}_U]$ for the location shift $\Delta$ whose endpoints coincide with ordered pairwise differences $X_i - Y_j$ when applied to Wilcoxon rank sums (the **Hodges–Lehmann estimator** and Moses confidence interval). ### Unbiased Monte Carlo Randomization $P$-Values (Phipson and Smyth, 2010) When $m$ and $n$ are moderate (e.g., $m = n = 20 \implies \binom{40}{20} \approx 1.378 \times 10^{11}$), enumerating all $\binom{N}{m}$ partitions is computationally intractable. Instead, we draw $B$ independent uniform random permutations $\mathbf{z}_1^*, \dots, \mathbf{z}_B^*$ (typically $B = 9{,}999$ or $99{,}999$) and count how many simulated values satisfy $T(\mathbf{z}_b^*) \ge t_{\text{obs}}$. > [!IMPORTANT] > **Why the $+1$ Pseudocount Is Mandatory**: If $K = \sum_{b=1}^B \mathbf{1}_{\{T(\mathbf{z}_b^*) \ge t_{\text{obs}}\}}$, the naive Monte Carlo estimator $\hat{p}_{\text{naive}} = K / B$ can equal $0$ with positive probability, violating the exact validity requirement $\mathbb{P}_{H_0}(\hat{p} \le \alpha) \le \alpha$. Including the observed data configuration itself as one of the $B + 1$ exchangeable realizations under $H_0$ yields the **exact valid Monte Carlo $p$-value** (Barnard, 1963; Phipson and Smyth, 2010): > $$\hat{p}_{\text{MC}} = \frac{1 + \sum_{b=1}^B \mathbf{1}_{\left\{ T(\mathbf{z}_b^*) \ge t_{\text{obs}} \right\}}}{B + 1}$$ > Under $H_0$, the rank of $t_{\text{obs}}$ among $\{t_{\text{obs}}, T(\mathbf{z}_1^*), \dots, T(\mathbf{z}_B^*)\}$ is discrete uniform on $\{1, 2, \dots, B + 1\}$, guaranteeing $\mathbb{P}_{H_0}\left(\hat{p}_{\text{MC}} \le \frac{j}{B+1}\right) = \frac{j}{B+1}$ exactly! ### Bootstrap Hypothesis Testing vs. Permutation Testing: Shifting to Impose $H_0$ and Pivotal Studentization Unlike permutation testing—which samples **without replacement** from the pooled data and strictly requires total exchangeability of the distributions ($F_X = F_Y$) under $H_0$—**Bootstrap Hypothesis Testing** (Efron and Tibshirani, 1993; Hall and Wilson, 1991) can test equality of a specific parameter such as means $H_0: \mu_X = \mu_Y$ even when the population variances or shapes differ ($\sigma_X^2 \neq \sigma_Y^2$, known as the nonparametric Behrens–Fisher problem)! However, executing a valid bootstrap hypothesis test requires obeying **Hall and Wilson's Two Golden Rules of Bootstrap Testing**: 1. **Resample From an Empirical Distribution That Explicitly Satisfies $H_0$**: If we naively resampled $X_1^*, \dots, X_m^*$ from $\{X_1, \dots, X_m\}$ (which has sample mean $\bar{X}$) and $Y_1^*, \dots, Y_n^*$ from $\{Y_1, \dots, Y_n\}$ (which has sample mean $\bar{Y} \neq \bar{X}$), the bootstrap world would have true mean difference $\bar{X} - \bar{Y} \neq 0$, mimicking $H_1$ rather than $H_0$! Instead, we must recenter each sample to share the common pooled mean $\bar{Z} = \frac{m\bar{X} + n\bar{Y}}{m+n}$ (or recenter both to $0$): $$\tilde{X}_i = X_i - \bar{X} + \bar{Z} \quad (i = 1, \dots, m), \qquad \tilde{Y}_j = Y_j - \bar{Y} + \bar{Z} \quad (j = 1, \dots, n)$$ and draw $X_1^*, \dots, X_m^*$ with replacement from $\{\tilde{X}_1, \dots, \tilde{X}_m\}$ and $Y_1^*, \dots, Y_n^*$ with replacement from $\{\tilde{Y}_1, \dots, \tilde{Y}_n\}$. 2. **Pivot / Studentize the Bootstrap Test Statistic**: Compute the Welch-studentized statistic on the observed sample and on each bootstrap resample $b = 1, \dots, B$: $$t_{\text{obs}} = \frac{\bar{X} - \bar{Y}}{\sqrt{\frac{s_X^2}{m} + \frac{s_Y^2}{n}}}, \qquad t_b^* = \frac{\bar{X}_b^* - \bar{Y}_b^*}{\sqrt{\frac{(s_{X,b}^*)^2}{m} + \frac{(s_{Y,b}^*)^2}{n}}}$$ Because $t_b^*$ is an **asymptotically pivotal quantity** (its limiting $\mathcal{N}(0, 1)$ distribution is free of unknown nuisance variances $\sigma_X^2, \sigma_Y^2$), Edgeworth expansion theory proves that the studentized bootstrap achieves **second-order asymptotic accuracy** (Type I error discrepancy $O(n^{-1})$ instead of $O(n^{-1/2})$). ---