r - 定义一个计算相关矩阵的协方差矩阵的函数

Question

我对矩阵的转换以及行和列的名称有一些问题。

我的问题如下：

作为输入矩阵，我有一个（对称）相关矩阵，如下所示：

相关向量由下三角矩阵的值给出：

现在，我想计算这些相关性的方差-协方差矩阵，它们与方差-协方差矩阵大致呈正态分布：

方差可以近似为

-> N 是样本大小（在本例中 N = 66）

协方差可以近似为

例如，r_02 和 r_13 之间的协方差由下式给出

现在，我想在 R 中定义一个函数，它将相关矩阵作为输入并返回方差-协方差矩阵。但是，我在执行协方差计算时遇到了问题。我的想法是为相关向量的元素命名，如上所示（r_01，r_02 ...）。然后我想创建一个空的方差-协方差矩阵，它的长度是correlation_vector。行和列应该与correlation_vector 具有相同的名称，因此我可以通过例如[01][03] 来调用它们。然后我想实现一个 for 循环，它设置 i 和 j 以及 k 和 l 的值，如协方差公式所示，我需要作为协方差公式的输入的相关性的列和行。这些必须始终是六个不同的值（ij；ik；il；jk；jl；lk）。这是我的想法，但我不

这是我的代码（不计算协方差）：

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); 
}

有谁知道如何使用上面显示的公式实现协方差的计算？

对于这个问题的任何帮助，我将不胜感激！！！

Answer 1

如果您保留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

Answer 2

我编写了一种方法，使用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

我希望，我做的一切都是正确的。

Answer 3

萨拉姆/你好

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])

}