时间戳

首页 > 解决方案 > 使用 Bootstrap 生成两个变量的相关性并计算置信区间

r - 使用 Bootstrap 生成两个变量的相关性并计算置信区间

问题描述

我正在用 R 编写循环或函数，但我还没有真正理解如何做到这一点。目前，我需要编写一个循环/函数（不确定哪个更好）来在同一数据帧中创建多个 Bootstrap 结果。

示例数据集如下所示：

"ID A_d B_d C_d D_d E_D f_D chkgp
M1  10  20  60  30  54  33  Treatment
M1  20  50  40  33  31  44  Placebo
M2  40  80  40  23  15  66  Placebo
M2  30  90  40  67  67  66  Treatment
M3  30  10  20  22  89  77  Treatment
M3  40  50  30  44  50  88  Placebo
M4  40  30  40  42  34  99  Treatment
M4  30  40  50  33  60  80  Placebo",header = TRUE, stringsAsFactors = FALSE)

我写了一个函数来找到斯皮尔曼相关性

k=cor(df$A_d,df$E_D,method="spearman")
k

结果是-0.325407

现在我必须运行引导方法，通过对两个变量中的数据进行混洗来获得 5000 次相关值

所以使用下面的代码

fc <- function(d, i){
    d2 <- d[i,]
    return(cor(df$A_d,df$E_D,method="spearman"))
}

定义函数 fc 后，我们可以使用 boot 命令，提供我们的数据集名称、我们的函数以及要绘制的引导样本数。

根据 5000 次引导重复计算 BOOTSTRAP 置信区间计算。

#turn off set.seed() if you want the results to vary
set.seed(626)
bootcorr <- boot(hsb2, fc, R=500)
bootcorr

我从 5000 次重复中找出置信区间

boot.ci(boot.out = bootcorr, type =c( "perc"))

结果

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 500 bootstrap replicates

CALL : 
boot.ci(boot.out = bootcorr, type = c("perc"))

Intervals : 
Level     Percentile     
95%   (-0.3254, -0.3254 )  
Calculations and Intervals on Original Scale

我需要写一个循环条件来得到结果如下

Variable1 Variable2 confidence interval
A_d       E_D        (-0.3254, -0.3254 )  
A_d       f_D
B_d       E_D
B_d       f_D
C_d       E_D
C_d       f_D
D_d       E_D
d_d       f_D

因为我有一个包含 100 多个变量的数据集，所以每次都很难做到，所以我需要自动化部分来做到这一点。

标签： rloopsbootstrap-4bootstrap-modalconfidence-interval

解决方案

我们可以创建一个矢量化函数并使用outer()：

corpij <- function(i,j,df) {cor(df[,i],df[,j],method="spearman")}
corp <- Vectorize(corpij, vectorize.args=list("i","j"))

outer(2:(ncol(df1)-1),2:(ncol(df1)-1),corp,df1)

#>            [,1]         [,2]         [,3]       [,4]        [,5]
#> [1,]  1.0000000  0.289588955 -0.480042672 0.22663483 -0.32540701
#> [2,]  0.2895890  1.000000000 -0.006379918 0.53614458 -0.35928788
#> [3,] -0.4800427 -0.006379918  1.000000000 0.01913975 -0.13952023
#> [4,]  0.2266348  0.536144578  0.019139754 1.00000000  0.02395253
#> [5,] -0.3254070 -0.359287879 -0.139520230 0.02395253  1.00000000
#> [6,]  0.7680403 -0.120481928 -0.421074589 0.33734940  0.07185758
#>             [,6]
#> [1,]  0.76804027
#> [2,] -0.12048193
#> [3,] -0.42107459
#> [4,]  0.33734940
#> [5,]  0.07185758
#> [6,]  1.00000000

另一种方法是使用psych::corr.test()：

library(psych)

corr.test(df1[,-c(1,ncol(df1))], method = "spearman")$r

数据：

df1 <- read.table(text="ID A_d B_d C_d D_d E_D f_D chkgp
                        M1  10  20  60  30  54  33  Treatment
                        M1  20  50  40  33  31  44  Placebo
                        M2  40  80  40  23  15  66  Placebo
                        M2  30  90  40  67  67  66  Treatment
                        M3  30  10  20  22  89  77  Treatment
                        M3  40  50  30  44  50  88  Placebo
                        M4  40  30  40  42  34  99  Treatment
                        M4  30  40  50  33  60  80  Placebo",
header = TRUE,stringsAsFactors = FALSE)

推荐阅读