Fit a mixed model by AI-REML

Fit a mixed model by AI-REML

Usage

model(
  formula,
  data,
  pedigree = NULL,
  genotypes = NULL,
  blend = 0.05,
  apy_core = NULL,
  missing_code = NULL,
  maxiter = 300L,
  tol = 1e-08,
  n_em = 4L,
  metafounders = NULL,
  gamma = NULL,
  verbose = interactive(),
  weights = NULL,
  start = NULL,
  residual_by = NULL
)

Arguments

formula

for example peso ~ cg + sexo + animal(id). An unmarked term is a fixed class effect, and fixed(cg) is the same term (it takes traits = in model_mt()); cov(x) is a fixed covariate; animal(id), maternal(dam) and sire(sire) are random with relationship, and sire(sire, mgs = "mgs") is the sire and maternal-grandsire model, 1 on the sire and 1/2 on the maternal grandsire of the record in the same effect (an unknown grandsire, "0", leaves the sire only, and an unknown sire the grandsire only; a missing value in that column is refused, see below); pe(id) and random(lote) are random without relationship. group = "nome" puts two random terms in the SAME covariance matrix, with the correlation estimated. The terms of a group index one set of levels, and level l of one term covaries with level l of the other: the pedigree animals in a relationship group, and in a group of terms without a relationship matrix (random(), pe()) the union of the level names of their columns, sorted (integer labels in numeric order, then the rest). The pairing is by name, so the order of the rows does not matter, and a level that appears in only one column is still an effect of the other term, with no record there. Their covariance is estimated only through the levels that have records in both columns: with none (sire and dam ids that never repeat between the sexes) it does not enter the likelihood, and every fitter that estimates components stops before fitting and names the two terms. At given components (start = with maxiter = 0 and n_em = 0, or theta_fixed = in gibbs()) the group is accepted. traits = "y1" on a random term is for model_mt() only: see there. The other fitters fit one trait and stop on it. ONLY THE OBSERVATION MAY BE MISSING. A column that gives the level of a term (a fixed class, the id of a random term, the maternal grandsire of sire(mgs =), the class of a nested covariate) may not hold a missing label on ANY row, including a row whose observation is missing: NA, NaN, Inf or -Inf in a numeric column, and NA, blank text (what read.csv() gives an empty cell of a text column; spaces, tabs and line breaks alone count as blank) or the text "NaN" (what factor() makes of a numeric NaN) in a text or factor column. Those rows would form one shared level. The texts "NA" and "Inf" are accepted, since they may be real names, and in the maternal-grandsire column blank text is the unknown grandsire, like "0". A column that enters as a VALUE, the covariate of cov() (nested or not) and the columns of base =, may not hold NA, NaN, Inf or -Inf on any row either: the covariate would leave the model without a word. Every case stops the fit with the term, the column, the number of such rows and the first of them (a text NA with the column and the first row); drop those rows, or fill them, before the fit. The rule is checked on every row, also where the observation is missing, so one data frame is accepted or refused alike by model(), model_mt(), gibbs(), snp_blup(), the eval_internal family and pegs() (id and environment =). Refusing is the conservative one of the field's conventions: ASReml stops on a missing explanatory value unless !MVINCLUDE (the missing value counts as zero) or !MVREMOVE (the record is dropped) is declared, and asreml-R defaults to na.method(y = "include", x = "fail"); BLUPF90 reads a 0 in a class column as the effect being absent from that record, so a BLUPF90 data file with 0 as its missing class code gives here one shared level "0", which cannot be refused because "0" may be a real name: drop those records, or give them their levels, before the fit. A missing observation (NA, NaN, Inf or -Inf in the trait, or missing_code) with its levels and covariates present only drops its record. missing_code is the code of the OBSERVATION only, as OPTION missing is in BLUPF90, where a data file with -999 for an unknown age in a cov column regresses those records on an age of -999. Here a covariate cell equal to missing_code is refused on every row, like NA, with the term, the count and the first row, when the code lies more than 10 robust SDs (the median absolute deviation, scaled to the SD) from the median of the other cells of the column, the signature of a code. Closer to them it may be a real value (the 0 of a 0/1 covariate under missing_code = 0, the BLUPF90 default), and the fit goes on with a warning that the cells enter as that value; WOMBAT likewise stops on a covariable of zero until told whether it is a value (COVZER FIT) or a missing record (IGN) (Meyer, WOMBAT manual, 2017, section 4.10.6.2). A level column holding the code is one more level "-999". Turn those cells into NA, which stops the fit with the column and the first row, and then drop or fill those records. VALUES THAT LOOK LIKE A CODE. Before the fit, every fitter that reads a numeric trait (model(), model_mt(), gibbs(), snp_blup(), the eval_internal family and pegs()) warns about two patterns among the observed records, naming the trait, the count, the first row and the value: a common missing-value code (-9, -99, -999, -9999, -99999, 99, 999, 9999, 99999, or -123456789, the WOMBAT default) held by two records or more, outside the range of the other values and at least 3 robust SDs from their median; and records more than 10 robust SDs from the median. The same two patterns in a covariate, on the records with the trait observed, give a warning that names the covariate. On simulated data, ten records at -999, or three records 1e4 above the rest, took the heritability from 0.35 to 1e-10 without a word before this check. A NUMBER STAYS A NUMBER: a factor or text column read as a number (the trait, a covariate, the traits of pegs()) is refused with the column, and not read as its codes 1, 2, ... The marker argument base = takes a value and not a column name. The fit evaluates it in the environment where you wrote the formula, so a formula that one function builds and another fits finds its basis. The fit keeps in formula the formula with the value of base = written in. So h2() and accuracy() read the value of the fit after the variable changes or goes away. The arguments that NAME something, group =, nome =, nested = and mgs =, are taken literally and never evaluated: group = g is the group called "g", the same as group = "g", and not the value of a variable g. Because a group meant by value would silently become a group of its own, a term that is alone in a group whose unquoted name is a variable of the formula's environment holding a different text is refused, with both spellings in the message; write group names in quotes.

data

data.frame with the columns referenced

pedigree

data.frame animal, sire, dam; required with a relationship term. Ids may be numeric or text, and each column may have its own type: a numeric id is written as its full integer everywhere (the data, the pedigree, the genotypes, the matrices this package builds and the names of the results), so 100000 as a double, as an integer and as the text "100000" are one animal. A numeric id that is not an integer keeps the digits that tell it apart from every other number (1234567.5 is "1234567.5", and 123456.7 never meets 123457). Text that R wrote in scientific notation ("1e+05" is what as.character(), factor() and rownames<- give for the double 100000) is a different label, and when the other side holds that same number the fit stops with an error rather than dropping the animal's records.

THE IDS OF THE DATA AGAINST THE PEDIGREE follow one rule in every fitter (model(), model_mt(), gibbs(), snp_blup() and the eval_internal family), and accuracy(), solutions() and snp_effects() redo it from the pedigree given, so one data frame gets one answer everywhere:

• an animal of the data (in any relationship term, on any row, a maternal grandsire included) with no line in the pedigree enters as a FOUNDER, with unknown parents, and a warning counts them by column and names the first. renumf90 numbers every animal with a record the same way; ASReml asks for every animal of the data in the pedigree file. When NO id of a column is in the pedigree the two use different ids, and the fit stops with the first id of each side;

• the unknown code 0 (or the text "NA") in the column of a relationship term is an unknown animal, as in BLUPF90, whose renumf90 readme warns that "missing animals have code 0": the effect is absent from that record (a level 0 in a BLUPF90 data file adds nothing to that effect), and the record stays in the other terms with zero incidence in this one. So a record with an unknown dam keeps its direct effect, and a message counts those records. The same 0 is absent too in a pe() or random() term that reads the same column (the maternal permanent environment of maternal(dam) + pe(dam), the pe(id) of animal(id) + pe(id)): there it is never one shared level of unknown dams. In a column that no relationship term reads, 0 stays an ordinary level: pe(dam) with no maternal(dam) puts every record with dam 0 in one level "0" (give those dams ids of their own, or add the relationship term, to avoid it). A term whose every row holds 0 stops the fit; in sire(mgs =) a record with an unknown sire keeps only the grandsire contribution, and none when the grandsire is unknown too, and the message counts each. The fit refuses a missing value in a relationship column (NA, blank text) as in any level column, and tells you to write an unknown parent as 0;

• an id with no line that matches a line written another way (spaces, tabs, carriage returns or no-break spaces at the ends, upper against lower case, leading zeros, a trailing .0, accents composed against decomposed, latin1 bytes read as UTF-8) stops the fit with both spellings and the difference: it is a mismatch, not a new animal. The same text in latin1 and in UTF-8 is one id;

• a numeric column that gives levels (any id column, and also a class such as a numeric contemporary-group or herd code, a nested class, a pen or a maternal grandsire), of the data, of the pedigree or of genotypes$ids, holding values at or above 2^53 (the 16- and 17-digit ids of electronic ear tags read by read.csv()) has already lost digits: it is refused, with the column and how to read it as text. A covariate is a number and is not checked. An integer64 column (data.table::fread()) is read as its exact digits;

• the pedigree is PRUNED to the animals of the data (every row, the trait observed or not), the genotyped and all their ancestors, the informative set renumf90 extracts ("genotyped and phenotyped animals + their ancestors traced back") with a PED_DEPTH deep enough to reach every ancestor (its default, 3, keeps three generations; 0 loads the whole pedigree), and a message counts the lines left out. The pruning is exact for the likelihood and for every breeding value kept. An animal left out (no record, no genotype, no recorded descendant: a young candidate) has NO breeding value in the fit: give it a row of the data with the trait missing, its other columns filled, to keep it. A national pedigree far larger than the data no longer reaches the factorization;

• when the records form blocks with no genetic link between them (no relationship through the pedigree nor the genotypes) and the fixed effects give a block levels of its own, a warning says so: any genetic difference between the blocks is absorbed by those effects, and the breeding values compare only within a block (Fernando, Gianola and Grossman 1983; Kennedy and Trus 1993). solutions() then has a block column

genotypes

list with ids and m for single-step: animals x markers coded 0/1/2, as a double, integer or raw matrix. NA (5 in a raw matrix, the BLUPF90 code) is imputed with the marker mean, never converted to zero. The matrix is read where it is, without a copy, so a raw one costs 1 byte per genotype (read_blupf90_snp() and read_plink(storage = "raw") return it). Every function that takes genotypes = checks the same things and refuses, naming the animals: genotypes without animals or without markers (pass genotypes = NULL for a fit on the pedigree alone), a missing or repeated id in ids (one genotype per animal), rownames of m that contradict ids (the genotype goes to the animal at its position in ids), and, everywhere a missing genotype is imputed by the mean, an animal with a call rate below 0.90: such a row is mostly marker means and looks like an average, accurately known animal. The 0.90 is the callrateAnim default of preGSf90, but the rate is counted over the markers that vary (the ones that enter G), where preGSf90 counts every SNP: on a monomorphic marker the mean imputes the genotype every animal carries. The output of qc_genotypes() always passes this check; options(uaitrem.min_call_rate_animal =) moves the line, and 0 turns it off. With a pedigree, a genotyped parent and offspring with opposite homozygotes on more than 1% of the markers called in both (the exclusion test of preGSf90, exclusion_threshold 1) are a warning naming the pairs, and the fit goes on with the data as given; see parentage_check(), which also says how to turn the check off and why sex-chromosome markers must be left out of that test. An optional element genotypes$sex_markers (column indices, names or one logical per column, as qc_genotypes() returns it when given sex_markers =) names them; they stay in G

blend

weight of A22 in the adjusted G, the usual 0.05

apy_core

ids of the genotyped animals that form the APY core; with it the inverse of G* is the APY approximation (cost in the size of the core, not cubic in the genotyped) and the result message SAYS it is an approximation and with which core (its size is part of the result). "auto" chooses the core with apy_core_select(): as many animals as eigenvalues of G explaining 98% of its trace, drawn at random

missing_code

missing-value code for observations, for example -999. It applies to the trait only: in a level column the code is an ordinary level, and in a covariate a cell at the code is refused when it lies far from the other values and warned about otherwise (see the missing-value rule under formula)

maxiter

maximum number of iterations of the damped step. The default 300 was raised from 100 after a measured case: a model warm-started from the reduced fit still had relDelta 1.6e-4 at iteration 100, no defect, a model that walks slowly along a covariance boundary. A fit that hits the ceiling says so in message and reports converged = FALSE, with one exception: when the Newton decrement at the last point is below its tolerance, the fit is certified on a flat direction of the likelihood. Then converged = TRUE at the ceiling, and message says that the components were still moving. The walk also stops before the ceiling when 25 iterations in a row give no gain above 1e-8 of the -2logL and no EM step moves, with the same verdict as at the ceiling. maxiter = 0 evaluates the start (after the n_em EM rounds) and always reports converged = FALSE: an evaluation is not a convergence, whatever the Newton decrement at that point, which message gives

tol

RELATIVE tolerance on the components, sqrt(sum delta^2 / sum theta^2). BLUPF90 note: airemlf90/blupf90+ VCE test the SQUARED quantity, so their conv_crit equals this tol squared (their 1e-10 is tol = 1e-5 here; this 1e-8 default is 1e-16 on their scale). A small step alone never certifies convergence: converged = TRUE additionally requires the Newton decrement of the free components, g' AI^-1 g restricted to the components not held at a boundary, to fall under 2e-4, near the optimum the decrement is about twice the -2logL gap to it, so the certificate bounds that gap by ~1e-4. The value is reported in newton_dec. Measured motivation: a fit that stalled against the singularity boundary with relDelta 5.1e-9 and the score far from zero sat 6.9 -2logL units above the optimum, and the step criterion alone declared it converged. A component held at a boundary is left out of the decrement only when the boundary passes the Karush-Kuhn-Tucker condition of a constrained optimum. Here that is the sign of the score: a variance at its floor is excluded only with the score pushing it out of the admissible cone, and the step never freezes a component whose score points inward. model_mt() checks it on the likelihood instead, because its step also freezes a block whose curvature has collapsed, whatever the sign (see tol there). The variance floor of the step is 1e-10 of the variance of y, fixed for the whole walk, so the walk does not depend on the unit of y (it was 1e-10 in the user's unit whenever the residual variance was below 1: with y times 1e-5 the permanent-environment variance stuck on it, 27% above its optimum, and the fit was declared converged with the additive variance 14% low)

n_em

EM iterations before the AI, to land in the right basin

metafounders

labels of unknown-parent groups; a parent with one of these labels needs no line of its own (any OTHER cited-without-line parent is still a declared error). Metafounders enter as virtual base rows of A(Gamma) after Legarra et al. (2015)

gamma

the base relationship matrix of the metafounders. Either a vector of length n_metafounders, read as the DIAGONAL, or a full symmetric n_metafounders x n_metafounders matrix. The off-diagonal gamma_jk is the ancestral relationship BETWEEN two base populations, which is what a multibreed analysis turns on; it may be NEGATIVE, for bases pulled apart by selection in opposite directions. Admissibility is positive SEMI-definiteness, tested by a spectral decomposition and not by a Cholesky, plus a diagonal below 2 so that a metafounder's offspring keeps a positive Mendelian variance. A fit needs gamma positive-definite: with a SINGULAR gamma (gamma = 0, or two metafounders whose rows are identical) a metafounder, or a combination of them, has no prior variance and leaves the mixed model equations, and the fit refuses it with the cause. For the limit of unknown parent groups give a small diagonal (gamma = 1e-6); for two labels of one population, use one label. a_inverse() takes a singular gamma, through the Moore-Penrose pseudo-inverse. An INDEFINITE gamma, a negative eigenvalue, is refused everywhere

verbose

print the fit as it walks: one line per AI iteration with the -2logL and the relative step, so a long fit is a progress report instead of silence. The relative step is half of the convergence criterion; the Newton decrement, reported in newton_dec, is the other half. Defaults to interactive(), live in a session, quiet in scripts and checks. Every fitter also honors Ctrl+C now

weights

a column of data, or a numeric vector: a record of weight w has residual variance s2e / w. Weights enter as a row scaling by sqrt(w), so the normal equations solved are the weighted ones. Use them when records are means of different sizes, or estimates that carry their own precision, a two-step analysis, a de-regressed proof. The -2logL leaves out the term -sum(log(w)) of log|R|. So the same model written with weights c * w and residual variance c * s2e gives a -2logL larger by n * log(c), with n the records (measured: 27.73 for c = 2 and 40 records). Compare the -2logL of fits that have the same weights only

start

starting values for the components, in the order the fit reports them. Use it to warm-start from a submodel, or to check that the optimum does not depend on where the search began. Without it the start comes from var(y)

residual_by

name of a column of data whose classes (an age class, a period) get their own residual variance, var(residual|class) for each class, in the order of the sorted labels; every record with a response needs a class, and there must be 2 classes or more. The variances are estimated with the others, and h2() reads the residual at the mean over the records of the class variances (the rule of the residual of a fit with weights, which this argument does not take)

Value

an object of class breeding_fit. Besides the fields below it carries used, one logical per row of data saying whether that record entered the equations: n_used counts them and used says WHICH, which is what anything that has to sum on the same base as the likelihood needs. dense_block is c(dense = k, columns = n): the last k of the n columns of the Cholesky factor are completely full, and each factorization costs about k^3 / 3 there. In a single step without apy_core =, k is at least the number of genotyped animals; with it, at least the core size plus one. The fill of the pedigree itself can add to it, and in a small herd it dominates. In a single step it carries h_prior, the diagonal of G* of each genotyped animal, named by animal, and h_prior_row, the row of that animal in the pedigree the fit built: the prior variance accuracy() uses for a genotyped animal instead of 1 + F. It also carries genomic, list(blend), which with apy, metafounders and gamma is what snp_effects() reads to rebuild the same G*. It carries, one value per row of data (NA where the record did not enter), the conditional residual residual, the fitted values fitted (X b + Z u) and fitted_fixed (X b) and the leverage leverage, which residuals.breeding_fit() reads, with missing_code and weights as given. With a relationship term it carries pedigree_used, list(lines, founders, kept): the lines of the pedigree given, the ids of the data it added as founders and the animals kept by the pruning (see pedigree), which accuracy() and snp_effects() read to rebuild the same pedigree. When the records form blocks whose genetic difference the fixed effects absorb it carries blocks, list(n, records, absorbed_by, animal): the number of blocks, the records of each, the fixed terms that give each block its own levels and the block of each animal. The object holds the components theta with their se, the fixed-effect solutions b (named term=level, in the order the columns of X entered), ebv and pev per covariance group, score, vcov, the convergence fields (converged, iters, reldelta, and newton_dec, the Newton decrement of the free components at the final point) and message. vcov is 2 AI^-1, the sampling covariance of the components. A component held at a boundary (a variance at zero, kernel(fixed =)) has no standard error, and the others come from the information conditional on it. When a covariance group stops at a wall (a correlation at +/-1), vcov is the covariance conditional on that wall: 2 J_F (J_F' AI J_F)^-1 J_F'. Here J is the derivative of theta in the coordinates of the step. F is the set of coordinates off the wall, so vcov has one rank less. PARAMETRIZATION OF b, and it matters when comparing against a book or another program: the model carries an implicit intercept and drops linearly dependent columns (listed in dropped_x; a dropped level has solution zero). COVARIATES ARE STANDARDIZED INSIDE: a cov() column enters the equations as (x - mean) / SD of the used rows, as WOMBAT deviates covariables from their means (Meyer, WOMBAT manual, 2017, section 4.8.1.1), and b comes back in the unit of the data, the intercept taking the centre, so it equals the fit of hand-centred data back-transformed, and -2logL keeps its value. cov(x, nested = "c") is centred too when c is also a fixed class of the model (the levels of c, or through the reference level the intercept and the other levels, take the centre). A date, a POSIXct moving by minutes or a timestamp in milliseconds is thereby fitted at full precision: before, age + 1e9 left the starting theta inadmissible, age + 3e8 gave a regression of 5.90 against 0.0102, and age + 1e11 or age x 1e-12 went to dropped_x. When c is not a fixed class of the model, the nested covariate is only scaled, and a warning names it when its values sit more than 1e3 SDs from zero. When c is a fixed class aliased with other fixed effects, so that the rank test drops more than one level of c that has used records (for example, a herd in which c is nested, fitted before c; a level with no used record is in dropped_x too, but it does not count), the nested covariate is also only scaled, with no warning, and a very large offset (around 1e7 SDs from zero) can still make the fit fail as it did before (an inadmissible or singular starting theta). Fit c before the class it is nested in, or leave that class out (its columns are already in the span of c). A covariate column with no variation beyond rounding (SD up to 100 machine epsilons of its largest value) is left as it is. The columns are taken in the order they entered, over the records that enter the fit (with weights, scaled by the square root of the weight), on the standardized scale of the covariates, and a column goes when its distance to the span of the columns kept before it is below 1e-9 of its norm (1e-9 itself when the norm is below 1), so of a dependent set the LAST column goes: the last level of a class against the intercept, a level with no used record, a covariate that is a combination of earlier ones. No more columns are kept than there are records. The design X is sparse, one entry per record and fixed term, but finding its rank keeps a dense triangular factor of p(p + 1)/2 doubles for p kept columns, with O(p^3) work in one thread: of the order of 0.4 GB and minutes at 10 000 fixed columns. The factor is sized before that work for p = the candidate columns with an entry on the used records (or the number of records, if smaller), dependent ones included, and the fit stops with a message when that passes 4 GB (fit the largest class as random). A column that nearly collinear columns before it span only with very large coefficients (x1, a near copy x2 = x1 + 1e-8 z, then x1 + 1e8 (x2 - x1)) sends the rank to a dense check of up to 8 x records x columns bytes and O(n p^2) work in one thread, which also stops with a message above 4 GB; dropping the near copy avoids it. Unscaled polynomial covariates, or near copies on their own, stay on the sparse route. A program that instead zeroes some other level – Mrode's examples zero one level per factor and fit no intercept – agrees with b only on CONTRASTS, differences between levels of the same factor, never on the raw values.

Structure checks

Before fitting, every fitter that estimates components (model(), model_mt() and gibbs()) reads the design and STOPS on the ones the data cannot identify, naming the term, the levels and what to do:

• a random term whose levels are a function of the levels of a fixed class effect (each level of the fixed effect holds the records of a single level of the term: one sire per contemporary group, a fixed litter with maternal(dam) or pe(dam) whenever each litter has one dam, whatever the number of litters per dam, cg = id), a random term whose records all carry the same single level (the intercept absorbs it; a pe() or random() term with one known level and the other records on the unknown 0 of a relationship column is not absorbed and only warns, see below), or a term whose records all carry the unknown code 0 in its column. The incidence of the term is then a sum of columns of X (or zero), the fixed effect absorbs the whole term and its variance does not enter the REML likelihood. ASReml lists a random term confounded with a fixed one among the causes of a singular average information and stops unless !AISINGULARITIES is given. The unknown 0 of a relationship column has no incidence, so a level of the fixed effect that holds only unknown parents does not count against the nesting. In model_mt() the check reads the rows with some trait observed and then the rows of each trait, since the fixed effects are per trait. The check takes one fixed effect at a time; a confounding by a combination of fixed effects is left to the SINGULAR verdict of the fit (see h2() for how the extractors read it).

At given components (maxiter = 0 with n_em = 0 here, maxiter = 0 in model_mt(), theta_fixed = in gibbs(), and snp_blup()) every design above is a warning instead: the BLUP exists, and the data say nothing about the term beyond the components given. maxiter = 0 with n_em > 0 moves the start by EM, which estimates, so it stops.

The fit goes on with a warning, also written in message, when:

• an animal() term has repeated records of some animal and no pe() or random() term reads the same column, or a column that groups the records the same way: the permanent-environment covariance between records of the same animal goes into the additive variance and inflates h2, by more the larger the share of repeated records. The warning counts the repeated animals and their records, and the advice is to add pe() of that column, the repeatability model; when less than 1% of the records repeat an animal it says that h2 moves little and asks to check for duplicated rows;

• a pe() or random() term has fewer than 5 known levels among the records: its variance comes from that many levels and is barely determined, and labels that number something inside another unit (the parity within each dam used as the litter) merge units. A term with one known level whose other records carry the unknown 0 of a relationship term that reads the same column (maternal(dam) + pe(dam) with a single known dam) warns too, counting the records of that level: the unknown 0 has no incidence, so the intercept does not absorb the term.

These checks read the rows with an observed response that enter the analysis.

References

Patterson, H.D. & Thompson, R. (1971). Recovery of inter-block information when block sizes are unequal. Biometrika 58:545-554.

Gilmour, A.R., Thompson, R. & Cullis, B.R. (1995). Average information REML. Biometrics 51:1440-1450.

VanRaden, P.M. (2008). Efficient methods to compute genomic predictions. Journal of Dairy Science 91:4414-4423.

Aguilar, I. et al. (2010). A unified approach... Journal of Dairy Science 93:743-752; Christensen, O.F. & Lund, M.S. (2010). Genetics Selection Evolution 42:2.

Misztal, I., Legarra, A. & Aguilar, I. (2014). Using recursion to compute the inverse of the genomic relationship matrix. Journal of Dairy Science 97:3943-3952.

Fernando, R.L., Gianola, D. & Grossman, M. (1983). Identifying all connected subsets in a two-way classification without interaction. Journal of Dairy Science 66:1399-1402.

Kennedy, B.W. & Trus, D. (1993). Considerations on genetic connectedness between management units under an animal model. Journal of Animal Science 71:2341-2352.

Misztal, I. et al. The BLUPF90 family of programs, readme of renumf90 (PED_DEPTH: the genotyped and phenotyped animals and their ancestors, default 3, 0 for the whole pedigree; Warnings: missing animals have code 0). University of Georgia, https://nce.ads.uga.edu/wiki/doku.php?id=readme.renumf90