Randomization Tutorial

I. Sampling Techniques in the Mosaic Package

We have three main sampling functions in the Mosaic package for R that help produce randomized versions of statistical tests and confidence intervals.

1. Shuffle. A permutation, size equal to number of observations.
2. Sample. Draws without replacement, size less than or equal to number of observations.
3. Resample. Draws with replacement, size unlimited.

Let’s create a sample vector Prime so we can see how these functions work.

Prime = c(2,3,5,7,11,13,17,19,23,29)
Prime

 [1]  2  3  5  7 11 13 17 19 23 29

1. Shuffle

The shuffle function permutes the values in the vector.

shuffle(Prime)

 [1] 23  3  5 13  7  2 11 19 29 17

Draws are without replacement, and the number of draws must be equal to the number of observations.

2. Sample

The sample function draws without replacement so the number of draws must be less than or equal to the number of observations.

sample(Prime, size = 6)

[1] 23 13  3  7  5 29

Left on its default setting for the size of draw, sample is equivalent to shuffle because the draws will equal the number of observations.

sample(Prime)

 [1] 29 23 11  7  2 19 17 13  5  3

3. Resample

The resample function draws with replacement, so the number of draws is unlimited.

resample(Prime, size = 20)

 [1]  7  7  2 17 11 19 23 17  2 29  2 19  2  2 19  7 13 17  7  2

If left to its default setting, the draws will equal the number of observations. Since a single observation may repeat multiple times, resample is not equivalent to either shuffle or sample, except by accident.

resample(Prime)

 [1] 29 29 23  7 11 13 23 11 11 13

II.Toy Data Frames

Along with the Primes vector, we will need a couple of toy data frames to demonstrate the randomization techniques for \(t\)-tests, correlation and proportion tests.

Zoo Data Frame

Zoo = as.data.frame(matrix(c(5:12,1:8),nrow=16))
colnames(Zoo) = "Number"
Animal = c("Lion","Lion","Lion","Lion","Lion","Lion","Lion","Lion","Zebra","Zebra","Zebra","Zebra","Zebra","Zebra","Zebra","Zebra")
Zoo$Animal = Animal
Zoo

We can use the t.test function on the Zoo data frame.

t.test(Number ~ Animal, data = Dinner)

    Welch Two Sample t-test

data:  Number by Animal
t = 3.266, df = 14, p-value = 0.005631
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.373184 6.626816
sample estimates:
 mean in group Lion mean in group Zebra 
                8.5                 4.5 

Dinner Data Frame

Dinner = as.data.frame((matrix(c(3, 4, 5, 8, 9, 11,150, 175, 325, 375, 475, 525),nrow=6)))
colnames(Dinner) = c("Predator", "Prey")
Dinner

We can now create a linear model using the lm function on the Dinner data frame.

lm(Predator ~ Prey, data = Dinner)

Call:
lm(formula = Predator ~ Prey, data = Dinner)

Coefficients:
(Intercept)         Prey  
   -0.02933      0.01984  

III. Bootstrapping

For one-sample \(t\)-tests and \(t\)-interval estimates, we use a randomization method called bootstrapping. We repeatedly resample the observations from the data frame to generate a distribution.

Confidence Interval

To create a randomized version of the \(t\)-interval, we draw 500 resamples and load their means into a data frame called bootstrap. Once we find the percentiles corresponding to the endpoints of the interval, we’re done.

bootstrap = do(500) * mean(resample(Prime))
qdata(~mean, p=c(0.025, 0.975), data = bootstrap)

   2.5%   97.5% 
 7.7475 18.1525 

\(t\)-test

Consider a one-sample \(t\)-test on Prime using the hypothesis \[H_0 : \mu_p = 16\\H_a : \mu_p > 16\] To create a randomized version, we input a logical operation into the sum function so that it will count for us. Again, we use the bootstrap distribution, and count the number of times the mean is greater than 16.

sum(bootstrap$mean > 16)

[1] 60

So our randomized \(p\)-value is \[p=\frac{56}{500}=\frac{14}{125}= .112\] and we fail to reject the null. Note the \(p-value\) from the \(t\)-test would be

t.test(~Prime, mu = 16, alternative = "less")

    One Sample t-test

data:  Prime
t = -1.0863, df = 9, p-value = 0.1528
alternative hypothesis: true mean is less than 16
95 percent confidence interval:
     -Inf 18.13107
sample estimates:
mean of x 
     12.9 

For some reason, the \(t\)-statistic is negative which forces us to reverse the inequality from our hypothesis. Still, the bootstrapping \(p\)-value (0.112) appears to be reasonable given \(p\)-value from the \(t\)-test \((p=0.153)\).

IV. Permutation Tests

For permutation tests, we shuffle the observations of the dependent variable.

Test of Significant Correlation

Consider the the linear model \[\text{Predator} \sim \text{Prey}\] which has a significant correlation.

mod = lm(Predator ~ Prey, data = Dinner)
summary(mod)

Call:
lm(formula = Predator ~ Prey, data = Dinner)

Residuals:
       1        2        3        4        5        6 
 0.05333  0.55733 -1.41867  0.58933 -0.39467  0.61333 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.029333   0.955223  -0.031  0.97697   
Prey         0.019840   0.002615   7.587  0.00162 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8952 on 4 degrees of freedom
Multiple R-squared:  0.935, Adjusted R-squared:  0.9188 
F-statistic: 57.56 on 1 and 4 DF,  p-value: 0.001619

To create a permutation test, we choose a random permutation of the dependent variable and pair those values with the independent variable. In general, we need several thousand permutations to generate accurate results.

perm = do(10000) * coef(lm(shuffle(Predator) ~ Prey, data = Dinner))
perm

We are testing the hypothesis that \[H_0 : \beta = 0\\H_a : \beta \neq 0\] By entering a logical expression into the function sum, we can count the number of randomized models that had slope greater than 0.0198.

sum(abs(perm$Prey) > 0.0198)

The estimated significance of the slope coefficient is \[p=\frac{28}{10,000} = 0.0028\] which is reasonable given what the linear regression \(t\)-test results were (\(p = 0.0016\)). We used a two-tailed test by taking the absolute value of the shuffled model slopes. We can randomize a one-tailed test by removing the absolute value and adjusting the inequality. Either way, we have strong evidence that the slope (or correlation) is non-zero.

Estimates of the Slope Coefficient

We can use a bootstrapping approach to estimate the slope coefficient with 95% confidence.

boots = do(500) * coef(lm(Predator ~ Prey, data=resample(Dinner)))
qdata(boots$Prey, c(0.025, 0.975))

      2.5%      97.5% 
0.01518438 0.02830419 

Independent Samples \(t\)-test

We can permute the dependent variable so that its association with the grouping variable is random.

perm = do(1000) * mean(shuffle(Number) ~ Animal, data = Zoo)
perm$diff = perm$Lion - perm$Zebra
perm

The observed mean difference in the between Lions and Zerbra was 4. To estimate the \(p\)-value for the difference in means, we need to count how many of the randomized differences were greater than or equal to 4.

sum(perm$diff >= 4)

[1] 6

Thus, we have an estimated \(p\)-value of \[p=\frac{4}{1000}=0.004\] and, therefore, strong evidence for a difference in group means.