It's easier if you keep r
as a matrix and use this helper function to make things clearer:
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
}
Now define this second function:
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)))
}
And voila:
> 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
EDIT:
For the transformed Z values, as in your comment, you can use this:
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))
}
And simply replace it in 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)))
}
New result:
> 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