# Equality of Means (Student, Bootstrap Permutation) Using R

I am getting acquainted with resampling methods. Here is a problem I tried to tackle: I have two normal independent samples of size $$60$$ of equal variance but possibly different means.

The two samples are i.i.d. with distribution $$\mathcal{N}(\mu_{X},\sigma{{}^2})$$ and $$\mathcal{N}(\mu_{Y},\sigma{{}^2})$$ respectively. I just wrote code to test equality of means using Student's t-test, a permutation test and a bootstrap-based test.

While the first two methods yield similar p-values (around 0.2) for the null hypothesis, the bootstrap-based one gets the value wrong by a mile (around 0.5). Is there a problem with my code?

    set.seed(2018)
n <- 60 # sample size
mean <- 20 # expected value of the underlying distribution of the samples
sigma <- 5 # standard deviation of the underlying distribution of the samples
B <- 2000
### Generating the sample in the form of a matrix with 2 rows and n (n=60) columns
### same underlying mean and variance for the samples
Sample <- rbind(rnorm(n=n, mean=mean, sd=sigma), rnorm(n=n, mean=mean, sd=sigma))
indexSet <- 1:n   # index set for the items of the sample

### Hypothesis testing
### Null hypothesis: the ratio of the means is equal to 1
### Alternative hypothesis: the ratio of the means is not equal to 1

x_mean <- mean(Sample[1,])
y_mean <- mean(Sample[2,])
sx2 <- sd(Sample[1,])**2
sy2 <- sd(Sample[2,])**2
sp2 <- ((n-1)*sx2+(n-1)*sy2)/(2*n-2)
t_obs <- (y_mean-x_mean)/(sqrt(sp2)*sqrt(2/n))
### The parametric approach: p-value = P(t_(n+m-2)>t_obs)
pv <- 1-pt(t_obs, df = 2*n-2)

############## Permutations ####################

library(gtools) # for permutations
tVector <- rep(0, B) # Vector with stat
for (j in 1:B){
cSample <- gtools::permute(Sample)
cSample1 <- cSample[1:n]
cSample2 <- cSample[(n+1):(2*n)]
xmean <- mean(cSample1)
ymean <- mean(cSample2)
sx2 <- sd(cSample1)**2
sy2 <- sd(cSample2)**2
sp2 <- ((n-1)*sx2+(n-1)*sy2)/(2*n-2)
tVector[j] <- (ymean-xmean)/(sqrt(sp2)*sqrt(2/n))
}
pVEstPerm<- (1 + length(which(tVector>t_obs)))/(B+1)
pVEstPerm # Approximately 0.2 (very close to the true value)

#############   Bootstrap ##############

tBVector <- rep(0, B) # Vector with stat
for (j in 1:B){
indexSet_star <- sample(indexSet, replace = TRUE)
cSample1 <- Sample[1,c(indexSet_star)]
cSample2 <- Sample[2,c(indexSet_star)]
xmean <- mean(cSample1)
ymean <- mean(cSample2)
sx2 <- sd(cSample1)**2
sy2 <- sd(cSample2)**2
sp2 <- ((n-1)*sx2+(n-1)*sy2)/(2*n-2)
tBVector[j] <- (ymean-xmean)/(sqrt(sp2)*sqrt(2/n))
}

pVEstBoot<- (1 + length(which(tBVector>t_obs)))/(B+1)
pVEstBoot # Approximately 0.5 (for B = 2000)far from the true value)


EDIT: I resample by taking columns of the (2*n) data matrix with replacement.

• I suspect there is at least one problem with the method you use to resample the data, but I'm unsure because you haven't described your bootstrapping method. Could you please include a description in your post? (As a test of your code, see what happens when you apply it to compare samples of different sizes.) – whuber Dec 27 '18 at 22:58
• 1. First explain what your code is supposed to be implementing, so that errors of comprehension (not correctly understanding the precise algorithm to implement) are separated from errors of implementation (simply failing to correctly code the algorithm - usually off topic here). $\quad$ 2. While something may well be wrong with your code, I am concerned by the statement that because one p-value differed from two others, that it's necessarily wrong. You might be suspicious of an error, but p-values of different tests being different is not usually surprising, unless the tests are equivalent. – Glen_b Dec 28 '18 at 0:28
• Please register &/or merge your accounts (you can find information on how to do this in the My Account section of our help center), then you will be able to edit & comment on your own question. – gung - Reinstate Monica Dec 28 '18 at 2:54

For the bootstrap, it is not an equivalent test. The answer being questioned by this test is essentially, "is the true $$\delta$$ between the two populations larger than the one we observed?"

You are resampling without replacement for each of the two populations and counting the times that the bootstrapped $$\hat{\delta}_b$$ is greater than $$\hat{\delta}_{obs}$$. It should not be surprising that this is approx. $$0.5$$.

If you you want the bootstrapped equivalent to the permutation test, you could do basically the same thing as the permutation test, only sample the combined data with replacement.

• (+1) What re-sampling would you do if you suspected $H_a$ might be true and were trying to use a bootstrap to estimate power against a particular alternative. – BruceET Dec 28 '18 at 3:35
• Please avoid posting new questions in the comments. I'm not sure what you gain by using the bootstrap for a power calculation. You could recenter the empirical distribution and sample from it, but that seems specious. – HStamper Dec 28 '18 at 17:27

(a) This seems to be a non-standard formulation of the t test. Ordinarily, it's $$H_0: \mu_1 - \mu_2 = 0$$ vs $$H_a: \mu_1 - \mu_2 \ne 0.$$

### Null hypothesis: the ratio of the means is equal to 1
### Alternative hypothesis: the ratio of the means is not equal to 1


(b) Your annotations fizzle out when you start the bootstrap, making it difficult (for me, anyhow) to make sense of your code.

(c) If you are going to ask people to check a simulation, it is a good idea to set a seed before you generate the fake data and another seed before you do the random re-sampling.

(d) I'm not sure, but it seems your bootstrapping code may be treating the data as paired.

(e) In order to get reproducible results, $$B = 2000$$ seems to me to be a very small number of re-samples.

Some explorations of my own.

# Fake data
set.seed(1227)
x1 = rnorm(60, 20, 5);  x2 = rnorm(60, 20, 5)
t.test(x1, x2, var.eq=T)

t.test(x1, x2, var.eq=T)

Two Sample t-test

data:  x1 and x2
t = -0.48065, df = 118, p-value = 0.6317
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.261090  1.377852
sample estimates:
mean of x mean of y
20.11781  20.55943


Here is my attempt at a bootstrap, which seems to give a P-value between 0.6 and 0.7. At least it's not off by an order of magnitude from the P-value of the actual pooled t test.

I am a fan of permutation tests, so my experience with bootstrapping has been mainly to make confidence intervals. Consequently, my method of doing a test via bootstrapping might not be optimal; criticism welcome. (For one thing, it is convenient, but very wasteful of running time, to use R's t.test procedure to get the re-sampled t statistics, via $-notation.) Just trying to explore some of the issues in my list above. # Bootstrap set.seed(1218); B = 10^5; t.re = numeric(B) for (i in 1:B) { x1.re = sample(x1, 60, rep=T); x2.re = sample(x2, 60, rep=T) t.re[i] = t.test(x1.re, x2.re, var.eq=T)$stat
}
mean(abs(t.re) > abs(-0.48065))
[1] 0.67019