Resampling Techniques
Permutation Tests
U G A R T E , M . D . ( * )
( * ) D E P A R T A M E N T O D E E S T A D Í S T I C A E I . O . , U N I V E R S I D A D P U B L I C A D E N A V A R R A , P A M P L O N A , S P A I N
E-MAIL: LOLA@UNAVARRA.ES
Material from the book Probability and Statistics with R
by Ugarte, Militino, and Arnholt. Chapman and Hall/CRC, 2008
Brno, December 2007 Lola Ugarte
Introduction
Permutation tests are computationally intensive techniques that actually predate
computers.
Until recently, permutation tests were more of a theoretical ideal than a useful
technique. With the advent of high-powered computers, permutation tests have
moved out of the abstract into the world of the practical.
The permutation test is examined in only one context here: the two-sample prob-
lem.
Brno, December 2007
1
Lola Ugarte
Introduction -Cont.
The fundamental idea behind the permutation test is that if there are no differences
between two treatments, all data sets obtained by randomly assigning the
data to the two treatments have an equal chance of being observed.
Permutation tests are especially advantageous when working with small samples
where verification of assumptions required for tests such as the pooled
Mest are difficult.
Brno, December 2007
2
Lola Ugarte
Introduction -Cont.
To test a hypothesis with a permutation test:
1. Choose a test statistic 0 that measures the effect under study. Note that
certain statistics will have more power to detect the effect of interest than
others.
2. Create the sampling distribution that the test statistic in Step 1 would have if
the effect is not present in the population (the null hypothesis is true).
Brno, December 2007
3
Lola Ugarte
3. Find the "observed test statistic" in the sampling distribution from Step 2.
Observed values in the extremes of the sampling distribution suggest that the
effect under study is "real". In contrast, observed values in the main body of
the sampling distribution imply that the effect is likely to have occurred by
chance.
4. Calculate the p-value based on the observed test statistic. This may be
^(|0| > \0obs\), in a two-sided test
P{0 < 9obs), in a left one-sided test
P0>9obs), in a right one-sided test
Brno, December 2007
4
Lola Ugarte
The two-sample problem
Suppose two independent random samples z = {zh..., zn} and y = {yh... ,ym}
are drawn from possibly different probability distributions F and G.
The question of interest is whether F = G.
Assuming H0 : F = G is true, it is possible to create the permutation sampling
distribution for some statistic of interest, 9.
Let N equal the combined sample size n + m.
Brno, December 2007
5
Lola Ugarte
The two-sample problem
* Let v = {vi,v2, ˇ ˇ ˇ ,VN}, be the combined and ordered vector of values, (the
first n values correspond to the first sample and the rest to the second).
* Let g = {</i, </2 ˇ ˇ ˇ I9N} be a vector indicating the group membership for v.
Since there are n ZÍ'S and m y/s in g, there are (^) possible ways of partitioning
N elements into two subsets of sizes n and m.
* Consequently, under the null hypothesis that F = G, the vector g has probability
1 / Q of equalling any one of its possible values. That is, all combinations
of Zi's and y/s are equally likely if F = G.
Brno, December 2007
6
Lola Ugarte
The two-sample problem
Suppose HA: F > G and 0 = z -y, then the exact p-value is found by finding
#{# > Oobs}
For all but relatively trivial sized samples n and m, the number (^) will be huge,
making the enumeration of all possible samples of the statistic of interest a
monumental task.
For example, consider that if n = 10 and m = 10 that complete enumeration would
require listing (f0) = 184,756 possible outcomes.
Brno, December 2007
7
Lola Ugarte
The two-sample problem
Consequently, an approximation to the exact p-value is often obtained by resampling
WITHOUT replacement the original data some "large" number of
times, B, which is usually at least 999, and approximating the p-value with
n vaiMP - Ľ1
+ #{^ > <M1 _ [i + Ef=i He > eobs}]p
~value
- (äTTj - (äTij
When the sampling is done with replacement, a bootstrap test is performed. Bootstrap
tests are somewhat more general than permutation tests since they apply
to a wider class of problems; however, they do not return "exact" p-values.
Brno, December 2007
8
Lola Ugarte
Example (Original data: Ott and
Mendenhall, 1985)
The data for this problem will be stored in a data frame called Ratbp. Researchers
wanted to know whether a drug was able to reduce blood pressure of rats. Twelve
rats were chosen and the drug was administered to six rats, the treatment group,
chosen at random. The other six rats, the control group, received a placebo. The
drops in blood pressure (mmHg) for the treatment group (with probability distribution
F) and the control group (with probability distribution G) are stored in the
variables Treat (z) and Cont (y), respectively. Note that positive numbers indicate
blood pressure decreased while negative numbers indicate that it rose.
Brno, December 2007
9
Lola Ugarte
Example
Under the null hypothesis, H0 : F = G, the data come from a single population
(i.e., we are not able to prove the drug efficacy). The question of interest is: How
likely are differences as extreme as those observed between the treatment and
control groups to be seen if the null hypothesis is correct? Use 0 = ž - y as the
statistic of interest and compute
* (a) the exact permutation p-value
* (b) an estimated permutation p-value based on 499 permutation replications
* (c) an estimated bootstrap p-value based on 499 bootstrap replications
Brno, December 2007
10
Lola Ugarte
Solution (a) Exact permutation p-value
The test statistic of interest (Step 1) has been specified to be 9 = z - y. Finding
the p - value requires creation of sampling distributions with different methods.
First, create Ratbp data frame.
Combine the treatment (Treat) and control (Cont) data in a single variable called
Blood.
Brno, December 2007
11
Lola Ugarte
## To introduce the data in R
Treat<-c(6 9,24,63,87.5,77.5,40)
Cont<-c(9,12,36,77.5,-7.5,32.5)
Ratbp<-data.frame(Treat, Cont) # to create a data frame
> Ratbp
Treat Cont
1 69.0 9.0
2 24.0 12.0
3 63.0 36.0
4 87.5 77.5
5 77.5-7.5
6 40.0 32.5
Blood<-c(Treat,Cont)
> Blood
[1] 69.0 24.0 63.0 87.5 77.5 40.0 9.0 12.0 36.0 77.5 -7.5 32.5
Brno, December 2007
12
Lola Ugarte
Solution (a) Exact permutation p-value
* We compute all of the possible combinations of the indices of Blood: (g2
) =
924 and put them in a vector called pdT6.
* To do that we may install in R the package BSDA. Alternatively, we may load
the function Combinations.R.
* Now we construct a matrix called Comb of dimension 924 x 12 where each
row contains the values of Blood.
* We also construct a vector Theta with all of the possible values of 0 = ž - y,
i.e., a value for each possible combination.
Brno, December 2007
13
Lola Ugarte
library(BSDA)
pdT6<-t(Combinations(12, 6) )
#Let us show the first three rows of pdT6:
pdT6[1:3,]
> pdT6[1:3, ]
N N N N N
[1, ] 1 2 3 4 5 6
[2, ] 1 2 3 4 5 7
[3, ] 1 2 3 4 6 7
## Now we construct the matrix with 924 rows -all the same-
Comb<-matrix(rep(c(Treat,Cont),924),ncol=12,byrow=T)
Comb[1:3,] #three first rows of Comb
> Comb[1:3,]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 69 24 63 87.5 77.5 40 9 12 36 77.5 -7.5 32.5
[2,] 69 24 63 87.5 77.5 40 9 12 36 77.5 -7.5 32.5
[3,] 69 24 63 87.5 77.5 40 9 12 36 77.5 -7.5 32.5
Brno, December 2007
14
Lola Ugarte
# This is for understanding what we are doing
Comb[2,pdT6[2, ]]
#we extract from row 2 of matrix Comb the columns
#specified by the second row of matrix pdT6
> pdT6[2,]
[1] 1 2 3 4 5 7
# In words, we take values of row 2 of Comb
#corresponding to columns 1, 2, 3, 4, 5 y 7
> Comb[2,pdT6 [2, ] ]
[1] 69.0 24.0 63.0 87.5 77.5 9.0
#Could you imagine the result of the following command?
Comb[2, -pdT6[2,]]
Brno, December 2007
15
Lola Ugarte
RESULT:
Comb[2, - p d T 6 [ 2 , ] ]
> Comb[2, - p d T 6 [ 2 , ] ]
[1] 4 0 . 0 1 2 . 0 3 6 . 0 7 7 . 5 - 7 . 5 3 2 . 5
#Note that it gives the data from row 2 that DO NOT correspond
#to columns 1, 2, 3, 4, 5 y 7
Leťs remember the contents of Comb[2,]
> Comb[2,]
[1] 69.0 24.0 63.0 87.5 77.5 40.0 9.0 12.0 36.0 77.5 -7.5 32.5
Brno, December 2007
16
Lola Ugarte
Leťs create now the vector Theta that contains all of the possible values of 9 = z - y.
Theta<-array(0,924)
for (i in 1:924) {
Theta[i]<-mean(Comb[i,pdT6[i,]])-mean(Comb[i, -pdT6[i,]])}
Theta.obs<-mean(Treat)-mean(Cont)
pvaK-sum (Theta >=Theta. obs) /choose (12, 6) # exact permutation p-value
> pval
[1] 0.03138528
Brno, December 2007
17
Lola Ugarte
Solution (a) exact permutation p-values
using the package coin
library(coin)
GR<-as.factor(c(rep("Treat",6), rep("Cont",6)))
oneway_test(Blood~GR, distribution^'exact",alternative="less")
Exact 2-Sample Permutation Test
data: Blood by groups Cont, Treat
Z = -1.871, p-value = 0.03139
alternative hypothesis: true mu is less than 0
18
Brno, December 2007 Lola Ugarte
Example (Solution (b)). Estimated
permutation p-value
To compute the estimated permutation p-value we draw randomly and without replacement a
certain amount of samples. From a total of 924 permutations we draw randomly B = 499.
set.seed(13)
theta.obs<-mean(Ratbp$Treat)-mean(Ratbp$Cont)
boot.blood<-sample(l:924,size=4 99, replace=F)
B<-499
theta<-array(0,B)
for (i in 1:length(boot.blood)){
j<-boot.blood[i]
mediaK-mean (Comb [ j, pdT6 [ j, ] ] )
media2<-mean(Comb[j,-pdT6[j,]])
theta[i]<-medial-media2}
pval.aprox<-(sum(theta >=theta.obs)+1)/(B+1) ##pval.aprox
pval.aprox
0.032
19
Brno, December 2007 Lola Ugarte
Example (Solution (b)). -Using package
boot*
For that we need to load package boot
* Let us construct a function called blood.funQ that computes the difference among the six
first values and the six last values of a vector of length 12.
* To obtain samples without replacement, we need to add the argument sim= "permutation"
in the function boot()
Let us see the R solution:
Brno, December 2007
20
Lola Ugarte
l i b r a r y (boot)
s e t . s e e d ( 1 3 )
B<-499
b l o o d . f u n < - f u n c t i o n ( d a t o s , i ) {
d < - d a t o s [ i ]
M D < - m e a n ( d [ l : 6 ] ) - m e a n ( d [ 7 : 1 2 ] )
MD}
b o o t . b l o o d . b < - b o o t ( B l o o d , b l o o d . f u n , R=B, s i m = " p e r m u t a t i o n " )
p l o t ( b o o t . b l o o d . b ) # T h i s g r a p h s t h e s a m p l i n g
# d i s t r i b u t i o n of t h e s t a t i s t i c
p v a l . b o o t . b < - ( s u m ( b o o t . b l o o d . b $ t >= b o o t . b l o o d . b $ t O ) + 1 ) / ( B + 1 )
> p v a l . b o o t . b
[1] 0 . 0 3
Brno, December 2007
21
Lola Ugarte
Histogram oft
i 1 1 1--h 1
-40 0 20 40 60
Brno, December 2007
22
- 3 - 2 - 1 0 1 2 3
Quantiles of Standard Normal
Lola Ugarte
Example (Solución (c)). Estimated bootstrap
p-value using boot
If we select the 499 samples with replacement, we are doing a BOOTSTRAP TEST
* Load package boot
* Let us construct a function called blood.funQ that computes the difference among the six
first values and the six last values of a vector of length 12.
* To obtain samples with replacement, we need to add the argument sim= "ordinary" in the
function boot()
Let us see the R solution:
Brno, December 2007
23
Lola Ugarte
l i b r a r y (boot)
s e t . s e e d ( 1 3 )
B<-499
b l o o d . f u n < - f u n c t i o n ( d a t o s , i ) {
d < - d a t o s [ i ]
M D < - m e a n ( d [ l : 6 ] ) - m e a n ( d [ 7 : 1 2 ] )
MD}
b o o t . b l o o d . b < - b o o t ( B l o o d , b l o o d . f u n , R=B, s i m = " o r d i n a r y " )
p l o t ( b o o t . b l o o d . b ) #Podemos r e p r e s e n t a r l a d i s t r i b u c i \ ' o n
t e n e l m u e s t r e o d e l e s t a d V i s t i c o
p v a l . b o o t . b < - ( s u m ( b o o t . b l o o d . b $ t >= b o o t . b l o o d . b $ t O ) + 1 ) / ( B + 1 )
> p v a l . b o o t . b
[1] 0 . 0 3 8
Brno, December 2007
24
Lola Ugarte
Example (Solution(c)). Estimated bootstrap
p-value without using boot
set.seed(13)
theta.obs<-mean(Ratbp$Treat)-mean(Ratbp$Cont)
boot.blood<-sample(l:924,size=4 99, replace=T)
B<-499
theta<-array(0,B)
for (i in 1:length(boot.blood)){
j<-boot.blood[i]
mediaK-mean (Comb [ j, pdT6 [ j, ] ] )
media2<-mean(Comb[j,-pdT6[j,]])
theta[i]<-medial-media2}
pval.boot<-(sum(theta >=theta.obs)+1)/(B+1) ##pval.boot
#plot(boot.blood)
pval.boot
0.04
25
Brno, December 2007 Lola Ugarte
Example Conclusions
* Using a permutation test the exact p-value is 0.0314.
* The estimated permutation p-value is 0.03.
* The estimated bootstrap p-value is 0.038.
* The three procedures lead to the same conclusion. We reject the null hypothesis,
showing an effect of the drug. It seems to decrease the blood pressure
in rats.
Brno, December 2007
26
Lola Ugarte
Exercise for the students
A company wants to know if it is efficient to teach some new tools to its workers
using internet courses. The company randomly select 7 workers and randomly
assign them to two groups of sizes 4 and 3. The first group attended traditional
courses and the second internet courses. After the courses a test was administered
to the workers whose results were:
Internet courses: 37, 49, 55, 57
Traditional courses: 23, 31, 46
Show if internet courses are preferred over traditional courses using a permutation
test and a bootstrap test. Compute the p-value if the company is interested
in checking if both learning methods are different.
Brno, December 2007
27
Lola Ugarte
Exercise for the students
A Japanese company and an American company claim that they have developed
new technology to increase network transmission speeds. The marketing managers
of both companies simultaneously announce that they can transmit 1 terabyte
per second. To substantiate their claims, each company submits trial data
(in seconds) to transmit one terabyte with their new technologies.
Japanese company
American company
0.98 0.85 0.9 1
0.95 0.94 0.8 1 0.99
Is there evidence to suggest the transmission speed using the technology developed
by the American company is superior to the transmission speed using
the technology developed by the Japanese company? Compute the p - value to
answer the question with the following techniques:
1. Enumerate all possible combinations to find the p - value for a permutation
test.
28
Brno, December 2007 Lola Ugarte
2. Use the function onewaytest from the coin package to calculate an appropriate
p - value. Does this p - value match the one in part 1 ?
3. Obtain an estimated permutation p - value using the boot function from the
boot package.
4. What conclusion do the p - value support?
Brno, December 2007
29
Lola Ugarte