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.
remotes::install_github("phyllype/BreedingR") # requires Rtools on Windows
# or, from a clone:
install.packages(".", repos = NULL, type = "source")
library(BreedingR)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.
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 by Laplace and EM (Foulley, Im,
Gianola and Hoeschele, 1987); 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 = ), 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.
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=) |
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() |
| G, dominance, epistasis to any order | g_matrix(), g_dominance(), g_epistasis_ad/dd/order() |
| 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(), qc_phenotypes(), read_plink(), describe(), thi(), heat_load(),
fst(), roh(), simulate_breeding(), mc_study(), suggest_model(),
benchmark_fit(), 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 54 of the 64 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.
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$m, min_maf = 0.01) # call rate, MAF, HWE, with counts
fit <- model(y ~ cg + animal(id), q$data, s$pedigree,
genotypes = list(ids = s$genotypes$ids, m = g$m))
fit # components, SEs, share of the phenotypic variance
h2(fit) # the same ratio as the var(animal) share
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()
cor(ebv(fit)[names(s$tbv)], s$tbv) # against the simulator's own truthA 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.
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 markers by call rate, minor allele frequency and
Hardy-Weinberg equilibrium, and reports how many each filter removed. describe() shows
the data before a model touches it, and suggest_model() reads its shape and names the
terms it calls for.
# 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")
# 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)
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 by Laplace and the EM
# step of Foulley et al. (1987), which runs low for a binary trait with few records
# per level (Tempelman, 1998). 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)
# 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 the single-step fit
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")$m, min_maf = 0.01, hwe_p = 1e-7)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 =, 1 for an iid term), 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 64 functions, grouped by kinship, is in
FUNCTIONS.md; the hands-on that works through 54 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 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.
Nothing here is checked against itself. Each piece answers to an independent path:
| 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 |
| 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 |
| 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() |
| 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 | 30 replicates per design: pens of 2 to 8 recover the components (REML and Gibbs); pens of one size from two full-sib families flag SINGULAR and give the same -2logL from different starts |
| 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, estimated components | the EM fixed point against the minimum of the Laplace -2logL found without the EM step; 80 sires with 50 daughters each, binary, planted 0.15: mean 0.148 over 10 replicates |
| 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); rhat() near 1 on iid chains, above 1.01 on shifted or rescaled ones; 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 |
| 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, fit and ssSNPBLUP, bit for bit, with missing values in the matrix; read_blupf90_snp() returns the matrix that was written, subsets by id |
survival_split() |
the subject and change tables reproduce hand-built elementary records exactly, with the same fit; S(t | e) S(e) = S(t) |
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.
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.
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.
Two endings get their own verdict in model(), model_mt() and model_ar1(). When
maxiter runs out with the Newton decrement already under tolerance, 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.
Fits start from var(y) by default, so a run does not inherit a neighbour's answer by
accident. Two adjustments in model_mt() and model_ar1() keep that default 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.
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. A pen mate missing from the pedigree is an error rather than a silent discard, because dropping him would quietly change who competed with whom. A parent cited without a line of its own is an error rather than a new founder.
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.
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.
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.