relfactor: Relationship factor from a pedigree

In Rpedigree/pedigreeTools: Versatile Functions for Working with Pedigrees

Relationship factor from a pedigree

Description

Determine the right Cholesky factor of the relationship matrix for the pedigree ped, possibly restricted to the specific labels that occur in labs.

Usage

relfactor(ped, labs = NULL)

getL(ped, labs = NULL)

Arguments

ped	pedigree
labs	a character vector or a factor giving individual labels to which to restrict the relationship matrix and corresponding factor using Colleau et al. (2002) algorithm. If labs is a factor then the levels of the factor are used as the labels. Default is the complete set of individuals in the pedigree.

Details

Note that the right Cholesky factor is returned, which is upper triangular, that is from A = LL' = R'R (lower (upper triangular) and not L (lower triangular) as the function name might suggest.

Value

matrix (dtCMatrix - upper triangular sparse)

Functions

• getL(): Relationship factor from a pedigree

References

Colleau, J.-J. An indirect approach to the extensive calculation of relationship coefficients. Genet Sel Evol 34, 409 (2002). https://doi.org/10.1186/1297-9686-34-4-409

Examples

ped <- pedigree(sire = c(NA, NA, 1,  1, 4, 5),
                dam =  c(NA, NA, 2, NA, 3, 2),
                label = 1:6)
(L <- getL(ped))
chol(getA(ped))

# Test for correctness
LExp <- matrix(data = c(1.0000, 0.0000, 0.5000, 0.5000, 0.5000, 0.2500,
                        0.0000, 1.0000, 0.5000, 0.0000, 0.2500, 0.6250,
                        0.0000, 0.0000, 0.7071, 0.0000, 0.3536, 0.1768,
                        0.0000, 0.0000, 0.0000, 0.8660, 0.4330, 0.2165,
                        0.0000, 0.0000, 0.0000, 0.0000, 0.7071, 0.3536,
                        0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.6847),
               byrow = TRUE, nrow = 6)
stopifnot(!any(abs(round(L, digits = 4) - LExp) > .Machine$double.eps))
LExp <- chol(getA(ped))
stopifnot(!any(abs(L - LExp) > .Machine$double.eps))

(L <- getL(ped, labs = 4:6))
(LExp <- chol(getA(ped)[4:6, 4:6]))
stopifnot(!any(abs(L - LExp) > .Machine$double.eps))

Rpedigree/pedigreeTools documentation built on Oct. 13, 2023, 9:49 p.m.