Multivariate Survival Mixed Models for Genetic Analysis of Longevity Traits

Longevity is an important trait often considered in animal breeding programs [31, 30, 33, 37, 27, 35, 17, 9, 6, 32, 1]. Even small changes in the longevity of a population under production might have remarkable economic, welfare and ethics consequences [30, 12]. Since the study of longevity involves several types of incomplete observation (e.g. censoring, truncation, late entry and competing risks), survival and event-history-analysis techniques are typically used [8, 19, 3]. However, the use of those techniques in the context of quantitative genetics of longevity involves several non-trivial challenges. We will present a model framework here that allows to overcome these challenges. These models will have a structure of means and covariances similar to the gaussian linear mixed models classically used in quantitative genetics. They will allow for a proper representation of quantitative genetic phenomena and for efficient implementations required in practice. We will show that these models extend the class of models currently in use for studying the genetic of longevity.

The first challenge in the use of survival models for characterizing the genetic aspects of longevity is the high complexity of the models. Typically it is necessary to adjust simultaneously for the effects of several explanatory variables, some continuous with linear effects (e.g. breed composition and heterosis [21]), some factors with many levels (e.g. herd) or even time-dependent effects (e.g. year and season). Further complexity is added by the necessity of representing complex genetic effects (several generations deep pedigrees) under different genetic models (e.g. sire model, sire-dam-model etc). A typical scenario of applications in quantitative genetics involves a very large number of observations (individuals); indeed often the analyses involve several hundred thousand of observations (we will present two relative small illustrative examples involving 142,133 and 200,084 individuals). This leads to inference problems of the complexity equivalent to solve systems of linear equations with several hundred of thousand, or even several million, simultaneous linear equations (in our examples between 8,248 and 118,432 simultaneous linear equations) when using classic gaussian linear mixed models. Now, the complexity of survival and event-history models is even larger since these models require to keep track of the individuals’ history. Moreover, as it will be apparent from the description of the methods available and from our discussion, some of those models might present very flat likelihood functions making the inference problem hard and numerically unstable. This rules out the naive use of the standard methods of survival analysis and makes the use of specially constructed models and approximations mandatory.

A second challenge is the proper representation of the genetic phenomena in play. The interpretation of the gaussian linear mixed models in genetic terms requires a linear decomposition of the so called total variance (i.e. the phenotypic variance after having corrected for the effects of a range of non-genetic related aspects with known effects on the trait of interest)[18, 34]. This decomposition of the total variance occurs naturally in the gaussian linear mixed models but not in all the adaptations of survival models specially constructed for genetic evaluations of longevity. We will show that a way to circumvent this problem is to consider models containing a dispersion parameter which plays the role of the residual variance (i.e. the part of the total variance that is neither attributed to genetic components nor to identifiable sources of variation) and base the inference on the quasi-likelihood theory [36, 4]. The use of generalised linear mixed models with a dispersion parameter occurs only occasionally in the literature of quantitative genetics [14, 5], and the genetics and the statistics consequences of using a free dispersion parameter has not being systematically explored yet. Moreover, survival models with dispersion parameter were, to the best of our knowledge, never considered in the literature of quantitative genetics.

A third challenge is the presence of incomplete observations following complex patterns. Variants of the so called Cox proportional model [8] were used for univariate traits describing longevity [9, 13, 10]. These techniques succeeded to implement models of approximations that were operational for some of the purposes of animal breeding. However, all those models were of univariate nature (i.e. they consider one trait at time) [30, 11], which makes it difficult, if not impossible, to address some key questions related to the issues of competing risks, informative censoring and inhomogeneity of the failure mechanism. A clear example is the study of time to death were the individuals on the study may die from one among several different causes of death. Moreover, the use of multivariate models is necessary to study the time death in a competing risk scenario, by evaluating the risk of death for each specific cause simultaneously. This allows to form the basis to study further aspects as the presence of common genetic determining factors and/or independent or specific genetic factors for each of the death causes. In the framework presented here, it will be possible to treat these issues by using suitable multivariate models.

Two types of models will be considered: models based on continuous time and models based on discrete time. The continuous time models (CTM) will be suitable approximations of the semi-parametric Cox proportional model in which the hazard function will be assumed to be piece-wise constant. These models allow for regular parametric likelihood-based inference by exploring a coincidence of their likelihood functions and the likelihood function of a Poisson model applied to a specially constructed data [2]. When including gaussian random components (log gaussian frailties) the inference can be based on a Laplace approximation to a high dimensional integral [4, 24, 28]. Here, we will include also a dispersion parameter that will play the role of the residual variance which is an essential element for the genetic interpretation of the models. The CTM will depend on an arbitrary choice of time cut points used to define the intervals at which the hazard function is assumed to be constant.

The discrete time models (DTM) used are multivariate variants of the classic discrete relative risk models [19]. These models present advantages in terms of computational and statistical complexity as compared to the CTM: the algorithms run faster relatively to the analogous algorithm for CTMs, are numerically more stable, and the statistical inference is more efficient. The loss of genetic information occurred when switching from the CTM to the DTM is minimal in typical applications in quantitative genetics and this loss is fairly compensated by the avoidance of capturing a large amount of noise. Here, we will also introduce the use of a dispersion parameter, which will be essential for well representing genetic scenarios.

The aim of this article is to introduce, characterize and discuss a class of multivariate mixed survival models for continuous and discrete time with a complex covariance structure in a context of quantitative genetics applications. The methods introduced here can be used in many other applications in quantitative genetics although the discussion presented concentrates on longevity.

The paper is organized as follows. Section 2 describes the basic set-up and the genetic scenario discussed. There, a suitable multivariate version of the proportional hazard model is introduced in general terms. Those models will encompass models for competing risks possibly defined with different types of time scale (continuous and discrete time). The techniques for statistical inference under those models are presented in section 3, including a general discussion on the calculations involved in the likelihood based inference and some connections with multivariate generalized linear mixed models (section 3.1). A Poisson approximation useful for efficiently implementing likelihood based inference for continuous time models is presented in section 3.2. Section 3.3 presents an extension of the models to a situation where there is a stratification variable. A discussion comparing the amount of statistical information between models constructed with discrete time and models based on continuous time is given in section 3.4. The quantitative genetic theory behind the models considered here requires a special decomposition of the phenotypic variance in terms of the variance of the random components and a scale parameter present in the models. This will be necessary for the calculation of the so called heritability, which is crucial in quantitative genetic applications. Section 4 will discuss methods for that and some technical details involving counting processes and martingale theory are presented in the appendix. Two illustrative examples involving longevity of sows and dairy cattle will be presented in section 5.1 and 5.2 respectively. Some discussion will be given in section 6.

The class of models we will discuss are thought to be applied in the following general scenario. The life spans of a population of $n$ individuals are observed. Genetic information (typically in the form of a reasonably deep pedigree or comprehensive genomic data) and the values of $k$ explanatory variables are available for each individual. The interest lies in modeling the longevity, i.e. the length of the life span or the length of the productive life of the individuals, with particular interest on some forms of genetic determination typically used in quantitative genetics.

The longevity can be operationally measured using a continuous time, by determining the time elapsed between two life events in the life of the individuals (e.g. the time elapsed from the first parity the to the culling of cows under production) or a discrete time, by counting a number of certain events during the life of the individual (e.g. by counting the number of survived parities of cows). We will present two illustrative examples: one involving the longevity of sows and the other studying the longevity of dairy cattle. Many similar examples can be found in the literature with other husbandry animals as ewes [25], salmon [20] among others.

An additional issue that might occur in longevity studies is the presence of competing risks, which arises when the individuals could be culled for one of two or more types of culling. For instance, in the example of dairy cattle longevity, the cows might die or be slaughtered, and the interest is in studying the culling rates for both causes. We assume in this scenario that there are two causes of death or culling, which will be labeled by the index $j$. It is straightforward to extend the setup described here to the general case with more than two culling causes or to the case with only one culling type.

The longevity will be characterized by the time development of the rate at which the individuals die or are culled. This is measured differently when using continuous or discrete time as we define below. In order to describe the models we have in mind, consider two random variables: $T$ representing the time (continuous or discrete) of culling and $J$ indicating the cause (or type) of culling. In the case of the continuous time, we define the cause-specific hazard function [19] for the $j^{\rm th}$ culling cause (for $j=1,2$) by

In contrast, when the time is discrete, we define the cause-specific hazard probability function [19] for the $j^{\rm th}$ culling cause (for $j=1,2$) by

Although these two characterisations of the time development of the culling rate are of different nature, they will essentially play the same rule in the models discussed and will be generically referred simply as the cause-specific hazard function. Note that in the case where only one cause of culling is studied the cause-specific hazard function and the cause-specific hazard probability function defined above coincide with the hazard function and the hazard probability function used in the literature of survival analysis (see [19]).

The models presented below will assume a scenario where there are two causes of culling (indexed by $j$). Here, the cause-specific hazard functions (or the hazard probability functions) are specified in terms of $k$ explanatory variables (the fixed effects) and a set of gaussian random components. Without loss of generality, we consider only two random components for each culling cause: $\mbox{\bf U}_{1}$ and $\mbox{\bf U}_{2}$ that will represent additive genetic effects (usually determined using information on the pedigree of the animals in play) for cause $1$ and $2$ respectively and $\mbox{\bf V}_{1}$ and $\mbox{\bf V}_{2}$, representing an environment effect (e.g. the effect of herd where an animal is kept) for cause $1$ and $2$, respectively. The random components related to genetics (i.e. $\mbox{\bf U}_{1}$ and $\mbox{\bf U}_{2}$) are independent of the components related to environment ($\mbox{\bf V}_{1}$ and $\mbox{\bf V}_{2}$) by construction. This constraint can be relaxed for describing possible interactions between environmental and genetic effects, but we will not pursue this project here. The extension of the present definition to a situation with a different number of random components is straightforward. According to the general model we describe, the vector of cause specific hazard functions of the $i^{\mbox{th}}$ individual (for $i=1,\dots,n$), conditionally on $\mbox{\bf U}_{1}=\mbox{\bf u}_{1}$, $\mbox{\bf U}_{2}=\mbox{\bf u}_{2}$, $\mbox{\bf V}_{1}=\mbox{\bf v}_{1}$ and $\mbox{\bf V}_{2}=\mbox{\bf v}_{2}$ (for the realizations $\mbox{\bf u}_{1}$ , $\mbox{\bf u}_{2}$, $\mbox{\bf v}_{1}$ and $\mbox{\bf v}_{2}$ of $\mbox{\bf U}_{1}$ , $\mbox{\bf U}_{2}$, $\mbox{\bf V}_{2}$ and $\mbox{\bf V}_{2}$, respectively), is given by

Here the times $t_{1}$ and $t_{2}$ might be continuous or discrete; moreover, $t_{1}$ and $t_{2}$ are not necessarily of the same type (as in the example of sows longevity considered in section 5.1). Equation (1) holds for $t_{j}\geq 0$ when $t_{j}$ is continuous and for $t_{j}=0,1,2,\dots$ if $t_{j}$ is discrete ($j=1,2$). Furthermore, $\lambda_{1}(\cdot)$ and $\lambda_{2}(\cdot)$ are baseline hazard functions describing a common time development of the cause specific hazard functions (typically considered as a nuisance parameter), $\mbox{\boldmath$\beta$}_{1},\mbox{\boldmath$\beta$}_{2}\in\mathbb{R}^{k}$ are (finite dimensional) parameters representing the fixed effects (also viewed as a nuisance parameter), $\mbox{\bf X}_{i}$, $\mbox{\bf Z}_{ji}$ and $\mbox{\bf W}_{ji}$ are incidence matrices for the fixed, genetic and environment effects for the $j^{\mbox{th}}$ cause and the $i^{\rm th}$ individual. When the time is continuous $\lambda_{[j],i}(\cdot)$ in (1) is a hazard function and the marginal model described will correspond to a variant of the Cox proportional hazard model with frailties [19, 3]. If the time is discrete, $\lambda_{[1],i}$ in (1) is a hazard probability function and the marginal model described will be a variant of the discrete relative risk model with frailties [19].

In the case of the models based on continuous time, we assume that the baseline function, $\lambda_{j}(\cdot)$ is piece-wise constant. That is, we assume that there is a set of disjoint intervals, $I_{1},\dots,I_{K}$ covering the positive real line, in which $\lambda_{j}(\cdot)$ is constant. This constraint in the class of possible baseline functions will allow us to connect the models described here with some generalized linear mixed models: Moreover, this will allow us to perform efficient inference with complex data using existent software.

The structure of covariance of the random components is a crucial part of the construction of the models for longevity described since this is the part of the model where the additive genetics and the environment contributions to the cause specific hazard functions are considered. We assume the random components to be multivariate normally distributed, i.e. $(\mbox{\bf U}_{1},\;\mbox{\bf U}_{2},\;\mbox{\bf V}_{1},\;\mbox{\bf V}_{2})^{{}^{\prime}}\sim\mathcal{N}\left({\bf 0},{\bf\Sigma}\right)$ where

Here, ${\bf I}$ is an identity matrix and ${\bf A}$ is the $n\times n$ additive relationship matrix constructed as follows. Each entry of ${\bf A}$ represents a pair of individuals in the population. $A_{ii}=1$ for $i=1,\dots,n$, and $A_{ij}=0$ if the pair of individuals $(i,j)$ are not related. In the case were the pair of individuals $(i,j)$ are related, we consider the degree of relationship $r$ given by 1 for relatives of first-degree, 2 for relatives of second-degree, 3 for relatives of third-degree and so on. In that case $A_{ij}=(1/2)^{r}$. The additive relationship between two individuals measures the proportion of identical by descent genes expected to be shared by pairs of individuals in the population. The methods introduced here could also be applied to genomic analyses where the relationship matrices are inferred from genomic data.

The covariance structure between the genetic random components is given by

In particular the genetic correlation between the two culling causes is given by $\rho_{g}=\frac{\sigma_{g.12}}{\sqrt{\sigma_{g.1}^{2}\sigma_{g.2}^{2}}}$. Finally, the covariance structure between the environment random components is given by

The notion of heritability is crucial in quantitative genetic applications; it measures the magnitude of the detected additive genetic signal relative to the total variance of a trait of interest and is typically used to quantify the potential response to selection. Heritability (in the narrow sense) is defined operationally by the ratio

where $\sigma_{g}^{2}$ is the variance of an unobserved random component representing all the additive genetic variation (termed the genetic variance) and $\sigma_{p}^{2}$ is the variance of the trait of interest (called the phenotypic variance). For details see [22, 18]. In the classic scenario, where gaussian linear mixed models are used, the calculation of the heritability is straightforward because the phenotypic variance can be expressed as the sum of the genetic variance $\sigma_{g}^{2}$ and the variances of the other random components present in the model. However, the calculation of the heritability is much more complicated in the context of survival analysis considered here and are, therefore, developed in full details below. It will be necessary to introduce some details of the counting process behind the models used and to define the trait of interest carefully. Here, we present the main results informally and expose the precise formal definitions and the technical details using the mathematical machinery of counting processes in appendix A. For simplicity of exposition, we discuss only the one-dimensional case were only one cause of culling is present (referred as death) and assume the same model structure presented in section 3 but with dimension one.

Let $T$ be a non-negative random variable representing the observed survival time and $D$ be an indicator variable taking the value 1 if a death is observed and 0 otherwise (censure). The random variable $T$ will take values in $\mathbb{R}_{+}$ if the time is continuous or in $\mathbb{Z}_{+}$ (i.e. in $\{0,1,2,\dots\}$) if the time is discrete. We assume that for each time $t$ a vector $X(t)$ of explanatory variables (possibly changing with time) is known. Denote by $\mbox{\bf X}(t)=\{X(s):0\leq s<t\}$ the trajectory (or path) function until the time just before $t$ (denoted by $t-$). Additionally, we consider two independent gaussian random components, U and V, representing the additive genetic effects and some environmental effects, respectively. Here, we follow the same convention as in section 2, so $Var(\mbox{\bf U})=\mbox{\bf A}\sigma_{g}^{2}$ and $Var(\mbox{\bf V})=\mbox{\bf I}\sigma_{e}^{2}$, where A is the relationship matrix and I is an identity matrix with suitable dimensions. The hazard function for the $i^{\rm th}$ individual takes then the form

In order to describe the models presented in section 3 properly, we need to introduce two stochastic processes: the death (or event) counting process and the at risk process. To define the death process we consider, for each $t$, the random variable $N(t)$, which indicates whether a death occurred until the time $t$, i.e. $N(t)=\mathbf{1}(T\leq t,D=1)$. The stochastic process $\mathbf{N}=\left\{N(t):t\in\mathbb{R}_{+}\mbox{ or }\mathbb{Z}_{+}\right\}$ is the death counting process. It is convenient to introduce the notation $dN(t)$ to denote the increments of the counting process $\mathbf{N}$ at each time $t$, i.e. $dN(t)=N(t)-N(t-)$, where $N(t-)=\lim_{\Delta\downarrow 0}N(t-\Delta)$. Furthermore, the at risk process given by $\mathbf{Y}=\left\{Y(t):t\in\mathbb{R}_{+}\mbox{ or }\mathbb{Z}_{+}\right\}$, where, for each time $t$, the random variable $Y(t)$ takes the value $1$ when the individual is at risk at time $t$ or $0$ otherwise.

We will calculate the phenotypic variance and the heritability for two characteristics of interest: the hazard function evaluated at a given time and the accumulated hazard function up to a given time. These two characteristics are associated to the intensity and the cumulative intensity of the counting process of death $\mathbf{N}$. Due to the structure of the models used here (see equation (1) and (11) for instance), the entire hazard function of an individual is multiplied by the factor $\exp(\mbox{\bf u})$, where u is the realization of the genetic random component U for that individual; so, according to this model, when selecting individuals by the predicted values of U (as it is the usual practice in animal breeding), the entire hazard curve is affected and therefore, it is indifferent if one measures the potential response to selection (i.e. the heritability) of the hazard evaluated at a given time-point, say $t$, or even of the accumulated hazard up to $t$.

In order to study the heritability related to the hazard function evaluated at a given time $t$ (in $\mathbb{R}_{+}$ or in $\mathbb{Z}_{+}$), we define the conditional random variable

Note that $E\left({\cal{Y}}(t)\,|\,\mbox{\bf U}=\mbox{\bf u},\,\mbox{\bf V}=\mbox{\bf v}\right)=\lambda\left(t\,|\,\mbox{\bf U}=\mbox{\bf u},\,\mbox{\bf V}=\mbox{\bf v}\right)$, therefore ${\cal{Y}}(t)$ will be treated as the trait of interest, for which we calculate the phenotypic variance and the heritability. A more precise definition of ${\cal{Y}}(t)$ using counting processes theory is given in appendix A. There we show, taking advantage of the basic theory of martingales and using a suitable Taylor expansion, that in the continuous time case the variance of ${\cal{Y}}(t)$ is approximated by

where $\lambda^{*}(t)$ is the hazard function evaluated at $t$ with vanishing random components, i.e.

Note that according to (13) the variance of ${\cal{Y}}(t)$ is additively decomposed in three components: one depending on $\sigma^{2}_{g}$, one depending on $\sigma^{2}_{e}$ and one depending on $\phi$. This is a situation analogous to the classic gaussian linear mixed models. Here, the component of the variance of ${\cal{Y}}(t)$ associated to $\sigma^{2}_{g}$, namely $\left[Y(t)\lambda^{*}(t)\right]^{2}\sigma_{g}^{2}$ , plays the role of the genetic variance and therefore the heritability for the hazard evaluated at the time $t$ is given by

Analogous calculations (see appendix A) yield for the discrete time case that the variance of ${\cal{Y}}(t)$ is approximated by

Here, as in the case of the continuous time, the variance of ${\cal{Y}}(t)$ is additively decomposed in three separate a components depending on $\sigma^{2}_{g}$, on $\sigma^{2}_{e}$ and on $\phi$, respectively. Therefore, the heritability is approximated by

Moreover, when a Poisson approximation is used, the approximated heritability takes the form given in (15).

We turn now to the calculation of the heritability of the accumulated hazard function evaluated at a time point $t$. In the continuous time case we define the accumulated hazard function, for each $t\in\mathbb{R}_{+}$, by

and in the discrete time case the accumulated hazard function is given, for each $t\in\mathbb{Z}_{+}$, by

where $\lambda(\cdot)$ is the hazard function in play. Furthermore, the conditional random variable

is such that $E\left({\cal{Z}}(t)\,|\,\mbox{\bf U}=\mbox{\bf u},\,\mbox{\bf V}=\mbox{\bf v}\right)=\Lambda\left(t\,|\,\mbox{\bf U}=\mbox{\bf u},\,\mbox{\bf V}=\mbox{\bf v}\right)$. Therefore, ${\cal{Z}}(t)$ will be treated as the trait of interest, for which we calculate the phenotypic variance and the heritability. A more precise definition of ${\cal{Z}}(t)$ using counting processes theory is given in appendix A. It can be shown (see appendix A) that in the case of the continuous time the variance of ${\cal{Z}}(\cdot)$ evaluated at a given $t\in\mathbb{R}_{+}$ is approximated by

where $\Lambda^{*}(t)=\int_{0}^{t}Y(s)\lambda^{*}(s)ds$. Therefore, the heritability is approximated by

In the case of the discrete time the variance of ${\cal{Z}}(\cdot)$ evaluated at a given $t\in\mathbb{Z}_{+}$ is approximated by

where $\gamma(t)=\sum_{i=0}^{t}{\phi Y(t)\lambda^{*}(t)[1-\lambda^{*}(t)]}$ and the cumulative hazard is $\Lambda^{*}(t)=\sum_{s=0}^{t}Y(s){\lambda^{*}(s)}$ .Then, the heritability is approximated by

When the Poisson approximation is used the heritability is approximated by the right side of (20).

The calculations above are performed for one individual. We propose that given a population of $n$ individuals, the phenotypic variances and the heritabilities for the hazard function or the cumulative hazard (evaluated at a time $t$) should be calculated using the estimated variance components and an estimated hazard or an estimated cumulative hazard as described below. In the continuous time case, where the PCHM is used, the observed time is split in intervals $I_{k}=[t_{k},t_{k+1})$, for $k=0,\dots,K$ where $0=t_{0}<t_{1}<\dots<t_{K}<t_{K+1}=\tau$, where $\tau$ is the largest observed time. Given that $\mbox{\bf U}=\hat{\mbox{\bf u}}$ and $\mbox{\bf V}=\hat{\mbox{\bf v}}$, where $\hat{\mbox{\bf u}}$ and $\hat{\mbox{\bf v}}$ are the BLUP of the random components, the conditional hazard is defined, for each $t\in I_{k}$, by

where $\eta_{i,k}(t)=\log(\lambda_{k})+\mbox{\bf X}_{i,k}^{{}^{\prime}}\mbox{\boldmath$\beta$}+\mbox{\bf Z}_{i}^{{}^{\prime}}\hat{\mbox{\bf u}}+\mbox{\bf W}_{i}^{{}^{\prime}}\hat{\mbox{\bf v}}$ and $\mbox{\bf X}_{i,k}=\mbox{\bf X}_{i}(t)$ for $t\in I_{k}$. We obtain $\hat{\eta}_{i,k}(t)=\widehat{\log(\lambda_{k})}+\mbox{\bf X}_{i,k}^{{}^{\prime}}\hat{\mbox{\boldmath$\beta$}}+\mbox{\bf Z}_{i}^{{}^{\prime}}\hat{\mbox{\bf u}}+\mbox{\bf W}_{i}^{{}^{\prime}}\hat{\mbox{\bf v}}$ from the adjusted model and then calculate, for each time $t\in I_{k}$,

where $\bar{\eta}_{k}=n^{-1}\sum_{i=1}^{n}{\hat{\eta}_{i,k}}$. The cumulative hazard is estimated, for $t\in I_{k}$, by

where $\Delta_{j}=t_{j+1}-t_{j}$ is the length of the interval $I_{j}$.

In the discrete time case the hazard, conditionally that $\mbox{\bf U}=\hat{\mbox{\bf u}}$ and $\mbox{\bf V}=\hat{\mbox{\bf v}}$, is given, for each $t\in\mathbb{Z}_{+}$, by

where $\eta_{i}(t)=\log(\lambda_{t})+\mbox{\bf X}_{i}^{{}^{\prime}}(t)\mbox{\boldmath$\beta$}+\mbox{\bf Z}_{i}^{{}^{\prime}}\hat{\mbox{\bf u}}+\mbox{\bf W}_{i}^{{}^{\prime}}\hat{\mbox{\bf v}}$. We obtain $\hat{\eta}_{i}(t)=\widehat{\log(\lambda_{t})}+\mbox{\bf X}_{i}^{{}^{\prime}}(t)\hat{\mbox{\boldmath$\beta$}}+\mbox{\bf Z}_{i}^{{}^{\prime}}\hat{\mbox{\bf u}}+\mbox{\bf W}_{i}^{{}^{\prime}}\hat{\mbox{\bf v}}$ from the adjusted model and then calculate for each time $t$ ($t=0,1,2,\dots,$)

where $\bar{\eta}(t)=n^{-1}\sum_{i=1}^{n}{\hat{\eta}_{i}(t)}$. The cumulative hazard is estimated by $\tilde{\Lambda}^{*}(t)=\sum_{j=0}^{t}\tilde{\lambda}^{*}(j)$. The heritabilities will be estimated at $\tilde{\lambda}^{*}(t_{m})$ and $\tilde{\Lambda}^{*}(t_{m})$ where $t_{m}$ represents the median survival time, i.e. $P(T\geq t_{m})=0.5$ in the illustrative examples.

The methods presented above will be illustrated by two examples: one involving the longevity of sows (section 5.1) and another characterizing the longevity of dairy cattle (section 5.2). In both examples, we will characterize the longevity of female animals in two different ways: first, by the length of the productive life of the animals (sows or cows) measured by the number of days from the first parity to the culling day (ND); second, by the number of survived parities (NP). Models based on continuous time and models based on discrete time were applied to study ND and NP, respectively.

Incomplete observations are present in both examples because some animals were still alive at the end of the study, some animals were moved, sold or exported during the study and due to intentional right truncation. In the study of sow longevity the proportion of censure was relatively low ($9.6\%$), while in the study of longevity of cows a much higher proportion of censure was observed ($23.5\%$). Furthermore, in the study of longevity of cows it was possible to distinguish two culling causes: slaughter ($67.7\%$ of the observations) and death (only $8.8\%$ of the observations), which characterizes a typical scenario of competing risks. Therefore, it is natural to use bivariate models in this example to model these two culling causes simultaneously. Bivariate models will be used in the example of sows longevity to model ND and NP simultaneously and use that to discuss how much the information of these two characterizations of the longevity overlap.

The several variants of models for studying longevity in quantitative genetics can be described using the following general approach. The models are one-dimensional and describe the conditional hazard function, given the random components, as taking the form

for $t$ in $\mathbb{R}_{+}$ or in $\mathbb{Z}_{+}$. The models based on continuous time used in the literature differ essentially in the form of the baseline hazard function $\lambda_{0}(\cdot)$: the classical Cox model (with frailties) assumes the hazard function to be a smooth function (almost everywhere) imposing in this way almost no restriction on $\lambda_{0}(\cdot)$ (see [19, 26, 28]); the piece-wise Weibull models (see [9]) assume $\lambda_{0}(\cdot)$ to be a saw function (i.e. a piecewise linear function); while the piece-wise constant hazard models (as described here) assume $\lambda_{0}(\cdot)$ to be a step function (i.e. a piece-wise constant function). These models are nested, in the sense that saw functions are smooth almost everywhere, and step functions are particular cases of saw functions. Moreover, smooth functions can be arbitrarily approximated (point-wisely) by saw functions or by step functions. A fourth category of models, which was introduced here, are the piece-wise constant hazard models with free dispersion parameter. These models are an instance of quasi-likelihood based models (see [36, 4]). They contain the traditional piece-wise constant hazard models as a particular case (by setting the dispersion parameter $\phi$ to be constant equal to 1) and might arbitrarily approximate the piece-wise Weibull models and the Cox frailty models. Although in many practical situations those models yield equivalent representations of the data and produce similar results, this is not always the case for the quantitative genetics applications discussed here.

One of the interesting features of the piece-wise constant hazard model with (free) dispersion parameter presented here is that the phenotypic variance can be approximately decomposed as the sum of parcels representing the additive genetic effects, environmental effects and unspecified sources of variability. This is analogous to the classic decomposition of the phenotypic variance obtained in gaussian linear mixed genetic models. Here, the component of the phenotypic variance involving the dispersion parameter $\phi$ plays the rule of the residual variance in the classic models. Note that, if we fix the dispersion parameter $\phi$, say by setting it equal to 1 as in the piece-wise constant hazard model, it might be that we obtain models that describe the data well and the referred decomposition phenotypic variance would still exists, but this would be equivalent to set the residual variance to be constant, which clearly causes problems in the representation of quantitative genetic phenomena. For instance, in a sire model the genetic variance transmitted by the dam is not captured by the random component representing the sire, therefore, this part of the variability should be represented in the residual variance (i.e. should compose part of the residual variance); but this is not possible to occur if the residual variance is set to be constant. This issue occurs also in the piece-wise Weibull models (see [9]). Furthermore, in the most general case where $\lambda_{0}$ is only assumed to be smooth (the Cox frailty model) it might be argued that $\lambda_{0}$ could be point wisely approximated (almost every where) by a sequence of step functions obtained even from distributions in a model with dispersion parameter, which would have a proper decomposition of the phenotypic variance. However, in this case it would not be possible to identify the dispersion parameter $\phi$ (as essentially argued in [13]); therefore, one would still have the representation issues referred above.

The piece-wise constant hazard model and the piece-wise Weibull model are parametric models, making it possible to implement relatively efficient inference techniques. This explains why they are implementable for complex and large models as in quantitative genetics typical applications. Moreover, when the baseline hazard function $\lambda_{0}(\cdot)$ is assumed to be piece-wise linear (a saw function), as in the piece-wise Weilbull models, the rate of approximation of $\lambda_{0}(\cdot)$ to a smooth function (as in the frailty model) is improved (specially if there are regions where the smooth baseline is steep) as well argued in [9]; however, in this case the inference via generalized linear mixed models would not be possible without using profile likelihood techniques (to estimate the shape parameter), which rules out the direct use of extensions containing a dispersion parameter or multivariate extensions already implemented. Finally, a disadvantage of the piece-wise constant hazard models and the piece-wise Weibull models is that both depend on an arbitrary choice of time cutting points.

Models based on discrete time (i.e. variants of the discrete relative risk models) has been occasionally used in quantitative genetic previously; however, without incorporating a dispersion parameter. These models do not assume any pre-specified form for the baseline hazard function and also not require any arbitrary definition of time cutting points; therefore, they are similar to the Cox frailty models. Due to the coincidence with the Bernoulli models with random components (with or without dispersion parameter), techniques for parameter inference are at hand, allowing for efficient implementations capable to handle complex large problems as illustrated here. Moreover, these models allow to incorporate time-dependent variables (by splitting data with record for each survived time for each individual) both as fixed effects and as random components. The Poisson approximation presented here might represent a substantial save in computational resources and yields essentially the same results as the exact inference, as illustrated in the example of sows longevity presented in section 5.1.

The continuous time models and the discrete time models yielded similar results in the example of sows longevity. However, the continuous time models failed in the example of dairy cattle longevity (section 5.2), which can be explained by the fact that, in this case, both rare events and high frequencies of censored observations are observed. This is in line with the limitations of the continuous time models discussed in section 3.4 and shows that those limitations do occur in practical situations.

The lack of multivariate techniques for analyzing several traits simultaneously has been a serious limitation in the application of survival models in quantitative genetics [30, 11]. This lack not only makes the study of the relation between longevity and other traits impossible, but also limits very much the possibility of performing time-to-event-analysis in the presence of competing risks. Indeed, if there are two competing culling reasons, one might use standard survival models for modeling the rate of occurrence of one of the culling reasons by considering the observations where the other culling reason occurred as censored. However, this naive approach implies in the implicit use of the assumption that the competing culling reasons are independent. The example of dairy cattle longevity presented here illustrates that this assumption might be not reasonable. Indeed, in this example we detected a significant negative correlation between the herd-year-season random components; moreover, a range of fixed effects were found to be significant for both culling reasons, implying that a naive use of marginal models would violate the assumption of absence of informative censoring.

Rafael Pimentel Maia was financed by the project ”Svineavl: Developing New Methods for Genetic Selection of Sow Durability” , Ministry of Food, Agriculture and Fisheries of Denmark.