Introduction

In many problems, we would like to test if two samples are coming from the same population. This is a common scenario in scientific studies. For instance, if a research lab invented a new drug for curing a disease, they will often do a clinical trial to test this. In a clinical trial, the researchers will recruit a couple of patients and then randomized them into two groups: a control group and a treatment group. For individuals in the control group, they will receive a placebo whereas the individuals in the treatment group will receive the actual drug. After a couple of weeks, the patients will be asked to come back and the researchers will examine disease characteristics of each individual. Finally, they will compare the responses from the control group versus the treatment group and test if the responses are different.

For simplicity, we assume that we have equal sample sizes in both the control and treatment groups. The above problem can be formulated as having two sets of observations \[\begin{align*} \mathcal{S}_X &= X_1,\cdots,X_n\sim P_X\\ \mathcal{S}_Y &= Y_1,\cdots,Y_m\sim P_Y, \end{align*}\] where \(X_i\)’s denote the response of individuals in the control group and \(Y_i\)’s are the response of individuals in the treatment group. The quantities \(n,m\) are the sample sizes. Namely, we assume that these responses are from two populations: the population of those taking placebo \(P_X\) and the population of those taking the drug \(P_Y\).

If the drug has no influence on the diseases, then there should be no difference between the two populations. Namely, \(P_X=P_Y\) if the drug has no effect. Thus, the two-sample test problem is to test \[\begin{equation} H_0: P_X = P_Y. \label{eq::h0} \end{equation}\]

Testing equation \(\eqref{eq::h0}\) might be non-trivial because both \(P_X\) and \(P_Y\) are functions (they are cumulative distribution functions). Thus, many methods will test certain characteristics of the cumulative distribution function. For instance, the \(t\)-test¹ (or \(Z\)-test) will test \[\begin{equation} H_0: \mu_X = \mu_Y, \label{eq::h1} \end{equation}\] where \(\mu_X ,\mu_Y\) are the mean of \(P_X\) and \(P_Y\), respectively. Two-sample tests based on comparing the mean is often called a mean test.

However, testing the mean has a potential problem – when there are outliers, the power will be low. Namely, even when the two samples are clearly very different from each other, we often are still unable to reject them. Consider the following very simple example: \[\begin{equation} \begin{aligned} \mathcal{S}_X & = -2,-1,-1,0,-2,-1,0,-1,-2,100\\ \mathcal{S}_Y & = 7,13,11,5,14,9,8,10,12,11. \end{aligned} \label{eq::SXY} \end{equation}\] The sample average of both samples is \(9\) but you can clearly see that except for the last observation of \(\mathcal{S}_X\), all other observations of \(\mathcal{S}_X\) are smaller than \(\mathcal{S}_Y\). Thus, the t-test or other mean tests fail in this case!

The mean tests are not ideal when the outliers are present. In the next few sections, we are going to introduce a few new methods that are much more robust to outliers.

Equal sample size: Sign Test

We first introduce a new set of variables \(Z_1,\cdots, Z_n,\) where \(Z_i = Y_i-X_i\) for each \(i=1,\cdots,n\). Under the \(H_0\) in equation \(\eqref{eq::h0}\), \(P_X=P_Y\) so the difference \(Z_i\) is coming from a distribution with \(0\) mean and \(0\) median.

Sign test utilizes the fact that the distribution of \(Z_i\) has \(0\) median. This implies that the chance of having positive value and the chance of having negative value of \(Z_i\) are the same: \(0.5\) (we assume that the two samples are from a continuous distribution). Thus, if we count the number of \(Z_i\)’s being positive, this number should follow a binomial distribution \({\sf Bin}(n,0.5)\) under \(H_0\). Let this number be \(T\).

In the example at the beginning, we have \[\mathcal{S}_Z = -9,-14,-12,-5,-16,-10,-8,-11,-14,89.\] Only one of them has a positive sign. Thus, \(T=1\) in our case.

If \(H_0\) is correct, the number of positive difference should be around \(n/2\) (half of the sample size). The number being very small (\(\approx 0\)) or very large (\(\approx n\)) are both evidences against \(H_0\). A simple way to compute a \(p\)-value using \(T\) is to calculate the probability of observing a more extreme event against \(H_0\) than what is observed. In our case, we see \(T=1\) so a more extreme event is \(T=0\) and each of them corresponds to a probability \[P(T=1) = \frac{{10\choose 1}}{2^{10}},\quad P(T=0) = \frac{{10\choose 0}}{2^{10}}.\] However, \(T\approx 10\) is also against \(H_0\) so here we need to do a two-sided test. Thus, the p-value will be two times the sum of the above two probabilities, which is \[{\sf pvalue}_1 = 2(P(T=1)+P(T=0)) = \frac{22}{1024} \approx2.15\%.\]

We can use asymptotic normality to carry this test as well. If \(H_0\) is correct, then \(T\) will follow \({\sf Bin}(n,0.5)\), which is approximately normal with \(N(n/2, n/4)\). Thus, we can use \[T_* = \frac{T - n/2}{\sqrt{n/4}}\] as our test statistic. \(T_*\) will follow a standard normal distribution so the \(p\)-value can be computed using \[{\sf pvalue}_2 = 2P(N(0,1)\geq |T_*|) \approx 1.14\%.\]

You may notice that p-values from the two methods are slightly different. The former one is the exact p-value whereas the latter one (using normal distribution) is an asymptotic approximation. When the sample size is small, the asymptotic approximation may not be very accurate.

Equal sample size: Signed-Rank Test

The sign test has a powerful extension that incorporates the information of ‘rank’. This extension is called the signed-rank test, also known as Wilcoxon signed-rank test. It incorporates a few more assumptions about the distribution of the differences \(Z_i\): symmetry. The signed-rank test requires that the distribution of \(Z_i\) is symmetric.

In our case, we have \(Z_1,\cdots,Z_n\), \(n\) differences. Let \({\sf sgn}(Z_i)\) denote the sign of \(Z_i\) (\(+1\) or \(-1\)). Now we define a new variable \(R_i = {\sf rank}(|Z_i|)\), where \(R_i=1\) means that \(Z_i\) is the smallest value of \(|Z_1|,\cdots,|Z_n|\) and \(R_i=n\) implies that \(Z_i\) is the largest value.

The signed-rank test uses a test statistic \[W = \sum_{i=1}^n {\sf sgn}(Z_i) R_i.\] Namely, \(W\) is the signed-rank sum. Under \(H_0\), the median of \(Z_i\) is \(0\), so \(W\) should be close to \(0\). We can then use the limiting distribution of \(W\) to test if the median is \(0\) or not. Note that because we assume that the distribution of \(Z\) is symmetric, the mean and median are the same.

Under \(H_0\), the variance of \(W\) is \[{\sf Var}(W)= \frac{n(n+1)(2n+1)}{6}\] and \[\frac{W }{\sqrt{{\sf Var}(W)}} \overset{D}{\rightarrow} N(0,1).\] Intuition of the variance: \(W= \sum_{i=1}^n {\sf sgn}(Z_i) R_i = \sum_{i=1}^n Y_i\) with \(\mathbb{E}(Y_i)=0\). Although \(Y_i\)’s are not independent, they are exchangeable. Then \({\sf Var}(W) = \mathbb{E}(W^2) = \sum_{i=1}^n \mathbb{E}(Y_i^2) + \sum_{j\neq \ell=1}^n {\sf Cov}(Y_j,Y_\ell)\). The covariance is \(0\) by symmetry and \(\mathbb{E}(Y_i^2) = \frac{1}{n}\sum_{i=1}^n i^2 = \frac{1}{6}(n+1)(2n+1)\).

Using the asymptotic normality, the p-value can be computed using \[{\sf pvalue} = 2P\left(N(0,1)\geq \left|\frac{(W-0) }{\sqrt{{\sf Var}(W)}}\right|\right) = 2P\left(N(0,1)\geq \sqrt{\frac{6 W^2}{n(n+1)(2n+1)}}\right) .\]

Now, going back to our example:

Thus, \(W= -35\), leading to a p-value \(2P\left(N(0,1)\geq \sqrt{\frac{6 W^2}{n(n+1)(2n+1)}}\right)\approx 7.45\%\).

Although the signed-rank test has a better power of testing \(H_0\) compared to the sign test, the sign test works for general ordinal data whereas the signed-rank test will not (think about why).

Non-equal sample size: Rank-Sum Test

The previous two methods assume that the two samples have equal sizes, which may not be true in reality. Here we introduce a novel generalization called the rank-sum test, which does not require the two samples to be equal in size. It is also called the Mann-Whitney U test, Mann-Whitney-Wilcoxon (MWW), Wilcoxon rank-sum test, and Wilcoxon-Mann-Whitney test. The rank-sum test directly tests the null hypothesis \(H_0: P_X=P_Y\).

The idea of the rank-sum test is similar to the signed-rank test but now we do not use the sign. Instead, we pull both samples together and compute the rank of each observation. Then we compute the summation of rank in one sample (the summation of rank in the other sample will also be determined). The summation of ranks will be our test statistics.

Recall that our original two samples are \(X_1,\cdots,X_n\) and \(Y_1,\cdots, Y_m\). We pull them together to form a new data \(D_1,\cdots,D_{n+m}\) such that \(D_i = X_i\) for \(i=1,\cdots, n\) and \(D_{i} = Y_{i-n}\) for \(i=n+1,\cdots, n+m\). Essentially, we are just concatenating the two samples to form a new one. Then let \(S_1,\cdots, S_{n+m}\) be the rank of \(D_1,\cdots,D_{n+m}\). Again, a smaller rank means a lower value whereas a higher rank indicates a higher value.

Before we formulate our test statistic, consider the following two statistics: \[Q_1 = \sum_{i=1}^n S_i,\quad Q_2 =\sum_{i=n+1}^{n+m} S_i.\] \(Q_1\) is the sum of ranks in \(\mathcal{S}_X\) while \(Q_2\) is the sum of ranks in \(\mathcal{S}_Y\). Note that \(Q_1+Q_2 = \frac{(n+m)(n+m+1)}{2}\).

Now we define two statistics: \[U_1 = Q_1 - \frac{n(n+1)}{2}, \quad U_2 = Q_2 - \frac{m(m+1)}{2}.\] We can use any of them as our test statistic because \(U_1+U_2 = n\cdot m\). You can view \(U_\ell\) as the “joint rank sum" (two sample jointly) minus individual group rank-sum (using only one sample). \(U_\ell\) has several nice properties: \[\mathbb{E}(U_\ell) = \frac{nm}{2}, \quad {\sf Var}(U_\ell) = \frac{nm(n+m+1)}{12}\footnote{Note that this result is assuming no ties in the data. When there are ties, there will be some modification for the variance.},\] and \[\frac{U_\ell - \mathbb{E}(U_\ell)}{\sqrt{{\sf Var}(U_\ell)}}\overset{D}{\rightarrow} N(0,1).\] Thus, we can use this asymptotic distribution to test \(H_0\). Essentially, the \(p\)-value will be \[{\sf pvalue} = 2P\left(N(0,1)\geq \left|\frac{U_1 - \mathbb{E}(U_1)}{\sqrt{{\sf Var}(U_1)}}\right|\right) =2P\left(N(0,1)\geq \left|\frac{U_\ell - \frac{nm}{2}}{\sqrt{\frac{nm(n+m+1)}{12}}}\right|\right).\]

Again, we turn to our example in equation \(\eqref{eq::SXY}\):

Using the above table, \(Q_1 = 65\) so \(U_1 = 65 - 55=10\). Thus, the \(p\)-value is \[{\sf pvalue} =2P\left(N(0,1)\geq \left|\frac{10 -50}{\sqrt{175}}\right|\right) \approx 0.25\%.\]

Assumptions (Rank-Sum Test).

• All observations are independent of each other.

• Observations within the same sample are IID.

The rank-sum test requires much weaker assumptions compared to the signed-rank test or sign test. And it also works for ordinal data. However, in general, the rank-sum test is a consistent test (power goes to infinity when the sample size goes to infinity) if \(P(X>Y)\neq P(Y>X)\). Namely, the distribution of one dataset is stochastically greater than the distribution of the other.

Median Test: One Sample Case

We end this section with a simple application of the median. We will consider testing the median in a one-sample problem. Namely, we only have one sample and our goal is to test if the median equals a certain value. Assume that we observe a sample \[\mathcal{S}_X = -2,-1,-1,0,-2,-1,0,-1,-2,100.\] Now we want to test \[H_0: m_X = 9.\] How should we test this statement?

If \(H_0\) is true, then the median is \(9\), so any observation having a value above \(9\) is about \(50\%\). We can simply compute the ratio of observations whose value is above \(9\). If \(H_0\) is correct, the ratio should be around \(50\%\). In our case, we see that there is only \(1\) out of \(10\) observations whose value is above \(9\).

Well, this could be caused purely by sampling error. So to see how unlikely our observation is, we compute its \(p\)-value. The \(p\)-value is the chance of observing a more extreme event than what we have observed, which corresponds to \(0\) out of \(10\) and \(1\) out of \(10\) observations whose value is above \(9\). \[P(0 \mbox{ out of } 10) = \frac{1}{2^{10}}, \quad P(1 \mbox{ out of } 10) =\frac{10}{2^{10}}.\] Notice that we need to do a two-sided test here (think about why), so the \(p\)-value will be 2 times the sum of the above two probabilities, which is \[{\sf pvalue} = 2\times \left(P(0 \mbox{ out of } 10)+P(1 \mbox{ out of } 10)\right) = \frac{22}{1024}\approx 2.15\%.\]

Essentially, the median test is based on the powerful property that under \(H_0\), the probability of any observation being greater than or less than the median value is \(0.5\). Then we just need to apply combinatorial methods to compute the \(p\)-value.

Here is a description of a general form of the median test. Assume we observe \(Z_1,\cdots,Z_n\sim P_Z\) and we want to test \[H_0: m_Z = m_0.\] Then we first count the number of observations such that \[T = \sum_{i=1}^n I(Z_i\geq m_0),\] where \(I(\cdot)\) is an indicator function such that if the input is a true statement, it returns \(1\) and \(0\) otherwise. Namely, \(T\) is the number of observations whose value is above or equal to the median. Because both large \(T\) and small \(T\) are the case where \(H_0\) is unlikely to be true, we use a symmetrized test statistic \[T_* = \begin{cases} T& \mbox{ if $T\leq n/2$}\\ n-T& \mbox{ if $T>n/2$}. \end{cases}\] The \(p\)-value will then be \[{\sf pvalue} = \frac{2\times \sum_{\ell=0}^{T_*} {n \choose \ell}}{2^n}.\]