Fits a model with AR(1)/CAR(1) residual by AI-REML

The structure checks of model() (see its section "Structure checks") apply here, with a record that has no finite time read as a record with no observation (it still counts as a pen mate of indirect()); the AR(1) residual with the animal as subject models the repetition, so repeated records of animal() give no warning then.

Usage

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

Arguments

formula

as in model(), including reaction norm (rn(id, base = ...)); a cbind(t1, t2) left-hand side fits the multi-trait AR(1)/CAR(1) with the separable residual Gamma (x) R0 (full R0 between traits, one rho in time). A record with some traits missing stays in the series: each missing cell is carried as its own fixed effect (the mv device of ASReml), which keeps the residual separable and gives exactly the marginal likelihood of the observed cells; a record with no trait observed leaves. A fixed level with no observation of one trait is refused by name. Gates: V form over the observed cells, a REML written from the raw data, collapse at rho = 0 onto the multi-trait fitter, finite differences on every parameter

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

subject

column identifying the subject (typically the animal). A missing subject, NA, NaN, Inf, -Inf, blank text or the text "NaN", is refused on every row with the column, the number of such rows and the first of them, as a level column in model(): those rows would form one series. So is the unknown code 0 when the subject column is also the column of a relationship term (animal(id) with subject = "id"): there 0 is an unknown animal, and distinct unknown animals would be chained as one series

time

numeric time column (a Date enters as days, a POSIXct as seconds); two records of the SAME subject at the SAME time are a declared error: with AR(1) the time identifies the record, and simultaneous repetition calls for a permanent environment effect. The field reads the axis in one of two ways: the discrete AR(1) on a grid (corAR1 of nlme, an integer covariate with the correlation phi^|lag| in its units; AR1 of ASReml, for equally spaced levels) and the continuous CAR(1) (corCAR1 of nlme, with the parameter in (0, 1); EXP of ASReml). Here the axis is a GRID when every gap between consecutive records of a subject is a whole number of one step: the step is the largest number that divides every gap, sought among the smallest gap divided by 1 to 1000. The rule reads only the ratios between gaps, so it does not depend on the unit or the origin of the time. On a grid the time is divided by the step, rho is the correlation PER STEP, in (-1, 1), and rho^k the correlation of records k steps apart. Dividing by the step is BreedingR's own parametrization (in nlme the step is the unit of the covariate); the step is returned in time_step, and when it is not 1 a message names it (also in message). So the same records give the same rho written in days, in seconds (POSIXct), in weeks or as days x 2. Otherwise the time is CONTINUOUS: rho is per time unit, the correlation is rho^dt, and rho lies in [0, 1), because a continuous-time Markov process has no negative correlation.

On a grid a negative correlation exists only between records an ODD number of steps apart. Where at least 10% of the gaps are odd and at least 10% even numbers of steps (records on Monday, Wednesday and Friday: gaps of 2 and 3 days) rho and -rho give the same correlation at the even gaps and the opposite one at the odd gaps, and the likelihood can have a mode of each sign. There the -2logL is also evaluated at the final point with rho of the other sign; when it is lower, a second walk starts from it and the better fit is returned (start_from "sign-flipped rho", the first fit in message).

The fit is then CHECKED in two situations. (1) A few gaps (at most 20% of them, or one, and each gap length among them in at most 5% of the gaps) are not whole numbers of a step H that every other gap is a whole number of, H larger than the step of the grid (on a continuous axis, any such H that leaves some gap out; H is the largest step that leaves at most those few out): the trace of a few records with a wrong time. One record at t = 4.5 among daily records makes the grid half days, where consecutive days are two steps apart and cannot correlate negatively; one Wednesday record at noon among records on Monday, Wednesday and Friday does the same. The data are fitted once more with those records moved to the nearest free point of the grid of H (a run of records joined by gaps that are whole numbers of H moves as a whole, a tie goes to the earlier point, and when a record has no free point between its neighbours the order of the records within subject is the time instead). (2) Some gap is not a whole number of the typical (median) gap and the correlation at the typical gap is near 0, |rho^k| < 0.05, on an axis that cannot make it negative for nearly all consecutive records (continuous, or fewer than 90% of the gaps an odd number of steps), or on one that can but where no gap is a single step (there rho = 0 is a stationary point the walk can stop at): the data are fitted once more with the order of the records within subject as time (equally spaced levels, as AR1 of ASReml reads them). Either way the new time goes in a new column, and the time column stays in the data for the formula, which may read it as a class or a covariate. The data are also fitted with no serial correlation (model(), or model_mt() under cbind()), on the records of the fit (a record with no time is out of both). The gain is the -2logL of the better of the time-axis fit and the fit with no serial correlation minus the -2logL of the refit: the likelihood ratio of rho = 0 in the refit (chi-square 1 df). With a NEGATIVE refit rho that the axis cannot represent and a gain above 10.83 (0.1%) the fit stops with an error that says what to do and names the gaps off the grid (how many, in how many subjects, the first ones in the order of the data with their length); with a gain from 3.84 (5%) to 10.83 a warning says the same and the fit on the time axis is kept. With a refit rho that the axis can represent (positive; negative at a step H that is an odd number of steps; or, in the second situation, negative with 90% of the gaps odd) and a gain above 3.84, the axis is fitted again from the correlation the refit gives at its step and the better fit is kept (start_from "refit from the dominant grid" or "refit from the order of the records"); a gain still above 10.83 means the spacing contradicts the correlation (a duplicated or mistyped time) and stops the fit with an error that names the gaps. The high bar is because with no serial correlation the refit is better by 3.84 in about 5% of data sets, half of them with a negative rho. A record at t = 2.01 among daily ones with a positive correlation is the second case. With one trait the check costs two more fits (three when the axis is fitted again), and only when it runs, which includes ordinary data: weekly records with an occasional gap of 8 or 10 days (the first situation) and irregular integer-day records with a weak correlation at the typical gap (the second); under cbind() with t traits it costs 2(t + 1) more fits (the refit is a multi-trait AR(1) fit with its own t univariate AR(1) starts, and the fit with no serial correlation a model_mt() fit with its own t univariate starts), one more when the axis is fitted again. The fit on such an axis can itself take every iteration of maxiter before the check, with rho flipping sign around 0. With POSIXct clock times in a time zone with daylight saving, records at the same clock time are 23 or 25 h apart across the change, so the step becomes 1 h; the messages say so. On a continuous axis a negative rho in start (or in the theta of eval_internal_ar1()) is refused, with that reason

genotypes

list with ids and m (0/1/2 matrix) for single-step; NA is imputed with 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

vecchia_k

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

missing_code

missing-value code for observations

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, under cbind(), the start is the univariate AR(1) optimum of each trait (model_ar1() with the same right-hand side, one trait at a time, run to its own convergence), every covariance between traits at zero and rho the mean of the univariate ones; a variance that a univariate fit left at the zero boundary starts at 1% of the var(y) start, the rule of model_mt(), where the evidence for this default and its measured cost are. With one trait, when a univariate fit fails, or with maxiter = 0, the start comes from var(y), divided by the geometric mean of a declared kernel's eigenvalues so that the same model with K and with c * K starts at equivalent points. When the fit from the univariate start stops short of maxiter without certifying an optimum, the var(y) start is fitted as well and the better of the two is returned, as in model_mt(). 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; the univariate fits of the default start are not bound by it, as in model_mt(). maxiter = 0 evaluates the start and reports converged = FALSE, as model() does: 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, after the checks on the likelihood, the second walk from inside the cone and the polish at a boundary described under tol in model_mt()), and, as in model(), converged requires it under 2e-4: the step walks in log-Cholesky coordinates, so it cannot jam against a boundary. At maxiter with the decrement already under the tolerance and those 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. The checks cover covariance blocks only: with one trait, a fit where rho goes to 1 and the residual variance to its floor (records repeated within a subject) can still be certified on an unbounded ridge, with the SINGULAR warning

metafounders

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

gamma

the base relationship matrix of the metafounders. Either a vector of length n_metafounders, read as the DIAGONAL, or a full symmetric n_metafounders x n_metafounders matrix. The off-diagonal gamma_jk is the ancestral relationship BETWEEN two base populations, which is what a multibreed analysis turns on; it may be NEGATIVE, for bases pulled apart by selection in opposite directions. Admissibility is positive SEMI-definiteness, tested by a spectral decomposition and not by a Cholesky, plus a diagonal below 2 so that a metafounder's offspring keeps a positive Mendelian variance. A SINGULAR gamma is accepted, through the Moore-Penrose pseudo-inverse: that covers gamma = 0, the unknown-parent-group limit, where the pseudo-inverse reproduces the A-inverse of unknown parent groups exactly, and two metafounders standing for one population, whose rows are identical. What is still refused is an INDEFINITE gamma, a negative eigenvalue, which does not generate a covariance matrix at all

verbose

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

Value

besides the components, rho(residual) with its standard error; |rho| >= 1 is never a result: a step that leaves the interval is rejected like any inadmissible theta. The fixed-effect solutions come in b, named term=level (with |trait appended under cbind()); the parametrization note of model() applies – dropped columns are in dropped_x and only contrasts compare against a reference-level convention. 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 where the record did not enter, with used, subject_index and time (in the unit of the data), and time_step, the step of the grid rho is per (NA on a continuous axis); residuals.breeding_fit() reads them. 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. The pattern of the equations holds only the pairs of columns of one record and of neighbouring records (R^-1 is tridiagonal within subject), so a long series with a time contemporary group keeps a sparse factor, and dense_block stays small. dense_block, h_prior, h_prior_row, k_prior and k_level_type are as in model(): in a single step, accuracy() divides a genotyped animal by its diagonal of G* and not by 1 + F, and a kernel() level by its diagonal of K. start_from is "univariate fits", "start=", "var(y)", "sign-flipped rho", "refit from the dominant grid" or "refit from the order of the records" (see time), and start the start the walk left from, named like theta (a start beyond the floors of the step is raised to them first, as in model_mt()). dropped_x_reason, parallel to dropped_x, says "no record" or "aliased" (a reference level or a collinear column); a level with no record of one trait stops the fit with an error that names it. The walk runs on each trait divided by the power of two just above its standard deviation, an exact rescaling in floating point (the variances have units and rho has none, so even with one trait the step criterion and the singularity check used to depend on the unit of y), the singularity verdict is taken on the information matrix in those units, and the result is evaluated back in the original units: y and 100 * y give the same verdict, standard errors scaled by 100^2 and the same standard error of rho

References

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.

Pinheiro, J.C. & Bates, D.M. (2000). Mixed-Effects Models in S and S-PLUS. Springer, section 5.3 (corAR1, integer time; corCAR1, continuous time with the parameter in (0, 1)).

Butler, D.G., Cullis, B.R., Gilmour, A.R., Gogel, B.J. & Thompson, R. (2017). ASReml-R Reference Manual, version 4. VSN International (ar1() for equally spaced levels, exp() for continuous distance).