运行 R 中的 1000 个排列测试相关性

Question

我想在 R 中的 "law" 数据集上进行运行 1000 次排列测试，以测试 LSAT 分数和 GPA 之间相关性的显着性。我有以下代码：

nperm <- 1000
law.perm <- rep(0,nperm)
for (i in 1:nperm) {
   ind <- sample(law)
   law <- ind
   Group1 <- law$LSAT[law==1]
   Group2 <- law$GPA[law==2]
   law.perm[i] <- cor(Group1,Group2)
}
law.perm

但是，运行使用上述代码生成相关性的所有 NA 值。谁能帮忙找出问题所在？

下面是一些示例输出：

str(law)
'data.frame':   15 obs. of  2 variables:
$ LSAT: num  576 635 558 578 666 580 555 661 651 605 ...
$ GPA : num  3.39 3.3 2.81 3.03 3.44 3.07 3 3.43 3.36 3.13 ...

Answer 1

数据集 law 在包 bootstrap 中。而你所做的似乎是一个非参数bootstrap。这里有两种不同的方式，使用 for 循环和使用函数 bootstrap::bootstrap.

在运行代码之前，加载数据集。

library(bootstrap)

data(law)

首先，您在问题中尝试的方式已更正。

set.seed(1234)    # Make the results reproducible

nperm <- 1000
law.perm <- numeric(nperm)
n <- nrow(law)
for (i in 1:nperm) {
  ind <- sample(n, replace = TRUE)
  law.perm[i] <- cor(law[ind, "LSAT"], law[ind, "GPA"])
}

第二种方式，使用bootstrap函数。这遵循函数帮助页面中的最后一个示例。

theta <- function(x, xdata){ 
  cor(xdata[x, 1], xdata[x, 2]) 
}

set.seed(1234)
res <- bootstrap(seq_len(n), nperm, theta = theta, law)

比较两个结果。

mean(law.perm)
#[1] 0.769645

mean(res$thetastar)
#[1] 0.7702782

中位数的差异较小。

median(law.perm)
#[1] 0.7938093

median(res$thetastar)
#[1] 0.7911014

并绘制两个结果的图表。

op <- par(mfrow = c(1, 2))
hist(law.perm, prob = TRUE)
hist(res$thetastar, prob = TRUE)
par(op)

运行 R 中的 1000 个排列测试相关性

Running 1000 permutation tests in R for correlation

r

resampling