Bootstrapping confidence intervals in randomization model and population model

Question

Ernst 2004 shows that the permutation tests in the random assignment scheme (e.g., controlled experiment) and in the random sampling scheme (e.g., observational study) share the same constructing mechanics but different reasoning (pages 677, 681). The scopes of inference in two cases are inherently different.

I then wonder about the construction of confidence interval of a statistic by bootstrapping. It seems sensible under the random sampling scheme. Does it make sense under the random assignment scheme?

Update: I continue with the code from @BruceET in which I use difference between means as the metric in the permutation test, p-val is similar. My concern is about the 95% bootstrap CI that I calculate below for this difference in the improvement scores between two groups of student, whether it makes sense in the random assignment setting.

set.seed(1234)
x1 = rnorm(15, 100, 14)
x2 = rnorm(15, 110, 14)

d.obs <- mean(x2) - mean(x1); d.obs
[1] 11.14442

x <- c(x1, x2)
d.per <- rep(0, 99999)
for (i in 1:99999){
    newx <- split(sample(x), rep(c("x1", "x2")), c(15,15))
    d.per[i] <- mean(newx$x2) - mean(newx$x1)
}

p <- (sum(abs(d.per)>=abs(d.obs))+1)/1e5; p
[1] 0.02118

d.boot <- rep(0, 99999)
for (i in 1:99999){
    newx1 <- sample(x1, replace = T)
    newx2 <- sample(x2, replace = T)
    d.boot[i] <- mean(newx2) - mean(newx1)
}
quantile(d.boot, c(0.025,0.975))
2.5%     97.5% 
 2.640416 20.407527 

Answer 1

I will show both a bootstrap CI and a permutation test, with comments along the way. In both cases it happens to be possible to get similar results without computationally intensive methods.

Bootstrap. A nonparametric bootstrap CI for a population standard deviation $\sigma.$

Let x be a sample of size $n = 100$ from a population with unknown SD $\sigma.$ Find a 95% bootstrap CI for $\sigma.$

set.seed(2022)
x = rnorm(100, 50, 7)
sd(x)
[1] 7.140652

The bootstrap procedure helps to estimate the variability of the sample SD $S.$

If we knew the distribution of $S=\sigma,$ then we could find bounds $L$ and $U$ such that $P(L < S-\sigma < U) = 0.95$ so that a 95% CI for $\sigma$ would be of the form $S-U < \sigma < S-U.$

Not knowing $L$ and $U,$ we repeatedly take re-samples, which are samples of size $n$ with replacement from x, find the sample standard deviation $S^*$ of each and use distances $S^* - S_{obs}$ to get an idea of the distribution of $S - \sigma,$ and thus estimates of $L$ and $U.$

The resulting 95% CI for $\sigma$ is $(6.13,\, 8.24).$

s.obs = sd(x);  s.obs
[1] 7.140652

set.seed(417)
d = replicate(2000, 
       sd(sample(x,100,rep=T))-s.obs)
LU = quantile(d, c(.975,.025))
s.obs - LU
   97.5%     2.5% 
6.125899 8.241970

In this particular example, we are using simulated normal data, for which $\sigma = 7.$ and for which an exact 95% CI $(6.27,\, 8.30)$ for $\sigma,$ based on a chi-squared distribution, is available. [In an application with real data, none of this would be known.]

sqrt(99*var(x)/qchisq(c(.975,.025),99))
[1] 6.269541 8.295119

Of course, the exact CI is better. Moreover, we could have improved the bootstrap CI in various ways; perhaps most notably by recognizing that $\sigma$ is a scale parameter and bootstrapping ratios instead of differences.

The main point is that even one of the simplest bootstrap methods gives a useful CI if we know nothing about the population distribution, except that it has a standard deviation, as does its estimate $S.$

Permutation test. Consider two independent samples. We wish to test whether their population means are equal.

set.seed(1234)
x1 = rnorm(15, 100, 14)
x2 = rnorm(15, 110, 14)
x = c(x1, x2)
g = rep(1:2, each=15)
boxplot(x ~ g, horizontal=T, col="skyblue2")

Suppose the 30 randomly chosen students are randomized into two groups, and taught how to do a task using two different methods. Data x1 and x2 are improvement scores derived from before and after tests, the same tests used for both groups.

Not having done this kind of training before, nor used these tests before, we wonder whether data are normal. Particularly because of the small sample sizes we wonder whether it is appropriate to use a pooled 2-sample test to judge whether the group means differ.

If the null hypothesis (no difference) is true, then we believe that the values in x1 and x2 are exchangeable. That is, under $H_0$ we believe the two groups are indistinguishable.

We may have doubts about the validity of a pooled 2-sample t test, because it may not be clear that the t statistic has a t distribution. Nevertheless, we are willing to use the pooled t statistic a reasonable way to measure any difference between the two groups. Accordingly, we will use the pooled t statistic as a 'metric' for a permutation test. And we will simulate the distribution of the metric under $H_0.$

The permutation distribution of the t statistic can be approximated by scrambling the improvement scores in x1 and x2 at random at each iteration and finding the t statistic for each scrambled set of samples. We use $m =100\,000$ iterations.

The following R code finds the $m$ pooled t statistics and the approximate P-value of the the permutation test. Note that 'sample(g)` randomly scrambles the group assignments at each iteration.

t.obs = t.test(x~g,var.eq=T)$stat;  t.obs
        t 
-2.374767 
set.seed(418)
t.stat = replicate(10^5, t.test(x~sample(g),var.eq=T)$stat)
mean(abs(t.stat) >= abs(t.obs))
[1] 0.0211  # aprx P-val of perm test

The histogram below shows the simulated permutation distribution of the t statistic. The P-value of the test is the total probability outside the two vertical bars.

hist(t.stat, prob=T, col="skyblue2")
 abline(v=c(t.obs,-t.obs), col="brown")

The P-value $0.025$ of the pooled t test is nearly the same as the P-value of the permutation test. [So, the pooled t test would have been OK.]

t.test(x~g, var.eq=T)$p.val
[1] 0.02464866