R 上相关系数的自举 p 值
Bootstrapped p-value for a correlation coefficient on R
在 R
,我使用了 boostrap 方法来获得相关系数估计和置信区间。
我想,为了获得 p 值,我可以计算不包含零的置信区间的比例。但这不是解决方案。
在这种情况下如何获得 p 值?
我正在使用 cor.test
进行系数估计。 cor.test
也可能会给我每个测试的 p 值。但是我怎样才能得到自举的 p 值?
非常感谢!
下面是一个例子:
n=30
data = matrix (data = c (rnorm (n), rnorm (n),rnorm (n), rpois(n,1),
rbinom(n,1,0.6)), nrow = n, byrow = F)
data= as.data.frame(data)
z1 = replicate( Brep, sample(1:dim(data)[1], dim(data)[1], replace = T))
res = do.call ( rbind, apply(z1, 2, function(x){ res=cor.test(data$V1[x], data$V2[x]) ; return ((list(res$p.value,res$estimate))) }))
coeffcorr = mean(unlist(res[,2]), na.rm = T) #bootstrapped coefficient
confInter1 = quantile(unlist(res[,2]), c(0.025, 0.975), na.rm = T)[1] #confidence interval 1
confInter2 = quantile(unlist(res[,2]), c(0.025, 0.975), na.rm = T)[2] #confidence interval 2
p.value = mean (unlist(res[,1]), na.rm = T ) # pvalue
R 中 bootstrapping 的标准方式是使用基础包 boot
。您首先定义 bootstrap 函数,该函数采用两个参数,即数据集和数据集的索引。这是下面的函数 bootCorTest
。在函数中,您对数据集进行子集化,仅选择由索引定义的行。
剩下的就简单了。
library(boot)
bootCorTest <- function(data, i){
d <- data[i, ]
cor.test(d$x, d$y)$p.value
}
# First dataset in help("cor.test")
x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
y <- c( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8)
dat <- data.frame(x, y)
b <- boot(dat, bootCorTest, R = 1000)
b$t0
#[1] 0.10817
mean(b$t)
#[1] 0.134634
boot.ci(b)
有关函数 boot
和 boot.ci
的结果的更多信息,请参阅它们各自的帮助页面。
编辑。
如果您想 return 来自引导统计函数 bootCorTest
的几个值,您应该 return 一个向量。在以下情况下,它 return 是一个具有所需值的命名向量。
请注意,我设置了 RNG 种子,以使结果可重现。我上面应该已经做了。
set.seed(7612) # Make the results reproducible
bootCorTest2 <- function(data, i){
d <- data[i, ]
res <- cor.test(d$x, d$y)
c(stat = res$statistic, p.value = res$p.value)
}
b2 <- boot(dat, bootCorTest, R = 1000)
b2$t0
# stat.t p.value
#1.841083 0.108173
colMeans(b2$t)
#[1] 2.869479 0.133857
在 R
,我使用了 boostrap 方法来获得相关系数估计和置信区间。
我想,为了获得 p 值,我可以计算不包含零的置信区间的比例。但这不是解决方案。
在这种情况下如何获得 p 值?
我正在使用 cor.test
进行系数估计。 cor.test
也可能会给我每个测试的 p 值。但是我怎样才能得到自举的 p 值?
非常感谢!
下面是一个例子:
n=30
data = matrix (data = c (rnorm (n), rnorm (n),rnorm (n), rpois(n,1),
rbinom(n,1,0.6)), nrow = n, byrow = F)
data= as.data.frame(data)
z1 = replicate( Brep, sample(1:dim(data)[1], dim(data)[1], replace = T))
res = do.call ( rbind, apply(z1, 2, function(x){ res=cor.test(data$V1[x], data$V2[x]) ; return ((list(res$p.value,res$estimate))) }))
coeffcorr = mean(unlist(res[,2]), na.rm = T) #bootstrapped coefficient
confInter1 = quantile(unlist(res[,2]), c(0.025, 0.975), na.rm = T)[1] #confidence interval 1
confInter2 = quantile(unlist(res[,2]), c(0.025, 0.975), na.rm = T)[2] #confidence interval 2
p.value = mean (unlist(res[,1]), na.rm = T ) # pvalue
R 中 bootstrapping 的标准方式是使用基础包 boot
。您首先定义 bootstrap 函数,该函数采用两个参数,即数据集和数据集的索引。这是下面的函数 bootCorTest
。在函数中,您对数据集进行子集化,仅选择由索引定义的行。
剩下的就简单了。
library(boot)
bootCorTest <- function(data, i){
d <- data[i, ]
cor.test(d$x, d$y)$p.value
}
# First dataset in help("cor.test")
x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
y <- c( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8)
dat <- data.frame(x, y)
b <- boot(dat, bootCorTest, R = 1000)
b$t0
#[1] 0.10817
mean(b$t)
#[1] 0.134634
boot.ci(b)
有关函数 boot
和 boot.ci
的结果的更多信息,请参阅它们各自的帮助页面。
编辑。
如果您想 return 来自引导统计函数 bootCorTest
的几个值,您应该 return 一个向量。在以下情况下,它 return 是一个具有所需值的命名向量。
请注意,我设置了 RNG 种子,以使结果可重现。我上面应该已经做了。
set.seed(7612) # Make the results reproducible
bootCorTest2 <- function(data, i){
d <- data[i, ]
res <- cor.test(d$x, d$y)
c(stat = res$statistic, p.value = res$p.value)
}
b2 <- boot(dat, bootCorTest, R = 1000)
b2$t0
# stat.t p.value
#1.841083 0.108173
colMeans(b2$t)
#[1] 2.869479 0.133857