Корреляционная петля Спирмена в R

Question

В предыдущем посте объяснялось, как сделать цикл хи-квадрат в R для всех ваших пар данных: Анализ хи-квадрат с использованием цикла for в R . Я хотел использовать этот код, чтобы сделать то же самое для корреляции Спирмена.

Я уже пытался изменить несколько переменных и смог вычислить корреляционные переменные Пирсона, используя этот код:

library(plyr)

combos <- combn(ncol(fullngodata),2)

adply(combos, 2, function(x) {
  test <- cor.test(fullngodata[, x[1]], fullngodata[, x[2]])

  out <- data.frame("Row" = colnames(fullngodata)[x[1]]
                , "Column" = colnames(fullngodata[x[2]])
                , "cor" = round(test$statistic,3)
                ,  "df"= test$parameter
                ,  "p.value" = round(test$p.value, 3)
                )
 return(out)

})  

Но поскольку я работаю с данными по порядковой шкале, мне нужно использовать корреляцию Спирмена.

Я думал, что смогу получить эти данные, просто добавив команду method = "spearman", но это, похоже, не работает. Если я использую код:

library(plyr)

combos <- combn(ncol(fullngodata),2)

adply(combos, 2, function(x) {
  test <- cor.test(fullngodata[, x[1]], fullngodata[, x[2]], method="spearman")

  out <- data.frame("Row" = colnames(fullngodata)[x[1]]
                , "Column" = colnames(fullngodata[x[2]])
                , "Chi.Square" = round(test$statistic,3)
                ,  "df"= test$parameter
                ,  "p.value" = round(test$p.value, 3)
                )
  return(out)

})  

Я получаю ответ:

Error in data.frame(Row = colnames(fullngodata)[x[1]], Column =    
colnames(fullngodata[x[2]]),  : 
arguments imply differing number of rows: 1, 0
In addition: Warning message:
In cor.test.default(fullngodata[, x[1]], fullngodata[, x[2]], method = "spearman") :
Cannot compute exact p-values with ties

что я делаю не так?

Answer 1

Попробуйте rcor.test функцию в пакете ltm.

mat <- matrix(rnorm(1000), 100, 10, dimnames = list(NULL, LETTERS[1:10]))
rcor.test(mat, method = "spearman")

  A      B      C      D      E      F      G      H      I      J     
A  ***** -0.035  0.072  0.238 -0.097  0.007 -0.010 -0.031  0.039 -0.090
B  0.726  ***** -0.042 -0.166  0.005  0.025  0.007 -0.231  0.005  0.006
C  0.473  0.679  *****  0.046  0.074 -0.020  0.091 -0.183 -0.040 -0.084
D  0.017  0.098  0.647  ***** -0.060 -0.151 -0.175 -0.068  0.039  0.181
E  0.338  0.960  0.466  0.553  *****  0.254  0.055 -0.031  0.072 -0.059
F  0.948  0.805  0.843  0.133  0.011  ***** -0.014 -0.121  0.153  0.048
G  0.923  0.941  0.370  0.081  0.588  0.892  ***** -0.060 -0.050  0.011
H  0.759  0.021  0.069  0.501  0.756  0.230  0.555  ***** -0.053 -0.193
I  0.700  0.963  0.690  0.701  0.476  0.130  0.621  0.597  ***** -0.034
J  0.373  0.955  0.406  0.072  0.561  0.633  0.910  0.055  0.736  *****

upper diagonal part contains correlation coefficient estimates 
lower diagonal part contains corresponding p-values

Answer 2

Проблема в том, что cor.test возвращает значение NULL для параметра, когда вы проводите тест на копье.Из параметра ?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)

str(cor.test(x, y, method = "spearman"))
List of 8
 $ statistic  : Named num 48
  ..- attr(*, "names")= chr "S"
 $ parameter  : NULL    
 $ p.value    : num 0.0968
 $ estimate   : Named num 0.6
  ..- attr(*, "names")= chr "rho"
 $ null.value : Named num 0
  ..- attr(*, "names")= chr "rho"
 $ alternative: chr "two.sided"
 $ method     : chr "Spearman's rank correlation rho"
 $ data.name  : chr "x and y"
 - attr(*, "class")= chr "htest"

Решение: если вы удалите из кода следующую строку, она должна работать:

,  "df"= test$parameter