The trait is a set of m ordered categories (binary is the m = 2 case). The model
lives on an unobserved normal liability with m - 1 thresholds: the probability of
category k is F(t_k - a) - F(t_(k-1) - a) with a = x'b + z'u (probit link,
Eqn 15.5 of Mrode & Pocrnic). The system of Gianola & Foulley (1983) is solved by
Fisher scoring, with the same A^-1 penalty on the genetic term the Gaussian
fitters use.
Usage
model_threshold(
formula,
data,
pedigree = NULL,
start = NULL,
k_inverse = NULL,
missing_code = NULL,
thresholds_start = NULL,
maxiter = 50L,
tol = 1e-08,
verbose = interactive(),
estimate = FALSE,
maxiter_em = NULL,
tol_em = NULL,
genotypes = NULL,
blend = 0.05,
apy_core = NULL,
vecchia_k = NULL,
max_evals = 500L,
tol_estimate = 1e-06,
profile = TRUE
)Arguments
- formula
as in
model(): fixed class effects,cov()covariates, and random terms amonganimal(),sire(),random()andindirect(), withgroup =to put several terms in one covariance matrix (direct and indirect, typically). The terms of a group share one set of levels: the relationship levels foranimal()/sire()/indirect(), the union of the column labels forrandom().rn(),kernel(),nested =andsire(mgs =)are not available here.cbind(quant, bin)on the left-hand side switches to the joint analysis, which takes one relationship term and neitherindirect()norgroup =. A missing level (NA,NaN,Inf,-Inf, blank text or the text "NaN") or a missing covariate value is refused on every row, with the same rule and message as inmodel(); a factor or text column incov(), or as the quantitative trait of the joint mode, is refused rather than read as its codes. A covariate enters standardized, as inmodel(): the thresholds take its centre in the ordinal mode, and the levels of the first class (or the intercept) in the joint mode;b,se_b, the thresholds and their standard errors come back in the unit of the data. The missing-value rule and the warnings about values that look like a code are those ofmodel().- data
data.frame with the columns referenced. The categories of the trait run from the lowest liability to the highest: numeric codes (or logical) by value; a factor by its levels (an
ordered()factor, or a factor whose levels were set in the liability order, as MASS::polr and ordinal::clm read a factor response); text that reads as numbers by value. Text that is not a number is REFUSED, and the message shows the order sorting would give: text carries no order, and sorted "low", "mid", "high" run high < low < mid, which inverts the scale. A factor whose labels are numbers out of numeric order is refused (factor()of a numeric text column gives "1" < "10" < "2"; a deliberate order such aslevels = 3:1must be declared withordered()), and so is anordered()factor whose levels are that text order, whatordered()gives withoutlevels =; labels that write one number twice ("1" and "1.0") are refused too. A factor that is not ordered and whose levels are in alphabetical order, whatfactor()andread.csv(stringsAsFactors = TRUE)give text whenlevels =is not set, runs with a WARNING that shows the order used and the count in each category (ordered()declares it and silences the warning; brms refuses every unordered factor for an ordinal family, MASS::polr and ordinal::clm take the levels silently). Anordered()factor is taken as the declared order, so give itlevels =: without them it sorts text alphabetically too, and only numeric labels in text order are caught. When the order does not come from numbers the fit message says which order was used, and for a binary trait which category is the upper one. A level of a fixed class effect with every record in the lowest or in the highest category (the extreme category problem; Misztal, Gianola and Foulley 1989) has no finite effect, and the fit is refused naming the level and its count: merge the level, drop its records, or fit the effect as random. In the joint mode the second trait must have exactly two values, ordered by the same rule (the upper one is the "difficulty"), and the same check runs on it. A missing observation (NA,NaN,Inf,-Inf, blank text, the text "NaN" ormissing_code) drops its record and is never a category, in either trait of the joint mode.- pedigree
data.frame animal, sire, dam; required with
animal(),sire()orindirect()unlessk_inverseis given. The ids are matched to the data by the rule ofmodel()(seepedigreethere): an absent animal enters as a founder, an unknown 0 is no effect, the same id written two ways stops the fit, and the pedigree is pruned to the data, the genotyped and their ancestors (withk_inverse =the levels are the ids of that K, an id outside it stops the fit, and the unknown 0 is no effect, counted by the same message). The joint mode refuses an unknown 0 in its relationship term on the records it fits, those with both traits observed (a 0 on a record that misses a trait leaves with that record), and counts the founders' records and the blocks on those records too- start
MANDATORY. Ordinal mode: the components on the liability scale (the residual is fixed at 1 and is not in
start), in the order ofmodel(): the covariance groups declared withgroup =in order of first appearance, then each term without a group in formula order, and inside a group the lower triangle by columns, asfit$thetanames them (var(animal),cov(indirect,animal),var(indirect)). Without groups that is one variance per random term in formula order. Joint mode:list(G = 2x2 genetic matrix, R = 2x2 residual matrix).- k_inverse
the inverse of the relationship matrix for the relationship terms, replacing the
A^-1built frompedigree(every relationship term uses it, so the terms of a group share it): either tripletslist(i, j, x, n, id)asa_inverse()returns, or a dense symmetric matrix with the ids asdimnames. This is the door for a relationship the pedigree path cannot build – the book's sire / maternal- grandsire matrix of Example 15.2 (p.278) enters here.- missing_code
missing-value code for the trait(s); records missing the trait (or either trait, in the joint mode) are dropped and counted. A record dropped for a missing trait still counts as a pen mate of
indirect().- thresholds_start
starting values for the m - 1 thresholds; by default the normal quantiles of the cumulative category proportions, as the book does (p.267)
- maxiter
maximum number of Fisher scoring iterations. With
estimate = TRUEit bounds each evaluation of the -2logL, and an evaluation that does not converge within it counts as an inadmissible point (the message counts them)- tol
RELATIVE tolerance on the full solution vector,
sqrt(sum(delta^2) / sum(sol^2)), same convention as the other fitters- verbose
print one line per scoring iteration with the relative step (the convergence criterion itself), and with
estimate = TRUEone line per evaluation of the Laplace -2logL instead- estimate
ordinal mode only: estimate the components as the minimum of the Laplace -2logL, with
startas the starting point, instead of taking them as given- maxiter_em, tol_em
no longer used: they controlled the EM step that the minimum of the Laplace -2logL replaced. Passing either (by name or in their old positions) gives a warning; use
max_evalsandtol_estimate- genotypes, blend, apy_core, vecchia_k
single step, as in
model(): the relationship terms get theH^-1ofh_inverse()instead ofA^-1- max_evals
with
estimate = TRUEand more than one component, the maximum number of Nelder-Mead evaluations of the Laplace -2logL, and again at each point of a correlation profile. It does not bound the Newton polish (up to 10 steps of 2k^2 + 1 evaluations for k components), the boundary check (one per variance), the standard errors (2k^2 + 1) or the number of profile points;n_evalsreports the total. Refused with a single component, whichstats::optimize()finds without such a bound- tol_estimate
with
estimate = TRUE, the tolerance of the estimation: the interval width ofstats::optimize()on the log of the variance for one component, and the largest Newton step (log variances and atanh correlations) that still moves the estimate for several- profile
with
estimate = TRUE:TRUE(default) profiles a correlation near +-1,FALSEskips the profile (each profile point re-optimizes every other component, which can cost more than the estimate) and only withholds its standard error. Refused when the model estimates no correlation
Value
an object of class breeding_fit_thr. Ordinal mode: thresholds (with
se_thresholds from the generalized inverse, the book's p.271 column), theta
(the components plus var(residual) = 1), se, b and se_b (named
term=level), ebv and pev per covariance group as in model() (one vector per
group, the blocks of its terms in sequence, each named by level; liability scale;
ebv() and accuracy() work), categories, and the convergence fields of every
fitter. With estimate = TRUE also n_evals, the evaluations of the -2logL,
vcov, the delta-method covariance of theta (zero for the fixed residual, NA for a
component without a standard error; absent when no component has one), and
profile, the correlation profiles when one
was run.
With k_inverse = it carries k_prior, the diagonal of the declared K named by
level (read from the selective inverse of K^-1), which accuracy() divides the
PEV by in place of 1 + F; with genotypes = or k_inverse = h_inverse(...) it
carries h_prior and h_prior_row as in model() instead. A fit with
k_inverse = made before k_prior existed has neither, and accuracy() then
divides it by 1 + F of the pedigree given: refit it.
predict() on the fit returns the per-category probabilities, the number the
book actually delivers (p.272-273). Joint mode: b and ebv come named
...|trait; nu holds the corrected liability solutions, b_regression the
residual regression, and G, R, Gc the matrices used; pev comes from the
generalized inverse, for u1 and for the ranking value u2 = nu + b1 u1, and
predict() gives the probability of Eqn 15.25.
Details
TWO CONVENTIONS, both forced by identifiability and both different from model():
the residual variance on the liability scale is FIXED at 1. The liability is not observed, so it is only defined up to scale; s2e = 1 sets the scale, as the book does throughout chapter 15.
there is NO implicit intercept: the thresholds take its place, and the FIRST level of every fixed class factor is dropped (solution zero), which sets the origin. Levels come in factor order when the column is a factor, alphabetical order otherwise. To compare against a table that zeroes some other level, compare CONTRASTS, or pass the column as a factor with the reference level first.
THE COMPONENTS ARE GIVEN BY DEFAULT, as the book uses the model (both examples of
chapter 15 fix them): start= holds the components on the liability scale. With
estimate = TRUE (ordinal mode) start is the starting point and the components are
ESTIMATED by approximate marginal maximum likelihood: the locations (thresholds,
fixed and random effects) are integrated by the Laplace approximation around the
mode, and the estimate is the MINIMUM of that Laplace -2 log likelihood, the very
number the fit reports in neg2logl: stats::optimize() in the log of the variance
when there is one component, Nelder-Mead followed by Newton steps on a numerical
Hessian (in log variances and the inverse hyperbolic tangent of the correlations)
when there are several. Earlier versions iterated the EM-type step of Foulley, Im,
Gianola and Hoeschele (1987) instead and reported the -2logL at its fixed point. That
point is NOT the minimum (the step ignores how the weights W move with the variance),
so the reported value was not a maximized likelihood, a likelihood-ratio test compared
values that were not maxima, and the standard error was read off the curvature away
from the optimum; validation/indirect_threshold_survival_recovery.R em_antigo measures
the gap against the code of that version.
neg2logl is comparable between fits of the same data and fixed effects; the standard
errors come from its numerical Hessian (Var = 2 H^-1, delta method to the components),
and vcov carries that covariance for se_function() and t2(). A correlation
between two terms of a group at or near +-1 is PROFILED instead (every other component
re-optimized on a grid of r): its delta-method standard error means nothing at the
boundary and is withheld, and the message gives the profile-likelihood 95% interval,
with the profile in profile (profile = FALSE skips the profile and still withholds
the standard error). A variance the data do not tell apart from zero (moved to 1e-6
with the rest fixed, the -2logL rises by less than 0.01) is held at its estimate with
the correlations it links, as ASReml holds a boundary component, and has no standard
error; the others come from the curvature in the free components. Where that
curvature is singular or nearly so (a combination of the components the data do not
identify), only the standard errors that depend on the combination are withheld, and
the message names them. An
indirect() term whose pens all hold a single animal has no incidence, and its
estimation is refused.
Read the estimate knowing its known bias: with a binary trait and few records per
level of the random effect the Laplace approximation UNDERESTIMATES the variance
(Tempelman 1998); many daughters per sire is where it is reliable. Gibbs sampling with
data augmentation (gibbs() with family = "probit" for a binary trait, or
THRGIBBS1F90) does not use the approximation. With indirect() on a binary trait
neither recovers the components: on pens of 2 to 8 with 3200 records and 4000 animals
(30 data sets; validation/indirect_gibbs_probit_recovery.R), where REML on the
liability itself gives 1.011 / 0.203 / 0.300 against 1 / 0.2 / 0.3, the Laplace
minimum sat at |r| > 0.99 in 29 of 30, averaging 0.204 / 0.173 / 0.150 over the 30
(bias of var(animal) -0.796, SE 0.010), and the gibbs() probit chains reached a bulk
ESS of 400 in 0 of 30.
INDIRECT GENETIC EFFECTS. indirect(id, pen = "pen", dilution = d) enters the
liability exactly as it enters the phenotype in model(): record i receives
(n_i - 1)^(-d) sum_j a_S,j over its distinct pen mates j, with n_i the number of
distinct animals of its pen counted over EVERY row of data (a mate without a
phenotype is still a mate). Declared in one group = with animal(), the direct
and indirect values share the 2x2 matrix G0 and the penalty kron(G0^-1, A^-1).
The pen column follows the engine's labelling (a numeric pen 10 and a text pen "10"
are the same pen), and a pen that is NA, NaN, Inf, blank text or the text "NaN" is
refused with its row: such rows have no known mates. Add random(pen) for the
environmental part shared by pen mates; leaving it out inflates the indirect
variance, which then absorbs it. t2() gives T2 and the direct heritability on the
liability scale, with the residual 1 in the denominator.
DECLARED LIMIT: the probit fixes the residual of every record at 1, and with an
environmental effect of the mates (variance s2_ES) the residual of a record in a pen
of n animals is s2_ED + (n - 1)^(1 - 2d) s2_ES, of which
(n - 1)^(-2d) (n - 2) s2_ES is shared with each mate (taking the direct and social
environmental deviations of one animal as uncorrelated; with a covariance between
them the shared part can be negative, and no random(pen) variance represents it).
One scale for every record is exact only when all pens have the same size (the
shared part is then random(pen)); with pens of unequal size the latent scale
changes with the pen size, which needs a heteroscedastic probit that this fitter does
not have, and the residual 1 that t2() divides by is then a convention, not the
variance of every record.
MEASURED BIAS, at prototype scale only (2329 records, one per animal, three
categories, 8 replicates; validation/indirect_threshold_survival_recovery.R): the
direct variance estimated with indirect() ran 0.031 (SE 0.007) below the direct-only
fit of the same direct values without an indirect effect, where REML on the same
liability observed gives 0.007 (SE 0.008). The extra is consistent with the Laplace
approximation (it shrank with ten categories) but has not been isolated from it.
JOINT QUANTITATIVE + BINARY ANALYSIS (Foulley et al. 1983; section 15.3 of the
book). With cbind(quant, bin) on the left-hand side the fit is the joint one: a
linear model for the quantitative trait and a threshold model for the binary one,
tied by the residual regression b1 = r12 / (sqrt(r11) sqrt(1 - r12^2))
(Eqn 15.22) and the reparametrized genetic covariance Gc = C G C' (Eqn 15.16).
In this mode start = list(G =, R =) with the 2x2 genetic and residual matrices
(R[2,2] on the liability scale), the FIRST fixed class factor keeps all its
levels (it carries the quantitative intercept and absorbs the threshold, which the
book says "has no interest", p.275), exactly one relationship term is accepted, and
the EBV of the binary trait is u2 = nu + b1 * u1 – the value the book recommends
for ranking (p.281), already using the quantitative information. The genetic
parameters reported in theta are read from G, NOT from Gc: reading h2 off Gc is
the documented trap (g22 of Gc gives 0.117 where the book prints 0.18).
The structure checks of model() (see its section "Structure checks") apply here: with
estimate = TRUE a design that does not identify a component stops (all the records in
one pen of indirect() gave components the data do not identify, reported as converged,
before this version),
with the components given it is a warning, and the warnings go to message.
References
Gianola, D. & Foulley, J.L. (1983) Sire evaluation for ordered categorical data with a threshold model. Genet. Sel. Evol. 15, 201-224. 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. Genet. Sel. Evol. 15, 401-424. Mrode, R.A. & Pocrnic, I. (2023) Linear Models for the Prediction of the Genetic Merit of Animals, 4th ed., chapter 15. Foulley, J.L., Im, S., Gianola, D. & Hoeschele, I. (1987) Empirical Bayes estimation of parameters for n polygenic binary traits. Genet. Sel. Evol. 19, 197-224. Tempelman, R.J. (1998) Generalized linear mixed models in dairy cattle breeding. J. Dairy Sci. 81, 1428-1444. Misztal, I., Gianola, D. & Foulley, J.L. (1989) Computing aspects of a nonlinear method of sire evaluation for categorical data. J. Dairy Sci. 72, 1557-1568. Bijma, P. (2010) Multilevel selection 4: modeling the relationship of indirect genetic effects and group size. Genetics 186, 1029-1031.
Examples
# Example 15.1 of Mrode & Pocrnic: calving ease in three categories, sire model,
# var(sire) = 1/19 (h2 = 0.20 on the liability scale)
sub <- data.frame(
herd = c(1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2),
sex = c("M","F","M","F","M","F","M","F","M","F",
"M","M","F","M","F","M","M","F","F","M"),
sire = c(1,1,1,2,2,2,3,3,3,1,1,1,2,2,3,3,4,4,4,4),
c1 = c(1,1,1,0,1,3,1,0,1,2,1,0,1,1,0,0,0,1,2,2),
c2 = c(0,0,0,1,0,0,1,1,0,0,0,0,0,0,1,0,1,0,0,0),
c3 = c(0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,1,0,0,0,0))
ind <- do.call(rbind, lapply(seq_len(nrow(sub)), function(i)
data.frame(herd = as.character(sub$herd[i]),
sex = factor(sub$sex[i], levels = c("M", "F")),
sire = as.character(sub$sire[i]),
score = rep(1:3, times = as.integer(sub[i, c("c1", "c2", "c3")])))))
ped <- data.frame(id = as.character(1:4), sire = c("0", "0", "1", "3"),
dam = rep("0", 4))
fit <- model_threshold(score ~ herd + sex + sire(sire), data = ind,
pedigree = ped, start = 1/19, verbose = FALSE)
fit$thresholds # 0.4378 and 1.0675 on p.271
#> t1 t2
#> 0.4377702 1.0675346
ebv(fit) # the four sires on the liability scale
#> 1 2 3 4
#> -0.04342738 0.05920912 0.04123085 -0.06600173