R中的自举相关答案

【问题标题】：Bootstrapped correlation in RR中的自举相关
【发布时间】：2019-10-15 11:26:39
【问题描述】：

我正在尝试在 R 中进行自举相关。我有两个变量 Var1 和 Var2，我想获得 Pearson 相关性的引导 p.value。

my variables look like this:
      x            y
1   .6080522    1.707642
2   1.4307273   1.772616
3   0.8226198   1.768537
4   1.7714221   1.265276
5   1.5986213   1.855719
6   1.0000000   1.606106
7   1.1678940   1.671457
8   0.6630012   1.608428
9   1.0842423   1.670619
10  0.5592512   1.107783
11  1.6442616   1.492832
12  0.8326965   1.643923
13  1.1696954   1.763181
14  0.7484543   1.762921
15  1.0842423   1.591566
16  0.9014748   1.718669
17  0.7604917   1.782863
18  0.8566499   1.796216
19  1.4307273   1.913675
20  1.7579695   1.903155

到目前为止，我有这个：

data = as.data.frame(data)
x = data$Var1
y = data$Var2
dat = data.frame(x,y)

library(boot)
set.seed(1)
bootCorTest3 <- function(data, i){
  d <- data[i, ]
  results  <- cor.test(d$x, d$y, method='pearson')
  c(est = results$estimate, stat = results$statistic, param = results$parameter, p.value = results$p.value, CI = results$conf.int)
}


b3 <- boot(dat, bootCorTest3, R = 1000)
b3

# Original (non-bootstrap) statistics with label
b3$t0
colMeans(b3$t)
boot.ci(b3, type = c("norm", "basic", "perc", "bca")) #bootstrapped CI.

引导的 p 值应该是我使用 colMeans(b3$t) 得到的值，对吧？

colMeans(b3$t) 给了我这个：

est.cor      stat.t    param.df     p.value         CI1         CI2
 0.28495324  2.13981008 48.00000000  0.14418623  0.01438146  0.51726022

似乎一切正常。问题是我在不同的软件上运行了相同的统计数据，结果大相径庭。我在这里得到的 p 值比另一个要高得多。我认为我在这里可能做错了什么，因为我在 R 方面并不强。

谁能给我一些关于这段代码的反馈？难道我做错了什么？你会得到 Pearson 相关性的引导 p.value 吗？

感谢您的宝贵时间。

【问题讨论】：

自举复制在 t 矩阵的行中。
是的，好的。所以如果我这样做 (mean(b3$t) (or colMeans(b3$t)) 我应该得到引导的 p.value，对吧？
我会想到一些使用 rowMeans 或仅仅意味着的东西。
rowMeans 应该是逐行的平均值。在 b3$t 我得到 1 列 1000 行。每行一个p.value（对吗？）。所以，如果我做 rowMeans(b3$t) 我会得到 1000 个结果。我认为 colMeans 或只是意味着有效，因为我只得到一个值，如果我在 excel 中导出并执行平均值，我会得到相同的值。问题是这个值和其他软件还是有很大区别的……
“大不相同”。或许你可以更具体一点？

标签： r correlation statistics-bootstrap

【解决方案1】：

如果你想引导你的相关性检验，你只需要从你的引导统计函数中返回相关系数。在这种情况下，自举相关检验的 p 值是不合适的，因为您忽略了相关检验的方向性。

在 CrossValidated 上查看这个问题，以获得有关执行引导假设检验的一些不错的答案：https://stats.stackexchange.com/questions/20701/computing-p-value-using-bootstrap-with-r

library("boot")
data <- read.csv("~/Documents/stack/tmp.csv", header = FALSE)
colnames(data) <- c("x", "y")

data <- as.data.frame(data)
x <- data$Var1
y <- data$Var2
dat <- data.frame(x,y)

set.seed(1)

b3 <- boot(data, 
  statistic = function(data, i) {
    cor(data[i, "x"], data[i, "y"], method='pearson')
  },
  R = 1000
)
b3
#> 
#> ORDINARY NONPARAMETRIC BOOTSTRAP
#> 
#> 
#> Call:
#> boot(data = data, statistic = function(data, i) {
#>     cor(data[i, "x"], data[i, "y"], method = "pearson")
#> }, R = 1000)
#> 
#> 
#> Bootstrap Statistics :
#>      original        bias    std. error
#> t1* 0.1279691 -0.0004316781    0.314056
boot.ci(b3, type = c("norm", "basic", "perc", "bca")) #bootstrapped CI. 
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 1000 bootstrap replicates
#> 
#> CALL : 
#> boot.ci(boot.out = b3, type = c("norm", "basic", "perc", "bca"))
#> 
#> Intervals : 
#> Level      Normal              Basic         
#> 95%   (-0.4871,  0.7439 )   (-0.4216,  0.7784 )  
#> 
#> Level     Percentile            BCa          
#> 95%   (-0.5225,  0.6775 )   (-0.5559,  0.6484 )  
#> Calculations and Intervals on Original Scale

plot(density(b3$t))
abline(v = 0, lty = "dashed", col = "grey60")

在这种没有 p 值的情况下，可以说大部分抽样分布的质量非常接近于零。

【讨论】：