在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中。你正在做的似乎是一个非参数的引导程序。这里有两种不同的方式，有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)