Paired Comparisons and Multiple Treatments

Reference: Lehmann, E.L. (1998). Nonparametrics: Statistical Methods Based on Ranks, Prentice Hall, Upper Saddle River, NJ.

The sign test

The rank sum test has the most power when the population is relatively homogeneous. If the population is extremely heterogeneous, the variability of responses across individuals can make it hard to detect a real treatment effect. When the population is heterogeneous, sometimes we can have more power to detect a real treatment effect by dividing the population into subgroups that are less heterogeneous than the population, making comparisons within subgroups, and then combining the results across subgroups. The farthest this idea can go is to divide the population into pairs of subjects as similar as possible to each other in the values of some measured variables.

Suppose we form pairs of subjects according to the values of some variables that might affect the outcome, then randomly assign one member of each pair to treatment and one to control, by a fair coin toss, for example. Let there be n pairs of subjects. Suppose the assignment to treatment or control is independent across pairs. We will look at two approaches to testing whether treatment affects the response. The first method, the sign test, works well even when different pairs are heterogeneous. (However, see Freedman, 2005. FIX ME!) The second method, the Wilcoxon signed rank test, works better when the different pairs are homogeneous.

Let S denote the number of pairs for which the response of the treated subject is larger than the response of the subject assigned to control. Suppose for the moment that there are no ties within pairs. Under the null hypothesis that treatment makes no difference, S has a binomial probability distribution with parameters n and p=½. If treatment has a positive effect, we would expect S to be larger than n/2, so we would reject for large values of S. Critical values can be calculated from the binomial distribution. This gives the sign test based on paired samples.

Suppose now that there are ties within pairs. Consider the test statistic S' = S + ½S₀, where S₀ is the number of ties. Under the strong null hypothesis, pairs are not separated, so the number of ties S₀ is always the same: the probability distribution of S' is that of a constant plus a binomial random variable with parameters n−S₀ and p=½. Critical values for the test thus can be computed using the binomial distribution (but for n−S₀ trials). (Ties might instead be broken at random independently; it turns out that discarding ties leads to the more powerful test.)

A nice feature of the test is that within pairs, one only needs to be able to judge which subject's response is better—that is, one only needs to rank pairs of observations. That is equivalent to using just the sign of the difference between the two responses, which is why this is called the sign test. However, because the test does not use the magnitude of the difference, typically there are tests with better power when quantitative differences are available; see below. Thus the primary application of the sign test is to situations where only the sign of the difference is observed.

Like the Wilcoxon rank sum test, the sign test also makes sense in some population models, not just the randomization model. With suitable restrictions, the critical values are the same, and one can calculate the power against various alternatives.

Suppose we draw n pairs of subjects from a population of pairs, and assign one from each pair to treatment and one to control, at random with probability ½, independently across pairs. Let X_j denote the response of the subject assigned to control from the jth pair, and let Y_j denote the response of the subject assigned to treatment from the jth pair. The pairs {(X_j, Y_j): 1 ≤ j ≤ n} are identically distributed. Following Lehmann's notation, define

M(x,y) = P(X_j ≤ x and Y_j ≤ y).

Let Z_j = Y_j−X_j. The random sampling assumption implies that {Z_j} are identically disributed. Let F_Z denote their common cdf. If there is no treatment effect, X_j and Y_j are exchangeable, so the distribution of Z_j is the same as the distribution of −Z_j. That is, under the hypothesis that treatment has no effect, F_Z(z) = 1 − F_Z(−z): the distribution of Z_j is symmetric about zero. If treatment tends to increase the response, the distribution of Z_j will favor positive values.

Suppose the population of pairs is large enough compared with n that we can neglect the fact that the sample of n pairs is drawn without replacement, so that the pairs are effectively iid, not just identically distributed. Then {Z_j} are also iid. This is the population model; the null hypothesis is that the distribution of the Z_j is symmetric about zero; the alternative hypothesis is that it favors positive values. Suppose that F_Z is continuous so that the probability that Z_j=0 is 0. Then, under the null hypothesis,

P(Z_j > 0) = P(Z_j < 0 ) = ½.

The Z_j's are independent, so the null probability distribution of the number of strictly positive Z_j's is binomial with parameters n and ½. Under the alternative that the probability that Z_j>0 is larger than ½, the number of positive Z_j's will tend to be larger; its probability distribution is binomial with parameters n and p=1−F_Z(0)>½. Thus the sign test derived above is applicable; the power of the test is clearly monotone in p, so the test is unbiased against alternatives p>½. In the shift model F_Z(z) = F₀(z−d), where F₀ is the cdf of a continuous distribution that is symmetric about zero, the level a sign test of the hypothesis that the distribution of the Z's is symmetric about zero against the alternative that positive values are favored is also a level a test of the null hypothesis d≤0 against the alternative d>0.

The number of positive Z_j also has a binomial distribution with parameters n and ½ in other situations, for example, when the Z_j's are independent, not necessarily identically distributed, but each has probability ½ of being positive or negative (for example, each has a probability distribution that is symmetric about zero). Then the sign test is a level a test of the null hypothesis against the alternative that the distributions of the Z_j's favor positive values. See Divenyi, Stark and Haupt, 2005. Decline of Speech Understanding and Auditory Thresholds in the Elderly, J. Acoustical Soc. Am., 118, 1089–1100, for examples.

Note that the kinds of violations of the null hypothesis that the sign test is sensitive to are those for which the median is positive—other violations of symmetry result in power equal to the significance level.

As before, when ties are possible, the null distribution of the number of strictly positive Z_j's given the number of zero Z_j's is binomial with parameters n−(#ties) and ½, so we can use the same approach we used for the Wilcoxon rank sum test to get conditional tests in the presence of ties.

The sign test also can be used to test hypotheses about the median of a continuous distribution from iid observations {Z'_j: 1 ≤ j ≤ n}, because then the variables Z_j = Z'_j − d are iid with P(Z_j < 0) = P(Z_j >0) if d is their median. This is essentially a special case of the approach we used to find confidence intervals for quantiles from iid observations.

The Wilcoxon Signed-Rank Test

The sign test ignores all information but the signs of the differences between treatment and control. That allows the pairs to be completely different from each other. When there is some homogeneity across pairs, it can help to retain some quantitative information about the response: the difference between the control and treatment responses for each pair, or at least the rank of the difference. Suppose there are n pairs of subjects, as before, with one member of each pair assigned to treatment and one to control, at random, with probability 50%, independently across pairs. We are interested in the alternative hypothesis that treatment increases the response. Let t be the number of pairs for which the treatment response "wins"—exceeds the control response—and let c = n−t denote the number of wins for control. For the moment, we assume that there are no ties among the differences, and that none of the differences is zero. We rank the response differences by absolute value, then give each rank the sign of the corresponding difference. Let R₁ < R₂ < …R_c be the ranks with negative signs, and let S₁ < S₂ < …< S_t be the ranks with positive signs. Then

{R₁, R₂, …R_c, S₁, S₂, …, S_t } = {1, 2, …, n}.

There are 2ⁿ possible sign combinations for the ranks, ±1, ±2, …, ±n. Specifying a combination of signs is equivalent to listing the positive-signed ranks {S₁, S₂, …, S_t }. If treatment has no effect, the rank of the response difference of each pair is fixed, but the signs of the response differences are iid, positive with probability ½ and negative with probability ½. Under the null hypothesis, each of the 2ⁿ assignments of signs is equally likely; each corresponds to a number N₊ = x ≤n of positive differences, and a subset {S₁, …, S_x} of x of the n ranks. Thus, under the strong null hypothesis,

P(N₊ = x; S₁=s₁, S₂=s₂, …, S_x=s_x ) = 1/2ⁿ.

On the other hand, if treatment tends to increase the response, we would expect more and larger treatment differences to be positive. That is, we would expect

V_s = S₁ + S₂ + … + S_t

to be large. Note that under the strong null hypothesis, t is random. The test statistic V_s leads to the Wilcoxon signed rank test: reject the strong null hypothesis if V_s ≥ C, where C is chosen so that the test has the right significance level. The critical values for the Wilcoxon signed rank test are not the same as the critical values for the Wilcoxon rank sum test.

The Wilcoxon signed rank test tends to be more powerful than the sign test. Critical values are in Lehmann (1998); they can also be estimated by simulation. There is a normal approximation as well: under the strong null hypothesis,

E(V) = N(N+1)/4; Var(V) = N(N+1)(2N+1)/24.

See Lehmann (1998). The Wilcoxon signed rank test can be extended to allow for ties by using midranks rather than ranks; the null distribution of the test statistic then depends on the configuration of ties, as it does for the Wilcoxon rank sum test and the sign test. Simulation is straightforward. Here is sample code in R:

        midRanks <- function(x) {
        # P.B. Stark, 10/5/05.  www.stat.berkeley.edu/~stark
        #
        # calculates vector of midranks
            sapply(x, function(z){(sum(x <= z) + sum(x < z) + 1)/2});
        }

        signedRankSum <- function(x,y) {
        # P.B. Stark, 10/5/05.  www.stat.berkeley.edu/~stark
        #
        # signed rank sum for control responses x and treatment
        # responses y.
        # Calls midRanks().
            r <- midRanks(abs(y-x));
            sum(r[(y-x)>0]);
        }

        wilcoxonSignedRankSim <- function(x, y, iter) {
        # P.B. Stark, 10/5/05.  www.stat.berkeley.edu/~stark
        #
        # simulated one-sided P-value for the Wilcoxon signed rank
        # test using midranks.
        #
        # x: vector of control responses
        # y: vector of treatment responses
        # iter: number of replications
        #
        # Calls midRanks()
            r <- midRanks(abs(y-x));
            v <- sum(r[(y-x) > 0]);
            sum(replicate(iter, sum(r[runif(x) > 0.5])) >= v)/iter
        }

Combining dissimilar blocks of homogeneous subjects

We have been working with pairs, which are the finest-grained blocks possible. Pairing is helpful when the population is extremely heterogeneous; it lets us control for differences that might influence the effect of treatment, provided we can find pairs of individuals who are similar with respect to those differences. Earlier (in the discussion that led to the Wilcoxon rank sum test), we assumed that individuals were homogeneous. Between these two extremes, we might have dissimilar blocks of homogeneous individuals, which we consider here.

Suppose we have b blocks. The jth block contains N_j subjects. Of those, n_j are assigned at random to treatment, and the remaining m_j = N_j − n_j are assigned to control. There are

N = N₁ + N₁ + … N_b

subjects in all, of which

n = n₁ + n₁ + … n_b

are assigned at random (within blocks) to treatment and

m = m₁ + m₁ + … m_b

are assigned to control. In this scheme, there are

N_assign = _N1C_n1 × _N2C_n2 × … × _NbC_nb

equally likely possible assignments of subjects to treatment and control. To allow for heterogeneity across blocks, but combine evidence from different blocks, we could calculate the Wilcoxon rank sum statistic separately for each block, then form a weighted sum of the statistics across blocks. Let W_j be the Wilcoxon rank sum calculated for block j by itself. One choice of test statistics for the collection of blocks is

W = W₁/(N₁+1) + W₂/(N₂+1) + … + W_b/(N_b+1).

This choice downweights blocks with larger numbers of subjects. The null distribution of this statistic can be calculated using the fact that under the null hypothesis, the N_assign assignments of ranks are equally likely; alternatively, the null distribution can be simulated using random sampling.

More than two treatments: the Kruskal-Wallis Test

Suppose there are s mutually exclusive categories, and let

n = n₁ + n₂ + … + n_s.

How many ways can the n objects be placed into the s categories, so that n₁ are assigned to cetegory 1, n₂ are assigned to category 2, etc.? There are _nC_n1 ways of selecting n₁ of the objects to assign to the first category; once that is done, there are _n−n1C_n2 ways of selecting n₂ of the remaining n−n₁ objects to assign to the second category, and so on. Finally, only n_s objects remain for the sth category. The total number of possible assignments is thus

_nC_{n1,n2,…, ns} = _nC_n1 × _n−n1C_n2 × … × _{n−n1−n2 − …− ns−1}C_ns = n!/(n₁!n₂! … n_s!)

This is the multinomial coefficient.

Suppose we wish to compare s different treatments using n subjects. We select n₁ at random and assign them to the first treatment; we select n₂ at random and assign them to the second treatment, and so on. Suppose that there are no ties among the responses. Under the strong null hypothesis that all the treatments are the same, all _nC_{n1,n2, … , ns} assignments of the n ranks among the treatments are equally likely.

Consider the ranks of the pooled responses, and suppose for the moment that there are no ties (we will generalize to midranks presently). Let R^j denote the mean of the ranks of the responses to treatment j. Under the null hypothesis, the expected value of R^j is the average of all the ranks (if there are no tied observations, this is (n+1)/2; more generally, the expected mean rank for each treatment is the average of the midranks of the pooled responses). Let r denote this expected value. We can test the null hypothesis using a weighted sum of squares of the deviations of the treatment mean ranks from their expected value:

K = (12/(n(n+1))) ( n₁(R¹−r)² + n₂(R²−r)² + …; + n_s(R^s−r)² ).

We reject the null hypothesis if K≥C, where C is chosen to get the right significance level. If there are no ties, under the null hypothesis K converges in distribution to a chi-square random variable with s−1 degrees of freedom as n₁,…,n_s grow. If there are ties, a different definition is usually used; see Lehmann (1998).

Here is R code to simulate the null distribution of K computed using midranks.

        kruskalWallisSim <- function(x, iter) {
        # P.B. Stark, 10/5/05.  www.stat.berkeley.edu/~stark
        #
        # Simulated P-value for the Kruskal-Wallis test using midranks.
        # Does not normalize the Kruskal-Wallis statistic.
        #
        # Calls midRanks()
        #
        # x: list of arrays of observed responses. The number
        #    of arrays is the number of treatments.
        # iter: number of iterations in the simulation
        #
            names(x) <- 1:length(x);
            resp <- stack(x);
            resp[,1] <- midRanks(resp[,1]);
            resp[,2] <- factor(resp[,2]);
            exRank <- mean(resp[,1]);          # expected mean rank for each treatment
            sqComp <- function(z){length(z)*(mean(z)-exRank)^2}
            ts <- sum(tapply(resp[,1],resp[,2],sqComp))
            sum(replicate(iter,(sum(tapply(sample(resp[,1]),resp[,2],sqComp))>=ts)))/iter;
        }

Other topics

2 × t contingency tables: See Lehmann (1998).