# Bootstrap bias of the sample correlation

I am working on a problem asking "Provide a Bootstrap bias of $$\hat{r}$$ using R = 1000 bootstrap estimators." However, I am stuck on what even to do. As far as I know, the formula for finding bootstrap bias is

$$Bias_{boot}(\hat{r}) = E[\hat{r}^\ast] - \hat{r}$$

Please correct me if I am wrong. So from here, I used simple bootstrap function as

boot.fn = function (data, index) {
X = data[index]
Y = data[index]
sample_cor_hat = cor(X, Y)
return(sample_cor_hat)
}

b = boot(sample_distribution, boot.fn, R = 1000)
b


Then I would get

Bootstrap Statistics :
original       bias     std. error
t1*        1 1.110223e-16 7.719711e-17


Is my formula for bias correct? am I using the correct way for the function?

update: here's the given data

sample_size = 20
sample_meanvector = c(2, 3)
sample_covariance_matrix = matrix(c(2, 0.4, 0.4, 1), ncol = 2)
sample_distribution = mvrnorm(n = sample_size, mu = sample_meanvector, Sigma = sample_covariance_matrix)
X = matrix(sample_distribution[,1])
Y = matrix(sample_distribution[,2])

• X and Y are the same in your boot,fn function (aren't they?) Jun 8, 2022 at 4:38
• I am pretty new to stats so I might misunderstand something. I forgot to give more context, but I updated my post. @JeremyMiles Jun 8, 2022 at 4:52
• you can pretty much ignore my codes, sorry about the confusion. I guess my question is about the formula. Jun 8, 2022 at 4:56
• if your question is about the formula then the formula is correct. take the difference between the average of the 1000 bootstrap correlations and the true correlation. Jun 8, 2022 at 7:38

## 1 Answer

It seems you are using a simple quantile style of bootstrap. Perhaps it would be better to use bias-corrected quantile bootstrap as below.

Roughly, each bootstrap re-sample of the sample correlation is compared to the observed sample correlation, and the distances d.re between them are found. Then the boostrap CI is based on quantiles .025 and .975 of the $$B$$ d.res.

For my fictitious normal data x and y, this bootstrap CI and the classical CI for population $$\rho$$ are essentially both $$(-0.14, 0.25.).$$

# Fictitious normal data and classical CI
set.seed(2022)
x = rnorm(100, 50, 6)
y = rnorm(100, 50, 6)
r.obs = cor(x,y)
cor.test(x,y)\$conf.int
[1] -0.1409677  0.2506400
attr(,"conf.level")
[1] 0.95
plot(x, y, pch=19)


# 95% Nonparametric bootstrap CI
set.seed(1234)
B= 2000;  seq = 1:100;  d.re=numeric(B)
for (i in 1:B)  {
seq.re = sample(seq, 100,rep=T)
d.re[i] = cor(x[seq.re],y[seq.re]) - r.obs
}
LR = quantile(d.re, c(.975,.025))
r.obs - LR
97.5%       2.5%
-0.1387333  0.2473008


Then for non normal data, the same bootstrap method provides a reasonable 95% CI for $$\rho.$$

set.seed(608)
x = rnorm(100, 50, 6) + 1:100
y = rnorm(100, 50, 6) + 1:100
r.obs = cor(x,y);  r.obs
[1] 0.9559053
plot(x, y, pch=19)
shapiro.test(x)$$p.val [1] 0.002321988 # Not normal shapiro.test(y)$$p.val
[1] 0.01022498        # Not normal


Because x and y are not normal one cannot rely on the CI from cor.test; see end note. However, my bootstrap procedure should give useful CIs for samples as large as 100.

set.seed(1234)
B= 2000;  seq = 1:100;  d.re=numeric(B)
for (i in 1:B)  {
seq.re = sample(seq, 100,rep=T)
d.re[i] = cor(x[seq.re],y[seq.re]) - r.obs
}
LR = quantile(d.re, c(.975,.025))
r.obs - LR
97.5%      2.5%
0.9412820 0.9751636


Note from R documentation for cor.test: "If method is 'pearson', the test statistic is based on Pearson's product moment correlation coefficient cor(x, y) and follows a t distribution with length(x)-2 degrees of freedom if the samples follow independent normal distributions. [Emphasis added.]