Search code examples
rmatrixinverse

Construct and inverse matrix from list using R


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

Solution

  • 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