Introduction to Resampling Methods Using R

Contents

1 Sampling from known distributions and simulation

1.1 Sampling from normal distributions
1.2 Specifying seeds
1.3 Sampling from exponential distributions
2 Bootstrapping
2.1 Bootstrap distributions
2.2 Bootstrap confidence intervals
2.2.1 Percentile method
2.2.2 Pivot method
2.2.3 Standard bootstrap
3 Randomization tests
3.1 Creating random permutations
3.2 Comparing groups
3.2.1 Exact randomization distribution
3.2.2 Random sampling the randomization distribution
3.2.3 Choice of test statistic
3.3 Wilcoxon rank-sum test
3.4 Selecting among two-sample tests
3.5 More than two groups
3.6 Contingency tables
4 Methods for correlation and regression
4.1 Randomization test for linear relation
4.1.1 Pearson correlation or slope of regression line
4.1.2 Rank correlation
4.2 Bootstrap intervals for correlation and slope
4.2.1 Bivariate bootstrap sampling
4.2.2 Confidence intervals
4.2.3 Fixed-X sampling for the slope
5 Two-sample bootstrap intervals

1
1. Sampling from known distributions and simulation

In introductory statistics courses we are told that the t-test is “robust” to departures from
normality, especially if the sample size is large. What this means is if we specify a particular
Type I error rate, then the actual proportion of false rejections will be close to the Type I error
rate. Let’s create and run a simulation to explore this.

Steps.
1. Generate a random sample from some population distribution
2. Calculate sample mean, standard deviation and t test statistic
3. Decide if the null hypothesis is rejected
4. Repeat 1-3, counting the number of rejections

1.1. Sampling from normal distributions

counter <‐ 0 # set counter to 0

[Link] <‐ qt(0.95,14) #5% critical value

for (i in 1:1000)
{
x <‐ rnorm(15, 25, 4) # draw a random sample of size 15
from a N(25,4) distribution
t <‐ (mean(x)‐25)*sqrt(15)/sd(x)
if (t >= [Link]) # check to see if result is significant
counter <‐ counter + 1 # increase counter by 1
}

counter/1000 #compute estimate of Type I error rate

## [1] 0.06

If we execute this code again, a different set of random samples will be selected, and a different
estimate obtained

counter <‐ 0 # set counter to 0

[Link] <‐ qt(0.95,14) #5% critical value

for (i in 1:1000)
{
x <‐ rnorm(15, 25, 4) # draw a random sample of size 15
from a N(25,4) distribution
t <‐ (mean(x)‐25)*sqrt(15)/sd(x)
if (t >= [Link]) # check to see if result is significant

2
 counter <‐ counter + 1 # increase counter by 1
}

counter/1000 #compute estimate of Type I error rate

## [1] 0.043

1.2 Specify a seed to get identical results each time.

[Link](4123)
counter <‐ 0 # set counter to 0

[Link] <‐ qt(0.95,14) #5% critical value

for (i in 1:1000)
{
x <‐ rnorm(15, 25, 4) # draw a random sample of size 15
from a N(25,4) distribution
t <‐ (mean(x)‐25)*sqrt(15)/sd(x)
if (t >= [Link]) # check to see if result is significant
counter <‐ counter + 1 # increase counter by 1
}

counter/1000 #compute estimate of Type I error rate

## [1] 0.043

Execute the same code again:

[Link](4123)
counter <‐ 0 # set counter to 0
[Link] <‐ qt(0.95,14) #5% critical value

for (i in 1:1000)
{
x <‐ rnorm(15, 25, 4) # draw a random sample of size 15
from a N(25,4) distribution
t <‐ (mean(x)‐25)*sqrt(15)/sd(x)
if (t >= [Link]) # check to see if result is significant
counter <‐ counter + 1 # increase counter by 1
}

counter/1000 #compute estimate of Type I error rate

## [1] 0.043

3
Instead of using a counter, we may want to store the results so they can be explored later. In the
code below, a vector is created and used to store the calculated t-statistics.

[Link](4123)
nsims <‐ 1000
[Link] <‐ qt(0.95,14) #5% critical value
results <‐ numeric(nsims) #Vector to store t statistics

for (i in 1:nsims)
{
x <‐ rnorm(15, mean=0, sd=1) # draw a random sample of size 15

from a N(25,4) distribution

results[i] <‐ (mean(x)‐0)*sqrt(15)/sd(x)
}

sum(results >= [Link])/nsims #compute estimate of error rate

## [1] 0.043

Having the results saved in a vector allows us to explore the actual sampling distribution. Below
we graphically assess agreement between theoretical and actual distributions.

hist(results, freq = F, ylim=c(0,0.4)) # Plot histogram of t

statistics
curve(dt(x,14), add = TRUE) # superimpose t(14) density

4
1.3 Sampling from an exponential distributions

[Link](4123)
nsims <‐ 1000
[Link] <‐ qt(0.95,14) #5% critical value
results <‐ numeric(nsims) #Vector to store t statistics

for (i in 1:nsims)
{
x <‐ rexp(15, rate=1/25) # draw a random sample of size 15 from an
Exp(mean=25) distribution
results[i] <‐ (mean(x)‐25)*sqrt(15)/sd(x)
}

sum(results >= [Link])/nsims #compute estimate of error rate

## [1] 0.015

Graphically assess agreement between theoretical and actual distributions.

hist(results, freq = F, xlim=c(‐6,4), ylim=c(0,0.4)) # Plot histogram of

t statistics
curve(dt(x,14), add = TRUE) # superimpose t(14) density

5
Available distributions ([Link]

Distribution Functions
Beta pbeta qbeta dbeta rbeta
Binomial pbinom qbinom dbinom rbinom
Cauchy pcauchy qcauchy dcauchy rcauchy
Chi-Square pchisq qchisq dchisq rchisq
Exponential pexp qexp dexp rexp
F pf qf df rf
Gamma pgamma qgamma dgamma rgamma
Geometric pgeom qgeom dgeom rgeom
Hypergeometric phyper qhyper dhyper rhyper
Logistic plogis qlogis dlogis rlogis
Log Normal plnorm qlnorm dlnorm rlnorm
Negative Binomial pnbinom qnbinom dnbinom rnbinom
Normal pnorm qnorm dnorm rnorm
Poisson ppois qpois dpois rpois
Student t pt qt dt rt
Studentized Range ptukey qtukey dtukey rtukey
Uniform punif qunif dunif runif
6
Weibull pweibull qweibull dweibull rweibull
Wilcoxon Rank Sum Statistic pwilcox qwilcox dwilcox rwilcox
Wilcoxon Signed Rank Statistic psignrank qsignrank dsignrank rsignrank

7
2. Bootstrap Confidence intervals

Suppose we want to estimate a population parameter, based on random sample.

- Classical World: Observe one sample and the value of the sample statistic. Sampling
distribution is determined by considering all possible (unobserved) samples from the same
assumed population. Cannot directly observe the sampling distribution.

- Bootstrap World: Rather than assume a population, consider the observed sample to be the
best estimate of the population. In fact, we will assume that it represents the probability
distribution for the population. We can then generate all (or at least very many) possible
samples by taking bootstrap samples (with replacement), from this estimated population and
thus observe the sampling distribution of the sample estimator.

2.1 Drawing bootstrap samples using R.

We start with a very small data set, a set of new employee test scores:

23, 31, 37, 46, 49, 55, 57

First select a sample of size 7, with replacement and compute the mean of the bootstrap sample.

score <‐ c(37,49,55,57,23,31,46)

mean <‐ mean(score)
mean

## [1] 42.57143

boot <‐ sample(score, size=7, replace=TRUE)

boot

## [1] 31 37 31 31 31 37 31

[Link] <‐ mean(boot)

[Link]

## [1] 32.71429

8
We need to do this many times to estimate the sampling distribution of the mean.

score <‐ c(37,49,55,57,23,31,46)

mean <‐ mean(score)
mean

## [1] 42.57143

N <‐ length(score)

nboots <‐ 10000

[Link] <‐ numeric(nboots)
for(i in 1:nboots)
{
[Link] <‐ sample(score, N, replace=TRUE)
[Link][i] <‐ mean([Link])
}

hist([Link])

9
2.2 Bootstrap confidence intervals

Example. Suppose we have a random sample of size 30 from an exponential distribution with
mean 25. We want to use the sample mean to estimate the population mean.

We will discuss three ways to construct confidence intervals using bootstrapping.

2.2.1 Percentile method

If the estimator of the population parameter is the statistic used to create the
distribution (e.g., the sample mean), then the confidence interval is simply the equal-tail
quantiles that correspond to the confidence level.

Bootstrap 95% percentile confidence interval

[Link](4123)
[Link] <‐ rexp(30, rate=1/25)
[Link]

## [1] 2.528853 2.235845 40.423011 5.255557 3.355874 2.724010

## [7] 10.787030 9.792154 31.882324 15.816492 13.925713 38.726646
## [13] 19.283214 70.833128 30.556819 31.620638 64.401698 34.713028
## [19] 6.153141 81.498008 16.034828 80.867315 41.204157 77.066567
## [25] 14.237094 47.647705 35.505526 13.203980 104.795832 11.897545

boxplot([Link])

10
n <‐ length([Link])
[Link] <‐ mean([Link])

nboots <‐ 10000

[Link] <‐ numeric(nboots)
for(i in 1:nboots)
{
[Link] <‐ sample([Link], n, replace=TRUE)
[Link][i] <‐ mean([Link])
}

hist([Link])

11
[Link]

## [1] 31.96579

quantile([Link], c(0.025,0.975))

## 2.5% 97.5%
## 22.49355 42.19327

2.2.2 Pivot method

A pivot quantity is a function of the estimator whose distribution does not depend on the
parameter being estimated.

Example: Estimating the population mean,  based on the sample mean, Y . Then the
Y 
statistic T  ~ t (n  1) has a Student’s t distribution with n-1 degrees of freedom.
S/ n
Because the distribution of T does not depend on  , T is a pivot quantity.

12
 When such a quantity exists, we can then use bootstrapping to estimate the distribution of
the pivot quantity—essentially a custom table—and use quantiles from the table to create
the confidence interval.

Y 
For the statistic T  , the confidence interval is given by
S/ n
 S   S 
Y  tb ,0.975      Y  tb ,0.025  ,
 n  n

where Y , S are the values from the original sample.

Bootstrap 95% T‐pivot confidence intervals

[Link](4123)
[Link] <‐ rexp(30, rate=1/25)

n <‐ length([Link])
[Link] <‐ mean([Link])
[Link] <‐ sd([Link])

nboots <‐ 10000

[Link] <‐ numeric(nboots)
for(i in 1:nboots)
{
[Link] <‐ sample([Link], n, replace=TRUE)
[Link][i] <‐ (mean([Link])‐[Link])*sqrt(n)/sd([Link])
}

hist([Link])

13
[Link] <‐ quantile([Link], 0.975)
[Link] <‐ quantile([Link], 0.025)

[Link] <‐ [Link] ‐ [Link]*[Link]/sqrt(n)

[Link] <‐ [Link] ‐ [Link]*[Link]/sqrt(n)
[Link]

## 97.5%
## 22.80743

[Link]

## 2.5%
## 44.39896

Quantiles of t(29) distribution compared to bootstrap distribution

quantile([Link], c(0.005,0.025,0.05,0.95,0.975,0.995))

## 0.5% 2.5% 5% 95% 97.5% 99.5%

## ‐3.515286 ‐2.416519 ‐1.977130 1.496785 1.780026 2.367328

qt(c(0.005,0.025,0.05,0.95,0.975,0.995),n‐1)

## ‐2.756386 ‐2.045230 ‐1.699127 1.699127 2.045230 2.756386

14
“Theory” behind T-pivot interval
Suppose we wish to estimate the mean of a population. The pivotal method can be used,
assuming we can find a statistic whose distribution does not depend on the parameters to be
estimated.
X 
Ex: t  . If X  N (  , ) , then t  t (n  1) . The distribution of t does not depend on
2
s
n
either  or  2 , and thus t is a pivotal quantity. Then since
 
 X  
P  t   t1   1   , we have
 2 s 2

 n 
 s   s 
t    X    t1  
2 n 2 n
 s   s 
 X  t        X  t1  
2 n 2 n
 s   s 
X  t      X  t1  ,
2 n 2 n
 s   s 
 X  t1      X  t  
2 n 2 n
is a 100(1   )% confidence interval for  .

Now, since t  t (n  1) where t (n  1) is symmetric, we have t  t1 . We can write

2 2

X  t1 (s )    X  t1 (s ),
2 n 2 n
where if   .05,1   2  .975 .

If the sample does not come from a normal population, t is still a pivotal quantity so we can still
write
 
 X  
P  t   t1   1  
 2 s 2

 n 
15
  s   s 
 X  t1      X  t   (1)
2 n 2 n

But we can not say that t  t (n  1) !

One Solution: Estimate the distribution of t using bootstrap sampling.

Xb  X
*For each bootstrap sample of size n from the data, compute tb  , then find tb , and
 sb  2
 
 n
tb,1 and substitute into (1) above.
2

16
2.2.3 Using bootstrap samples to estimate standard error.

The “original” or often called “standard” bootstrap method. This method is motivated by
the “Wald” interval which assumes that many statistics are approximately
normally distributed for large sample sizes, and creates the interval as

Estimator  Z / 2 * SE .

If a formula for the SE is not available, then SE can be estimated using bootstrapping. For
example, to estimate the standard error of the mean:

  X b ,i  X  , where X is the
1 n 2
1. Compute the mean squared error, MSE 
n i 1
mean of the original sample and X b ,i is the mean of the ith bootstrap sample.

2. Compute SE  MSE

3. The confidence interval is X  Z / 2 * SE

[Link](4123)
[Link] <‐ rexp(30, rate=1/25)

n <‐ length([Link])
[Link] <‐ mean([Link])

nboots <‐ 10000

[Link] <‐ numeric(nboots)
for(i in 1:nboots)
{
[Link] <‐ sample([Link], n, replace=TRUE)
[Link][i] <‐ (mean([Link])‐[Link])^2
}

SE <‐ sqrt(mean([Link]))

[Link] <‐ [Link] ‐ qnorm(0.975)*SE

[Link] <‐ [Link] + qnorm(0.975)*SE
[Link]

## [1] 22.07738

[Link]

17
## [1] 41.85421

18
3. Randomization/Permutation tests--Comparing two or more groups

A company is trying to decide whether to augment its traditional instruction for new employees
with computer assisted instruction. Seven new employees are selected. Four are assigned at
random to the new method and the remaining three to the traditional method. A test is given at
the end of instruction for all employees, and the scores are given below.

New method Smith (37), Lin (49), O’Neal (55), Zedillo (57)
Traditional method Johnson (23), Green (31), Zook (46)

Suppose we would like to test for evidence that the new method tends to produce higher scores.
We might test H 0 : μN  μT vs. H a : μN  μT , using a t-test.

 t  2.08 , p-value = 0.046. Statistically significant?

 Independent random samples from normally distributed populations?

Permutation/Randomization test.

 Sampling distribution based upon all possible assignments of the experimental units to
treatments.
 Important assumption: Each assignment is equally likely under the null hypothesis—
guaranteed by random assignment.

3.1 Creating a random permutation using R.

[Link](4123)
score <‐ c(37,49,55,57,23,31,46)
perm <‐ sample(score, replace=F)
perm

## [1] 46 55 37 31 23 57 49

Now, compute the mean of the first 4 entries, the last 3 and compute the difference
[Link] <‐ mean(perm[1:4])
[Link]

## [1] 42.25

[Link] <‐ mean(perm[5:7])

[Link]

19
## [1] 43

diff <‐ [Link]‐[Link]

diff

## [1] ‐0.75

Another way is to create an index vector, and sample for just one group.
[Link](4123)
N <‐ length(score)
index <‐ sample(N, size=4,replace=F)
index

## [1] 7 3 1 6

score[index]

## [1] 46 55 37 31

score[‐index]

## [1] 49 57 23

diff2 <‐ mean(score[index])‐mean(score[‐index])

diff2

## [1] ‐0.75

3.2 Comparing groups

Steps of the permutation test using mean difference as the test statistic:

1. Compute the test statistic: Dobs  X N  X T on the observed data

2. Randomly assign units to treatments, and recompute D

3. Repeat #2 for all possible random assignments of units to treatments

4. The p-value is P  D  Dobs   (# D values at least as large as Dobs )/(# randomizations).

3.2.1. Exact randomization distribution.

The table below lists all 35 possible randomizations and their corresponding mean differences.
Note that the permutation we created is that listed in row 19.

20
Randomization New1 New2 New3 New4 Trad5 Trad6 Trad7 X New  X Trad Sum New
1 46 49 55 57 23 31 37 21.4167 207
*2 37 49 55 57 46 23 31 16.1667 198
3 37 46 55 57 31 49 23 14.4167 195
4 31 49 55 57 46 23 37 12.6667 192
5 31 46 55 57 37 49 23 10.9167 189
6 37 46 49 57 31 23 55 10.9167 189
7 37 46 49 55 31 23 57 9.7500 187
8 23 49 55 57 46 31 37 8.0000 184
9 31 46 49 57 23 37 55 7.4167 183
10 23 46 55 57 37 49 31 6.2500 181
11 31 46 49 55 23 37 57 6.2500 181
12 31 37 55 57 46 49 23 5.6667 180
13 23 46 49 57 31 37 55 2.7500 175
14 31 37 49 57 46 23 55 2.1667 174
15 23 46 49 55 31 37 57 1.5833 173
16 23 37 55 57 46 49 31 1.0000 172
17 31 37 49 55 46 23 57 1.0000 172
18 31 37 46 57 23 49 55 0.4167 171
19 31 37 46 55 23 49 57 -0.7500 169
20 23 31 55 57 46 49 37 -2.5000 166
21 23 37 49 57 46 31 55 -2.5000 166
22 23 37 49 55 46 31 57 -3.6667 164
23 23 37 46 57 31 49 55 -4.2500 163
24 31 37 46 49 23 55 57 -4.2500 163
25 23 37 46 55 31 49 57 -5.4167 161
26 23 31 49 57 46 37 55 -6.0000 160
27 23 31 49 55 46 37 57 -7.1667 158
28 23 31 46 57 37 49 55 -7.7500 157
29 23 31 46 55 37 49 57 -8.9167 155
30 23 37 46 49 31 55 57 -8.9167 155
31 23 31 46 49 37 55 57 -12.4167 149
32 23 31 37 57 46 49 55 -13.0000 148
33 23 31 37 55 46 49 57 -14.1667 146
34 23 31 37 49 46 55 57 -17.6667 140
35 23 31 37 46 49 55 57 -19.4167 137

21
 Two of 35 possible assignments of units to observations were as large or larger than the
observed value of Dobs  16.2 (one of these is Dobs  16.2 ). Thus the p-value is
P  D  Dobs   2 / 35  0.057 .

 Exact p-value—does not depend upon unverifiable assumptions

 The p-value (0.046) from the t-test can be viewed as an approximation to the exact
permutation p-value.

3.2.2 Using R to generate many randomizations and compute p-value.

Now, we use a for loop as before to create many random permutations and corresponding
mean differences, store all the mean differences and then compute the p-value.

score <‐ c(37,49,55,57,23,31,46)

[Link] <‐ mean(score[1:4] ‐ mean(score[5:7]))
N <‐ length(score)
[Link](4132)
nperms <‐ 9999
[Link] <‐ numeric(nperms) # vector to save the random differences
for(i in 1:nperms)
{
index <‐ sample(N, size=4, replace=FALSE) #Choose 4 values w/o replacement
[Link][i] <‐ mean(score[index]) ‐ mean(score[‐index])
}

hist([Link])

22
(sum([Link] >= [Link])+1)/(nperms + 1) #P‐value

## [1] 0.0562

Note that this p-value is not exact.

When sample sizes are moderate to large, enumerating all possible arrangements may be very
time consuming at best and practically impossible at worst. As the table below illustrates, with
two samples of 25, there are over 100 trillion arrangements to consider!

Total sample size Sample 1 Sample 2 Permutations

10 5 5 252
20 10 10 184,756
30 15 15 155,117,520
40 20 20 137,846,528,820
50 25 25 126,410,606,437,752

Solution:

Randomly sample the population of permutations. While enumerating 1 trillion

permutations may be computationally time-prohibitive, in many cases enumerating
10,000 or even 100,000 is not. When estimating an exact p-value based on a random
23
 p(1  p)
sample of R permutations, the exact p-value would be expected to be within 2
R
with 95% confidence. For example, if the true p-value is p  0.05 , the estimate would be
expected to be within the margins of error below:

95% margin
R of error
1000 0.013784
5000 0.006164
10000 0.004359
100000 0.001378

24
Example. This dataset contains results from an experiment in visual perception using random dot
stereograms, such as that shown below. Both images appear to be composed entirely of random
dots. However, they are constructed so that a 3D image (of a diamond) will be seen, if the
images are viewed with a stereo viewer, causing the separate images to fuse. Another way to fuse
the images is to fixate on a point between them and defocus they eyes, but this technique takes
some effort and practice. An experiment was performed to determine whether knowledge of the
form of the embedded image affected the time required for subjects to fuse the images. One
group of subjects (group NV-43 subjects) received either no information or just verbal
information about the shape of the embedded object. A second group (group VV-35 subjects)
received both verbal information and visual information (e.g., a drawing of the object).

[Cleveland, W. S. (1993). Visualizing Data. Original source: Frisby, J. P. and Clatworthy, J.L.,
"Learning to see complex random-dot stereograms," Perception, 4, (1975), pp. 173-178]

 78 
A randomization test on these data involves over    1.8  10 22 permutations. In the following
 43 
code, the data are read and a boxplot created that shows several outliers in each group. Thus we
consider using the median difference as the test statistic instead of the mean difference. The p-
value is estimated based on 10,000 randomly sampled randomizations.
fusion=[Link]("C:/Users/sjricht2/Documents/DataSets/Independent samples T
‐test/Fusion_data.txt", header=TRUE)
boxplot(fusion$time~fusion$treatment)

25
names(fusion)

## [1] "time" "treatment"

N <‐ length(fusion$time)
Treat1 <‐ subset(fusion, Select=time, treatment=="NV", drop=T)
Treat2 <‐ subset(fusion, Select=time, treatment=="VV", drop=T)
N1 <‐ length(Treat1)

head(Treat1,5)

## time treatment
## 1 47.20001 NV
## 2 21.99998 NV
## 3 20.39999 NV
## 4 19.70001 NV
## 5 17.40000 NV

head(Treat2,5)

## time treatment
## 44 19.70001 VV
## 45 16.19998 VV
## 46 15.90000 VV
## 47 15.40002 VV
## 48 9.70000 VV

26
[Link] <‐ median(Treat1$time)‐median(Treat2$time)
[Link]

## [1] 3.3

nperms <‐ 9999

[Link](4132)

result <‐ numeric(nperms)

for(i in 1:nperms)
{
index <‐ sample(N, size=N1, replace = FALSE)
result[i] <‐ median(fusion$time[index]) ‐ median(fusion$time[‐index])
}

(sum(result >= [Link])+1)/(nperms + 1)

## [1] 0.3213

27
What hypotheses are being tested by the permutation test?

No population distribution is assumed, and thus it does not make sense to test parameters
(e.g., equality of means).

Null hypothesis: Population distributions are identical

Alternative hypothesis: Observations tend to be larger in one of the populations

3.2.3 Choice of test statistic

Another advantage of the permutation test is that the function of the sample that is best suited
to address the research question may be used. For a test of location difference, we may use

 Difference of means
 Difference of medians
 others (e.g., trimmed means, ratios)

28
3.3 Wilcoxon Rank sum test: Permutation test on rank transformed data

 Requires only ordinal level data

 Reduces effect of outliers for interval level data
 Hypotheses are the same as permutation test on raw data

Rank transformation: replace data with their respective ranks

 For WRS test, data from both samples are combined, then assigned ranks
 Tied observations generally all receive the “average rank”

Steps for WRS test.

1. Combine observations and assign ranks, with tied observations receiving the average rank
2. Perform permutation test on ranks (mean difference or sum of ranks in sample 1 can be
used as test statistic)

fusion=[Link]("C:/Users/sjricht2/Documents/DataSets/Independent samples T
‐test/Fusion_data.txt", header=TRUE)

fusion$[Link] <‐ rank(fusion$time)

head(fusion,5)

## time treatment [Link]

## 1 47.20001 NV 78.0
## 2 21.99998 NV 77.0
## 3 20.39999 NV 76.0
## 4 19.70001 NV 74.5
## 5 17.40000 NV 73.0

N <‐ length(fusion$time)
Treat1 <‐ subset(fusion, Select=[Link], treatment=="NV", drop=T)
Treat2 <‐ subset(fusion, Select=[Link], treatment=="VV", drop=T)
N1 <‐ length(Treat1)

[Link] <‐ mean(Treat1$[Link])‐mean(Treat2$[Link])

[Link]

## [1] 11.42791

nperms <‐ 9999

[Link](4132)

result <‐ numeric(nperms)

29
for(i in 1:nperms)
{
index <‐ sample(N, size=N1, replace = FALSE)
result[i] <‐ mean(fusion$[Link][index]) ‐
mean(fusion$[Link][‐index])
}

(sum(result >= [Link])+1)/(nperms + 1)

## [1] 0.2064

3.4 Selecting Among Two-sample tests

1) t-test—If selecting independent random sample from normal populations, is optimal for
detecting location difference

2) Permutation test using means—Gives exact p-value regardless of distribution of

population. Power will be similar to t-test

3) Permutation test using medians—Can have higher power than tests on means, especially
for skewed and heavy-tailed distributions

4) Permutation test using ranks (Wilcoxon rank-sum test)-- Can have higher power than
tests on means, especially for skewed and heavy-tailed distributions

The WRS test has been studied extensively in relation to the t-test. The t-test tends to have higher
power for symmetric distributions, especially for lighter tailed distributions and smaller sample
sizes. The WRS test generally has higher power for heavier-tailed distributions and moderate to
large sample sizes.

3.5 More than two groups

The methods of this section can be extended to more than two groups. The ANOVA F-statistic
can be used as the test statistic if using the raw, numeric data. The randomization test based on
ranks is known as the Kruskal-Wallis test.

30
3.6 Randomization tests for contingency tables

Example. Seven patients are included in a study to compare two methods of relieving
postoperative pain. Three are allowed to control the amount of pain-relief medicine themselves,
while the other four are given a physician prescribed level of medicine according to standard
practice. Afterward, the patients evaluate their satisfaction as either “not satisfied (NS)”,
“somewhat satisfied (SS)” or “very satisfied (VS)”.

NS SS VS
Physician prescribed 2 2 0
Self-administered 0 1 2

Is there an association between the method of pain-relief and satisfaction?

Null hypothesis—Method and satisfaction are independent (or, the distributions of patients
among the satisfaction categories are identical)
Alternative hypothesis—There is an association between method and satisfaction (or, the
distributions of patients among the satisfaction categories are different).

Typically a chi-squared test would be considered to assess these hypotheses, using the test
statistic

 observed  exp ected 

2

X  
2
,
all cells exp ected

which has an approximate  2   rows  1 cols  1  distribution for large sample sizes.
However, the approximation can be poor when there are small expected counts.

31
The randomization test can be viewed as identical to the two-group permutation test with a
quantitative response, i.e.,
Physician prescribed NS1, NS2, SS3, SS4
Self-administered SS5, VS6, VS7

Then a randomization test can be performed by randomly assigning observations to groups, as

follows:

 observed  exp ected 

2

1. Compute the test statistic: X obs

2
  on the observed data
all cells exp ected

2. Randomly assign units to treatments, construct the corresponding table and compute
X perm
2

3. Repeat #2 for all possible random assignments of observations to groups

4. The p-value is P  X perm   (# X perm values at least as large as Xobs

2
 X obs
2 2 2
)/(#
randomizations).

Example. The expected frequencies are computed below as (row total)*(column total)/n.
table <‐ matrix(
c(2, 2, 0,
0, 1, 2),
nrow = 2, byrow = TRUE
)
table

## [,1] [,2] [,3]

## [1,] 2 2 0
## [2,] 0 1 2

chitest <‐ [Link]([Link](table))

## Warning in [Link]([Link](table)): Chi‐squared approximation may be

## incorrect

chitest

32
##
## Pearson's Chi‐squared test
##
## data: [Link](table)
## X‐squared = 4.2778, df = 2, p‐value = 0.1178

chitest$expected

## A B C
## A 1.1428571 1.714286 1.1428571
## B 0.8571429 1.285714 0.8571429

For the observed table, the chi-squared statistic is 4.28.

One particular randomization would be the following

Physician prescribed SS5, NS2, SS3, SS4
Self-administered NS1, VS6, VS7

and the resulting table would be

NS SS VS
Physician prescribed 1 3 0
Self-administered 1 0 2

The chi-squared statistic for this table is 4.96.

This process would be somewhat cumbersome to code. Luckily, R has built-in functions that can
do this.

Approximate p-value based on a random sample of permutations

table <‐ matrix(

c(2, 2, 0,
0, 1, 2),
nrow = 2, byrow = TRUE
)

33
[Link]([Link](table),[Link]=T, B=10000)

##
## Pearson's Chi‐squared test with simulated p‐value (based on 10000
## replicates)
##
## data: [Link](table)
## X‐squared = 4.2778, df = NA, p‐value = 0.3247

Exact p-value based on all possible permutations

table <‐ matrix(

c(2, 2, 0,
0, 1, 2),
nrow = 2, byrow = TRUE
)

[Link]([Link](table))

##
## Fisher's Exact Test for Count Data
##
## data: [Link](table)
## p‐value = 0.3143
## alternative hypothesis: [Link]

Approximate p-value based on chi-squared distribution

table <‐ matrix(

c(2, 2, 0,
0, 1, 2),
nrow = 2, byrow = TRUE
)

[Link]([Link](table))

## Warning in [Link]([Link](table)): Chi‐squared approximation may be

## incorrect

##
## Pearson's Chi‐squared test
##
## data: [Link](table)
## X‐squared = 4.2778, df = 2, p‐value = 0.1178

34
35
4. Methods for correlation and regression

4.1. Randomization test for linear relation between quantitative variables

Data are either

(a) a random sample of bivariate pairs (Xi ,Yi ) , or
(b) values of Yi are obtained for fixed values of Xi.

1(a): There is usually no a priori reason to assign either X or Y as independent or dependent.

Thus use correlation as a measure of strength of linear relationship.

E[( X   X )(Y  Y )]
Correlation --   (population correlation)
 XY

(X i  X )(Yi  Y )
Estimate of  -- r 
i 1
n

(X
i 1
i  X ) 2  (Yi  Y ) 2

To test H0 :   0 :

 n2 
The test statistic t    r  t (n  2) can be used if the (X,Y) pairs are a random sample
2 
 1  r 
from a bivariate normal population.

36
 1(b): Fit the model Yi  0  1 Xi  i where  i ’s are iid with mean 0 and variance  2 . 1
is the slope of the regression line.

1 : ˆi 
 ( X  X )(Y  Y )
i i

( X  X )
Estimator of 2
i

To test H0 : 1  0 , the test statistic is t

 (X  X )2 ˆ
i
1  t ( n  2) , if  i ’s are normally
MSE
distributed.

 SY 
It can be shown that ˆ1  r   and also that tcorr  tslope .
 SX 

So, to test for a linear relation between X & Y, we can use either statistic.

4.1.1 Pearson correlation or slope of regression line

A permutation test may be used to obtain an exact p-value, regardless of the form of the
distribution of  i ’s. Under H0 : 1  0 or Ho :   0, X does not affect the value of Y, so an
observed Y is just as likely to occur with any X. Thus, the permutation distribution is derived
from all possible assignments of the observed Ys to the observed Xs.

The procedure is exactly the same as for comparing groups.

Steps of the permutation test using r or ˆ1 as the test statistic:

1. Compute the test statistic, e.g., r, on the observed data

2. Randomly assign observations (Y’s) to treatments (X’s), and recompute the test
statistic, rperm .

37
 3. Repeat #2 for all possible random assignments of observations to treatments

4. The p-value is
P  rperm  robs   (# rperm values at least as large as robs )/(# permutations).

Example. Lea (1965) discussed the relationship between mean annual temperature and the
mortality rate for a type of breast cancer in women. The subjects were residents of certain
regions of Great Britain, Norway, and Sweden.

Randomization test for linear association

[Link](4123)
cancer <‐ [Link]('C:/Users/sjricht2/Documents/DataSets/Regression/Breast
Cancer [Link]', header=T)
cancer

## Mortality Temperature
## 1 102.5 51.3
## 2 104.5 49.9
## 3 100.4 50.0
## 4 95.9 49.2
## 5 87.0 48.5
## 6 95.0 47.8
## 7 88.6 47.3
## 8 89.2 45.1
## 9 78.9 46.3
## 10 84.6 42.1
## 11 81.7 44.2
## 12 72.2 43.5
## 13 65.1 42.3
## 14 68.1 40.2
## 15 67.3 31.8
## 16 52.5 34.0

plot(cancer$Temperature,cancer$Mortality)

38
[Link] <‐ cor(cancer$Temperature,cancer$Mortality)
[Link]

## [1] 0.8748544

[Link] <‐ lm(cancer$Mortality~cancer$Temperature)$coeff[2]

[Link]

## cancer$Temperature
## 2.357695

n <‐ length(cancer$Mortality)
nperms <‐ 9999 #set number of times to repeat this process

result.r <‐ numeric(nperms)

[Link] <‐ numeric(nperms)
for(i in 1:nperms)
{
index <‐ sample(n, size=n, replace = FALSE)
result.r[i] <‐
cor(cancer$Temperature,cancer$Mortality[index])
[Link][i] <‐
lm(cancer$Mortality[index]~cancer$Temperature)$coeff[2]
}
'Permutation test p‐value'

39
## [1] "Permutation test p‐value"

(sum(result.r >= [Link])+1)/(nperms + 1)

## [1] 1e‐04

(sum([Link] >= [Link])+1)/(nperms + 1)

## [1] 1e‐04

4.1.2. Rank correlation

The randomization test can be carried out on rank transformed data. The X and Y values are
ranked separately, and Pearson correlation calculated on the rank-transformed data. Pearson
correlation calculated on rank transformed data is called Spearman correlation. The alternative
hypothesis for this test is that the rank-transformed data have a linear association, or that the
original data have a monotonic (strictly increasing or decreasing) relation.

Randomization test for rank correlation

[Link](4123)
cancer <‐ [Link]('C:/Users/sjricht2/Documents/DataSets/Regression/Breast
Cancer [Link]', header=T)

[Link] <‐ rank(cancer$Temperature)

[Link] <‐ rank(cancer$Mortality)

plot([Link],[Link])

40
[Link] <‐ cor([Link],[Link])
[Link]

## [1] 0.9029412

n <‐ length([Link])
nperms <‐ 9999 #set number of times to repeat this process

result.r <‐ numeric(nperms)

for(i in 1:nperms)
{
index <‐ sample(n, size=n, replace = FALSE)
result.r[i] <‐
cor([Link],[Link][index])
}
'Permutation test p‐value'

## [1] "Permutation test p‐value"

(sum(result.r >= [Link])+1)/(nperms + 1)

## [1] 1e‐04

41
4.2 Bootstrap confidence intervals for correlation and slope

Intervals for 

Suppose we have a random sample of ordered pairs, ( Xi ,Yi ), i  1,2,..., n . If the distribution of
(Xi ,Yi ) is bivariate normal, then it can be shown that:

1 1  r  1 1   1 
Z ln    N  ln  , ,
2 1  r   2  1   n  3

which can be used to construct an interval estimator for  . This interval is not robust to
departures from normality, however. Notice also that the distribution of Z depends on the
population correlation,  , and thus Z is not a pivot quantity.

Bootstrap interval.

1) Draw a specified number of bivariate bootstrap samples of size n, i.e., sample pairs
of observations.

2) Compute rb o o t , the sample correlation for each bootstrap sample.

3) Use the percentile method to construct confidence interval. (No pivot quantity exists).

4.2.1 Bivariate bootstrap sampling

[Link] <‐ sample(1:nrow(cancer), replace = TRUE)

[Link]

## [1] 8 1 16 13 13 16 1 3 15 3 6 14 7 2 14 12

cancer

## Mortality Temperature
## 1 102.5 51.3
## 2 104.5 49.9
## 3 100.4 50.0
## 4 95.9 49.2
## 5 87.0 48.5
## 6 95.0 47.8
## 7 88.6 47.3
42
## 8 89.2 45.1
## 9 78.9 46.3
## 10 84.6 42.1
## 11 81.7 44.2
## 12 72.2 43.5
## 13 65.1 42.3
## 14 68.1 40.2
## 15 67.3 31.8
## 16 52.5 34.0

[Link] <‐ cancer[[Link],]

[Link]

## Mortality Temperature
## 8 89.2 45.1
## 1 102.5 51.3
## 16 52.5 34.0
## 13 65.1 42.3
## 13.1 65.1 42.3
## 16.1 52.5 34.0
## 1.1 102.5 51.3
## 3 100.4 50.0
## 15 67.3 31.8
## 3.1 100.4 50.0
## 6 95.0 47.8
## 14 68.1 40.2
## 7 88.6 47.3
## 2 104.5 49.9
## 14.1 68.1 40.2
## 12 72.2 43.5

4.2.2 Bootstrap confidence interval for correlation, using bivariate sampling

[Link](4123)
cancer <‐ [Link]('C:/Users/sjricht2/Documents/DataSets/Regression/Breast
Cancer [Link]', header=T)

nboot <‐ 10000

[Link] <‐ numeric(nboot)
[Link] <‐ numeric(nboot)

for (i in 1:nboot) {
[Link] <‐ sample(1:nrow(cancer), replace = TRUE)
[Link] <‐ cancer[[Link],]
[Link][i]=cor([Link]$Mortality, [Link]$Temperature)
}
43
cor(cancer$Mortality,cancer$Temperature)

## [1] 0.8748544

quantile([Link],c(0.01,.025,.05,.10,.90,.95,.975,0.99))

## 1% 2.5% 5% 10% 90% 95% 97.5%

## 0.7428583 0.7683708 0.7915637 0.8150605 0.9461222 0.9569523 0.9647395
## 99%
## 0.9722582

The observed correlation is 0.875, and a 95% confidence interval is (0.768, 0.965).

44
Extra.

Bootstrap C.I for slope

Case 1: Data consists of a sample of bivariate pairs  Xi ,Yi  .

1) Draw a specified number of bivariate bootstrap samples of size n, i.e., sample pairs of
observations.

2) Compute ˆ1,b , the slope of the regression line, for each bootstrap sample.

3) Use the percentile method to construct confidence interval.

Since a pivot quantity exists for the sample slope, a t-pivot interval may also be computed.

ˆ1  ˆ MSE
The statistic t  , where SE ( ˆ1 )  , is a pivotal quantity.
SE ( ˆ1 ) ( n  1) S X2

Thus:

1) Draw a specified number of bivariate bootstrap samples of size n, i.e., sample pairs
of observations.
ˆ1,b  ˆ MSE b
2) Compute tb  , where SE ( ˆ1,b )  ;
SE(ˆ1,b ) ( n  1) S X2 ,b

3) For a given confidence level, 1   , determine the quantiles tb , and tb ,1 .

2 2

4) The confidence interval is given by:

ˆ1  tb,1 /2 SE (ˆ1 )  1  ˆ1  tb, /2 SE ( ˆ1 ) .

4.2.3 Fixed X sampling for the slope

 May be used for inferences on slope of regression line.

45
  Assumes Y  h ( X )   , where h ( X ) is some function, say a linear function, are
independent, identically distributed with mean 0, variance  2 .

 Idea: Sample errors and then add to the mean.

Steps:
1) Compute hˆ( X ) from the observed sample.

2) Compute the errors, ei  Yi  hˆ( X i ) .

3) Select a bootstrap sample of n errors ( ei ) and compute Yi  hˆ( X i )  ei ,b , i  1,..., n .

4) Repeat (3) Nb times.

This yields Nb bootstrap samples, each consisting of n pairs,  X i , Yi ,b  .

1) For each bootstrap sample compute

ˆ1,e
te  , where ˆ1, e is the slope of the pairs ( Xi , ei ) for each bootstrap sample, and
SE(ˆ1,e )

 (e  e
i i ,b )2
MSEe n2
SE(ˆ1,e )   ;
(n 1)S X2 (n 1)S X2

2) Then for a given confidence level, 1   , determine the quantiles te , and te ,1 .
2 2

3) The confidence interval is given by ˆ1  te,1 /2 SE (ˆ1 )  1  ˆ1  te, /2 SE (ˆ1 ) .

Fixed X sampling: Assumes

 Correct function being fit.

46
  X-values constants

 Constant error variance

If these assumptions are questionable, use bivariate sampling. Bivariate sampling will
usually result in a larger standard error if assumptions for fixed-X sampling appear
violated.

*See example 8.4.1

47
5. Two-Sample Confidence Intervals

Two cases, as in parametric case: Model: Yij  i  ij

1) ij  F( ) (common error distribution)

Pivot quantity:
(Y1  Y2 )  ( 1   2 )
t
1 1
sp 
n1 n2
SE(Y1  Y2 )

Bootstrap interval:

i) Compute eij  Yij  Yi for each observation.

ii) Select n1 errors, with replacement, from the set of all errors and assign to 1st
sample. Then select n2 in similar fashion and assign to second sample.

iii) Compute t 
e1,b  e2,b
, s p
2

 (e ij  ei )2
1 1 (n1  n2  2)
s p 
n1 n2

iv) Interval is (Y1  Y2 )  te,1 /2 SE(Y1  Y2 )  1  2  (Y1  Y2 )  te, /2 SE(Y1  Y2 ) .

2) Unequal error distributions: Select errors within samples.

e1  e2 (eij  ei ) 2
SE (Y1  Y2 ) 
s12 s 22
 . Compute t  , where s E2  
n1 n 2 s21 s22 ni  1

n1 n2
Interval is: (Y1  Y2 )  te,1 /2 SE(Y1  Y2 )  1  2  (Y1  Y2 )  te, /2 SE(Y1  Y2 ) .

48
49