I have a relationship matrix generated from GCTA, that I can import into R using the following function
ReadGRMBin=function(prefix, AllN=F, size=4){
sum_i=function(i){
return(sum(1:i))
}
BinFileName=paste(prefix,".grm.bin",sep="")
NFileName=paste(prefix,".grm.N.bin",sep="")
IDFileName=paste(prefix,".grm.id",sep="")
id = read.table(IDFileName)
n=dim(id)[1]
BinFile=file(BinFileName, "rb");
grm=readBin(BinFile, n=n*(n+1)/2, what=numeric(0), size=size)
NFile=file(NFileName, "rb");
if(AllN==T){
N=readBin(NFile, n=n*(n+1)/2, what=numeric(0), size=size)
}
else N=readBin(NFile, n=1, what=numeric(0), size=size)
i=sapply(1:n, sum_i)
return(list(diag=grm[i], off=grm[-i], id=id, N=N))
}
It then lists the off diagonal and diagonal.
$ diag: num [1:850] 0.878 0.815 1.11 1.161 1.062 ...
$ off : num [1:360825] 0.0181 -0.0304 -0.0663 -0.0211 -0.0583 ...
$ n : int 850
I wish to create a grm I can inverse from this and ideally in the output row, column, value
I have tried the following code but it doesn't read the off diagonal in the correct format
m <- matrix(NA, ncol = length(grm$diag), nrow = length(grm$diag))
m[lower.tri(m)] <- grm$off
m[upper.tri(m)] <- t(m)[upper.tri(t(m))]
diag(m) <- grm$diag
m
want=cbind(which(!is.na(m),arr.ind = TRUE),na.omit(as.vector(m)))
Instead of reading the diagonal values as
2 1, 3 1, 3 2, 4 1, 4 2 etc
It is reading the diagonal going length wise as
2 1, 3 1, 4 1, 5 1, 6 1 etc
So the resulting matrix (shortened) ends up like this
[,1] [,2] [,3] [,4] [,5]
[1,] 0.87798703 0.018129893 -0.03044302 -0.066282429 -0.02106927
[2,] 0.01812989 0.814602911 0.07577287 -0.004078172 -0.03182918
[3,] -0.03044302 0.075772874 1.10976517 -0.055698857 -0.03960679
[4,] -0.06628243 -0.004078172 -0.05569886 1.160611629 -0.01021352
[5,] -0.02106927 -0.031829182 -0.03960679 -0.010213521 1.06245303
When preference is this
[,1] [,2] [,3] [,4] [,5]
[1,] 0.87798703 0.018129893 -0.03044302 -0.02106927 -0.04011643
[2,] 0.01812989 0.814602911 -0.06628243 -0.00582625 -0.06237402
[3,] -0.03044302 -0.06628243 1.10976517 0.1315616 -0.1601102
[4,] -0.02106927 -0.00582625 0.1315616 1.160611629 -0.1388046
[5,] -0.04011643 -0.06237402 -0.1601102 -0.1388046 1.06245303
If you know how to amend the above code to give the wanted format it would be much appreciated. The end desired output would be the inverse of the matrix in long format if possible. Thanks
1 1 12456
1 2 78910
1 3 34568
1 4 68942
One simple solution is to adapt your code to fill the upper triangle before the lower (since it is the upper triangle that should be filled in column order):
grm = list(
diag = 1:5 * 11,
off = 0:9)
m <- diag(grm$diag)
m[upper.tri(m)] <- grm$off
m[lower.tri(m)] <- t(m)[lower.tri(t(m))]
# [,1] [,2] [,3] [,4] [,5]
# [1,] 11 0 1 3 6
# [2,] 0 22 2 4 7
# [3,] 1 2 33 5 8
# [4,] 3 4 5 44 9
# [5,] 6 7 8 9 55