Коэффициент корреляции от рандомизированных переменных в R

Моя цель - получить новую переменную коэффициентов корреляции (коэффициенты Спирмена), где каждое число соответствует корреляции между двумя рандомизированными переменными.

e.g.

var1=c(1, 2, 3, 0, 2)
var2=c(3, 6, 0, 1, 2)

я пробовал

set.seed(1)
f1=numeric(10000)
for (i in 1:10000) {rand <- replicate(10000, sample(var1))
 rand1 <- replicate(10000, sample(var2))
 f1[i]=cor(rand, rand1, use ="everything", method=c("spearman"))
 }

что дает мне это сообщение: Предупреждение: в f1 [i] = cor (rand, rand1, use = "all", method = c ("spearman")): количество заменяемых элементов не кратно длине замены

Я пробовал это:

cof <- cor((replicate(1000, sample(var1))), (replicate(1000, sample(var2))), use ="everything", method=c("spearman"))

который возвращает матрицу коэффициентов корреляции для каждого значения, а не для каждой переменной

В качестве альтернативы, если есть способ попросить R сопоставить, например строка 1 в одном кадре данных с строкой 1 в другом, затем строки 2, затем строки 3 и т. д. Я могу получить матрицы только моих рандомизированных переменных с помощью этого:

set.seed(1)
f1=numeric(10000)
for (i in 1:10000) {rand <- replicate(10000, sample(var1))
  rand1 <- replicate(10000, sample(var2))
  }

которые мне тогда пришлось бы соотносить друг с другом

есть ли способ вычислить коэффициент корреляции между каждой парой рандомизированных переменных по мере их создания, а затем создать новую переменную, состоящую из коэффициентов корреляции для каждой рандомизации?

Спасибо


r random correlation
person confused    schedule 03.05.2020    source источник

Ответы (3)


arrow_upward
0
arrow_downward

Я не совсем уверен, что понимаю, что вы пытались сделать. Возможно, это решит вашу проблему:

var1=c(1, 2, 3, 0, 2)
var2=c(3, 6, 0, 1, 2)

set.seed(1)
n=100
rand <- replicate(n, sample(var1))
rand1 <- replicate(n, sample(var2))

# That is maybe what you are searching for
f1 <- apply(rand,2,cor,rand1)

У вас будет матрица nxn с каждым (i, j), представляющим корреляцию между i -м столбцом rand и j th столбец rand1.

person Levon Ipdjian    schedule 03.05.2020

arrow_upward
0
arrow_downward

Я думаю, вам должно быть проще использовать фактическую формулу корреляции Спирмена без использования cor ().

Это выглядело бы так:

spearman<-function(x,y){
  X<-as.matrix(x)
  Y<-as.matrix(y)
  y<-rowSums(X)
  a<-rowSums(Y)
  spearman<-2*cor(y,a)/(1+cor(y,a))
  return(spearman)
}

После запуска вы можете использовать

spearman(data1$firstrow,data2$secondrow)

для расчета желаемых корреляций.

И тогда, я думаю, вы могли бы использовать что-то вроде этого цикла:

for (i in nrow(dat)) {
  for (i in nrow(dat)) {
  correlation<-spearman(datmat[i,],datmat2[i,])
  print(correlation[i])
  }
}
person dev120342    schedule 03.05.2020

arrow_upward
0
arrow_downward

Что касается вашего второго вопроса, кажется, ваши матрицы rand и rand1 имеют 5 строк и много столбцов, и вы хотели бы сопоставить каждый столбец из rand с эквивалентным столбцом из rand1? Если я правильно понял, вы можете использовать cor.test для получения корреляции рангов копейщика, например в петле. Поскольку это относительно медленно, вы также можете переписать формулу для ранговой корреляции Спирмена в векторизованной форме и использовать ее (см. Ниже). Если вас интересуют построчные корреляции, матрицы легко настроить или транспонировать.

var1=c(1, 2, 3, 0, 2)
var2=c(3, 6, 0, 1, 2)
set.seed(1)
n=10000
rand <- replicate(n, sample(var1))
rand1 <- replicate(n, sample(var2))

library(matrixStats)
colwiseSpearman <- function(m1, m2, correct=TRUE){
    require(matrixStats)
    n <- dim(m1)[2]
    l <- dim(m1)[1]
    if (correct){
        Txy <- t(sapply(seq_len(n), function(x){
            t0 <- tabulate(rand[,x])
            t1 <- tabulate(rand1[,x])
            return(c(Tx=sum(t0^3-t0)/12, Ty=sum(t1^3-t1)/12))
        }))
        return(((l^3-l)/6 - rowSums((colRanks(rand, ties.method="average")-colRanks(rand1, ties.method="average"))^2) - Txy[,1] - Txy[,2])/sqrt(((l^3-l)/6 - 2*Txy[,1])*((l^3-l)/6 - 2*Txy[,2]))) # Spearman cor.coeff. corrected for ties 
    } else {
        return(1-(6*rowSums((colRanks(rand, ties.method="average")-colRanks(rand1, ties.method="average"))^2) / (l^3-l)))}
}

library(microbenchmark)
microbenchmark(a=colwiseSpearman(rand, rand1),
               b=as.numeric(sapply(seq_len(n), function(x) cor.test(rand[,x], rand1[,x], method="spearman")$estimate)), times=10L )
#> Unit: milliseconds
#>  expr        min         lq       mean    median         uq       max neval cld
#>     a   65.47719   68.06543   74.83393   69.2682   72.90266  109.9133    10  a 
#>     b 2769.97084 2789.39907 2826.01399 2821.6867 2849.08012 2880.5115    10   b
a <- colwiseSpearman(rand, rand1)
b <- as.numeric(sapply(seq_len(n), function(x) cor.test(rand[,x], rand1[,x], method="spearman")$estimate))
all.equal(a, b)
#> [1] TRUE

Создано 03.05.2020 с помощью пакета REPEX (v0.3.0)

person user12728748    schedule 03.05.2020