Estimate the pen-size heterogeneous residual of the associative model

Fits var(e_i) = s2_ED + (n_i - 1)^(1 - 2d) s2_ES, with n_i the number of distinct animals in the pen of record i and d the dilution of the indirect() term in the formula (d = 0, the default, gives the book's (n_i - 1) s2_ES). The same d that dilutes the genetic side dilutes the social environmental one, so the residual and the genetic incidence describe one model. The ratio k = s2_ES / s2_ED is profiled: each candidate k is a model() fit with weights = 1 / (1 + (n_i - 1)^(1 - 2d) k) (exact REML at that k), the reported -2logL is brought to the standard scale by subtracting sum(log w) (the weight jacobian, constant in theta but not in k), and the minimum over k is located by a coarse grid followed by golden-section search (stats::optimize()) in the bracketing interval. The declared stopping rule is the golden-section interval width falling under tol_k.

Usage

indirect_residual(
  formula,
  data,
  pedigree = NULL,
  k_max = 5,
  n_grid = 6L,
  tol_k = 0.001,
  verbose = interactive(),
  ...
)

Arguments

formula: the same formula passed to model(); it must contain an indirect() term, which is where the pen column comes from
data: data.frame with the columns referenced
pedigree: data.frame animal, sire, dam, as in model()
k_max: upper end of the k search; a minimum found AT this boundary is reported in message and asks for a rerun with a larger k_max
n_grid: points of the coarse grid on [0, k_max] that brackets the minimum before the golden-section refinement
tol_k: absolute tolerance on k: the golden-section search stops when its interval is narrower than this
verbose: print one line per profiled k
...: passed on to every inner model() call (maxiter, tol, n_em, ...). weights and start are used by the profile itself and are refused here

Value

a list of class breeding_indirect_residual: k (the estimate), s2_ED (the var(residual) of the fit at k), s2_ES (k * s2_ED), fit (the model() object at the estimated k, whose -2logL is on the WEIGHTED scale as always), profile (data.frame of every k evaluated and its corrected -2logL), n (the pen size of each record) and message

Details

The pen and the animal columns are read from the indirect() term of the formula, so the pen sizes here are exactly the ones the incidence used. Declared limit: the profile handles the residual VARIANCE heterogeneity only. The residual COVARIANCE between pen mates is a different object, and this function does not carry it.

A previous version of this page said that covariance is absorbed by a random(pen) term. It is not, in the regime this function exists for. The pen-mate covariance is (n - 2) s2_ES, which GROWS with pen size, while a random(pen) term has one variance per pen: the two agree only when every pen is the same size, and equal pens are exactly the case where s2_ES does not separate from the intercept anyway. Use associative_matrix() instead, which builds the exact structure s2_ED * I + s2_ES * D as a declared kernel() covariance; s2_ES then comes back as an ordinary component with a standard error, estimated jointly, and neither the profile nor the weight jacobian is needed. This function stays as the cheap route when only the variance heterogeneity is wanted.

k is profiled and not walked by a gradient because a one-dimensional profile is cheap and cannot be trapped the way a joint update can.

Read the answer knowing what the data identifies. What unequal pens pin is the SLOPE s2_ES of the residual variance against n_i - 1; the intercept s2_ED separates from the genetic variance only through relatives when every animal has one record, so the RATIO k inherits that noise and its profile can be flat on the high side (measured on 780 simulated records with a planted k of 1: s2_ES came back between 5.8 and 6.7 for a truth of 6 across seeds, k between 0.77 and 1.77 – the gate in test-indirect-dilution.R encodes exactly this). profile is returned so the flatness can be seen instead of guessed; pens of size 1, whose residual is s2_ED alone, are the observations that pin the intercept.

References

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.

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

Mrode, R.A. & Pocrnic, I. (2023). Linear Models for the Prediction of the Genetic Merit of Animals, 4th ed. CABI, ch. 9.