|
|
|
A Manual for use of PyPedal
A software package for pedigree analysis |
|
|
|
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,
, may be written as
, where
is a lower triangular matrix and
is a diagonal matrix [HendersonHenderson1976,ThompsonThompson1977].
traces the flow of genes from one generation to another, while
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
and
. 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,
, 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
without accounting for inbreeding taking advantage of the fact that
, 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
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
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
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
is the same when calculating
regardless of the form
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
presented above.
Release 2.0.3, documentation updated on November 29, 2005
Revised May 15, 2012.
See About this document... for information on suggesting changes.