Определение функции, которая вычисляет ковариационную матрицу корреляционной матрицы

У меня проблемы с преобразованием матрицы и имен строк и столбцов.

Моя проблема в следующем:

В качестве входной матрицы у меня есть (симметричная) корреляционная матрица, подобная этой:

введите описание изображения здесь

вектор корреляции задается значениями нижней треугольной матрицы:

введите описание изображения здесь

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

введите описание изображения здесь

Отклонения можно приблизительно рассчитать как

введите описание изображения здесь

-> N - размер выборки (в этом примере N = 66)

Ковариации можно приблизительно определить как

введите описание изображения здесь

Например, ковариация между r_02 и r_13 определяется как

введите описание изображения здесь

Теперь я хочу определить функцию в R, которая получает на вход корреляционную матрицу и возвращает ковариационную матрицу. Однако у меня есть проблемы с выполнением вычисления ковариаций. Моя идея состоит в том, чтобы дать имена элементам correlation_vector, как показано выше (r_01, r_02 ...). Затем я хочу создать пустую матрицу дисперсии-кокариации, длина которой равна correlation_vector. Строки и столбцы должны иметь те же имена, что и correlation_vector, поэтому я могу назвать их, например, [01] [03]. Затем я хочу реализовать цикл for, который устанавливает значения i и j, а также k и l, как показано в формуле ковариации для столбцов и строк корреляций, которые мне нужны в качестве входных данных для формулы ковариации. Это всегда должно быть шесть разных значений (ij; ik; il; jk; jl; lk). Это моя идея, но сейчас я не знаю, как ее реализовать в R.

Это мой код (без расчета ковариаций):

require(corpcor)

correlation_matrix_input <- matrix(data=c(1.00,0.561,0.393,0.561,0.561,1.00,0.286,0.549,0.393,0.286,1.00,0.286,0.561,0.549,0.286,1.00),ncol=4,byrow=T)

N <- 66 # Sample Size

vector_of_correlations <- sm2vec(correlation_matrix_input, diag=F) # lower triangular matrix of correlation_matrix_input

variance_covariance_matrix <- matrix(nrow = length(vector_of_correlations), ncol = length(vector_of_correlations)) # creates the empty variance-covariance matrix


# function to fill the matrix by calculating the variance and the covariances

variances_covariances <- function(vector_of_correlations_input, sample_size) {

    for (i in (seq(along = vector_of_correlations_input))) {
        for (j in (seq(along = vector_of_correlations_input))) {

            # calculate the variances for the diagonale
            if (i == j) {
                variance_covariance_matrix[i,j] = ((1-vector_of_correlations_input[i]**2)**2)/sample_size 
            }

            # calculate the covariances
            if (i != j) {

                variance_covariance_matrix[i,j] = ???

            }
        }
    }

return(variance_covariance_matrix); 
}

Есть ли у кого-нибудь идеи, как реализовать вычисление ковариаций по приведенной выше формуле?

Буду признателен за любую помощь по данной проблеме !!!


r correlation covariance normal-distribution
person jeffrey    schedule 31.08.2013    source источник
comment
Что такое Korrelationsmatrix_Studie_i? Не вижу, где это определяется.   -  person Mark Miller    schedule 31.08.2013
comment
Это была ошибка. Это correlation_matrix_input.   -  person jeffrey    schedule 31.08.2013

Ответы (3)


arrow_upward
4
arrow_downward

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

covr <- function(r, i, j, k, l, n){
    if(i==k && j==l)
        return((1-r[i,j]^2)^2/n)
    ( 0.5 * r[i,j]*r[k,l]*(r[i,k]^2 + r[i,l]^2 + r[j,k]^2 + r[j,l]^2) +
      r[i,k]*r[j,l] + r[i,l]*r[j,k] - (r[i,j]*r[i,k]*r[i,l] +
      r[j,i]*r[j,k]*r[j,l] + r[k,i]*r[k,j]*r[k,l] + r[l,i]*r[l,j]*r[l,k]) )/n
}

Теперь определите вторую функцию:

vcovr <- function(r, n){
    p <- combn(nrow(r), 2)
    q <- seq(ncol(p))
    outer(q, q, Vectorize(function(x,y) covr(r, p[1,x], p[2,x], p[1,y], p[2,y], n)))
}

И вуаля:

> vcovr(correlation_matrix_input, 66)
            [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
[1,] 0.007115262 0.001550264 0.002917481 0.003047666 0.003101602 0.001705781
[2,] 0.001550264 0.010832674 0.001550264 0.006109565 0.001127916 0.006109565
[3,] 0.002917481 0.001550264 0.007115262 0.001705781 0.003101602 0.003047666
[4,] 0.003047666 0.006109565 0.001705781 0.012774221 0.002036422 0.006625868
[5,] 0.003101602 0.001127916 0.003101602 0.002036422 0.007394554 0.002036422
[6,] 0.001705781 0.006109565 0.003047666 0.006625868 0.002036422 0.012774221

РЕДАКТИРОВАТЬ:

Для преобразованных значений Z, как в вашем комментарии, вы можете использовать это:

covrZ <- function(r, i, j, k, l, n){
    if(i==k && j==l)
        return(1/(n-3))
    covr(r, i, j, k, l, n) / ((1-r[i,j]^2)*(1-r[k,l]^2))
}

И просто замените его на vcovr:

vcovrZ <- function(r, n){
    p <- combn(nrow(r), 2)
    q <- seq(ncol(p))
    outer(q, q, Vectorize(function(x,y) covrZ(r, p[1,x], p[2,x], p[1,y], p[2,y], n)))
}

Новый результат:

> vcovrZ(correlation_matrix_input,66)
            [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
[1,] 0.015873016 0.002675460 0.006212598 0.004843517 0.006478743 0.002710920
[2,] 0.002675460 0.015873016 0.002675460 0.007869213 0.001909452 0.007869213
[3,] 0.006212598 0.002675460 0.015873016 0.002710920 0.006478743 0.004843517
[4,] 0.004843517 0.007869213 0.002710920 0.015873016 0.003174685 0.007858948
[5,] 0.006478743 0.001909452 0.006478743 0.003174685 0.015873016 0.003174685
[6,] 0.002710920 0.007869213 0.004843517 0.007858948 0.003174685 0.015873016
person Ferdinand.kraft    schedule 31.08.2013
comment
У меня есть еще один короткий вопрос: в моем исследовании мне нужны z-преобразованные значения дисперсий и ковариаций. в этом случае формула дисперсии предназначена для всех элементов на диагонали: Var (z_ij) = 1 / (n-3), а знаменатель формулы ковариации - не n, а n [(1-r_ij ^ 2) * ( 1-р_кл ^ 2)]. Я попытался ввести это в ваш код, но мой результат не может быть правильным, потому что значения по диагонали имеют разные результаты. Не могли бы вы сказать мне, как изменить ваш код для вычисления z-преобразованной матрицы var.-cov? - person jeffrey; 31.08.2013
comment
@ Ferdinand.kraft, значения по диагонали должны быть (1 / (66-3)) = 1/63 = 0,01587302. 1/66 будет 0,01515152, как в матрице выше. Но я не могу понять, почему это так. Можешь мне помочь? - person jeffrey; 01.09.2013
comment
@jeffrey, спасибо, это опечатка, условие диагонали if(i==k && j==l). Исправлены обе функции covr и covrZ. Сейчас он работает. - person Ferdinand.kraft; 02.09.2013

arrow_upward
2
arrow_downward

Я написал подход с использованием combn и индексов строки / столбца для генерации различных комбинаций p.

variances_covariances <- function(m, n) {
  r <- m[lower.tri(m)]
  var <- (1-r^2)^2

  ## generate row/column indices
  rowIdx <- rep(1:nrow(m), times=colSums(lower.tri(m)))
  colIdx <- rep(1:ncol(m), times=rowSums(lower.tri(m)))

  ## generate combinations
  cov <- combn(length(r), 2, FUN=function(i) {
    ## current row/column indices
    cr <- rowIdx[i] ## i,k
    cc <- colIdx[i] ## j,l

    ## define 6 cases
    p.ij <- m[cr[1], cc[1]]
    p.ik <- m[cr[1], cr[2]]
    p.il <- m[cr[1], cc[2]]
    p.jk <- m[cc[1], cr[2]]
    p.jl <- m[cc[1], cc[2]]
    p.kl <- m[cr[2], cc[2]]

    ## calculate covariance
    co <- 0.5 * p.ij * p.kl * (p.ik^2 + p.il^2 + p.jk^2 + p.jl^2) +
          p.ik * p.jl + p.il * p.jk -
          (p.ij * p.ik * p.il + p.ij * p.jk * p.jl + p.ik * p.jk * p.kl + p.il * p.jl * p.kl)
    return(co)
  })

  ## create output matrix
  com <- matrix(NA, ncol=length(r), nrow=length(r))
  com[lower.tri(com)] <- cov
  com[upper.tri(com)] <- t(com)[upper.tri(com)]
  diag(com) <- var

  return(com/n)
}

Вывод:

m <- matrix(data=c(1.000, 0.561, 0.393, 0.561,
                   0.561, 1.000, 0.286, 0.549,
                   0.393, 0.286, 1.000, 0.286,
                   0.561, 0.549, 0.286, 1.00), ncol=4, byrow=T)

variances_covariances(m, 66)
#            [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
#[1,] 0.007115262 0.001550264 0.001550264 0.003101602 0.003101602 0.001705781
#[2,] 0.001550264 0.010832674 0.010832674 0.001127916 0.001127916 0.006109565
#[3,] 0.001550264 0.010832674 0.007115262 0.001127916 0.001127916 0.006109565
#[4,] 0.003101602 0.001127916 0.001127916 0.012774221 0.007394554 0.002036422
#[5,] 0.003101602 0.001127916 0.001127916 0.007394554 0.007394554 0.002036422
#[6,] 0.001705781 0.006109565 0.006109565 0.002036422 0.002036422 0.012774221

Надеюсь, я все сделал правильно.

person sgibb    schedule 31.08.2013

arrow_upward
-2
arrow_downward

салам / привет

variance_covariance_matrix<- diag (variance vector, length (r),length (r))
pcomb <- combn(length(r), 2)
for (k in 1:length(r)){
    i<- pcomb[1,k]
    j<- pcomb[2,k]
    variance_covariance_matrix[i,j]<- variance_covariance_matrix [j,i]<- genCorr[k] * sqrt (sig2g[i])  * sqrt (sig2g[j])

}
person safa    schedule 16.11.2013