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,
  vecchia_k = NULL,
  missing_code = NULL,
  maxiter = 300L,
  tol = 1e-08,
  n_em = 4L,
  metafounders = NULL,
  gamma = NULL,
  verbose = interactive(),
  weights = NULL,
  start = NULL
)

Arguments

formula

for example peso ~ cg + sexo + animal(id). An unmarked term is a fixed class effect; 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, the ids of K in a kernel() 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. 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, the pen of indirect()) 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(), model_ar1() (subject = included), gibbs(), snp_blup(), the eval_internal family, model_threshold(), model_survival() (subject = included) 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(), model_ar1(), gibbs(), snp_blup(), the eval_internal family, model_threshold(), model_survival() 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 (not for the survival time, whose tail is long by nature, nor for the categories of an ordinal threshold trait, which have no absurd value). 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. describe() reports the same findings. A NUMBER STAYS A NUMBER: a factor or text column read as a number (the trait, a covariate, the time and entry= of model_survival(), the quantitative trait of the joint threshold mode, the traits of pegs(), the time= and entry= of predict()) is refused with the column, and not read as its codes 1, 2, ... indirect(id, pen = "pen") is the indirect (associative) genetic effect of Mrode & Pocrnic (2023, ch. 9): the incidence of a record marks the animal's DISTINCT pen mates, and it shares a group with animal() so the direct-social covariance is estimated. Its dilution = d argument scales the entry of every mate to (n_i - 1)^(-d), where n_i is the number of distinct animals in the pen of record i: d = 0, the default, is the book's plain sum (coefficient 1 per mate, exactly the old behaviour); d = 1 is the mate mean; any d >= 0 in between is the dilution of Bijma (2010). With pens of unequal size the sum grows with n_i - 1 and the choice of d is an empirical question: fit a small grid (say 0, 0.5, 1) and compare -2logL, which is comparable across d because only the incidence changes. A pen of size 1 keeps its zero social row under every d. The pen column may be text, factor or numeric (codes 10, 20, 30 group exactly like "10", "20", "30"), but it may not hold a missing pen: NA, NaN, Inf, -Inf, blank text (what read.csv() gives an empty cell of a text column, spaces included) or the text "NaN". A record without a pen has pen mates nobody knows, so every fitter stops, names the column, the row count and the first row, and asks to drop those rows or assign them a pen. See indirect_residual() for the residual side of the same problem. kernel(id, K = D) is a random term with a user-supplied (user-defined) covariance matrix, called a DECLARED covariance throughout this package: K is a symmetric positive-definite matrix whose rownames are the level identifiers, and every row of K gets an equation, with or without a record, a dominance D (dominance_matrix(), g_dominance()), an epistatic G_AA (g_epistasis()), a partial multibreed matrix (partial_a()), or any relationship the pedigree and the markers do not already provide. A row of K that is ENTIRELY zero, diagonal included, declares a level with no contribution to this term: it gets no equation and its records stay in the analysis with zero incidence here, the generalized-inverse pattern of the multibreed partial matrices (Mrode & Pocrnic, 4th ed., p.243-244). An id absent from K altogether still excludes the record (an animal missing from the PEDIGREE enters as a founder instead, see pedigree): a zero row is a declaration, an absence is a gap, and when no record is left the fit stops with the first id of each side. Two kernel terms need nome= to tell their components apart. The inversion of K is dense, so K = is for matrices of moderate size, the size of a genotyped set, not of a national pedigree. kernel(id, Kinv = Q) declares the PRECISION instead, Q = K^-1, and nothing is inverted: Q is the lower triangle in triplets, list(i, j, x, n, id) as a_inverse() returns (i >= j), or a dense matrix with the level ids as dimnames. It is the route of dominance_inverse(), the sparse inverse of the dominance matrix by sire x dam subclasses: kernel(id, Kinv = dominance_inverse(ped)) has the component, -2logL, score and average information of kernel(id, K = dominance_matrix(ped)) (what the tests check, and in which fitter, is on the page of dominance_inverse()), with the subclass levels as latent effects, so ebv() of the term returns one value per level, animals and pairs; the fit carries k_level_type, the type of each level by term ("animal", "subclass", "ancestral"), to tell them apart, and solutions() adds a type column and lists the animals first. A record with an observed response of an animal that dominance_inverse(animals = ) left out is refused, not dropped, and so is a record whose id is the label of a pair level. A dense Kinv follows the rules of a dense K (square, named, symmetric, finite); the range, the lower triangle and the finiteness of triplets are checked by the engine in one pass. Exactly one of K and Kinv is given. The marker arguments that take a value and not a column name, base =, dilution =, K =, Kinv = and the fixed = of kernel() (which holds that term's variance at the value given), are evaluated in the environment of the formula, the one where it was written. A grid over d run as lapply(c(0, 0.5, 1), function(d) model(y ~ cg + animal(id, group = "g") + indirect(id, pen = "pen", group = "g", dilution = d), data, ped)) finds its d, and so does a formula built in one function and fitted in another. The fit keeps in formula the formula with the values of base =, dilution = and fixed = written in, so h2(), t2() and accuracy() read the value the fit was made with even after the variable changes or is removed. The arguments that NAME something, group =, pen =, 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(), model_ar1(), gibbs(), snp_blup(), the eval_internal family, model_threshold() and model_survival()), and accuracy(), solutions(), snp_effects() and residual_diagnostics() 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 pen mate of indirect() and 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(), random() or kernel() 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; in indirect() every pen mate must be a known animal; and the subject column of model_ar1() and of model_survival() cannot hold it when it is the column of a relationship term, since distinct unknown animals would become one series or one history. A missing value there (NA, blank text) is refused like in any level column, with the advice 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, a maternal grandsire or the subject of model_ar1()), 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, and describe() counts the blocks before the fit

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(BreedingR.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

vecchia_k

neighbors per animal in the Vecchia inverse of G*: each genotyped animal conditions on its k strongest previous relationships instead of a global core (the generalization of APY; Henderson's A^-1 is the pedigree case, with the parents as the conditioning set). k >= n - 1 reproduces the exact inverse; the result message says it is an approximation and with which k. Mutually exclusive with apy_core

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 direct-indirect 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. 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() and model_ar1() check it on the likelihood instead, because their 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 SINGULAR gamma is accepted, through the Moore-Penrose pseudo-inverse: that covers gamma = 0, the unknown-parent-group limit, where the pseudo-inverse reproduces the A-inverse of unknown parent groups exactly, and two metafounders standing for one population, whose rows are identical. What is still refused is an INDEFINITE gamma, a negative eigenvalue, which does not generate a covariance matrix at all

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

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)

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. With a kernel() term it carries k_prior, one vector per kernel term, named by term, holding the diagonal of the declared K by level: the prior variance accuracy() divides the PEV by. With a kernel(id, Kinv = dominance_inverse(ped)) term it also carries k_level_type, the type of each level by term ("animal", "subclass", "ancestral"). 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, vecchia_k) (vecchia_k = 0 without Vecchia), 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() and residual_diagnostics() read, 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. 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(), model_ar1(), gibbs(), model_threshold() with estimate = TRUE and model_survival() without sigma2 =) 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. A kernel(fixed = ) term, whose variance is given, is not checked. 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). model_survival() with lambda given warns instead: the first level of the fixed effect is held at the baseline, and the single level of the term in it is all that informs the variance;

• indirect() with every pen of one animal (the indirect incidence is all zero), with all the records in one pen, or with the pen (or anything finer) as a fixed effect: the fixed part absorbs the sum of each pen, and what is left of the indirect incidence of a record is minus its own indirect effect, which the direct effect of the same animal cannot be told apart from (Cantet and Cappa 2008, J. Anim. Breed. Genet. 125:371; Bijma 2010, Genetics 186:1013). The message suggests the pen as random(pen), with the caveat of Cantet and Cappa (2008) that a random pen does not always remove the collinearity: what separates the two effects is families spread over several pens. With dilution > 0 and pens of different sizes only the dilution weights separate the two, a very weak identification, and that is a warning.

At given components (maxiter = 0 with n_em = 0 here, maxiter = 0 in model_mt() and model_ar1(), theta_fixed = in gibbs(), estimate = FALSE in model_threshold(), sigma2 = in model_survival(), 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 (with every pen of one animal, the indirect breeding values come only from their covariance with the direct ones). 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(), random() or kernel() with a diagonal K and an estimated variance reads the same column, or a column that groups the records the same way (nor, in model_ar1(), is the animal the subject): 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 the one of suggest_model(); when less than 1% of the records repeat an animal it says that h2 moves little and asks to check for duplicated rows;

• fewer than 10 pens of indirect() hold more than one animal: the indirect components rest on that many pen sums (Bijma 2010 puts the detection of indirect effects at about 250 to 500 groups even with an optimum design);

• the pens of indirect() split by the levels of a fixed class effect in a way chance does not explain, which is how pens merged by a reused label look: when pen numbers restart within each batch, pens that never met are merged into one and every record gets mates it never had; give each pen a unique label, such as paste(batch, pen). The test: at least 90% of the pens that hold more than one level split into groups of two or more animals per level, and that count exceeds what animals placed in the pens regardless of the effect would give (each level at its overall proportion; the exact multinomial chance of a pen with no group of one) with P < 0.0001. Only a fixed effect that can be the batch is read: no animal has records in two of its levels (the period or age class of each record is not read), and it has three or more levels among the pens with mates, or another fixed class effect with more levels is nested in it (the contemporary group within the batch). A two-level effect of the animal, such as sex, is not read (unless a fixed class effect is nested in it), whether it is drawn at random within the pens or balanced by design in every pen. The data cannot tell merged pens from pens that truly mix levels when the groups are large, or when two batches have nothing nested in them: two batches of merged pens then look like pens balanced for a two-level effect and are not read. Measured on simulated data (validation/structure_cases.R, 50 labels, every animal present): merged pens of 2 or 3 over three or more batches are detected, also when only a contemporary group nested in the batch is in the formula (there from 30 labels); merged pens of 5 over three or more batches are detected only through a contemporary group nested in the batch; over two batches merged pens of 2 or 3 are detected only through such a group (pens of 3 from about 45 labels), and merged pens of 5 are not, since they look the same as a contemporary group that includes sex in pens that mix the sexes. Animals missing from the data leave groups of one: with 5% of the animals missing merged pens of 2 are not detected, and with 20% neither are merged pens of 3. An effect of three or more levels is read only when its levels meet in the mixed pens the way merged batches do (most pairs of levels share some pen): a contemporary group defined within sex in pens that mix the sexes (the sexes of one batch meet, never two batches) and a birth-week group in pens filled in birth order (neighbouring weeks meet) are not read. An effect of the animal with three or more levels balanced in every pen, with two or more animals of each level per pen (three sex classes in equal numbers in pens of 6), does look like merged pens and is reported (pens of 3 with one animal per class are not). options(BreedingR.warn_pen_labels = FALSE) silences this check;

• 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 (a row whose kernel() level has no line in K leaves it), and, for indirect(), the pen membership of every row, as the incidence does.

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.

Bijma, P. (2010). Multilevel selection 4: modeling the relationship of indirect genetic effects and group size. Genetics 186:1029-1031.

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