This section describes how to compute the inbreeding coefficients without having to set up the full L matrix that was described in a
previous section . The motivation for the algorithm is based off the Henderson (1976) paper that
outlined a method to generate A-1 directly from a list of sires and dams and the inbreeding coefficients for each animal. For
a non-inbred population, A-1 can be computed without having to compute either A or L and the entire matrix doesn't need to
be stored in memory. However, for an inbred population, the diagonal elements of L or A must first be found and when L is
too large to store in memory this becomes an issue. In order to make Henderson's algorithm scale to a larger number of animals, modifications to
generate the diagonal elements of L are needed. Shortly after the Henderson (1976) paper, the Quaas (1976) paper was published, which finds
the diagonal elements of L or A without having to store the full L matrix in memory.
As outlined in the section that describes the generation of L the diagonal
elements of L can be computed as functions of the inbreeding coefficients of the known parents. As outlined in Quaas (1976) the inbreeding
coefficients can be transformed into sums of squares of rows corresponding to the kth column of L. As a result, the kth
column or inbreeding coefficient for animal "k" can be generated as a summation of the kth row only. The formulas are outlined below
(p = parent; s = sire; d = dam):
One parent is known: \( \sqrt{1.0 - 0.25\sum_{k=1}^{p} {l_{p,k}}^2} \)
Both parents not known: 1.0
Similar to previous methods, the pedigree has to be sorted so that parents come before progeny. Lastly, if animals are numbered from 1 to the very last animal
then the sire and dam values can be used to index where the respective elements are located within the algorithm. Algorithms to sort and renumber a pedigree is
outlined in the following section . The algorithm outlined below generates the diagonal values of L
(i.e. D values) based on the method outlined above and in Quaas (1976) and then the remaining portion of the algorithm is the exact same as outlined in a
Recursive Method to Create A Inverse section.
Ainv_Quaas1976 <- function(animal,sire,dam)
{
if(is.vector(animal) == FALSE){stop("The animal input variable needs to be a vector!!")}
if(is.vector(sire) == FALSE){stop("The sire input variable needs to be a vector!!")}
if(is.vector(dam) == FALSE){stop("The dam input variable needs to be a vector!!")}
m <- length(animal)
## u: sums of squares each row of L are accumulated #
u = numeric(m)
## v: diagonal elements of L are stored as they are computed #
v = numeric(m)
Ainv <- array(0, dim = c(m, m))
for(k in 1:m)
{
######################################
## Part 1: Calculate Diagonals of L ##
######################################
## Part 1a
if (sire[k] == 0 & dam[k] == 0)
{
v[k] <- 1 # Both parents unknown #
}
if (sire[k] == 0 & dam[k] != 0)
{
v[k] = sqrt(1-0.25*u[dam[k]]) # Only Dam is known #
}
if (sire[k] != 0 & dam[k] == 0)
{
v[k] = sqrt(1-0.25*u[sire[k]]) # Only Sire is known #
}
if (sire[k] != 0 & dam[k] !=0 )
{
v[k] = sqrt(1-0.25*u[sire[k]]-0.25*u[dam[k]]) # Both parents known #
}
## Part 1b: Sums of Squares
if (k < m)
{
for(i in (k+1):m)
{
if(sire[i] >= k & dam[i] >= k){v[i] = 0.5*v[sire[i]]+0.5*v[dam[i]]}
if(sire[i] < k & dam[i] >= k){v[i] = 0.5*v[dam[i]]}
if(dam[i] < k & sire[i] >= k){v[i] = 0.5*v[sire[i]]}
if(sire[i] < k & dam[i] >= k){v[i] = 0.5*v[dam[i]]}
if(sire[i] < k & dam[i] < k){v[i] = 0}
}
}
## Part 1c: Square v and add to u
u[k:m] = u[k:m] + v[k:m]^2
############################
## Part 2: Calculate Ainv ##
############################
bi = v[k]^(-2)
if(sire[k] != 0 & dam[k] != 0) # Both parents known #
{
Ainv[animal[k],animal[k]] = Ainv[animal[k],animal[k]] + bi
Ainv[sire[k],animal[k]] = Ainv[sire[i],animal[k]] - bi/2
Ainv[animal[k], sire[k]] = Ainv[animal[k],sire[k]] - bi/2
Ainv[dam[k],animal[k]] = Ainv[dam[k],animal[k]] - bi/2
Ainv[animal[k],dam[k]] = Ainv[animal[k],dam[k]] - bi/2
Ainv[sire[k],sire[k]] = Ainv[sire[k],sire[k]] + bi/4
Ainv[sire[k],dam[k]] = Ainv[sire[k],dam[k]] + bi/4
Ainv[dam[k],sire[k]] = Ainv[dam[k],sire[k]] + bi/4
Ainv[dam[k],dam[k]] = Ainv[dam[k],dam[k]] + bi/4
}
if(sire[k] != 0 & dam[k] == 0) # Only Sire is known #
{
Ainv[animal[k],animal[k]] = Ainv[animal[k],animal[k]] + bi
Ainv[sire[k],animal[k]] = Ainv[sire[k],animal[k]] - bi/2
Ainv[animal[k],sire[k]] = Ainv[animal[k],sire[k]] - bi/2
Ainv[sire[k],sire[k]] = Ainv[sire[k],sire[k]] + bi/4
}
if(sire[k] == 0 & dam[k] != 0) # Only Dam is known #
{
Ainv[animal[k],animal[k]] = Ainv[animal[k],animal[k]] + bi
Ainv[dam[i],animal[k]] = Ainv[dam[k],animal[k]] - bi/2
Ainv[animal[k],dam[i]] = Ainv[animal[k],dam[k]] - bi/2
Ainv[dam[i],dam[i]] = Ainv[dam[k],dam[k]] + bi/4
}
if(sire[k] == 0 & dam[k] == 0) # Both parents unknown #
{
Ainv[animal[k],animal[k]] = Ainv[animal[k],animal[k]] + bi
}
}
return(Ainv)
}
The following pedigree file from Henderson (1976) can be utilized with the
R code above. The columns are animal, sire and dam and the pedigree is already ordered so that parents come
before progeny. Lastly the animals go from 1 to the total number of animals. Outlined below is what (A-1), u and v look like
at the end of each iteration of the for loop. Compare the final v vector to the diagonals of L
outlined in a previous section.
- Henderson, C. R. 1976. A simple method for computing the inverse of a numerator relationship matrix used in prediction of
breeding values. Biometrics 32(1):69-83.
- Quaas R. L. 1976. Computing the Diagonal Elements and Inverse of a Large Numerator Relationship Matrix. Biometrics, 32(4):949-953.
