Выполнение 1000 тестов перестановки в R для корреляции

Question

Я хотел бы выполнить 1000 тестов перестановки на наборе данных "law" в R, чтобы проверить значимость корреляции между оценками 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)