1
$\begingroup$

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])
$\endgroup$
4
  • $\begingroup$ X and Y are the same in your boot,fn function (aren't they?) $\endgroup$ Jun 8, 2022 at 4:38
  • $\begingroup$ I am pretty new to stats so I might misunderstand something. I forgot to give more context, but I updated my post. @JeremyMiles $\endgroup$
    – so0jong
    Jun 8, 2022 at 4:52
  • $\begingroup$ you can pretty much ignore my codes, sorry about the confusion. I guess my question is about the formula. $\endgroup$
    – so0jong
    Jun 8, 2022 at 4:56
  • 1
    $\begingroup$ 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. $\endgroup$
    – seanv507
    Jun 8, 2022 at 7:38

1 Answer 1

0
$\begingroup$

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)

enter image description here

# 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.

enter image description here

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.]

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.