Fits a multi-trait mixed model by AI-REML

The structure checks of model() (see its section "Structure checks") apply here, on the rows with some trait observed and then on the rows of each trait, since the fixed effects are per trait (a random term nested in a fixed effect among the records of one trait stops the fit, naming the trait); the warning about repeated records names the first trait that has them. The univariate fits of the start run with those warnings silenced, so each is given once.

Usage

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

Arguments

formula

with cbind() on the left-hand side: cbind(p1, p2) ~ cg + animal(id). Each effect applies to all the traits unless it has traits = (below). Each covariance group gains the GENETIC covariances between traits, and the residual becomes the full R0 matrix. A missing trait (NA or the missing_code code) does not drop the record: it enters with the R0 submatrix of the traits that it has.

traits = puts a random term on some of the traits only. The usual case is a trait with repeated records and a trait with one record per animal: cbind(y1, y2) ~ cg + animal(id) + pe(id, traits = "y1"). On a trait with one record per animal, the data do not separate pe() from the residual. The information matrix is then singular, and every standard error is NaN. With traits = "y1" the term has no component and no solution on y2: the names are var(pe@y1) and id|y1. This model has no non-genetic covariance between y2 and the permanent effect of y1: it is 0 by construction. To estimate that covariance, put pe() on both traits. The data then separate cov(pe@y2,pe@y1), but not var(pe@y2) from var(res@y2), and the fit says that the information matrix is singular along that pair. The terms of one covariance group (group =) can name different traits. An example is a random regression on y1 and a one-record trait y2: rn(id, base = c("phi0", "phi1"), traits = "y1", group = "g") + animal(id, traits = "y2", group = "g"). This G is over rn[0]@y1, rn[1]@y1 and animal@y2, and gives the covariance of each coefficient with the breeding value for y2. A fixed effect can be on some traits only too: fixed(cg, traits = "y1") for a class and cov(x, traits = "y1") for a covariate. On the other traits the effect has no equation. dropped_x does not list it, because no record is missing and no column is collinear with others. A column of a term may be missing (NA, blank) on a row that has none of the traits of the term. That cell does not enter the likelihood. The fit puts a value there that changes nothing: the first label of the column, or the mean of a covariate. On the other rows the rule of model() stands. The univariate fits of the default start leave the term out of the traits that it does not name.

A y2 record can be in the row of a y1 record, or in a row of its own with y1 missing. When no row has both traits, the data do not estimate the residual covariance between them. The fit then holds it at 0 and gives it no standard error, and message says so. The fit estimates the genetic covariance in both layouts, through the relationship. Gates: tests/testthat/test-mt-traits.R (the identity with the univariate fits, the dense V form, the score); recovery in validation/mt_traits_recovery.R (simulation).

data

data.frame

pedigree

data.frame animal, sire, dam. The ids are matched to the data by the rule of model() (see pedigree there): 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

genotypes

list with ids and m (0/1/2 matrix) for single-step; NA is imputed by the mean, as in model()

blend

weight of A22 in the G blend (0.05 by default)

apy_core

ids of the APY core, or "auto" as in model(); NULL uses the exact inverse of G

missing_code

missing-value code

start

starting values for the components, in the order the fit reports them. Use it to warm-start from a submodel, or to check that the optimum does not depend on where the search began. Without it the start is the univariate optimum of each trait (model() with the same right-hand side, one trait at a time, run to its own convergence whatever maxiter is) with every covariance between traits at zero, the start one would build by hand in BLUPF90 or ASReml. A variance that a univariate fit left at the zero boundary (below 1e-4 of the var(y) start) starts at 1% of the var(y) start, so the walk does not begin on the wall; every other component starts where the univariate fit left it. The univariate fits cost one model() each. When any of them fails, or with maxiter = 0, the start falls back to var(y). When the fit from the univariate start stops short of maxiter without certifying an optimum, the var(y) start is fitted as well, a second whole fit, and the better of the two is returned with the other named in message: the lower -2logL, unless only one of them certified and the other is within 2e-4 of it, and then the certified one. seconds counts both fits and iters is the one returned. Evidence for this default: on real data, with the engine of the time (before the walk was standardized per trait and before the second walk from inside the cone described under tol), a bivariate fit with an autoregressive residual from the var(y) start walked to the r = -1 and r = +1 boundaries and stopped 41.3 -2logL units above the optimum that the univariate start reached in 6 iterations. With the current engine no simulated case is known where the var(y) start fails: over 24 bivariate data sets in four designs (trait scales up to 100 apart, litter and permanent-environment correlations up to 0.95) both starts certified, at most 1.8e-4 -2logL units apart (the boundary spread described under tol), and that real case has not been rerun. So this default rests on that one real case. Measured cost where both starts reach the same optimum, 3 replicates on a machine shared with other jobs (validation/unit_shift_start_scale.R): an interior model_mt() with 24 000 records took 22.3 to 24.1 s against 14.8 to 21.8 s from var(y), one with the litter variance of a trait at zero 45.8 to 50.4 s against 41.8 to 48.8 s (both pay the second walk described under tol, 16 iterations instead of 8), medians 1.6 and 1.1 times. Which start was used is in start_from.

maxiter

maximum number of iterations; a fit that hits the ceiling says so in message and how to raise it (converged = TRUE at the ceiling only when the Newton decrement certifies a flat direction, as in model()). The univariate fits of the default start are not bound by it: the start is their optimum, so a small maxiter still pays for them; the second, var(y) fit runs only when the first stops short of maxiter. maxiter = 0 evaluates the start and reports converged = FALSE: an evaluation is not a convergence

tol

relative tolerance on the components, sqrt(sum delta^2 / sum theta^2); see the BLUPF90 scale note in model(). The Newton decrement g' AI^-1 g is computed at the final point and reported in newton_dec (components resting at a covariance boundary excluded), and, as in model(), converged requires it under 2e-4: the step walks in log-Cholesky coordinates, so it cannot jam against a boundary. The step freezes a boundary block whose curvature has collapsed, and the decrement leaves out every component of a block resting at a boundary, so such a block is excluded only after checks on the likelihood itself. First the Karush-Kuhn-Tucker condition: points inside the cone along the negative eigenvectors of the block's gradient matrix, with the length of the quadratic model of the AI (and, as a fallback, each Cholesky diagonal raised to 0.1%, 1% and 10% of its scale). Then the boundary face itself: the directions of the block that move along it, by central differences of the -2logL (step 1e-2 in log-Cholesky), with a joint Newton step together with the free components. A point that lowers the -2logL by more than 1e-4 sends the walk there, or refuses the certificate when the walk cannot move (at maxiter, after 10 re-entries into the cone or 20 moves along the face). Last, a fit that ends with a block at a boundary, certified or stopped short of maxiter, is walked again from inside the cone, within the same maxiter budget: each such block made diagonal with its variances at least 10% of the var(y) start, the rest where the first walk left it. The second walk is returned when it lowers the -2logL by more than 1e-4, or certifies within 2e-4 of a first walk that did not, and message says what it found either way; a second round follows only when the second walk was lower by more than 1e-4. The second walk costs a whole walk more in every fit that ends at a boundary (measured at 24 000 records with a litter variance at zero: 16 iterations instead of 8). Measured on real data before these checks: a covariance boundary 46.8 -2logL units above the interior optimum had been certified; on simulated starts at the boundary (48 fits) the KKT probe alone still certified 6 of them 0.0056 to 4.73 units above the optimum, and with the second walk none. A boundary face is flat, so the certified point then depended on the path: in sweeps of four starts per data set, up to 3.9e-4 -2logL units above the best of the four. A certified fit with a block at a boundary is therefore polished: for a trait whose variance sits at the floor, the same point with that trait's covariances at zero, and the Newton step along the face, each kept while it lowers the -2logL by more than 1e-5 and the new point still certifies, up to 10 times. It never turns a certified fit into an uncertified one. Measured after it, in the same sweeps (368 fits, the default, the var(y) start and two starts at the wall per data set), no certified point was more than 9.4e-5 above the best of its four starts (5 were 1.0e-4 to 3.9e-4 above before); over 24 other data sets (48 fits, the two default starts each) one pair still certified points 1.7e-4 apart, a model_mt() fit with two blocks at a boundary. These are searches, not a proof. At maxiter with the decrement already under the tolerance and the checks passed, the fit is declared converged and message says the likelihood is flat along some direction. A stop where no damped step and no EM step improves the likelihood is never certified: message gives the iteration and the decrement

metafounders

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

gamma

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

verbose

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

Value

an object of class breeding_fit_mt; besides the fields of model(), the fixed-effect solutions b come named term=level|trait, one entry per column of X and per trait. The parametrization note of model() applies: dropped columns are in dropped_x and only contrasts compare against a reference-level convention. A column dropped for all traits is named term=level; one dropped for a single trait is term=level|trait. The reference level of a class is chosen among the levels with records of every trait whenever one exists, so a |trait drop is a level with no record of that trait, and dropped_x_reason, parallel to dropped_x, says which case each drop is: "no record" or "aliased" (a reference level or a collinear column) for a column dropped for all traits, and "no record of y2" or "aliased within y2", with the trait's name, for one dropped for a single trait. start_from is "univariate fits", "start=" or "var(y)", and start the start the walk left from, named like theta: a start beyond the floors of the step (a Cholesky diagonal below 1e-3 of the largest in its block, a correlation within about 1e-6 of +/-1) is raised to them first, because the walk cannot represent it and the likelihood there is evaluated with lost digits (measured: 0.0028 and 51.3 -2logL units below the dense V form at 1 - r^2 = 1e-12 and 1e-15, and the fit stopped at its first iteration); with maxiter = 0 the start is evaluated as given. dense_block, h_prior, h_prior_row, k_prior and k_level_type are as in model(). The walk runs on each trait divided by the power of two just above its standard deviation (an exact rescaling in floating point, so y and 1024 * y give bit for bit the same fit) and the result is evaluated back in the original units, so a trait recorded in other units (y1 versus 0.01 y1) gives the same estimates up to the scale factors and the same -2logL up to the Jacobian 2 (n - p) log(0.01), with n the records and p the fixed equations of that trait. The singularity verdict is taken on the information matrix in those standardized units. With a factor that is not a power of two the data themselves round differently: an interior optimum agrees to about 1e-6, but components resting at a boundary sit on a flat face of the likelihood and can move along it (measured: 0.4% and 0.6% for y2 times 0.01 and 1000 in a two-term group with ten components at a covariance boundary, the -2logL equal to 2.5e-5 after the Jacobian). vcov is 2 AI^-1, the sampling covariance of the components, with the names of theta; a component with no standard error has its row and column NA. A block of the step can stop at a wall (a correlation at +/-1). Then vcov is conditional on that wall, as in model(), with one rank less. se_function() reads it for r_g and h2. residual, fitted, fitted_fixed and resid_var (the variance of each residual) are matrices with one row per row of data and one column per trait, NA in a missing cell, read by residuals.breeding_fit(). They come from one more assembly, factorization and selective inverse at the final components, at the end of every fit (about the cost of one more REML iteration); if that pass fails, the fit is kept with NA residuals and the reason in message.