|
|
||||||||
ANIMAL GENETICS |
Grup de Recerca en Remugants, Departament de Ciència Animals i dels Aliments, Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Spain
| Abstract |
|---|
|
|
|---|
Key Words: Beef Cattle Heritability Survival Analysis Survival to Weaning
| Introduction |
|---|
|
|
|---|
Despite its economic importance for suckling cow production systems, phenotypic and genetic studies on calf survival at weaning are scarce and little information is available (Cundiff et al., 1986
; Ray et al., 1989
; Goyache et al., 2003
). Until now, calf survival has usually been studied as a dichotomous variable at a given age; however, the use of a continuous measure as failure time is preferable to analyzing it with regression models because it makes use of all the information available more efficiently by not restricting observations to an arbitrarily defined point (Ducrocq, 1997
). The problem of working with these continuous measures is the presence of live animals at the completion of the study that have inexact information of their failure time. Records from these animals are classified as censored. Survival analysis is the method of choice to analyze the failure time data allowing the inclusion of censored and uncensored (dead animals) records. The aim of this paper was to analyze calf survival using survival analysis techniques. The general strategy for the analysis of survival data in beef calves was as follows: description of the observed survival function by the Kaplan-Meier method, estimation of fixed effects with a semiparametric proportional hazards model and validation of the baseline hazard function, and finally, the estimation of the sire variance component for hazard and the heritability of survival in the binomial scale.
| Materials and Methods |
|---|
|
|
|---|
Data were recorded between 1994 and 2002 in three breeding herds participating in the Yield Recording Scheme of the breed, and included data from 3,335 beef calves. The pedigree of the calves was only partly known, as 86 bulls sired 2,578 calves, with the rest of the sires being unknown. All dams were known. The survival time was estimated as the difference between the date of death and the date of birth. When the date of death was unknown, we assumed that the calf record was censored at the date of weaning (180 d), and any known dates of death after weaning were not used in the analysis. After editing and excluding twin births, our database included data of survival time of 2,504 calves, with 68 complete records (2.7%) and 2,436 censored records (97.3%).
Survival Analysis
Our survival analysis followed the general strategy proposed by Ducrocq (1997)
. In the preliminary analysis, the observed distribution of survival time was estimated empirically by the Kaplan-Meier method (Kaplan and Meier, 1958
). This observed distribution was the result of the influence of different factors involved in the survival process. To investigate the influence of these factors on the instantaneous death rate (h(t, X), the hazard), the following semiparametric proportional hazards model (Cox, 1972
) was assumed:
![]() |
The first term of the right hand side is the unknown baseline hazard function h0(t). The second term (between brackets) was dependent on a subset of explanatory variables, which we call survival factors. Herd x year of calving (HY), month of calving (MC), cows length of productive life at calving (LPL), calving difficulty (CD) and birth weight (BtW) were treated as time-independent variables and included sequentially from left to right in Models 1 to 4. Hypothesis testing was performed via likelihood ratio tests (LRT). The sequential test included the effects in the model in sequential order. The last test compared the full model with models excluding one effect at a time (this was done for each effect separately). These and the subsequent computations were performed by means of the Survival Kit package (Ducrocq and Sölkner, 1998
).
Because the form of the phenotypic relationship between survival factors and the hazard was not known, all variables included in the models were categorized (Table 1
). This way, no assumption about the nature of this relationship has to be made (Ducrocq, 1997
). The herd x year effect was treated as a fixed variable. The inclusion of an interaction term between herd and year of calving would permit precise modeling of changes in the mortality process within a farm over time, similar to that observed by Ducrocq (1997)
for longevity in dairy cattle. After a preliminary study, classes with a similar effect on the risk of mortality were grouped. Month of calving had two classes with a similar number of observations (September to February and March to August). Length of productive life of the cow at calving (i.e., the difference between the age of the cow at calving and the age at first calving) had two classes, 1 to 1,300 d (cows up to four calvings) and productive life longer than 1,300 d. Calving difficulty had four classes: calving without assistance, calving slightly assisted by the farmer, calving strongly assisted by the farmer or the veterinary practitioner, and missing values. Finally, calf birth weight was grouped into three classes: small weights (up to 42.9 kg), medium to large weights, and missing values.
|
Validation of the Baseline Hazard Function
The attractive feature of the semiparametric proportional hazard model is that it permits the estimation of the regression coefficients in ß, without making any assumption about the form of the baseline distribution function. However, the estimation of fixed and random effects and of genetic parameters is much less demanding of computing resources with a parametric model than with a semiparametric model (Yazdi et al., 2002
). A classical approach to choosing an adequate form for the distribution function of a random variable is to estimate its empirical survival function (Kaplan-Meyer estimate) and to compare its shape to that of known parametric distributions (Ducrocq, 1997
). In animal breeding, the Weibull model has been the parametric proportional hazards model of choice when longevity is analyzed. The Weibull model is a very flexible parametric model, in particular, when time-dependent covariates are included in the model (Yazdi et al., 2002
).
The Weibull hazard distribution is h0(t) =
(
t)
1, the Weibull baseline survival function being S0(t) = exp[(
t)
]. Then, a plot of ln(ln S0(t)) or ln(ln
KM(t)) against ln t gives a straight line: ln(ln S0(t)) =
ln
+
ln t (Ducrocq, 1997
). When
= 1, the Weibull hazard function is constant and is reduced to an exponential hazard function, h0(t) =
, the exponential baseline survival function being S0(t) = exp[
t]. Similarly, a plot of ln S0(t) against t gives a straight line as ln S0(t) =
t, thus the graphical test indicating whether an exponential model is suitable consists of checking that a plot of ln
KM(t) against t gives a straight line.
The inclusion of time-dependent covariates Pc converts the exponential hazard function into a piecewise exponential function:
![]() |
These formulas assume that the hazard of an individual is constant within each period of time but is different for the same individual between the c = 1, 2, ..., k + 1 different periods, corresponding to k cut points. Furthermore, the survival function associated with any hazard function is the exponential of the negative cumulative hazard function,
, and the piecewise exponential function becomes:
![]() |
There are k + 1 different formulas for piecewise exponential functions. These formulas can be summarized as follows:
![]() |
One property of the piecewise exponential function is that a plot of ln S0(t) against t yields c = 1, 2, ..., k + 1 different straight lines:
![]() |
This property suggests a graphical test for the piecewise exponential model that consists of checking that a plot of ln
KM(t) against t yields c = 1, 2, ..., k + 1 different straight lines. The slopes of these lines are the hazards in each period of time (
1,
2, ...,
k+1). To choose the cut points (c1, c2, ..., ck) of the piecewise function, we fit splines in the regression of the log of the KM survival function on time. To do that, all possible combinations of 1, 2, ..., k cut points in the time space (180 d) were explored, avoiding combinations of adjacent days.
As before, the inclusion of time-dependent covariates Pc converts the Weibull hazard function into a Weibull time-dependent hazard function:
![]() |
Again, the survival function associated is the exponential of the negative cumulative hazard function,
(formulas not shown). The Weibull time-dependent function has more degrees of freedom than its embedded distributions (Weibull, exponential and piecewise), and thus improves the likelihood of the model. This hypothesis was tested via LRT. In our analysis, a comparison between a Weibull time-dependent model and its embedded models was performed. The LRT under the null hypothesis was defined as:
![]() |
In this case, the LRT statistic has 1, k, and k + 1 degrees of freedom for the piecewise, Weibull, and exponential models, respectively.
Random Effects
The proportional hazards model can be extended to a frailty model with the addition of a random sire effect, si:
![]() |
The vector si of additive sire genetic effects, under polygenic inheritance, was assumed to follow a multivariate normal distribution s ~ MVN (0,A
S2), where A is the additive genetic relationship matrix among sires. The sire variance (
S2) was estimated using the Bayesian approach described in Ducrocq and Casella (1996)
. An obvious point estimate of the sire variance is the mode of this approximate marginal posterior density
^S2 . The general characteristics of the distribution through the computation of its first three moments by unidimensional numerical integration based on Gauss-Hermite quadrature have also been computed (Ducrocq and Casella, 1996
).
Heritability
The heritability on the binary scale (hbin2) is given by the formula
![]() |
which is similar to that reported by Yazdi et al. (2002)
. Here, the mode of the sire variance
^S2 is used. The heritability changes along time with the value of the survival function. Its derivation is presented in the Appendix.
| Results |
|---|
|
|
|---|
|
|
Hazard Ratios for Survival Factors
Herd-year effect had the strongest influence on calf risk of mortality (Table 2
). The mortality rate could change every year due to several sources of variation, such as natural forage availability, changes in herd management practices, and other environmental circumstances. The estimates of the corresponding hazard ratios were not shown, as they did not follow any consistent trend.
Calves born from September to February had the lowest mortality risk, increasing from March to August (P < 0.001; Table 3
). The calving season starts after the return of the cows from the mountains, in September or October, and by the beginning of March, half the cows had calved. That is, the risk of mortality increases, as does the cumulated distribution of calvings, probably due to a higher competition for resources and a higher risk of exposure to infectious diseases.
|
Unassisted calvings presented the least risk of mortality, and the risk increased as birth became more difficult (Table 3
). When assistance was very strong, calf mortality risk was more than five times that of calves born in an unassisted parturition (P < 0.001). Missing values probably corresponded to animals that were at a high risk of death before recording the variable whose value appears as missing, as happens when a calf dies just after birth. The inclusion of this category, however, might limit the bias in the estimation of the effects of the remaining categories.
Calves with medium and large birth weights had the lowest risk of death (Table 3
). This hazard tended to increase slightly when calf birth weights were small (P < 0.10). Future studies including a larger data set will likely help clarify the true effect of the birth weight of a calf on its future survival.
Validation of a Baseline Hazard Function
The Kaplan-Meier survival function may be approximated by a parametric function. Figure 2
presents plots for the Kaplan-Meier estimate of the survival function compared with plots of four parametric functions estimated in a model with no covariates from our data set. At first sight, both piecewise exponential and Weibull time-dependent functions showed a close fit to the observed (Kaplan-Meier estimate) survival function. The plots presented in Figures 3a and b
confirm this assertion and show that the observed function fits neither the Weibull or the exponential time independent functions. However, the observed curve was well approximated both by a piecewise exponential and a Weibull time-dependent distribution with cut points at 16 and 32 d. Furthermore, as the likelihoods of both piecewise and Weibull time-dependent survival functions were not statistically different for all models studied (results not shown in tables), the piecewise exponential model could be the function of choice to describe the mortality process of beef calves from birth to weaning.
|
|
1 = 0.001226 during the first 2 wk,
2 =0.000332 from 16 to 32 d, and
3 = 0.000046 after the first month until weaning (Figure 3b
Frailty Models
Similar estimates of hazard ratios for survival factors were obtained under a frailty model assuming different baseline hazard functions, with the exception of the exponential model (Table 4
). The estimates were also very similar to those obtained under the fixed model.
|
Conversely, the estimates of the random effect (sire variance) decreased with the number of survival factors included in the model (results not shown in tables). For Models 2 and 3, the sire variance estimates for the different baselines assumed were in the 0.52 to 0.54 and 0.46 to 0.51 ranges, respectively, whereas for Model 4, the estimates oscillated between 0.28 and 0.31.
Heritability
Heritability is not constant along the period studied, as it depends on the survival at each time point and the modal estimate of the sire variance. Assuming a sire variance of 0.3 and baseline values close to the KM estimates of the survival curve, heritability was always low, particularly during the first days of life. At 15 d, heritability is 0.023 and only reaches a value of 0.037 at weaning (180 d), assuming a piecewise exponential model (Figure 4
).
|
| Discussion |
|---|
|
|
|---|
According to Ray et al. (1989)
the two major components of percent calf crop weaned are birth rate and calf survival. Thus, calf survival is a parameter conditional on a calf born alive. In our case, evidence exists that some of the animals were born alive although an early death some hours after calving led the farmer to record this calf as dead. The consequence is that our survival would be an upward estimate of the true survival.
The survival experience was not constant over the lactation period. The mortality rate was more pronounced during the first 2 wk of life, diminished from 16 to 32 d, and after that was kept fairly constant until weaning. This observation is in agreement with several previous reports. Meijering (1984)
, Bellows et al. (1987)
, and Patterson et al. (1987)
reported in beef-calf surveys that 75% of mortality of all causes occurred by 7 d of age. Furthermore, Thurmond (1986)
found in dairy bull and heifer calves that the majority of deaths occurred between 9 and 16 d of age. These different mortality rates at different periods caused the Weibull distribution to fail in the validation of the baseline distribution, and a time-dependent distribution was needed. In this paper, the piecewise exponential function and a simple graphical method to validate its adequacy is proposed for calf survival analysis.
Several survival factors had an influence on calf survival. Herd-year effect had the strongest influence (Table 3
), as was reported previously (Ray et al., 1989
; Goyache et al., 2003
). Calf survival also depended on the season of calving. A similar phenomenon has been observed in Ripollesa lambs (J. Casellas, personal communication). Furthermore, Cubas et al. (1989)
and Goyache et al. (2003)
also included a calving season effect with two levels (from January 1 to June 30 and from July 1 to December 31) in the analysis of calf survival.
Calves born from cows younger than 1,300 d of productive life had a higher risk of death than calves born from the older cows. Greater death losses among calves from primiparous cows is a common observation that has been previously reported in Hereford, Angus, Charolais, Brahman, and crossbred cattle (see Ray et al., 1989
, and the review by Cundiff et al., 1986
). This fact could be explained, among other reasons, by an improvement of the maternal behavior as the cow matures.
Unassisted calvings showed the minimum instantaneous mortality rate, with the rate increasing as birth became more difficult. Cundiff et al. (1986)
reviewed several reports that have shown that calving difficulty has a significant effect on calf survival (Laster and Gregory, 1973
; Smith et al., 1976
; Gregory et al., 1978
; Menissier et al., 1981
). More recently, dystocia has been associated with calf mortality because of the possibility of uterine fluid inhalation (Bellows et al., 1987
). However, a proportion of this effect is genetic, as perinatal survival has a high positive genetic correlation with calving ease (Cubas et al., 1989
). In view of this high correlation, breeders interested in the improvement of calf survival rates normally have been recommended to use a correlated response to direct selection for calving ease (Cubas et al., 1989
).
Calves with medium and large birth weights had the lowest instantaneous mortality rate. This rate tended to increase slightly when calf birth weights were small. Similar results have been reported in the literature (Kulkarni et al., 1994
; review of Moore et al., 2002
), although there are some other problems associated with lighter birth weights. Bellows et al. (1987)
reported that 20% of abnormal beef calves (determined by necropsy) had birth weights at least 3.6 kg lighter than the average normal birth weight of 35 kg. Earlier onset of calf diarrhea was seen with the lightest weight dairy heifer calves (Pare et al., 1993
). Fallon et al. (1987)
reported that lighter Friesian bull-calves (<40 kg) were at higher risk for diarrhea and mortality than calves greater than 39 kg. A lower probability of survival of small birth weights was also reported in piglets. Casellas et al. (2004)
described a direct effect of birth weight on piglet survival justifying that small piglets need to spend more energy for thermoregulation purposes and generally suffer a delay in the first colostrum intake and a variable degree of physiological immaturity. On the other hand, Moore et al. (2002)
also reported that heavier calves had higher mortality because, as suggested by Meijering (1984)
, they may be responsible for dystocia due to feto-pelvic incompatibility. This effect was not detected in our study probably because of the inclusion of calving difficulty in the model. In future studies, calving difficulty might be analyzed jointly with other factors, such as birth weight, weaning weight, and survival rate in a multivariate analysis taking advantage of the genetic correlations among these traits.
The genetic random effect was estimated assuming a sire model. A sire-maternal grandsire model was also explored. Variance components estimated with both models were very similar. The sire-maternal grandsire model does not distinguish between direct and maternal effects that usually are very important in preweaning traits of beef cattle (see for example Meijering, 1984
; Cubas et al., 1989
; Goyache et al., 2003
). The ideal model of analysis would have been an animal model with maternal effects, but the structure of our data set did not allow us to develop such a model.
The sire variance estimated was around 0.3 for calf survival. This variance, estimated by assuming a proportional hazards model, could not be compared with the literature because published articles about calf survival estimated genetic variances assuming a threshold model. Until now, survival analysis techniques were only used to analyze cow longevity. Rogers et al. (2004)
estimated sire variances for beef cow longevity (between 0.028 and 0.037). In this sense, our sire variance for calf survival was very much higher than the cow longevity variances; however, cow longevity and calf survival are different characters and may not be comparable.
Despite this relatively high genetic variation, the corresponding heritability estimate for calf survival was low (0.037) at weaning (piecewise exponential model). This finding can be explained by the fact that heritability in the binary scale also depends on survival at a certain age, which in our case is very high at weaning (96.9%). For some other estimates of survival at weaning published in the literature (e.g., 0.95, 0.90 and 0.85), the heritabilities would be 0.059, 0.116, and 0.176, respectively. These estimates agree with those reported in the literature in which a threshold model was assumed. Cundiff et al. (1986)
, using a multibreed population, estimated a heritability for the direct effect of survival from birth to weaning of 0.07 within sire breeds and of 0.11 for the total population. Goyache et al. (2003)
estimated heritabilities considering the variables as calf traits ranging from 0.057 for late mortality (survival from 72 h to weaning) to 0.106 for weaning survival (from birth to weaning).
| Implications |
|---|
|
|
|---|
| Appendix |
|---|
|
|
|---|
![]() |
where the parameter p = S(t) is the probability that a progeny of sire i is still alive at time t. The analysis on the 0/1 scale gives an expectation of Y conditional on si:
![]() |
and a variance of Y conditional on si:
![]() |
The sire variance
s(b)2 is obtained using the delta method:
![]() |
where gradµ[S(t)] is the first derivative
evaluated at the mean value si = µs. Then, the heritability on the binary scale (hbin2) can be approximated as:
![]() |
where the suffix µ means that the function is evaluated at the mean breeding value si = µs. Because [gradµ (S(t))]2= [ln(Sµ (t))]2 Sµ (t)2, in the case of proportional hazard models, the heritability on the binary scale (hbin2) can be simplified to:
![]() |
This formula is similar to the one derived in a different way by Yazdi et al. (2002)
under the assumption of a Weibull distribution. Our derivation only supposes proportionality of hazards and can be applicable to any base distribution.
| Footnotes |
|---|
2 Correspondencephone: 34 93 5811399; e-mail:jesus.piedrafita{at}uab.es.
Received for publication September 27, 2004. Accepted for publication December 16, 2004.
| Literature Cited |
|---|
|
|
|---|
This article has been cited by other articles:
![]() |
A. Cecchinato, V. Bonfatti, L. Gallo, and P. Carnier Survival analysis of preweaning piglet survival in a dry-cured ham-producing crossbred line J Anim Sci, October 1, 2008; 86(10): 2486 - 2495. [Abstract] [Full Text] [PDF] |
||||
![]() |
X. F. de Sevilla, E. Fabrega, J. Tibau, and J. Casellas Effect of leg conformation on survivability of Duroc, Landrace, and Large White sows J Anim Sci, September 1, 2008; 86(9): 2392 - 2400. [Abstract] [Full Text] [PDF] |
||||
![]() |
J. Casellas, G. Caja, X. Such, and J. Piedrafita Survival analysis from birth to slaughter of Ripollesa lambs under semi-intensive management J Anim Sci, February 1, 2007; 85(2): 512 - 517. [Abstract] [Full Text] [PDF] |
||||
![]() |
J. Casellas, J. Tarres, J. Piedrafita, and L. Varona Parametric bootstrap for testing model fitting in the proportional hazards framework: An application to the survival analysis of Bruna dels Pirineus beef calves J Anim Sci, October 1, 2006; 84(10): 2609 - 2616. [Abstract] [Full Text] [PDF] |
||||
![]() |
J. J. Ghirardi, G. Caja, D. Garin, J. Casellas, and M. Hernandez-Jover Evaluation of the retention of electronic identification boluses in the forestomachs of cattle J Anim Sci, August 1, 2006; 84(8): 2260 - 2268. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |