Fit a Weibull proportional hazards frailty model — model

The trait is the time until failure (length of productive life, days to culling) and some records are RIGHT-CENSORED: the animal was still alive when recording stopped, so its time is a lower bound, not an observation. The hazard is h(t) = rho * lambda * (lambda t)^(rho-1) * exp(x'b + z'a) (Eqns 16.3 and 16.4 of Mrode & Pocrnic): a Weibull baseline scaled by the fixed risk factors and by the log-frailty a, which carries the usual A^-1 penalty. The mode of the joint log-likelihood (Eqn 16.5) is found by damped Newton, whose fixed-point form is exactly the iterative system of Eqn 16.6 (Kachman 1999).

Usage

model_survival(
  formula,
  data,
  pedigree = NULL,
  censor = NULL,
  rho = NULL,
  lambda = NULL,
  sigma2 = NULL,
  k_inverse = NULL,
  maxiter = 200L,
  tol = 1e-08,
  verbose = interactive(),
  entry = NULL,
  subject = NULL,
  gaps = c("error", "allow"),
  genotypes = NULL,
  blend = 0.05,
  apy_core = NULL,
  vecchia_k = NULL,
  mates = c("all", "present"),
  start = NULL,
  max_evals = 500L,
  tol_estimate = 1e-06,
  profile = TRUE
)

Arguments

formula: fixed class effects, cov() covariates, and one or more random terms among animal(), sire(), random() and indirect() – the frailties – with group = to put several in one covariance matrix. rn(), kernel(), nested = and sire(mgs =) are not available here (k_inverse= is the door for a relationship matrix the pedigree cannot build). The left-hand side is the time column, numeric and strictly positive: a factor or text time is refused, naming the column, rather than read as its codes 1, 2, ...
data: data.frame with the columns referenced. Records with a missing time or a missing censoring indicator are dropped and counted; they still count as pen mates of indirect(). 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 in model(); a factor or text column in cov(), in the time or in entry= is refused rather than read as its codes. There is no missing_code here: a missing time is NA, and a common missing-value code in the time (9999 held by two records or more, outside the range of the others) is warned about, with the first row (see the rule under formula in model()). A cov() covariate enters standardized, as in model(): centred and scaled when lambda is estimated (the intercept takes the centre), only scaled when lambda is given (there is no intercept, and a warning names a covariate whose values sit more than 1e3 SDs from zero); b, se_b and lambda come back in the unit of the data. A level of a fixed class effect with no failure among its records (all right-censored; with elementary records, among its pieces) has no finite relative risk, the monotone likelihood of Heinze and Schemper (2001): its solution runs to minus infinity, or, as a reference level held at the baseline when lambda is given, the failures it lacks are loaded onto the frailties of its animals. The fit is refused naming the level and its count: merge the level, drop its records, or fit the effect as a random frailty. When every record is censored there is no failure to estimate anything from, and the fit is refused unless rho, lambda and sigma2 are all given and there is no fixed effect; then it predicts the frailties from the given values alone, and the fit message says so. Separation by a cov() covariate (every record with a nonzero value censored) is not checked: it shows as converged = FALSE with a diverging solution.
pedigree: data.frame animal, sire, dam; required with animal(), sire() or indirect() unless k_inverse is given. 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 (with k_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 records with an unknown 0, the founders' records and the blocks are counted on the records the fit uses, those with the time, the censoring indicator and entry observed
censor: MANDATORY: the censoring indicator, as a column name or a vector – 1 (or TRUE) for a complete record, 0 (or FALSE) for a right-censored one. If no record is censored, say so explicitly with a column of ones. A factor or text indicator is read by its LABELS ("0" and "1", or "TRUE" and "FALSE"), never by the codes 1, 2, ... of the factor, and the fit message says so; any other label ("dead", "alive") is refused with the labels it found. This differs from survival::Surv(), which reads a factor event as a multi-state endpoint whose FIRST level is censoring: here factor(ev, levels = c(1, 0)) is still 1 complete.
rho: the Weibull shape: NULL (default) estimates it inside the Newton system; a positive number fixes it (rho = 1 is the exponential, the book's choice in Example 16.1). Risk rises with age when rho > 1 and falls when rho < 1.
lambda: the Weibull scale: NULL (default) estimates it through the intercept; a positive number fixes the baseline and drops the intercept (lambda = 1 reproduces the book's parametrization, reference classes at risk 1)
sigma2: the frailty components: NULL (default) estimates them by Laplace; numbers take them as GIVEN, in the order of model() (fit$theta names them: covariance groups first, the lower triangle of each by columns, then the terms without a group), which with a single frailty is its variance (Example 16.1 publishes its solutions under sigma2 = 0.4 – see the note in the gate file about the misprinted 20)
k_inverse: the inverse of the relationship matrix for the relationship terms, replacing the A^-1 built from pedigree (every relationship term uses it): either triplets list(i, j, x, n, id) as a_inverse() returns, or a dense symmetric matrix with the ids as dimnames
maxiter: maximum Newton iterations per inner fit
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 Newton iteration, and one per evaluation of the Laplace -2logL when the components are being estimated
entry: column with the start of each elementary record (0 for the first piece of a subject observed from time zero); NULL means every row starts at 0. It is numeric, as the time: a factor or text column is refused
subject: column with the subject of each elementary record; required with entry. The level of every frailty term other than indirect() must be the same in every piece of a subject (so random(pen) refuses an animal that changes pen). A missing subject (NA, NaN, Inf, -Inf, blank text or the text "NaN") is refused on every row, as the subject of model_ar1(): those pieces would form one history. So is the unknown code 0 on the records the fit uses (time, censoring indicator and entry observed) when the subject column is also the column of a relationship term (animal(id) with subject = "id"): there 0 is an unknown animal, and the pieces of distinct unknown animals would be read as the history of one subject. Give them ids of their own (an id with no line in the pedigree enters as a founder) or drop those rows
gaps: "error" (default) refuses a gap between two pieces of a subject; "allow" accepts it as time out of observation, in which no risk is counted
genotypes, blend, apy_core, vecchia_k: single step, as in model(): the frailty of a relationship term gets the H^-1 of h_inverse() instead of A^-1
mates: indirect() with elementary records: "all" (default) counts every animal that was in the pen as a mate on every piece; "present" counts on each piece only the mates with a piece overlapping it, and dilutes by their number
start: the starting point of the estimation when there is more than one component, in the order of sigma2; by default 0.1 for every variance and 0 for every covariance. A single component is searched on the whole interval [1e-6, 1e4] and takes none: start= is then refused, as it is with sigma2=
max_evals: with more than one component estimated, 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, the boundary check (one evaluation per variance), the standard errors or the number of profile points; n_evals reports the total. Refused with sigma2= and with a single component, which stats::optimize() finds without such a bound
tol_estimate: tolerance of the estimation: the interval width of stats::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. Refused with sigma2=
profile: TRUE (default) profiles a correlation near +-1; FALSE skips the profile and only withholds its standard error. Refused with sigma2= and when no correlation is estimated

Value

an object of class breeding_fit_surv: rho, lambda (with se_log_rho when rho was estimated), theta (the frailty components, with an se from the curvature of the Laplace -2logL when they were estimated), b and se_b (named term=level, log relative risks), ebv and pev per covariance group as in model() (log-frailty scale; ebv() works, and exp(ebv()) is the RRS of the book's table), n_censored, loglik_joint (the penalized joint log-likelihood at the mode), marginal_loglik (the Laplace value), the convergence fields of every fitter, vcov (the delta-method covariance of theta, NA for a component without a standard error) when the components were estimated and at least one has a standard error, profile when a correlation was profiled, and predict() for relative risks and survival probabilities S(t) – the book's p.292 numbers. With k_inverse = it carries k_prior, the diagonal of the declared K named by level, 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.

Details

CENSORING IS NOT MISSINGNESS, and that is why this fitter exists. A missing time says nothing; a censored time says "at least this much". Treating a censored record as complete biases the evaluation against the animals still alive – the best ones – and dropping it throws away exactly the selection candidates. The censor= indicator is therefore mandatory: 1 for a complete (uncensored) record, 0 for a censored one. Left- and interval-censoring are outside this fitter.

TIME-DEPENDENT COVARIATES enter as ELEMENTARY RECORDS, the device of the Survival Kit (Ducrocq and Solkner): a covariate that changes during a life is constant by pieces, and each piece is one row (entry, stop] of the same subject, with its own values of the covariates and censor = 1 only on the last piece when the subject failed. The cumulative hazard of a piece is Lambda0(stop) - Lambda0(entry), so a subject split into pieces with constant covariates has exactly the likelihood of the unsplit record. Build the pieces from the moment each covariate changes, and never freeze a covariate that is only known later: an indicator of "ever had the disease" set from the start gives the survivors the exposure for time in which they could not yet have had it (immortal time), and it can reverse the sign of the effect. A first entry above zero is left truncation, and the message says that with a frailty the conditional likelihood it uses is naive (van den Berg and Drepper 2016).

WHAT IS ESTIMATED AND HOW. The fixed effects and the log-frailties always. The Weibull rho joins the Newton system as one more coordinate unless it is given. lambda is the same parameter as the intercept (intercept = rho * log(lambda), Eqn 16.3): estimating it adds the intercept column, giving it fixes the baseline and removes the intercept, so the reference classes sit at risk 1 – the convention of Example 16.1, which fixes rho = 1 and lambda = 1. The frailty components are estimated as the MAXIMUM of the Laplace approximation of the marginal likelihood (the frailties are integrated at their conditional mode, the fixed effects and rho are profiled), the very marginal_loglik the fit reports, because the marginal has no closed form here and the package's AI-REML machinery needs a Gaussian residual this model does not have: stats::optimize() in the log of the variance when there is one component, Nelder-Mead and Newton steps on a numerical Hessian (log variances, inverse hyperbolic tangent of the correlations) when there are several. As in model_threshold(), a correlation at or near +-1 is profiled instead of given a delta-method standard error, a variance the data do not tell apart from zero (the -2logL at 1e-6 within 0.01 of its value at the estimate) is held at its estimate without a standard error, a curvature that is singular in some combination of the components withholds the standard errors that depend on it and keeps the others, an inner fit that does not converge counts as an inadmissible point, and an indirect() term whose pens all hold one animal is refused. Give sigma2= to skip the estimation (the book's route: its example takes the variance as known). There is no residual variance, so neither h2() nor t2() has a denominator here. MEASURED BIAS, at prototype scale only (2329 records, one per animal, 40% censored, 8 replicates; validation/indirect_threshold_survival_recovery.R): the direct frailty variance ran low, 0.195 (SE 0.016) against 0.25 without indirect() and 0.207 (SE 0.015) with it, so the indirect term added no bias of its own (paired +0.012, SE 0.014).

SEVERAL FRAILTY TERMS AND INDIRECT EFFECTS. The linear predictor takes any number of random terms, as in model(): a sire and a herd-year frailty, or indirect genetic effects on the log-hazard of animals kept in groups, animal(id, group = "g") + indirect(id, pen = "pen", group = "g", dilution = d), where record i receives (n_i - 1)^(-d) sum_j a_S,j over its distinct pen mates and the direct and indirect frailties share one 2x2 matrix. (Ellen et al. 2010 reached this model in two steps, a survival analysis and then a linear associative model, because the survival software of the time took no associative effect.) The pen is labelled as the engine labels it, and a pen that is NA, NaN, Inf, blank text or the text "NaN" is refused with its row. Add random(pen) for the frailty shared by the pen; left out, it is absorbed by the indirect variance. A frailty other than indirect() belongs to the subject, so with elementary records an animal that moves to another pen cannot carry a random(pen) frailty (refused), while its indirect effect follows it into both pens. With elementary records, mates = "all" (the default) keeps every animal that was ever in the pen as a mate on every piece, culled and dead ones included, with the dilution of the whole pen. That is the choice the literature supports: Ask et al. (2020; pigs, daily gain) found that omitting culled animals, or weighting their indirect effects by the time they spent in the pen, reduced predictive ability, and Brinker et al. (2015; laying hens) found that making the indirect effect of a mate stop at its death reduced EBV accuracy compared with an analysis of survival time. mates = "present" is that time-dependent alternative: each piece keeps only the mates with a piece of their own that overlaps it, and the dilution uses their number (this weight is the package's choice, the group-size dilution of Bijma 2010 applied to the mates present; it is not taken from Brinker et al.). With it, cut the records where the composition of the pen changes, so that each piece has one set of mates: a mate whose last piece ends at t is absent from a piece that starts at t, since pieces are (entry, stop].

NO IMPLICIT INTERCEPT unless lambda is estimated, and the FIRST level of every fixed class factor is dropped (solution zero), as the book does. Levels come in factor order when the column is a factor, alphabetical order otherwise. Solutions read as log relative risks: exp(b) is the risk ratio RRS against the reference level, and exp(a) the frailty of the animal – POSITIVE means MORE risk of failure, so a good animal has a negative solution.

The structure checks of model() (see its section "Structure checks") apply here, but for the warning about repeated records: the pieces of one animal are not repeated records. Without sigma2 = a design that does not identify a frailty component stops, with it the design is a warning. With lambda given there is no intercept and the first level of each fixed class effect is held at the baseline: a frailty term nested in a fixed effect is then informed by the single level of the term in that reference level alone, which does not determine its variance, and the fit goes on with a warning that says so.

References

Kachman, S.D. (1999) Applications in survival analysis. J. Anim. Sci. 77 (suppl. 2), 147-153. Ducrocq, V. (1997) Survival analysis, a statistical tool for longevity data. 48th Annual Meeting of the EAAP, Vienna. Mrode, R.A. & Pocrnic, I. (2023) Linear Models for the Prediction of the Genetic Merit of Animals, 4th ed., chapter 16. Ellen, E.D., Ducrocq, V., Ducro, B.J., Veerkamp, R.F. & Bijma, P. (2010) Genetic parameters for social effects on survival in cannibalistic layers. Genet. Sel. Evol. 42, 27. Brinker, T., Ellen, E.D., Veerkamp, R.F. & Bijma, P. (2015) Predicting direct and indirect breeding values for survival time in laying hens using repeated measures. Genet. Sel. Evol. 47, 75. Ask, B., Christensen, O.F., Heidaritabar, M., Madsen, P. & Nielsen, H.M. (2020) The predictive ability of indirect genetic models is reduced when culled animals are omitted from the data. Genet. Sel. Evol. 52, 8. Bijma, P. (2010) Multilevel selection 4: modeling the relationship of indirect genetic effects and group size. Genetics 186, 1029-1031. Heinze, G. & Schemper, M. (2001) A solution to the problem of monotone likelihood in Cox regression. Biometrics 57, 114-119.

Examples

# Example 16.1 of Mrode & Pocrnic (data on p.284, solutions on p.291): length of
# productive life of 12 cows in 2 herds, 4 of them censored, animal frailty with
# the full 19-animal pedigree, rho = 1, lambda = 1 and the frailty variance given.
cows <- data.frame(
  cow  = as.character(8:19),
  herd = as.character(c(1,1,1,1,1,1,2,2,2,2,2,2)),
  ysp  = as.character(c(3,4,1,2,3,1,4,1,2,3,4,2)),
  code = c(0,1,0,1,1,1,1,1,0,1,0,1),
  lpl  = c(40,47,22,28,50,33,49,29,23,37,35,30))
ped <- data.frame(
  animal = as.character(1:19),
  sire = c(rep("0",7), c("1","1","4","4","5","5","1","1","5","5","4","4")),
  dam  = c(rep("0",7), c("2","3","2","9","3","8","6","7","14","6","7","3")))
fit <- model_survival(lpl ~ herd + ysp + animal(cow), data = cows, pedigree = ped,
                      censor = "code", rho = 1, lambda = 1, sigma2 = 0.4,
                      verbose = FALSE)
coef(fit, "fixed")       # herd=2 -1.631; ysp 2:4 -2.346 -3.149 -2.982 (p.291)
#>    herd=2     ysp=2     ysp=3     ysp=4 
#> -1.631055 -2.346427 -3.149264 -2.982340 
exp(ebv(fit)[["1"]])     # the RRS of animal 1, 0.459 on p.291
#> [1] 0.458965
predict(fit, data.frame(herd = "1", ysp = "4", cow = "1"), time = 40,
        type = "survival")   # 0.394 on p.292
#> [1] 0.3944331