r - 相関行列の共分散行列を計算する関数の定義

Question

行列の変換と行と列の名前に問題があります。

私の問題は次のとおりです。

入力行列として、次のような(対称)相関行列があります。

相関ベクトルは、下三角行列の値によって与えられます。

ここで、これらの相関の分散共分散行列を計算したいと思います。これは、分散共分散行列でほぼ正規分布しています。

分散は次のように概算できます。

-> N はサンプルサイズです (この例では N = 66)

共分散は次のように近似できます。

たとえば、r_02 と r_13 の間の共分散は次のように与えられます。

ここで、相関行列を入力として取得し、分散共分散行列を返す R の関数を定義したいと思います。ただし、共分散の計算を実装するには問題があります。上記のように、correlation_vector の要素に名前を付けることを考えています (r_01、r_02...)。それから、correlation_vector の長さを持つ空の分散共カリアンス行列を作成したいと思います。行と列はcorrelation_vectorと同じ名前にする必要があるため、たとえば[01] [03]で呼び出すことができます。次に、共分散式の入力として必要な相関の列と行に対する共分散の式に示すように、i と j の値と k と l を設定する for ループを実装したいと思います。これらは常に 6 つの異なる値でなければなりません (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
}

次に、この 2 番目の関数を定義します。

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

サラム/こんにちは

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

}