BreedingR

Genetic parameters, breeding values and genomic prediction in R. One engine fits the animal model and a good deal beyond it: direct-maternal, reaction norms, indirect genetic effects, several traits at once, a residual with serial correlation, ordered categorical and censored traits, and any relationship matrix you supply yourself.

It is where I try out new ideas and put specific models into practice: reaction norms on an environmental gradient, indirect genetic effects in group housing, a residual that carries serial correlation. That is still what it is for. The numerics are written in the package itself, in src/, and R CMD INSTALL compiles them. No separate binary, no service, no run-time dependency.

Models	Genomics	Estimation and tools
animal model, repeatability direct-maternal; sire and maternal grandsire reaction norms on a gradient indirect genetic effects, with group-size dilution several traits; AR(1) and CAR(1) residuals ordered categories; Weibull survival with censoring dominance, epistasis, multibreed: any declared covariance	single step: exact, APY or Vecchia inverse of G ssSNPBLUP, which never forms G metafounders, Gamma estimated from the genotypes marker effects from a single-step fit pegs() for many environments of one trait Fst and runs of homozygosity genotype QC and parentage; PLINK and BLUPF90 files	AI-REML with a convergence certificate block Gibbs sampler, probit for binary traits breeding values, accuracy, h2, genetic correlations selection index, ranking drift residual and MCMC diagnostics gene-dropping simulator, Monte Carlo studies C++ engine; the same bits with any number of threads

The map below lists each model with the term that asks for it.

Installation

# straight from GitHub with base R only; it compiles, so Windows needs Rtools and macOS
# the Xcode command-line tools. Downloading first keeps RStudio from trying to read the
# URL as a local file (harmless "tar.exe: Error opening archive" warnings otherwise).
f <- file.path(tempdir(), "BreedingR.tar.gz")
download.file("https://github.com/phyllype/BreedingR/archive/refs/heads/main.tar.gz", f,
              mode = "wb")
install.packages(f, repos = NULL, type = "source")
# or with remotes, installed first with install.packages("remotes")
remotes::install_github("phyllype/BreedingR")
# or from a local clone, giving the folder that holds DESCRIPTION
install.packages("C:/path/to/BreedingR", repos = NULL, type = "source")

library(BreedingR)
packageVersion("BreedingR")   # the examples here need 0.4.0.9000

The CI installs the package both ways from GitHub on every push and runs the quick start below with the installed copy.

To replace an older version, first close every R session that has loaded BreedingR (RStudio included: quit it, or restart R and do not call library(BreedingR)), then run the install before anything else. Windows keeps the DLL of a loaded package locked, and the install stops with “cannot remove earlier installation, is it in use?” or “Permission denied”, leaving the old version in place.

What it does

The package covers the chain end to end. The pedigree becomes the relationship matrix A and its sparse inverse by Henderson’s (1976) rules, with inbreeding by the Meuwissen and Luo (1992) trace and, where the base population is not one homogeneous pool, metafounders (Legarra et al., 2015). A pedigree of sires and maternal grandsires is declared with sire_mgs(), never guessed. The model is written as a formula; the mixed model equations are assembled sparse, one record at a time, and factored by a sparse Cholesky after a minimum degree ordering (George and Liu, 1989) that sets dense nodes aside, as the dense-row rule of AMD does (Amestoy, Davis and Duff, 1996). In model(), model_mt(), model_ar1() and gibbs() the ordering and the symbolic analysis are computed once and reused; model_threshold() and model_survival(), written in R, factor again at every step but reuse the ordering and the symbolic analysis of the pattern. The dense tail of the factor, where the genotyped animals end up, is factored in tiles on br_threads() threads (OpenMP, default 1), and so are its inverse inside the selected inverse, the dense inverses of G* and A22 and the product Z Z’ of the G in the single step; each number has one owner thread and a fixed summation order, so the result is the same bit for bit with any number of threads. The variance components come from AI-REML (Gilmour, Thompson and Cullis, 1995): analytic score, average information, EM warm-up and an EM rescue (Dempster, Laird and Rubin, 1977), a damped step that walks in log-Cholesky coordinates so a covariance boundary is a limit rather than a wall, and convergence that takes both a small RELATIVE step and a Newton decrement under tolerance, so a stalled step cannot pass for an optimum. Breeding values fall out of the same solution, and their accuracy out of the selected inverse (Takahashi, Fagan and Chin, 1973), over the prior variance each animal has: 1 + F, or the diagonal of G* for a genotyped animal in a single step.

A-inverse of a simulated pedigree of 180 animals and its Cholesky factor, fill-in in amber. In pedigree order the factor fills in to 11,961 nonzeros; after the minimum degree ordering, 1,142.

Genotypes enter as an argument to the same fit: G by VanRaden (2008), brought to the scale of A22 by an affine adjustment and a blend, and the single-step H^-1 built as A^-1 plus a correction on the genotyped block. When the genotyped set is large enough that inverting G* hurts, the same fit accepts APY (Misztal, Legarra and Aguilar, 2014) with a core, chosen by hand or by apy_core = "auto", which takes as many animals as eigenvalues of G explain 98% of its trace (Pocrnic et al., 2016), and builds it without any matrix of the size of the genotyped set squared (G only on the core rows, A22 by Colleau’s algorithm, A22^-1 by a sparse Schur complement); the Vecchia (1988) recursion with per-animal conditioning sets; or snp_blup(), which never builds G at all and solves the marker equations by conjugate gradients.

Several further model families are formula terms here, because the unit of layout is the covariance group rather than the term: direct-maternal (Willham, 1972), reaction norms on an environmental gradient (Kirkpatrick, Lofsvold and Bulmer, 1990), indirect genetic effects among pen mates (Griffing, 1967; Muir and Schinckel, 2002), multi-trait with a full residual covariance, and an AR(1)/CAR(1) residual for repeated measures (Wade and Quaas, 1993). The same models can be sampled instead of maximized, through a block Gibbs sampler (Geman and Geman, 1984) over the same equations.

The trunk is Gaussian, but the trait does not have to be. An ordered categorical trait fits on a probit liability with model_threshold(), alone or jointly with a quantitative one, with the components given or estimated as the minimum of the Laplace approximation of the marginal likelihood, indirect genetic effects among pen mates included; gibbs(family = "probit") samples a binary trait by data augmentation (Albert and Chib, 1993; Sorensen et al., 1995). Time until failure fits with model_survival(), where a right-censored record enters as a lower bound rather than a missing value and a covariate that changes during a life enters as elementary records. model_threshold() and model_survival() take genotypes= for a single step, through h_inverse(). And a random term can carry any user-supplied covariance matrix through kernel(id, K = ), or its inverse through kernel(id, Kinv = ), which is how the sparse inverse of the dominance matrix from dominance_inverse() enters, called a DECLARED covariance throughout this package: the dominance and epistasis constructions of chapter 13 of Mrode and Pocrnic (2023), and the multibreed partial matrices of chapter 14, are matrices built by their own constructors and handed to the same engine.

The map

What the package fits, and the term that asks for it. Everything in the first table is the same engine and the same kron(C^-1, K^-1) penalty, which is why none of these has a fitter of its own.

To fit	Write	The fitter
animal model, repeatability	animal(id), pe(id)	model()
direct-maternal, correlation estimated	animal(id, group=) + maternal(dam, group=)	model()
reaction norm on a gradient	rn(id, base=)	model()
indirect (social) genetic effects	indirect(id, pen=, group=, dilution=)	model()
dominance, epistasis, multibreed, any declared K	kernel(id, K=), or its precision kernel(id, Kinv=)	model()
several traits, full residual matrix	cbind(y1, y2) ~ ...	model_mt()
serial correlation in the residual	subject=, time=	model_ar1()
the same models, sampled instead of maximised	same formula, prior=	gibbs()
binary trait, sampled on the liability	same formula, family = "probit"	gibbs()
ordered categorical trait	same formula, components given or estimate = TRUE	model_threshold()
time to failure, right-censored	censor=; entry=, subject= for covariates that change	model_survival(); survival_split() builds the pieces
competitive ability from grouped contests	contest table	competition_strength()

Relationships and genomics are arguments, not different programs.

For	Use
pedigree A-inverse, inbreeding	pedigree(), a_inverse(), a22_inverse()
a pedigree of sires and maternal grandsires	sire_mgs(), or pedigree(type = "sire_mgs"); dams where known and grandsires elsewhere with sire_mgs(ped, dam = "dam")
base populations that are not one pool	metafounders=, gamma= (full matrix, singular allowed); estimate_gamma() from the genotypes
single step	genotypes=, and apy_core= or vecchia_k= when G is large
the APY core, by the eigenvalues of G	apy_core_select(), or apy_core = "auto"
the single-step H-inverse, as triplets	h_inverse()
the single step without ever forming G	snp_blup()
many environments of one trait, fast, from markers	pegs() (residuals uncorrelated across environments)
marker effects from a single-step fit	snp_effects(), with the blend, APY core, Vecchia k and metafounders read from the fit
G, dominance, epistasis to any order	g_matrix(), g_dominance(), g_epistasis_ad/dd/order()
dominance by pedigree, dense or as a sparse inverse by sire x dam subclasses	dominance_matrix(), dominance_inverse()
multibreed partial matrices	partial_a()
the associative residual, exactly	associative_matrix()

And around the fit: solutions(), h2(), t2(), ebv(), accuracy(), rg(), h2_curve(), h2_observed(), h2_liability(), selection_index(), rank_drift(), profile_theta(), se_function(), qc_genotypes(), parentage_check(), qc_phenotypes(), read_plink(), describe(), thi(), heat_load(), fst(), roh(), simulate_breeding(), mc_study(), suggest_model(), benchmark_fit(), residual_diagnostics(), mcmc_diagnostics(), ess(), geweke_z(), rhat().

Where to read next. Start with Your first evaluation, which goes from two files on disk to breeding values you can act on, and assumes nothing about this package. After that: Theory and practice walks a full evaluation in order, explaining each matrix, each algorithm and each iteration alongside the code that runs it; Hands-on works through 58 of the 74 exported functions, step by step; Contest models derives the competitive-ability estimators from the group multinomial, one identity at a time. FUNCTIONS.md maps the whole surface.

Quick start

A complete run on data the package simulates itself, paste and go:

library(BreedingR)
s <- simulate_breeding(n_founders = 60, n_generations = 3,
                       offspring_per_generation = 150, h2 = 0.4,
                       n_markers = 500, seed = 1)      # pedigree + phenotypes + genotypes

q  <- qc_phenotypes(s$data, "y", classes = "cg")        # flag, count, never drop silently
g  <- qc_genotypes(s$genotypes, min_maf = 0.01)         # animals and markers, call rate, MAF

fit <- model(y ~ cg + animal(id), q$data, s$pedigree, genotypes = g)
fit                                # components, SEs, share of the phenotypic variance
h2(fit)                            # the var(animal) share; 4 var(sire) / Vp in a sire model
head(solutions(fit, s$pedigree))   # id, ebv, se, acc, sorted by breeding value
accuracy(fit, s$pedigree)[1:5]     # prior 1+F, G* diagonal if genotyped, K[i,i] for kernel(K=), the fit's prior for kernel(Kinv=)
cor(ebv(fit)[names(s$tbv)], s$tbv) # against the simulator's own truth

The fit above against the simulator’s truth: each point is an animal, colored by the accuracy the fit reports. The figures on this page are drawn by tools/figuras_readme.R from simulated data.

A long fit reports itself as it goes. With verbose = TRUE, the default in an interactive session, each AI iteration prints its -2logL, its relative step and the current value of every component being estimated, not just the residual, so a run that is drifting shows it while there is still time to stop; Ctrl+C interrupts any fitter. What converged means, and why the relative step alone is not enough, is under Choices worth knowing about.

Before the first iteration model() also prints the size of the dense block of the factor, and model(), model_mt(), model_ar1() and gibbs() keep it in fit$dense_block. That k is what each factorization pays, about k^3 / 3 flops: in a single step without apy_core= it is at least the number of genotyped animals, and the fill of the pedigree can add to it.

Quality control before the model

qc_phenotypes() marks the missing code, turns outliers beyond a Tukey (1977) fence into missing values, and reports class levels too small to estimate. It flags rather than deletes: removing a row would reshape contemporary groups and pen compositions without saying so. qc_genotypes() filters animals and markers by call rate (worst first, in passes, so what it keeps passes every fitter), markers by minor allele frequency and Hardy-Weinberg equilibrium, and, with pedigree =, the offspring whose genotype contradicts a genotyped parent; it reports what each filter removed. parentage_check() lists the genotyped parent-offspring pairs with their opposite homozygotes. Every function that takes genotypes = refuses a repeated or missing id and rownames that contradict the ids. The fitters and the functions that build G from the genotypes also refuse an animal below a 0.90 call rate, since its missing markers would be imputed by the mean. The single-step fitters warn when a genotyped pair contradicts the pedigree, leaving out the markers named in genotypes$sex_markers. Every fitter also reads the structure of the design before it fits: a random term nested in a fixed effect, indirect() without pen mates or with the pen as a fixed effect stop with what to do, and repeated records without pe(), few pens and pen numbers reused across batches are a warning (see “Structure checks” in ?model). describe() shows the data before a model touches it, and suggest_model() reads its shape and names the terms it calls for. Before a fit, a missing-value code or an absurd value in the trait or a covariate is named with its first row, and covariates such as dates or timestamps are standardized internally.

Only the observation may be missing. Some columns give a level: a class, the id of a random or relationship term, a nested class, the maternal grandsire, the pen of indirect(), the subject of model_ar1() and model_survival(), and the id and environment of pegs(). They may not hold NA, NaN, Inf, -Inf, blank text (an empty cell of a text column, as read.csv() reads it) or the text “NaN”; in the maternal-grandsire column of sire(mgs =) blank text is the unknown grandsire, like “0”, and only the other values stop the fit. A covariate of cov() or base = may not hold NA, NaN, Inf or -Inf. Every fitter checks every row, including rows whose trait is missing, and stops with the term, the column, the count and the first row, so missing rows are never grouped into one shared level and a covariate never leaves the model in silence. qc_phenotypes(classes =) counts the missing labels before the fit, in missing_labels. missing_code marks the trait, as OPTION missing does in BLUPF90. A covariate cell equal to the code (an unknown age of -999 in a BLUPF90 data file) is refused on every row, like NA, when it lies more than 10 robust SDs from the other cells, and warned about otherwise; a BLUPF90 class code 0 is one level “0” here. Turn those cells into NA before the fit, and the fit names them. The exception is the column of a relationship term (animal(), maternal(), sire(), rn(), the grandsire of sire(mgs =)) and a pe(), random() or kernel() on that same column: there 0 is an unknown animal, as in a BLUPF90 pedigree, and the record stays in the analysis with that effect absent. pe(dam) with no maternal(dam) reads 0 as one ordinary level.

What it fits

# animal model
model(weight ~ cg + cov(age) + animal(id), data, pedigree = ped)

# repeatability: the permanent environment of a subject with repeated records, what
# those records share and is not additive genetic, so it also carries the non-additive
# genetic effects; repeatability is share(animal) + share(pe) in the printed table
model(weight ~ cg + animal(id) + pe(id), data, ped)

# direct-maternal, with the correlation BETWEEN the two estimated
model(weight ~ cg + animal(id, group = "g") + maternal(dam, group = "g"), data, ped)

# the full maternal model (Mrode & Pocrnic, 2023, Eqn 8.1): ONE permanent environment, the
# DAM's, which carries her non-additive maternal genetics as well.
# A SECOND pe() on the animal itself is not part of Eqn 8.1, and with one record per
# animal it IS the residual: the likelihood cannot separate the two. On simulated data
# both models returned the same -2logL, 388.9537, and the extra term only split the
# 0.578 residual into 0.143 and 0.435. On a second simulated set the fit drifted
# instead: var(animal) 0.375 -> 0.279, the direct-maternal covariance -0.031 -> -0.001,
# and it stopped at the zero boundary. That second pe() belongs where the animal itself
# has REPEATED records, and there both pe() must be named (nome=), so that no component
# is renamed by the arrival of another term
model(weight ~ cg + animal(id, group = "g") + maternal(dam, group = "g") +
        pe(dam), data, ped)

# reaction norm on a Legendre basis; pe() needs repeated records, since with one record
# per animal it IS the residual and the fit says SINGULAR. h2() refuses a reaction norm:
# h2_curve() gives h2 along the gradient, with every random regression, a pe() on the
# same basis included, evaluated at the point
d <- cbind(d, legendre(d$thi, order = 1))
model(y ~ cg + rn(id, base = c("phi0", "phi1")) + pe(id), d, ped)

# indirect genetic effects (the associative model of Muir and Schinckel, 2002, and of
# Bijma et al., 2007): the fit returns
# var(animal), var(indirect) and the covariance between them, and the SIGN of that
# covariance is what separates heritable competition from heritable co-operation. The
# response, though, follows the TOTAL breeding value A_D + (n-1) A_S, so reading it
# takes the group size n as well (Bijma et al., 2007). With unequal pens, dilution=d
# scales every mate's entry to (n_i - 1)^(-d) (Bijma, 2010): d=0 is the book's plain
# sum and the default, d=1 the mate mean, and d is chosen by a small grid compared on
# -2logL. The choice is not cosmetic: on pens of 2 to 12 generated with d=1, forcing
# d=0 crushed var(indirect) to ~2% of its true value (0.046 against 2). The same d
# dilutes the social environmental deviation, so indirect_residual() and
# associative_matrix(dilution = d) describe the model the formula declares.
# Whether the three components separate is a property of the design (Cantet and Cappa,
# 2008): with pens of one size built from two full-sib families only two combinations
# are identified, and the fit says SINGULAR and names the components. The share column
# is blank here, because the phenotypic variance of a record depends on its pen size:
# h2(fit, n = , r = ) and t2(fit, n = , r = ) take the size and the relationship
# between mates (Bijma et al., 2007)
model(y ~ cg + animal(id, group = "g") + indirect(id, pen = "pen", group = "g"), d, ped)

# single step (ssGBLUP); the same genotypes=, apy_core= and vecchia_k= work in
# model_mt(), model_ar1() and gibbs(), and in model_threshold() and model_survival(),
# which get their H^-1 from h_inverse(). For the inverse of G*: exact, apy_core= (a
# global core; "auto" sizes it by the eigenvalues of G), or vecchia_k= (per-animal
# neighborhoods, the generalization of APY and of Henderson's own A^-1)
model(y ~ cg + animal(id), d, ped, genotypes = list(ids = gids, m = M))
gen <- read_blupf90_snp("snp.dat")   # the BLUPF90 SNP_FILE as a raw matrix, 1 byte/genotype
core <- apy_core_select(list(ids = gids, m = M))   # choose once, pass it to every fit

# a pedigree of sires and maternal grandsires is DECLARED, never guessed: the grandsire
# path weighs 1/4 (Mrode & Pocrnic, 2023, secs. 3.6 and 3.7), and a third column named
# mgs, mgsire or maternal_grandsire is refused unless declared; sire(sire, mgs = "mgs")
# is the sire and maternal-grandsire model, 1 on the sire and 1/2 on the grandsire
model(y ~ herd + sire(sire), d, sire_mgs(bulls))
model(y ~ herd + sire(sire, mgs = "mgs"), d, sire_mgs(bulls))

# multi-trait with full R0 and missingness by pattern; a fixed level with no record for
# one trait drops as the pair (column, trait), reported in dropped_x as 'CG=5|y1'
model_mt(cbind(t1, t2) ~ cg + animal(id), d, ped);  rg(fit, "animal", "t1", "t2")

# AR(1)/CAR(1) residual for longitudinal data; cbind() on the left fits the
# multi-trait version with the separable residual Gamma (x) R0
model_ar1(y ~ cg + animal(id), d, ped, subject = "id", time = "day")
# rho is per step of the time grid (fit$time_step); a continuous time is CAR(1), rho in [0, 1)

# the Bayesian half: block Gibbs with conjugate updates. prior = "jeffreys" (the
# default), "flat", "uniform_sd" (Gelman, 2006) or a proper c(df =, scale =);
# family = "probit" samples a binary trait on its liability, the residual fixed at 1
# (Albert and Chib, 1993; Sorensen et al., 1995). With a covariance group (direct and
# indirect, direct and maternal) "jeffreys" puts infinite posterior mass at a correlation
# of +-1 (Hobert and Casella, 1996): a chain that gets there stays, even when the data
# identify the group, and gibbs() warns. On a Gaussian trait "flat" or a proper prior
# avoids it; with family = "probit" only a proper c(df =, scale =) does, since the flat
# prior leaves that posterior improper as the variances grow. Fit model() first and read
# its SINGULAR verdict, which the chain does not give
gibbs(y ~ cg + animal(id), d, ped, n_iter = 20000)

# ordered categorical trait on a probit liability (Gianola & Foulley 1983);
# thresholds replace the intercept, predict() gives per-category probabilities.
# The components are GIVEN, or with estimate = TRUE estimated as the minimum of the
# Laplace -2logL, which runs low for a binary trait with few records per level
# (Tempelman, 1998), lower still with indirect(). indirect(), group= and random(pen)
# enter as in model(). cbind(quant, bin) fits the joint analysis of Foulley
# et al. (1983), with pev for u1 and for the ranking value u2, and predict() the
# probability of Eqn 15.25
model_threshold(score ~ herd + sex + sire(sire), d, ped, start = 1/19)

# time until failure with right-censoring: the Weibull frailty model of Kachman
# (1999); a censored record is a lower bound, not a missing value, and censor=
# is mandatory. Solutions are log relative risks; predict() gives RRS and S(t).
# A covariate that changes during a life enters as elementary records (entry, stop]
# of one subject, the device of the Survival Kit; gaps are refused unless
# gaps = "allow", and a late first entry is accepted and flagged
model_survival(lpl ~ herd + ysp + animal(cow), d, ped, censor = "code")

# a random term with a DECLARED covariance matrix: dominance beside the additive
# term (chapter 13), or any K that is neither A nor H. kernel(id, K = , fixed = v)
# holds the component at v in model() and gibbs(); model_mt() and model_ar1() refuse it
model(y ~ pen + animal(id) + kernel(id, K = dominance_matrix(ped)), d, ped)
# the same dominance on a large pedigree: the sparse inverse by sire x dam subclasses
# (Hoeschele & VanRaden 1991, generalized to inbreeding), for the animals with records
model(y ~ pen + animal(id) + kernel(id, Kinv = dominance_inverse(ped, animals = d$id)), d, ped)

# multibreed by pedigree: the partial relationship matrices of Garcia-Cortes &
# Toro (2006), one kernel() per founder breed and per segregating pair
pa <- partial_a(ped, breed = c("1" = "A", "2" = "A", "3" = "B", "4" = "B"))
model(y ~ herd + kernel(id, K = pa$K[["A"]], nome = "uA") +
        kernel(id, K = pa$K[["B"]], nome = "uB") +
        kernel(id, K = pa$K[["A:B"]], nome = "uAB"), d, ped)

# unknown-parent groups as metafounders (Legarra et al., 2015). gamma takes a vector for a
# diagonal, or a full MATRIX whose off-diagonal is the ancestral relationship BETWEEN two
# base populations, which is the parameter a multibreed analysis exists for. A SINGULAR
# gamma is accepted through the pseudo-inverse: gamma = 0 is the unknown-parent-group
# limit, and two metafounders standing for one population have identical rows
model(y ~ cg + animal(id), d, ped, metafounders = c("L1", "L2"),
      gamma = matrix(c(0.7, 0.2, 0.2, 0.6), 2, 2))

# marker effects backsolved from a single-step fit f (model(..., genotypes =)); blend,
# APY core, Vecchia k and metafounders are read from f
snp_effects(f, ped, genotypes = list(ids = gids, m = M))

# the single step WITHOUT G: markers as equations (ssSNPBLUP; Liu et al., 2014), conjugate gradients,
# A22^-1 applied matrix-free; theta is given, as in routine practice
snp_blup(y ~ cg + animal(id), d, ped, genotypes = list(ids = gids, m = M),
         theta = c(0.4, 0.6))

# PLINK .bed/.raw straight into genotypes=, with QC that reports what it removed
g <- qc_genotypes(read_plink("chip"), min_maf = 0.01, hwe_p = 1e-7, pedigree = ped)

Around the fit: pedigree() (topological order plus Meuwissen-Luo inbreeding), a_inverse(), a22_inverse(), solutions() (id, ebv, se and acc in one table, sorted by breeding value; a term column when a group has more than one effect per level), h2() (over the phenotypic variance of the trait, covariances included and never the rho(residual) of model_ar1(); a direct-maternal covariance enters with coefficient 1, Willham, 1972), t2() (the total heritable variance of an indirect-effect model, with delta-method standard errors), ebv(), accuracy() (prior per level: 1 + F, the diagonal of G* for a genotyped animal in a single step, K[i, i] for a kernel() term or a declared k_inverse =, the prior carried in the fit for kernel(Kinv =) (1 for an animal, F_cc / 4 for a subclass of dominance_inverse()), 1 for an iid term), residuals() and fitted() (conditional, standardized by the leverage, marginal, Pearson and the AR innovations), residual_diagnostics() (homoscedasticity by contemporary group, fitted fixed part and covariate, normality and within-animal serial correlation, each with an effect size and a parametric-bootstrap p-value, or a pivotal one for a gibbs() chain; next to them the largest standardized residuals with a Bonferroni p and a summary of the standardized breeding values), h2_curve() and plot() for the reaction norm, indirect_residual() for the pen-size residual of the associative model, var(e_i) = s2_ED + (n_i - 1)^(1 - 2d) s2_ES with d taken from the indirect() term (d = 0 is the book’s (n_i - 1) s2_ES), by profile REML over exact weighted fits (Bijma, 2010), associative_matrix(dilution = d) for the same residual as an exact covariance, describe() to look at the data before estimating, the selection-signature scans fst() (Weir and Cockerham, 1984) and roh() (the F_ROH of McQuillan et al., 2008, and islands), simulate_breeding(), a gene-dropping simulator so that examples and method studies share one honest generator, thi() and heat_load() for the heat-stress axis, selection_index() and rank_drift() for the selection side, mc_study() for repeated simulate-and-refit studies, and suggest_model(), which reads the shape of the data and names the term each shape asks for (and the trap it guards against). Timing claims go through benchmark_fit(), which replicates at least three times and checks the runs returned identical numbers: the package’s own timing rule as a tool.

The full map of the 74 functions, grouped by kinship, is in FUNCTIONS.md; the hands-on that works through 58 of them, step by step on data simulated in the document itself, is the vignette vignettes/hands-on.Rmd (every chunk runs at build time, so it cannot rot). The theory behind apy_core= (why APY works and what the Mendelian residual means) is in APY.md. The derivations behind competition_strength() and the contest estimators (exact pair conditioning, the aliasing of a uniform indirect effect, composite likelihood and its failed Bartlett identity, the Laplace variance components) are in vignettes/contest-models.Rmd.

The design decision that carries everything

The layout unit is the covariance group, not the term. group = "g" puts two random terms in the same covariance matrix with the correlation estimated, and that is why direct-maternal, the reaction norm and the associative model have no dedicated fitter: they are the same engine with different incidences and the same kron(C^-1, K^-1) penalty. 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, the ids of a declared K, or, for terms without a relationship matrix, the union of the level names of their columns.

The formula therefore departs from (1 | group) on purpose: that notation has nowhere to say that two different terms share a covariance matrix.

How this was validated

Nothing here is checked against itself. Each piece answers to an independent path:

The table: each piece and the path it answers to

what	against what
A^-1 and F	tabular A by the classic recursion; A^-1 A = I
sparse Cholesky, selected inverse	solve() and the package’s own dense path
-2logL of the sparse MME	dense V form, a path with nothing in common
analytic score	central finite differences, ALL parameters
full fit	recovery of the components used to simulate the data
A22^-1	inverse of the tabular-A block, plus the trap gate (block 22 of A^-1 != A22^-1)
single step	identity: blend 1 forces H^-1 == A^-1 through the whole fit
multi-trait	V form, finite differences on all parameters, and the collapse: both missing == row removed
AR(1) residual	V form, finite differences including rho, and the collapse: rho = 0 == identical iid path
APY	identity: the output is the exact inverse of the G that APY implies; core = everyone == exact
Vecchia	the bridge: k = 2 on a pedigree without full sibs IS Henderson’s A^-1 (1e-10); k = n-1 == exact fit
Fst and ROH	constructed references: alternate fixation gives exactly 1, a planted run is found
fixed-effect solutions	the published fixed effects of Examples 4.1, 5.1, 5.2, 8.1 and 9.1 (by contrast), plus GLS and dense-MME rebuilds in plain R at 1e-6
threshold model	Examples 15.1 and 15.2: published thresholds, solutions, standard errors, category probabilities, and the joint quantitative-binary analysis
kernel(K=), non-additive	Examples 13.1-13.5: the printed D and D^-1, solutions to 1e-3, and the MME-vs-V-form identity with a kernel term in the model
dominance_inverse(), kernel(Kinv=)	the subclass F printed in HV91 Tables 1 and 3; Example 13.1 through the sparse inverse (DV and BV of p.228 to 1e-3); Q^-1 = D to 1e-10 by both routes on random inbred pedigrees; at fixed components -2logL, score and AI equal to kernel(K = D) in model() and model_mt(), -2logL and score in model_ar1(); 50 columns of Q against Cockerham-from-Colleau columns on ~50,000 animals with pair keys past 2^31 (exactness, 4.6e-15); recovery of the base-population sigma2_D by gene dropping, 20 replicates of 6,300 animals (0.158, Monte Carlo SE 0.0075, target 0.15). A full fit at >= 20k animals was not run
multibreed partial matrices	Examples 14.1 and 14.2: the printed partial A’s and solutions, and the identity model 14.8 == variance-weighted 14.3
survival model	Example 16.1: the 23 published solutions, RRS and S(40); a censoring gate where treating censored as observed provably distorts the fit
survival, time-dependent covariates	a record split into two pieces with the same covariates changes nothing; the joint -2logL equals a Weibull likelihood written from scratch, with a late entry
categorical and survival data structure	text categories refused with the order sorting would give, and an unordered factor in alphabetical order warned; numeric text, and a factor or ordered() factor in liability order, give the numeric-code fit (thresholds and EBV to 1e-12), and the probit chain of a factor is the chain of its 0/1 column draw for draw; a 0/1 censoring factor in either level order equals the integer column in model_survival(), residuals() and survival_split(); a fixed level all in an end category, or with no failure, refused by name with its count, while a level all in a middle category and a one-record category fit; with every record censored and rho, lambda, sigma2 given, the frailties equal the penalized mode written in R (tests/testthat/test-estrutura-e5.R, 340 simulated records; refusal and identity gates, nothing at scale or on real data)
h2() and the share column	the denominators rebuilt by hand: Willham’s for direct-maternal, rho(residual) kept out in AR(1), the phenotypic variance of Bijma et al. (2007) with indirect(); 4 var(sire) in a sire model and f var(sire) / 4 of the known maternal grandsires with sire(mgs =); the residual at the harmonic mean weight; animal(id) and kernel(id, K = A), and metafounders and kernel(id, K = A_Gamma), give the same h2 on the same likelihood; NA with a warning when the SINGULAR verdict holds the numerator (tests/testthat/test-estrutura-e6.R, simulated data at given components)
structure of the data before the fit	a random term nested in a fixed effect, a single level, indirect() without mates, in one pen or with a fixed pen stop in every fitter that estimates and warn at given components; a relationship term whose records all carry the unknown code 0 stops in every fitter, in the ids rule; repeated records without pe(), few pens, few levels and reused pen labels warn, with the reused-label rates measured on correct and on merged data (tests/testthat/test-estrutura-e6.R, validation/structure_cases.R; simulated data only)
two terms in one group, several traits	the bivariate with between-trait covariances at zero equals the sum of the univariates, by component NAME, to 1.7e-10
multi-trait, a level missing for one trait	sparse and dense V routes agree to 7.5e-12 after the per-trait drop; the components land where the single-trait fits land
multi-trait AR(1) breeding values	EBV and PEV of the dense mixed-model equations built in R from the raw data
indirect effects, identifiability	full-sib families with one record per offspring, planted 1 / 0.2 / 0.3 / 1, 1600 records and 2000 animals per data set, 30 replicates per design (validation/indirect_gibbs_recovery.R; means over the replicates with their Monte Carlo SE, coverages with the binomial SE of 0.95 over 30, 0.040). Pens of 2 to 8 identify the components. REML gives 0.982 (0.033) / 0.210 (0.012) / 0.295 (0.009) / 1.029 (0.028), SINGULAR in 0 of 30, Wald coverage 1.00 / 0.90 / 0.83 / 0.93 on these 30; over 300 replicates it gives 0.998 / 0.206 / 0.300 / 0.999 with Wald coverage 0.96 / 0.93 / 0.92 / 0.93 (binomial SE 0.013, so cov, var(indirect) and var(residual) undercover by about 2 SE), and the 0.83 of var(indirect) is a draw of the first 30 seeds. gibbs() with the components estimated and prior = "flat", 2 chains of 50 000 to 200 000 iterations so that every component reached a bulk ESS of 400 (30 of 30): 1.056 (0.034) / 0.207 (0.012) / 0.306 (0.009) / 0.988 (0.029), coverage of the 95% interval 0.97 / 0.90 / 0.87 / 0.93 against the REML Wald 1.00 / 0.90 / 0.83 / 0.93 on the same replicates; across replicates the posterior mean correlates 0.999 to 1.000 with REML and sits -0.28 to +0.33 posterior SD from it, the shift a flat prior gives a skewed posterior, and the breeding values correlate with the true ones as the BLUP at the REML components does (0.622 direct, 0.675 indirect). Under the default prior = "jeffreys" a chain was held at a direct-indirect correlation of +-1 in 4 of 30 replicates (6, 9, 22, 30), with REML correlations of 0.53 to 0.71 and no SINGULAR; c99ad26 gave an empty message there, and gibbs() now warns in those 4 and in no other (gate in tests/testthat/test-indireto-identificabilidade.R). A bulk ESS of 400 was reached in 25 of 30: not in those 4, nor in replicate 19, where 12.6% of the draws of one chain had var(residual) below 0.1 (the same improper mass at zero, for a scalar variance) and the ESS stayed at 82 after 200 000 iterations. All 30 give 0.984 (0.042) / 0.220 (0.014) / 0.296 (0.010) / 1.046 (0.038), coverage 0.97 / 0.83 / 0.83 / 0.90, R-hat up to 1.38; the other 26, replicate 19 among them, give 1.033 / 0.203 / 0.305 / 0.995, coverage 1.00 / 0.92 / 0.88 / 0.96. Pens all of 4 holding 2 + 2 full sibs of two families: REML flags SINGULAR in 30 of 30 and gives the same -2logL from different starts. gibbs() is not refused there and has no SINGULAR verdict of its own: under Jeffreys its message was empty in 29 of 30, and in the 30th only the boundary warning came (1.1% of the draws of one chain), under flat empty in 30 of 30 (at c99ad26 empty in 30 of 30 under both). It recovers the two identified combinations, va + 6 cov + 9 vs 4.891 (Jeffreys) and 4.961 (flat) against 4.9 (REML 4.833) and va - 2 cov + vs 0.895 and 0.936 against 0.9 (REML 0.865), and wanders along the ridge for the separate components (R-hat above 1.01 in 30 of 30, up to 2.48 under Jeffreys and 1.89 under flat). Records on unrelated founders only (fundadores mode): REML SINGULAR on var(animal) and var(residual) in 30 of 30, and the message of the chain empty in 30 of 30
sire / maternal-grandsire pedigree	exact: the pedigree expanded with a dummy dam per animal; the A^-1 printed for Example 15.2; the mixed pedigree (dam =) against the same expansion only where the dam is missing, A to 1e-12 and the fit identical; at scale, 20 replicates of 36 000 records recover var(sire) = va / 4 within 1.6% (validation/sire_mgs_recovery.R)
APY core by eigenvalues	the count by two routes (eigenvalues of G, singular values of Z); "auto" == the same core passed by hand; the Lanczos estimate against the exact count on both sides of the Gram matrix, within 2% and 4 standard errors + 2 (20 probe seeds: at most 1.04% and 3.8 standard errors, validation/apy_core_lanczos_se.R); the quadrature error, paired with the same probes on the exact eigenvectors and averaged over v from 80 to 99.5%, under 0.1%; the standard error against the spread of 20 independent probe sets, ratio within 0.7 to 1.4 (0.82 to 1.31 over eight draws at the gated levels); a probe that stops at an invariant subspace gives the count of its probes on the exact eigenvectors, and the test fails with that branch removed; at 10 000 x 10 000 over ten seeds, mean 5957.3 against 5953 exact (z = 1.1) and a paired quadrature error of +2.3 counts, 0.04% (validation/apy_core_lanczos.R)
h_inverse()	the single-step formula rebuilt in R, exact and with APY; kernel(K = H) == genotypes= at the same theta
threshold and survival, estimated components	the estimate is the minimum of the Laplace -2logL the fit reports: no point of a grid around it in any optimizer coordinate is lower, a second start and an independent optimize() land on it; a variance the data do not tell apart from zero is held at its minimum without an SE while the other SEs stay exact; 80 sires with 50 daughters each, binary, planted 0.15: mean 0.161 (SE 0.007) over 10 replicates (validation/indirect_threshold_survival_recovery.R touro)
indirect() in model_threshold() and model_survival()	exactness on 105 to 296 records: the incidence built in R == the engine’s (Gaussian BLUP against model()); mode, system, PEV and Laplace with a direct-indirect group and random(pen) against dense references written from the definition with Z from a loop of their own (1e-8 or better); the correlation near +-1 profiled. Recovery at prototype scale only (2329 records, 8 replicates), and the declared gates FAIL: the threshold components run LOW, more with the indirect term than without (paired s2D -0.031, SE 0.007); the Weibull components run low with one record per animal, equally with and without it. Nothing at >= 20k records; a binary trait with the indirect components estimated is the next row
indirect effects, binary trait, components estimated	the families and pens of 2 to 8 of the identifiability row, on a probit liability (planted 1 / 0.2 / 0.3 and the residual 1, incidence 0.309 to 0.388), 3200 records and 4000 animals per data set, 30 replicates (validation/indirect_gibbs_probit_recovery.R). The design identifies the components: REML on the liability itself, which real data do not have, gives 1.011 / 0.203 / 0.300 with Wald coverage 0.933 / 1.000 / 0.933 (binomial SE 0.0455). From the binary records neither estimator recovers them. The Laplace minimum of model_threshold() averages 0.204 / 0.173 / 0.150 over the 30 (bias -0.796, -0.027 and -0.150; SE 0.010, 0.005 and 0.004) at |r| > 0.99 in 29 of 30, and its Wald interval covers var(animal) and var(indirect) in 0 of 30. gibbs(family = "probit") with the default prior and 2 chains of 6000 reaches a bulk ESS of 400, or an R-hat under 1.01, in 0 of 30 replicates (ESS 2.5 to 23.3 on the variances, R-hat up to 2.341), and its mean still rises from the first half of the chain to the second (var(animal) +0.346, SE 0.119): the average over replicates, 1.138 / 0.213 / 0.340, is where the chains had got to, not a posterior mean. On one data set with 2 chains of 40 000 the ESS of var(animal) is 3.4; under Jeffreys one chain stays near r = 1 at the Laplace point for its first six tenths (var(animal) 0.168 to 0.265, r 0.974 to 0.996), and under the flat prior var(animal) climbs to 33. mcmc_diagnostics() flags such chains: on 20 seeds of a test-sized design the verdict is ok in 0 of 20, with the indirect effect planted and without it
Gibbs, probit and kernel()	with the components held, the posterior mean tracks the threshold-model mode (probit) and the model() BLUP (kernel); K = I == random(id), the same chain
Gibbs, the whole chain	with flat priors the posterior of the components is the REML likelihood (Harville, 1974): four chains match a 2-D quadrature of the -2logL at z = 0.15 and 0.08 (validation/gibbs_harville.R); mcmc_diagnostics() equals posterior 1.7.0 to 1e-10 (R-hat, bulk/tail ESS, MCSE), and its verdict is calibrated on simulated chains: never not mixed on stationary ones; not mixed 83% and 100% of the time with one of 4 chains shifted 0.5 and 1 sd; burn-in too short every time for a 50 sd start with 4 chains (validation/mcmc_diagnostics_calibration.R). At 20 000 levels and 4 chains the diagnostics take 0.22 to 0.23 s (3 runs; validation/mcmc_diagnostics_scale.R, iid levels); serial and parallel chains identical
threads	the tiled tail, the dense inverses, the products and the selected inverse against chol(), solve() and the formulas; the same bits with 1 and 4 threads, and the whole suite passes at both
metafounders with genotypes	H(Gamma) rebuilt from the dense A(Gamma) and G05 at 1e-8; kernel(K = H) == genotypes=; snp_blup() == model() at the same theta; 20 simulated replicates of two bases reduce the trend bias from -0.141 to -0.081 (validation/hgamma_recovery.R)
estimate_gamma()	one pseudo-EM step is the metafounder block of H(Gamma) built densely; GLS shows its published bias; "ml" against the identity A(gamma) = (1 - gamma/2) A + gamma 11’ and the maximum of the dense likelihood; the engine version against the previous R version: pseudo-EM and GLS within 1e-10 with the same iteration count, with genotypes in the last generation only, with genotyped parents and with every animal genotyped; ML summaries within 1e-12 and gamma within 1e-6 with A22^-1 dense, sparse and with no non-genotyped animal; A22^-1 X with 70 columns against the dense A22; a metafounder with no genotyped descendant gives the Gamma of the problem without it
snp_effects()	G*^-1 u read back from the exported H^-1, ((H^-1 - A^-1) u)_g + A22^-1 u_g, with ungenotyped ancestors, on each route (exact, APY, Vecchia) with and without metafounders, and, on the exact route, on a sire / maternal-grandsire pedigree; APY also against the APY inverse built densely in R; with metafounders, the marker effects of snp_blup() at the same theta; the previous R version within 1e-9; 1 vs 4 threads identical; a different blend, other genotypes and reordered rows under Vecchia are refused
A22 by Colleau, the sparse Schur	H^-1 against the formula with A22 from the dense A, with inbreeding and with a sire-MGS pedigree; a22_inverse() against solve(); the APY route with core = everyone gives the exact H^-1 to 1e-12; with two metafounders and a full Gamma, core = everyone gives the exact H(Gamma) to 1e-10 and a partial core the APY inverse written in R on the dense G* to 1e-8
pegs()	at fixed variances the exact multivariate ridge, a dense mk x mk solve; the structures are identities where they must be; recovery with 2000 animals x 2000 markers x 3 traits; the Julia reference gives the same numbers on the same data
sire(sire, mgs =)	-2logL and BLUP against the dense GLS with the incidence built by hand
indirect(dilution =) in the siblings	the bivariate with no between-trait covariance == the sum of the univariate model() fits with the same d; AR(1) at rho = 0 == model(); the Gibbs chain with fixed components == the diluted BLUP
snp_blup(), any structure	the dense single-step solve built in R with H^-1 from G* without the affine step: breeding values and marker effects to 1e-6 of their SD with one animal term, direct-maternal in one group, direct and maternal in separate groups, a reaction norm and direct-indirect, the latter also with dilution = 0.7 in pens of 1 to 7, with mates without a phenotype and a repeated record; rpg near 1 falls back to the pedigree BLUP; a planted QTL comes out on top
genotype storage	double, integer and raw matrices give the same G, H^-1 (exact and APY), APY core (both routes), genomic F and D, Gamma (pseudo-EM, GLS, ML), fit, backsolved SNP effects (exact and APY) and ssSNPBLUP, bit for bit, with missing values in the matrix; read_blupf90_snp() returns the matrix that was written, subsets by id; F, D and G are also checked past one band of 1024 animals and two blocks of 2048 markers, D against the previous version and G against R’s tcrossprod, through the tiles and dsyrk
survival_split()	the subject and change tables reproduce hand-built elementary records exactly, with the same fit; S(t \| e) S(e) = S(t)
residuals, leverage, residual_diagnostics()	the dense R P y, GLS fixed part and diag(R P R) in plain R for model() (weights, pe, rn, mgs, indirect, kernel, single step), model_ar1() (rho of either sign) and model_mt(); the closed-form Cook distance against the refit without the record; Var(u^) = s2 H - PEV for genotyped and non-genotyped animals; type I of the bootstrap p measured in five simulated layouts of 600 to 2400 records, and of the pivotal p of a gibbs() chain in 640 records (the probit with the contemporary group random) (validation/residual_diagnostics_calibration.R). Not validated on real data
fixed-effect design (sparse)	the sparse rank keeps the columns of the dense Gram-Schmidt it replaced on 150 random designs (crossed, nested and disconnected classes, levels with no used record, collinear and nested covariates, weights), on edge columns (tiny, NaN, Inf, nearly collinear), and on unscaled polynomials and near copies. In near-singular designs every decision is certified or the rank falls back to that Gram-Schmidt, and it never keeps more columns than records. dropped_x of model() (with and without weights), model_mt() (per trait) and model_ar1() equals a dense reference rebuilt in R. The threshold and survival fitters drop what qr() drops and fit the reduced model. The bits are the same with 1 and 4 threads. With the rest of the engine unchanged, a check found the results bit-identical to the dense build over the 3 085 values returned by the fits of every test file but this change’s own (validation/fixed_design_identity.R). Scale (memory, not exactness): 200 000 records x 701 candidate columns, fixed part 145-153 MB of peak working set (232 MB committed), three sets of 3 replicates, with dropped_x checked against the table (validation/fixed_design_memory.R)
missing labels and values	the rule has two implementations, the engine’s and R’s. One battery of values (numeric NA, NaN, Inf, -Inf; blank, tab and line-break text, ” NaN “) gives the same types, the same count and the same first row through both, in every fitter. A missing trait drops only its record, with the -2logL of the data without that row
unit, shift and boundaries in model_mt() and model_ar1()	y, 100 y and 1e4 y give the same SINGULAR verdict and the same standard errors up to the factor in all three fitters; one trait in other units gives the same fit up to the scale factors and the Jacobian; a constant added to y leaves the REML -2logL alone; a covariance boundary that is not the optimum is not certified, and a certified point on a flat face is polished to the face optimum (tests/testthat/test-unidade-partida-fronteira.R). At 21 000 records the standard errors of a model_ar1() fit agree to 8.2e-7 between y and 100 y, and at 21 000 and 24 000 records the -2logL agrees to 6e-15 (relative) under a shift of up to 1000 standard deviations (validation/unit_shift_start_scale.R). In simulated sweeps of 368 fits from four starts per data set, run during development and not kept in validation/, no certified point sits more than 9.4e-5 above the best of its starts, and one known pair of starts still certifies points 1.7e-4 apart. The real fits that motivated the change were not rerun
ids of the data against the pedigree (founders, unknown 0, pruning, spellings, 2^53, blocks)	tests/testthat/test-estrutura-e1.R (34 gates; 33 fail on 12aca18, the other is a negative control); validation/pedigree_ids_scale.R: pruning and founders exact (difference 0 in -2logL, EBV and accuracy) at 20 200 records with a 125 200-line pedigree; the rule alone takes 0.03 s, and the fit with the national pedigree takes the time of the fit pruned by hand (median of 3: 11.97 s against 10.90 s, replicates 11.1-14.5 against 10.8-12.4 s on a shared machine); simulated only
AR(1) time axis (grid step, CAR(1), the axis check)	tests/testthat/test-estrutura-e3.R (32 gates; all fail on 12aca18) and test-ar1-unidade-tempo.R; validation/ar1_time_axis_cases.R: every case NEWS quotes, on both builds; validation/ar1_time_axis_scale.R: a time contemporary group of 715 to 2858 levels keeps dense_block at 8 (5 subjects x 5 000 and x 20 000 days: 0.7 s and 9.7 s, median of 3); 30 x 1 000 days with a daily group 1.0 s (bivariate 19.1 s); rho = 0 equals the iid model() to 1e-8 in -2logL; 4 000 x 6 days give the same fit in days and POSIXct; with no serial correlation, 0 of 20 data sets of 1 000 x 6 with a half day and 0 of 200 of 60 x 5 in continuous time stop (4 warn); with rho -0.4, 192 of 200 stop and 7 warn; simulated only

The tests in tests/testthat run these comparisons on every build, so a change that breaks one of the identities cannot pass quietly. simulate_breeding() is what they are built on: it generates the pedigree, the phenotypes and the genotypes together, so every check has the truth beside it.

References

The works cited on this page

Abdollahi-Arpanahi, R., Lourenco, D. & Misztal, I. (2022). A comprehensive study on size and definition of the core group in the proven and young algorithm for single-step GBLUP. Genetics Selection Evolution 54:34.

Aguilar, I., Misztal, I., Johnson, D.L., Legarra, A., Tsuruta, S. & Lawlor, T.J. (2010). A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score. Journal of Dairy Science 93:743-752.

Albert, J.H. & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association 88:669-679.

Amestoy, P.R., Davis, T.A. & Duff, I.S. (1996). An approximate minimum degree ordering algorithm. SIAM Journal on Matrix Analysis and Applications 17:886-905.

Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. & Sorensen, D. (1999). LAPACK Users’ Guide, 3rd ed. SIAM, Philadelphia.

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

Bijma, P., Muir, W.M. & Van Arendonk, J.A.M. (2007). Multilevel selection 1: quantitative genetics of inheritance and response to selection. Genetics 175:277-288.

Bradley, R.A. & Terry, M.E. (1952). Rank analysis of incomplete block designs: I. The method of paired comparisons. Biometrika 39:324-345.

Cantet, R.J.C. & Cappa, E.P. (2008). On identifiability of (co)variance components in animal models with competition effects. Journal of Animal Breeding and Genetics 125:371-381.

Christensen, O.F. & Lund, M.S. (2010). Genomic prediction when some animals are not genotyped. Genetics Selection Evolution 42:2.

Cockerham, C.C. (1954). An extension of the concept of partitioning hereditary variance for analysis of covariances among relatives when epistasis is present. Genetics 39:859-882.

Dempster, A.P., Laird, N.M. & Rubin, D.B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B 39:1-38.

Dempster, E.R. & Lerner, I.M. (1950). Heritability of threshold characters. Genetics 35:212-236.

Ducrocq, V. (1997). Survival analysis, a statistical tool for longevity data. 48th Annual Meeting of the European Association for Animal Production, Vienna.

Ford, L.R., Jr. (1957). Solution of a ranking problem from binary comparisons. American Mathematical Monthly 64(8, part 2):28-33.

Foulley, J.L., Gianola, D. & Thompson, R. (1983). Prediction of genetic merit from data on binary and quantitative variates with an application to calving difficulty, birth weight and pelvic opening. Genetics Selection Evolution 15:401-424.

Foulley, J.L., Im, S., Gianola, D. & Hoeschele, I. (1987). Empirical Bayes estimation of parameters for n polygenic binary traits. Genetics Selection Evolution 19:197-224.

Fragomeni, B.O., Lourenco, D.A.L., Tsuruta, S., Masuda, Y., Aguilar, I., Legarra, A., Lawlor, T.J. & Misztal, I. (2015). Use of genomic recursions in single-step genomic best linear unbiased predictor with a large number of genotypes. Journal of Dairy Science 98:4090-4094.

Garcia-Cortes, L.A. & Toro, M.A. (2006). Multibreed analysis by splitting the breeding values. Genetics Selection Evolution 38:601-615.

Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis 1:515-534.

Geman, S. & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence 6:721-741.

George, A. & Liu, J.W.H. (1989). The evolution of the minimum degree ordering algorithm. SIAM Review 31:1-19.

Geweke, J. (1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In Bernardo, J.M., Berger, J.O., Dawid, A.P. & Smith, A.F.M. (eds), Bayesian Statistics 4. Oxford University Press, Oxford.

Geyer, C.J. (1992). Practical Markov chain Monte Carlo. Statistical Science 7:473-483.

Gianola, D. & Foulley, J.L. (1983). Sire evaluation for ordered categorical data with a threshold model. Genetics Selection Evolution 15:201-224.

Gilmour, A.R., Thompson, R. & Cullis, B.R. (1995). Average information REML: an efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51:1440-1450.

Griffing, B. (1967). Selection in reference to biological groups. I. Individual and group selection applied to populations of unordered groups. Australian Journal of Biological Sciences 20:127-140.

Henderson, C.R. (1950). Estimation of genetic parameters (abstract). Annals of Mathematical Statistics 21:309-310.

Henderson, C.R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics 31:423-447.

Henderson, C.R. (1976). A simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values. Biometrics 32:69-83.

Hobert, J.P. & Casella, G. (1996). The effect of improper priors on Gibbs sampling in hierarchical linear mixed models. Journal of the American Statistical Association 91:1461-1473.

Hoeschele, I. & VanRaden, P.M. (1991). Rapid inversion of dominance relationship matrices for noninbred populations by including sire by dam subclass effects. Journal of Dairy Science 74:557-569.

Hunter, D.R. (2004). MM algorithms for generalized Bradley-Terry models. Annals of Statistics 32:384-406.

Kachman, S.D. (1999). Applications in survival analysis. Journal of Animal Science 77(suppl. 2):147-153.

Kirkpatrick, M., Lofsvold, D. & Bulmer, M. (1990). Analysis of the inheritance, selection and evolution of growth trajectories. Genetics 124:979-993.

Legarra, A., Christensen, O.F., Vitezica, Z.G., Aguilar, I. & Misztal, I. (2015). Ancestral relationships using metafounders: finite ancestral populations and across population relationships. Genetics 200:455-468.

Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares. Quarterly of Applied Mathematics 2:164-168.

Liu, Z., Goddard, M.E., Reinhardt, F. & Reents, R. (2014). A single-step genomic model with direct estimation of marker effects. Journal of Dairy Science 97:5833-5850.

Luce, R.D. (1959). Individual Choice Behavior: A Theoretical Analysis. Wiley, New York.

Marquardt, D.W. (1963). An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics 11:431-441.

Masuda, Y., Misztal, I., Legarra, A., Tsuruta, S., Lourenco, D.A.L., Fragomeni, B.O. & Aguilar, I. (2017). Technical note: avoiding the direct inversion of the numerator relationship matrix for genotyped animals in single-step genomic best linear unbiased prediction solved with the preconditioned conjugate gradient. Journal of Animal Science 95:49-52.

McQuillan, R., Leutenegger, A.-L., Abdel-Rahman, R., Franklin, C.S., Pericic, M., Barac-Lauc, L. et al. (2008). Runs of homozygosity in European populations. American Journal of Human Genetics 83:359-372.

Meuwissen, T.H.E. & Luo, Z. (1992). Computing inbreeding coefficients in large populations. Genetics Selection Evolution 24:305-313.

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.

Mrode, R.A. & Pocrnic, I. (2023). Linear Models for the Prediction of the Genetic Merit of Animals, 4th ed. CABI, Wallingford. doi:10.1079/9781800620506.0000.

Muir, W.M. (2005). Incorporation of competitive effects in forest tree or animal breeding programs. Genetics 170:1247-1259.

Muir, W.M. & Schinckel, A.P. (2002). Incorporation of competitive effects in breeding programs to improve productivity and animal well being. Proceedings of the 7th World Congress on Genetics Applied to Livestock Production, Montpellier.

National Research Council (1971). A Guide to Environmental Research on Animals. National Academy of Sciences, Washington DC.

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

Plackett, R.L. (1975). The analysis of permutations. Applied Statistics 24:193-202.

Pocrnic, I., Lourenco, D.A.L., Masuda, Y., Legarra, A. & Misztal, I. (2016). The dimensionality of genomic information and its effect on genomic prediction. Genetics 203:573-581.

Quaas, R.L. (1976). Computing the diagonal elements and inverse of a large numerator relationship matrix. Biometrics 32:949-953.

Schaeffer, L.R. & Dekkers, J.C.M. (1994). Random regressions in animal models for test-day production in dairy cattle. Proceedings of the 5th World Congress on Genetics Applied to Livestock Production, Guelph, 18:443-446.

Schafer, F., Katzfuss, M. & Owhadi, H. (2021). Sparse Cholesky factorization by Kullback-Leibler minimization. SIAM Journal on Scientific Computing 43:A2019-A2046.

Sorensen, D.A., Andersen, S., Gianola, D. & Korsgaard, I. (1995). Bayesian inference in threshold models using Gibbs sampling. Genetics Selection Evolution 27:229-249.

Takahashi, K., Fagan, J. & Chin, M.-S. (1973). Formation of a sparse bus impedance matrix and its application to short circuit study. Proceedings of the 8th PICA Conference, 63-69.

Tempelman, R.J. (1998). Generalized linear mixed models in dairy cattle breeding. Journal of Dairy Science 81:1428-1444.

Tukey, J.W. (1977). Exploratory Data Analysis. Addison-Wesley, Reading, MA.

Vandenplas, J., Calus, M.P.L., Eding, H. & Vuik, C. (2019). A second-level diagonal preconditioner for single-step SNPBLUP. Genetics Selection Evolution 51:30.

Vandenplas, J., Eding, H., Calus, M.P.L. & Vuik, C. (2018). Deflated preconditioned conjugate gradient method for solving single-step BLUP models efficiently. Genetics Selection Evolution 50:51.

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

Varadhan, R. & Roland, C. (2008). Simple and globally convergent methods for accelerating the convergence of any EM algorithm. Scandinavian Journal of Statistics 35:335-353.

Vecchia, A.V. (1988). Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society, Series B 50:297-312.

Vitezica, Z.G., Varona, L. & Legarra, A. (2013). On the additive and dominant variance and covariance of individuals within the genomic selection scope. Genetics 195:1223-1230.

Wade, K.M. & Quaas, R.L. (1993). Solutions to a system of equations involving a first-order autoregressive process. Journal of Dairy Science 76:3026-3032.

Wang, C.S., Rutledge, J.J. & Gianola, D. (1993). Marginal inferences about variance components in a mixed linear model using Gibbs sampling. Genetics Selection Evolution 25:41-62.

Wang, C.S., Rutledge, J.J. & Gianola, D. (1994). Bayesian analysis of mixed linear models via Gibbs sampling with an application to litter size in Iberian pigs. Genetics Selection Evolution 26:91-115.

Wang, H., Misztal, I., Aguilar, I., Legarra, A. & Muir, W.M. (2012). Genome-wide association mapping including phenotypes from relatives without genotypes. Genetics Research 94:73-83.

Weir, B.S. & Cockerham, C.C. (1984). Estimating F-statistics for the analysis of population structure. Evolution 38:1358-1370.

Willham, R.L. (1972). The role of maternal effects in animal breeding: III. Biometrical aspects of maternal effects in animals. Journal of Animal Science 35:1288-1293.

Wright, S. (1922). Coefficients of inbreeding and relationship. American Naturalist 56:330-338.

If a work that should be cited here is missing, please open an issue or write to the maintainer address in DESCRIPTION.

Choices worth knowing about

Convergence takes TWO things, not one. The step criterion is the RELATIVE change in the components, sqrt(sum(dtheta^2) / sum(theta^2)) < tol, with a default of 1e-8: never an absolute threshold on the score, which grows with the number of records. Coming from the BLUPF90 family, mind the scale: those programs test that quantity squared, so a card’s conv_crit is this tol squared. A 1e-12 there is tol = 1e-6 here, and the 1e-8 default here would be 1e-16 on that scale.

A small step is not an optimum, though. A damped step can accept a tiny move with the gradient still far from zero, so converged also requires the NEWTON DECREMENT, g' AI^-1 g, which is about twice the remaining gap in -2logL, to fall under 2e-4. It is reported as newton_dec, next to score, so the certificate can be read rather than trusted. Components resting at a covariance boundary are excluded from it, because a component pinned there points out of the cone by construction and its gradient never vanishes; when that happens the message says how many were excluded, and the certification is conditional on that pinning.

In model() a component at its floor is excluded only when the score pushes it out of the cone. In model_mt() and model_ar1() the step freezes a boundary block whatever the sign, so the certificate is checked on the likelihood: points inside the cone along the negative eigenvectors of the block’s gradient, steps along the boundary face by finite differences, and a second walk from inside the cone. A point that lowers the -2logL by more than 1e-4 sends the walk there, or refuses the certificate when the walk cannot move. Without those checks a real bivariate fit had been certified 46.8 units above its interior optimum. A boundary face is flat, so a certified point is then polished along it while each step gains more than 1e-5 and the point still certifies. In simulated sweeps of 368 fits from four starts per data set (run during development, not kept in validation/), every certified point then came within 9.4e-5 of the best of its starts, and in a separate search one pair of starts still certified points 1.7e-4 apart.

Two endings get their own verdict in model(), model_mt() and model_ar1(). When maxiter runs out with the Newton decrement already under tolerance (and, in model_mt() and model_ar1(), the boundary checks passed), the likelihood is flat along some direction: the fit returns converged = TRUE and says so, because more iterations would not choose a point on that ridge. When the average information is singular at the end, the message says SINGULAR and names the components the data do not separate; their standard errors are NaN, the reported point depends on start=, and only combinations of them are estimable. That is the design speaking, the identifiability question Cantet and Cappa (2008) raise for the indirect-effect model, and the answer is more information, a simpler model, or a component held fixed. In model_mt() and model_ar1(), a stop where no damped step and no EM step improves the likelihood is never a certificate: the message gives the iteration and the decrement. maxiter = 0 evaluates the start and reports converged = FALSE in all three fitters.

model() starts from var(y), so a run does not inherit a neighbour’s answer by accident. model_mt(), and model_ar1() with cbind(), start from the univariate optimum of each trait with every covariance between traits at zero, the start one builds by hand in BLUPF90 or ASReml. On real data a bivariate CAR(1) from var(y) had walked to r = -1 and r = +1 and stopped 41.3 -2logL units above the optimum the univariate start reached in 6 iterations. That was measured before the walk was standardized per trait and before the boundary checks above, and with the current engine no simulated case is known where the var(y) start fails. A variance that a univariate fit left at zero starts at 1% of the var(y) start, and the univariate fits run to their own optimum whatever maxiter is. When the fit from that start stalls without a certificate, the var(y) start is fitted too, and start_from says which start the returned fit used. The default rests on that one real case, and the univariate fits cost time: at 21 000 and 24 000 records, 1.1 to 1.8 times the var(y) start on the median of 3 runs when both reach the same optimum (validation/unit_shift_start_scale.R). The var(y) start, still used with one trait, with maxiter = 0 and as the fallback, carries two adjustments that keep it meaning the same thing when the model is written differently. The share of a term whose covariance is user-supplied is divided by the geometric mean of that matrix’s eigenvalues, so the same model with K and with c * K starts at equivalent points; without it, c = 25 used to land 146 units of -2logL away from c = 1, reporting convergence in both. And the AR(1) rho starts from the correlation at the typical spacing, 0.3^(1/dt), so the same series written in days or in weeks starts at the same place; rho = 0 is a stationary point of that parameterisation whenever no two times differ by exactly 1, which used to return rho = 0 with converged = TRUE. model() carries neither adjustment and does not need them: measured on the same cell, it reaches the same optimum from either scaling of K.

The multi-trait and AR(1) fitters walk on each trait divided by the power of two just above its standard deviation, and report in the original units. Neither the walk, nor the singularity verdict, nor the standard errors depend on the unit a trait is recorded in. y1 and 0.01 * y1 give the same estimates up to the scale factors and the same -2logL up to the Jacobian, and a power-of-two factor gives the same fit bit for bit. Components resting at a boundary are the exception: when the data round differently they can move along the flat face they rest on (up to 0.6% measured, with the -2logL equal within 2.5e-5 after the Jacobian). model() gets the same invariance from a variance floor proportional to the variance of y. The -2logL is formed from the residuals, so a constant added to y leaves it unchanged to 3e-14 (relative).

When you do want a different start, start= takes one in model(), model_mt() and model_ar1() alike: to warm-start from a submodel, or to check that the optimum does not depend on where the search began, which is the only direct evidence of a global optimum a non-convex likelihood offers.

A covariance that stops being positive-definite ends the fit instead of being nudged back into range, and a fit that did not converge says so in the print, in the message and in the object.

The data rules are equally deliberate. An unknown genotype code becomes NA and is imputed by the marker mean; it never becomes the zero genotype, which is a real observation. An animal of the data missing from the pedigree, a pen mate or a maternal grandsire included, enters as a founder with a warning that counts and names it, as renumf90 numbers every animal with a record; dropping a pen mate would quietly change who competed with whom. An id that matches a pedigree line written another way (a trailing space, the case, leading zeros) stops the fit with both spellings. A parent cited without a line of its own is an error rather than a new founder.

Questions and problems

Questions and suggestions are welcome at felipeoliveirafreitas@usp.br. If something looks like a defect, opening an issue keeps the answer where the next person who hits it will find it:

https://github.com/phyllype/BreedingR/issues

CONTRIBUTING.md says what makes a report easy to act on, and what a pull request is checked against.

License

GPL-3. Copyright (C) 2026 Felipe Andre Oliveira Freitas.

This program is free software: you can redistribute it and modify it under the terms of the GNU General Public License, version 3, as published by the Free Software Foundation. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY, without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. The full text is in LICENSE.md and at https://www.gnu.org/licenses/.

What that means day to day: use it, read it, change it and pass it on. What it asks in return is that anything built on it and distributed carries the same licence and ships its source. Running it for research and publishing the results is exactly what it is for, and the citation below is the only thing asked there.

Citation and funding

If this package contributed to published work, please cite it:

Freitas, F. A. O. (2026). BreedingR: genetic parameters, breeding values and genomic prediction in R. R package version 0.4.0.

Developed during doctoral research at the Escola Superior de Agricultura Luiz de Queiroz (ESALQ), Universidade de São Paulo, and at Purdue University, supported by

• the São Paulo Research Foundation (FAPESP), grants #2024/15502-6 and #2025/02949-5 (BEPE);
• the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Finance Code 001;
• the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).

The opinions, hypotheses and conclusions expressed here are the author’s own and do not necessarily reflect the views of the funding agencies.

Work that uses this package should carry the same acknowledgement, as the funding terms ask.