5.5 Decomposition and Direct Inverses of Numerator Relationship Matrices

PyPedal provides routines in the pyp_nrm module for decomposing and forming direct inverses of numerator relationship matrices (NRM). The NRM, $ A$ , may be written as $ A=TDT'$ , where $ T$ is a lower triangular matrix and $ D$ is a diagonal matrix [HendersonHenderson1976,ThompsonThompson1977]. $ T$ traces the flow of genes from one generation to another, while $ D$ provides varianced and covariances of Mendelian sampling. Mrode ref224 presents a more detailed discussion of this decomposition. The examples in this section all use the pedigree presented in Table 2.1 of that work, which is contained in the file "examples/new_decompose.ped":
# Pedigree from Table 2.1 of
# Mrode (1996)
1 0 0
2 0 0
3 1 2
4 1 0
5 4 3
6 5 2
The results presented below may be obtained by running the "new_decompose.py" example program. The pyp_nrm.a_decompose() function takes as its argument a pedigree object and returns the matrices $ T$ and $ D$ . The code
print 'Calling a_decompose()'
print '====================='
D, T  = pyp_nrm.a_decompose(example)
print 'D: ', D
print
print 'T: ', T
produces identical answers to those presented by Mrode ref224 (note that outputs may have been reformatted slightly to improve readability).
Calling a_decompose()
=====================
D:  [[ 1.       0.       0.       0.       0.       0.     ]
 [ 0.       1.       0.       0.       0.       0.     ]
 [ 0.       0.       0.5      0.       0.       0.     ]
 [ 0.       0.       0.       0.75     0.       0.     ]
 [ 0.       0.       0.       0.       0.5      0.     ]
 [ 0.       0.       0.       0.       0.       0.46875]]

T:  [[ 1.     0.     0.     0.     0.     0.   ]
 [ 0.     1.     0.     0.     0.     0.   ]
 [ 0.5    0.5    1.     0.     0.     0.   ]
 [ 0.5    0.     0.     1.     0.     0.   ]
 [ 0.5    0.25   0.5    0.5    1.     0.   ]
 [ 0.25   0.625  0.25   0.25   0.5    1.   ]]
It is not generally feasible to directly compute the inverse, $ A^{-1}$ , of large NRM. However, it is possible to formulate simple rules for computing the inverse directly from. Henderson ref143 presented a simple procedure for calculating $ A^{-1}$ without accounting for inbreeding taking advantage of the fact that $ A^{-1}=(T^{-1})'D^{-1}T^{-1}$ , which is implemented as pyp_nrm.form_d_nof():
print 'Calling a_inverse_dnf()'
print '======================='
Ainv = pyp_nrm.a_inverse_dnf(example)
print 'Ainv: ', Ainv
produces $ A^{-1}$ without accounting for inbreeding in the pedigree:
Calling a_inverse_dnf()
=======================
Ainv:  [[ 1.83333333  0.5        -1.         -0.66666667  0.          0.        ]
 [ 0.5         2.         -1.          0.          0.5        -1.        ]
 [-1.         -1.          2.5         0.5        -1.          0.        ]
 [-0.66666667  0.          0.5         1.83333333 -1.          0.        ]
 [ 0.          0.5        -1.         -1.          2.5        -1.        ]
 [ 0.         -1.          0.          0.         -1.          2.        ]]
Quaas ref235 built on Henderson's earlier work by providing an algorithm for the computation of $ A^{-1}$ accounting for inbreeding, which is implemented as the function pyp_nrm.a_inverse_df().
print 'Calling a_inverse_df()'
print '======================'
Ainv = pyp_nrm.a_inverse_df(example)
print 'Ainv: ', Ainv
You will note that the numbers in this matric are slightly different from those in the previous example. Those differences are due to the inbreeding in the pedigree, which was not accounted for by a_inverse_dnf().
Calling a_inverse_df()
======================
Ainv:  [[ 1.83333333  0.5        -1.         -0.66666667  0.          0.        ]
 [ 0.5         2.03333333 -1.          0.          0.53333333 -1.06666667]
 [-1.         -1.          2.5         0.5        -1.          0.        ]
 [-0.66666667  0.          0.5         1.83333333 -1.          0.        ]
 [ 0.          0.53333333 -1.         -1.          2.53333333 -1.06666667]
 [ 0.         -1.06666667  0.          0.         -1.06666667  2.13333333]]
The correctness of the results from can be verified by using some linear algebra to compute $ A^{-1}$ directly:
    print 'Calculating Ainv from D and T'
    print '============================='
    l = example.metadata.num_records
    Tinv = numpy.linalg.inv(T)
    print 'Tinv: ', Tinv
    Tpinv = numpy.linalg.inv(T.T)
    print 'Tpinv: ', Tpinv
    Dinv = numpy.linalg.inv(D)
    print 'Dinv: ', Dinv
    Ainvhalf = numpy.dot(Tpinv,Dinv)
    Ainv = numpy.dot(Ainvhalf,Tinv)
    print 'Ainv: ', Ainv
The output includes several intermediate matrices that have been verified against Mrode's book [MrodeMrode1996].
Calculating Ainv from D and T
=============================
Tinv:  [[ 1.   0.   0.   0.   0.   0. ]
 [ 0.   1.   0.   0.   0.   0. ]
 [-0.5 -0.5  1.   0.   0.   0. ]
 [-0.5  0.   0.   1.   0.   0. ]
 [ 0.   0.  -0.5 -0.5  1.   0. ]
 [ 0.  -0.5  0.   0.  -0.5  1. ]]

Tpinv:  [[ 1.   0.  -0.5 -0.5  0.   0. ]
 [ 0.   1.  -0.5  0.   0.  -0.5]
 [ 0.   0.   1.   0.  -0.5  0. ]
 [ 0.   0.   0.   1.  -0.5  0. ]
 [ 0.   0.   0.   0.   1.  -0.5]
 [ 0.   0.   0.   0.   0.   1. ]]

Dinv:  [[ 1.          0.          0.          0.          0.          0.        ]
 [ 0.          1.          0.          0.          0.          0.        ]
 [ 0.          0.          2.          0.          0.          0.        ]
 [ 0.          0.          0.          1.33333333  0.          0.        ]
 [ 0.          0.          0.          0.          2.          0.        ]
 [ 0.          0.          0.          0.          0.          2.13333333]]

Ainv:  [[ 1.83333333  0.5        -1.         -0.66666667  0.          0.        ]
 [ 0.5         2.03333333 -1.          0.          0.53333333 -1.06666667]
 [-1.         -1.          2.5         0.5        -1.          0.        ]
 [-0.66666667  0.          0.5         1.83333333 -1.          0.        ]
 [ 0.          0.53333333 -1.         -1.          2.53333333 -1.06666667]
 [ 0.         -1.06666667  0.          0.         -1.06666667  2.13333333]]
The function pyp_nrm.form_d_nof() is provided as a convenience function. The form of $ D$ is the same when calculating $ A^{-1}$ regardless of the form $ T$ takes (i.e., accounts for inbreeding or not). The code
print 'Calling form_d_nof()'
print '===================='
D  = pyp_nrm.form_d_nof(example)
print 'D: ', D
produces the result
Calling form_d_nof()
====================
D:  [[ 1.    0.    0.    0.    0.    0.  ]
 [ 0.    1.    0.    0.    0.    0.  ]
 [ 0.    0.    0.5   0.    0.    0.  ]
 [ 0.    0.    0.    0.75  0.    0.  ]
 [ 0.    0.    0.    0.    0.5   0.  ]
 [ 0.    0.    0.    0.    0.    0.5 ]]
which may be verified against the values of $ D$ presented above.
See About this document... for information on suggesting changes.