Vidakovic B. Statistics for Bioengineering Sciences: With Matlab and WinBugs Support

Подождите немного. Документ загружается.

12.5 Wilcoxon Sum Rank Test and Wilcoxon–Mann–Whitney Test 489

%-------------------------- alternatives -----------------

% alt == 0 for two sided;

% alt == -1 for mu1 < mu2; alt == 1 for mu1 > mu2.

if alt == 0

p = 2

normcdf(- abs( Z ) - cc);

elseif alt == -1

p = normcdf( Z - cc );

elseif alt == 1

p = normcdf( - Z - cc );

else

error(’Input "alt" should be either 0, -1, or 1.’)

end

%---------------------------------------------------------

probes against various biomolecules is very important in the furthering of ba-

sic studies of cellular biology and pathology. One way to improve the binding

properties of these probes is to have multiple binding domains or ligands on

the surface of the probes. The more ligands, the greater the chance of binding.

One issue that comes up with such probes is whether, due to the multiple lig-

ands, they will cause aggregation of target molecules within the cell. In order

to show that new probes do not induce aggregation, researchers in the lab of

Dr. Phil Santangelo at Georgia Tech compared the number of granules using a

monovalent and a tetravalent probe in a cell plated at the same time and un-

der the same biological conditions. The numbers of granules detected by each

probe were recorded, and the researchers were interested to see if there were

real differences between the numbers.

Using WSuRT we test the hypothesis of equality of distributions and, sub-

sequently, all theoretical moments at the signiﬁcance level

α =0.05.

monovalent = [117 92 84 213 89 76 96 104 114 142 ...

122 154 124 65 129 67 100 127 63 82 ...

114 93 117 83 82 83 111 78 92 91];

tetravalent = [103 78 155 107 113 75 74 80 120 112 ...

158 72 81 124 110 90 64 74 110 149 ...

97 70 105 94 110 93 115 114 110 95];

[sumranks1, tstat, pval] = wsurt(monovalent, tetravalent, 0)

% sumranks1 = 925.5000

% tstat = 0.1553

% pval = 0.8737

The statistical evidence that the two samples come from the same distribu-

tion is overwhelming since the p-value of 0.8737 exceeds the signiﬁcance level

α =0.05.

Example 12.4. Nanoscale Probes. The development of intracellular nanoscale

490 12 Distribution-Free Tests

12.6 Kruskal–Wallis Test

The Kruskal–Wallis (KW) test is a generalization of the WSuRT. It is a non-

parametric test used to compare three or more samples. It is used to test

the null hypothesis that all populations have identical distribution functions

against the alternative hypothesis that at least two of the samples differ only

with respect to location (median), if at all.

The KW test is an analog to the F-test used in one-way ANOVA. While

ANOVA tests depend on the assumption that all populations under compar-

ison are independent and normally distributed, the KW test places no such

restriction. Suppose the data consist of k independent random samples with

sample sizes n

,... , n

. Let n = n

+···+n

Sample 1 X

, X

, . . . X

1,n

Sample 2 X

, X

, . . . X

2,n

Sample k

−1 X

k−1,1

, X

k−1,2

, . . . X

k−1,n

k−1

Sample k X

, X

, . . . X

k,n

Under the null hypothesis, we can claim that all of the k samples are from

a common population.

The expected sum of ranks for sample i,

E(R

), would be n

times the ex-

pected rank for a single observation. That is, n

(n +1)/2, and the variance can

be calculated as

Var (R

) = n

(n +1)(n −n

)/12. One way to test H

is to calcu-

late R

j=1

r(X

i j

) – the total sum of ranks in sample i. The statistic

i=1

−

(n +1)

(12.1)

will be large if the samples differ, so the idea is to reject H

if (12.1) is “too

large.” However, its distribution is quite messy, even for small samples, so we

can use the normal approximation

−E(R

)

Var (R

)

appr

∼ N (0,1) ⇒

i=1

(

−E(R

)

Var (R

)

appr

∼ χ

−1

where the

-statistic has only k −1 degrees of freedom due to the fact that

only k

−1 ranks are unique.

Based on this idea, Kruskal and Wallis (1952) proposed the test statistic

i=1

−

n(n +1)

, (12.2)

where

12.6 Kruskal–Wallis Test 491

n −1

i=1

j=1

r(X

i j

)

−

n(n +1)

If there are no ties in the data, (12.2) simpliﬁes to

n(n +1)

i=1

−

(n +1)

. (12.3)

Kruskal and Wallis (Fig. 12.2a,b) showed that this statistic has an approxi-

mate

-distribution with k −1 degrees of freedom.

(a) (b) (c)

Fig. 12.2 (a) William Henry Kruskal (1919–2005), (b) Wilson Allen Wallis (1912–1998), and

The MATLAB routine kw.m implements the KW test using a vector to

represent the responses and another to identify the population from which the

response came. Suppose we have the following responses from three treatment

groups:

(1,3,4), (3,4,5), (4,4,4,6,5).

The code for testing the equality of locations of the three populations computes

a p-value of 0.1428.

data = [ 1 3 4 3 4 5 4 4 4 6 5 ];

belong = [ 1 1 1 2 2 2 3 3 3 3 3 ];

[H, p] = kw(data, belong)

% [H, p] = 3.8923 0.1428

Example 12.5. The following data are from a classic agricultural experiment

measuring crop yield in four different plots. For simplicity, we identify the

treatment (plot) using the integers {1,2,3,4}. The third treatment mean mea-

sures far above the rest, and the null hypothesis (the treatment means are

equal) is rejected with a p-value less than 0.0002.

data= [83 91 94 89 89 96 91 92 90 84 91 90 81 83 84 83 ...

88 91 89 101 100 91 93 96 95 94 81 78 82 81 77 79 81 80];

492 12 Distribution-Free Tests

belong = [1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 ...

3 3 3 3 3 3 3 3 4 4 4 4 4 4 4];

[H, p] = kw(data, belong)

% H = 20.3371

% p = 1.4451e-004

Kruskal–Wallis Pairwise Comparisons. If the KW test detects treatment

differences, we can determine if two particular treatment groups (say, i and j)

are different at level

α if

−

n−k,1−α/2

(n −1 −H

)

n −k

. (12.4)

Example 12.6. Since in Example 12.5 we found the four crop treatments sig-

niﬁcantly different, it would be natural to ﬁnd out which ones seem better and

which ones seem worse. In the table below, we compute the statistic

−

(n−1−H

)

n−k

for every combination of 1 ≤ i 6= j ≤4 and compare it to t

30,0.975

=2.042.

(i, j)

1 2 3 4

1 0 1.856 1.859 5.169

1.856 0 3.570 3.363

1.859 3.570 0 6.626

5.169 3.363 6.626 0

This shows that the third treatment is the best, but not signiﬁcantly differ-

ent from the ﬁrst treatment, which is second best. Treatment 2, which is third

best, is not signiﬁcantly different from treatment 1, but is different from treat-

ments 4 and 3.

12.7 Friedman’s Test

Friedman’s test is a nonparametric alternative to the randomized block de-

sign (RBD) in regular ANOVA. It replaces the RBD when the assumptions

of normality are in question or when the variances are possibly different from

population to population. This test uses the ranks of the data rather than their

12.7 Friedman’s Test 493

raw values to calculate the test statistic. Because the Friedman test does not

make distribution assumptions, it is not as powerful as the standard test if

the populations are indeed normal.

Milton Friedman (Fig. 12.2c) published the ﬁrst results for this test, which

was eventually named after him. He received the Nobel Prize in economics

in 1976, and one of the listed breakthrough publications was his article “The

use of ranks to avoid the assumption of normality implicit in the analysis of

variance,” published in 1937 (Friedman, 1937).

Recall that the RBD design requires repeated measures for each block at

each level of treatment. Let X

i j

represent the experimental outcome of subject

(or “block”) i with treatment j, where i

=1,..., b, and j = 1,..., k.

Treatment

Block 1 2 ... k

1 X

... X

2 X

... X

b X

... X

To form the test statistic, we assign ranks {1,2,..., k} to each row in the

table of observations. Thus the expected rank of any observation under H

is (k +1)/2. We next sum all the ranks by columns (by treatments) to obtain

r(X

i j

), 1 ≤ j ≤ k. If H

is true, the expected value for R

is E(R

) =

b(k +1)/2. The statistic

j=1

−

b(k +1)

is an intuitive formula to reveal treatment differences. It has expectation

bk(k

−1)/12 and variance k

b(b −1)(k −1)(k +1)

/72. Once normalized to

bk(k +1)

j=1

−

b(k +1)

, (12.5)

it has moments

E(S) = k −1 and Var (S) = 2(k −1)(b −1)/b ≈ 2(k −1), which

coincide with the ﬁrst two moments of

−1

. Higher moments of S also ap-

proximate well those of

−1

when b is large.

In the case of ties, a modiﬁcation to S is needed. Let C

= bk(k +1)

/4 and

∗

r(X

i j

)

. Then,

k −1

∗

−bk(k +1)

j=1

−bC

(12.6)

is also approximately distributed as

−1

494 12 Distribution-Free Tests

Although the Friedman statistic makes for a sensible, intuitive test, it

turns out there is a better one to use. As an alternative to S (or S

), the test

statistic

(b −1)S

b(k −1) −S

is approximately distributed as F

k−1,(b−1)(k−1)

, and tests based on this approx-

imation are generally superior to those based on chi-square tests that use S.

For details on the comparison between S and F, see Iman and Davenport

(1980)

Example 12.7. In an evaluation of vehicle performance, six professional drivers

(labeled I, II, III, IV, V, VI) evaluated three cars (A, B, and C) in a randomized

order. Their grades concern only the performance of the vehicles and suppos-

edly are not inﬂuenced by the vehicle brand name or similar exogenous infor-

mation. Here are their rankings on a scale of 1 to 10:

Car

I II III IV V VI

A 7 6 6 7 7 8

8 10 8 9 10 8

9 7 8 8 9 9

To use the MATLAB procedure

friedmanGT the input data should be format-

ted as a matrix in which the rows are blocks and the columns are treatments.

data = [7 8 9; 6 10 7; 6 8 8; ...

7 9 8; 7 10 9; 8 8 9];

[S,F,pS,pF] = friedmanGT(data)

%S = 8.2727

%F = 11.0976

%pS = 0.0160

%pF = 0.0029

% pF p-value is more reliable

Friedman’s Pairwise Comparisons. If the p-value is small enough to war-

rant multiple comparisons of treatments, we consider two treatments i and j

to be different at level

α if

−R

(b−1)(k−1),1−α/2

2 ·

∗

−

(b −1)(k −1)

. (12.7)

12.8 Walsh Nonparametric Test for Outliers* 495

Example 12.8. From Example 12.7, the three cars A, B, and C are considered

signiﬁcantly different at test level

α = 0.01 (if we use the F-statistic). We can

use the function

friedmanpairwise(x,i,j,alpha) to make a pairwise com-

parison between treatments i and j at level

alpha. The output = 1 if treat-

ments i and j are different; otherwise it is 0. The Friedman pairwise compari-

son reveals that car A is rated signiﬁcantly lower than both cars B and C, but

cars B and C are not considered to be different.

12.8 Walsh Nonparametric Test for Outliers*

Suppose that r outliers are suspect, where r ≥1 and ﬁxed. Order observations

,... , X

and obtain the order statistic X

(1)

≤ X

(2)

≤ ··· ≤ X

(n)

, and set the

signiﬁcance level

α. The test is nonparametric and no assumptions about the

underlying distribution are made.

We will explain the steps and provide a MATLAB implementation, but

readers interested in details are directed to Walsh (1962). The Walsh recipe

has the following steps, following Madansky (1988).

STEP 1: Calculate c = b

2nc, where bxc is the largest integer ≤ x,

=1/α, k = r +c, and a =(1 +b

c−b

c−1

)/(c −b

−1).

STEP 2: The smallest r values X

(1)

≤···≤ X

(r)

are outliers if

= X

(r)

−(1 +a)X

(r+1)

+aX

<0,

and the largest r values X

(n−r+1)

≤···≤ X

(n)

are outliers if

= X

(n−r+1)

−(1 +a)X

(n−r)

+aX

n−k+1

>0.

The sample size has to satisfy n

≥

(1 +

)

. To achieve α = 0.05, a sample

size of at least 221 is needed. As an outcome, one either rejects no outliers,

rejects the smallest r, the largest r, or even both the smallest and largest r,

thus potentially rejecting a total of 2r observations. Here is an example of

Walsh’s procedure implemented by m-function

walshnp.m.

dat = normrnd(0, 2, [300 1]);

data = [-11.26; dat; 13.12];

walshnp(data) %default r=1, alpha=0.05

%Lower outliers are: -11.26

%Upper outliers are: 13.12

%ans = -11.2600 13.1200

496 12 Distribution-Free Tests

12.9 Exercises

12.1. Friday the 13th. The following data set is part of the larger study from

Scanlon et al. (1993) titled “Is Friday the 13th bad for your health?” The

data analysis in this paper addresses the issues of how superstitions re-

garding Friday the 13th affect human behavior. The authors reported and

analyzed data on trafﬁc accidents for Friday the 6th and Friday the 13th

between October 1989 and November 1992. The data consist of the num-

ber of patients accepted in SWTRHA (South West Thames Regional Health

Authority, London) hospital on the dates of Friday the 6th and Friday the

13th.

# of accidents

Year, month Friday 6th Friday 13th Sign

1989, October 9 13 −

1990, July 6 12 −

1991, September 11 14 −

1991, December 11 10 +

1992, March 3 4 −

1992, November 5 12 −

Use the sign test at the level α =10% to test the hypothesis that the “Friday

the 13th effect” is present. The m-ﬁle

signtst.m could be applied.

12.2. Reaction Times. In Exercise 10.22 the paired t-test was used to assess

the differences between reaction times to red and green lights. Repeat this

analysis using both sign test and WSiRT.

12.3. Simulation. To compare the t-test with the Wilcoxon signed-rank test, set

up the following simulation in MATLAB:

(1) Generate n

=20 observations from N (0,1) as your ﬁrst sample X .

(2) Form

Y = X + randn(size(x)) + 0.5 as your second sample paired with

the ﬁrst.

(3) For the test of H

: µ

= µ

versus H

: µ

< µ

, perform the t-test at

α =0.05.

(4) Run the Wilcoxon signed-rank test.

(5) Repeat this simulation 1000 times and compare the powers of the tests

by counting the number of times H

was rejected.

12.4. Grippers. Measurements of the left- and right-hand gripping strengths of

ten left-handed writers are recorded.

Person 1 2 3 4 5 6 7 8 9 10

Left hand (X )

140 90 125 130 95 121 85 97 131 110

Right hand (Y )

138 87 110 132 96 120 86 90 129 100

12.9 Exercises 497

(a) Does the data provide strong evidence that people who write with their

left hand have a greater gripping strength in the left hand than they do in

their right hand? Use the Wilcoxon signed-rank test and

α =0.05.

(b) Would you change your opinion on the signiﬁcance if you used the paired

t-test?

12.5. Iodide and Serum Concentration of Thyroxine. The effect of iodide

administration on serum concentration of thyroxine (T

) was investigated

in Vagenakis et al. (1974). Twelve normal volunteers (9 male and 3 female)

were given 190 mg iodide for 10 days. The measurement X is an average of

in the last 3 days of administration while Y is the mean value in three

successive days after the administration stopped.

Subject 1 2 3 4 5 6 7 8 9 10 11 12

Iodide (X )

7.9 9.1 9.2 8.1 4.2 7.2 5.4 4.9 6.6 4.7 5.2 7.3

Control (Y )

10.2 10.2 11.5 8.0 6.6 7.4 7.7 7.2 8.2 6.2 6.0 8.7

Assume that the difference D = X −Y has a symmetric distribution. Com-

pare the p-values of WSiRT and the paired t-test in testing the hypothesis

that the mean of control measurements exceed that of the iodide measure-

ments.

12.6. Weightlifters. Blood lactate levels were determined in a group of ama-

teur weightlifters following a competition of 10 repetitions of 5 different

upper-body lifts, each at 70% of each lifter’s single maximum ability. A

sample of 17 randomly selected individuals was tested: 10 males and 7 fe-

males. The following table gives the blood lactate levels in female and male

weightlifters, in units of mg/100 ml of blood.

Gender N Blood lactate

Female 7 7.9 8.2 8.7 12.3 12.5 16.7 20.2

Male 10 5.2 5.4 6.7 6.9 8.2 8.7 14.2 14.2 17.4 20.3

Test an H

that blood lactate levels were not different between the two gen-

ders at the level

α =0.05. Assume that the data, although consisting of mea-

surements of continuous variables, are not distributed normally or vari-

ances are possibly heteroscedastic (nonequal variances). In other words,

the traditional t-test may not be appropriate.

12.7. Cartilage Thickness in Two Osteoarthritis Models. Osteoarthritis

(OA), characterized by gradual degradation of the cartilage extracellular

matrix, articular cartilage degradation, and subchondral bone remodeling,

is the most common degenerative joint disease in humans. One of the most

common ways for researchers to study the progression of OA in a controlled

manner is to use small animal models. In a study conducted in the labo-

ratory of Robert Guldberg at the Georgia Institute of Technology, a group

of 14 rats was randomly divided into two subgroups, and each subgroup

was subjected to an induced form of OA. One of the subgroups (n

=7) was

498 12 Distribution-Free Tests

subjected to a chemically induced form of OA via an intra-articular injec-

tion of monosodium iodoacetate (MIA), while the other subgroup (n

= 7)

was subjected to a surgically induced form of OA by transecting the me-

dial meniscus (MMT). In this study, all rats had OA induced by either MIA

or MMT on the left knee, with the right knee serving as a contralateral

control.

The objectives of this study were to quantify changes in cartilage thickness

of a selected area of the medial tibial plateau, compare thickness values

between the treatments and the controls, and determine if OA induced by

MIA produced different results from OA inducted by MMT. Cartilage thick-

ness values were measured 3 weeks after the treatments using an ex vivo

micro-CT scanner; the data are provided in the table below.

Rat MIA Treated Control Rat MMT Treated Control

1 0.1334 0.2194 8 0.2569 0.2726

2 0.1214 0.1929 9 0.2101 0.2234

3 0.1276 0.1833 10 0.1852 0.2216

4 0.1152 0.1879 11 0.1798 0.1905

5 0.1047 0.2529 12 0.1049 0.1444

6 0.1312 0.2527 13 0.2649 0.2841

7 0.1222 0.2595 14 0.2383 0.2731

(a) Using the Wilcoxon signed-rank test, test the hypothesis that the differ-

ence between MIA-treated rats and its control thickness is signiﬁcant.

(b) Find the difference in thickness (treatment–control) for MIA and MMT

and test the hypothesis that they are the same. Apply a nonparametric

version of a two-sample t-test on the differences.

Conduct both tests at a 5% signiﬁcance level.

12.8. A Claim. Professor Scott claims that 50% of his students in a big class get

a ﬁnal score of 90 and higher.

A suspicious student asks 17 randomly selected students from Professor

Scott’s class and they report the following scores:

80 81 87 94 79 78 89 90 92 88 81 79 82 79 77 89 90

Test the hypothesis that Professor Scott’s claim does not conform to the

evidence, i.e., that the 50th percentile (0.5-quantile, median) is different

than 90. Use

α =0.05.

12.9. Claustrophobia. Sixty subjects seeking treatment for claustrophobia are

independently sorted into two groups, the ﬁrst of size n

= 40 and the sec-

ond of size m

= 20. The members of the ﬁrst group each individually re-

ceive treatment A over a period of 15 weeks, while those of the second

group receive treatment B. The investigators’ directional hypothesis is that

treatment A will prove to be more effective. At the end of the experimental

treatment period, the subjects are individually placed in a series of claus-

trophobia test situations, knowing that their reactions to these situations