doi:10.2478/v10014-011-0025-5
COBISS: 1.01 Agris category code: L01
ANALYSIS OF LONGEVITY IN SLOVENIAN HOLSTEIN CATTLE
Klemen POTOČNIK Vesna GANTNER 2, Jurij KRSNIK Miran ŠTEPEC Betka LOGAR 3, Gregor GORJANC 1
Received August 20, 2011; accepted September 22, 2011. Delo je prispelo 20. avgusta 2011, sprejeto 22. septembra 2011.
Analysis of longevity in Slovenian holstein cattle
The longevity of Slovenian Holstein population was analysed using survival analysis with a Weibull proportional hazard model. Data spanned the period between January 1991 and January 2010 for 116,200 cows from 3,891 herds. Longevity was described as the length of productive life - from first calving till culling or censoring. Records above the sixth lactation were censored to partially avoid preferential treatment. Statistical model included the effect of age at first calving, stage of lactation within parity, yearly herd size deviation, season defined as year, herd, and sire-maternal grandsire (mgs). Some effects had time varying covariates, which lead to 1,839,307 or on average 16 elementary records per cow. Herd and sire-maternal grandsire effects were modelled hierarchically. Pedigree for sires and maternal grandsires included 2,284 entries. Estimated variance between herds was 0.12, while between sire variance was 0.04. Heritability was evaluated at 0.14. Genetic trend for sires was unfavourable, but not significant. A further research is needed to define the required number of daughters per sire and the dynamics of genetic evaluation for sires whose majority of daughters still have censored records.
Key words: cattle / breeds / Slovenian Holstein / longevity / Weibull proportional hazards model
1 INTRODUCTION
Longevity is a trait with great impact on dairy production economy and is, therefore, of considerable importance in dairy cattle breeding programmes (Charffed-
Analiza dolgoživosti pri črno-beli pasmi goveda v Sloveniji
Za analizo dolgoživosti smo pri slovenski črno-beli populaciji govedi uporabili metodologijo analize preživetja in Weibullov model sorazmernih ogroženosti. V analizo smo vključili podatke 116.200 krav iz 3.891 čred skozi obdobje od januarja 1991 do januarja 2010. Dolgoživost je bila predstavljena kot doba produktivnega življenja, ki je definirana kot število od prve telitve do izločitve ali do datuma zajema podatkov za živali, ki so na ta datum bile še žive. Šesto in kasnejše laktacije smo okrnili na konec šeste laktacije, da smo omilili precenjenost boljših živali. V statistični model smo vključili vpliv starosti ob prvi telitvi, stadija laktacije ločeno za vsako zaporedno laktacijo, spreminjanje velikosti črede med leti, leto telitve, čredo, očeta in materinega očeta. Ravni nekaterih vplivov so časovno spremenljivi, kar povzroči, da smo v analizi obravnavali 1.839.307 zapisov ali povprečno 16 osnovnih zapisov na kravo. Čreda in vpliv očeta z materinim očetom sta bila v model vključena hierarhično. Rodovnik za očete in materine očete je obsegal 2.284 zapisov. Ocenjena varianca za vpliv črede je znašala 0,12, medtem ko je ocena variance med očeti znašala 0,04. Dednostni delež je bil ocenjen na 0,14. Genetski trend ima negativno smer, a ni statistično značilen. Potrebne bodo nadaljnje raziskave, da bomo določili zadostno število hčera po biku in dinamiko obračunov plemenskih vrednosti za bike, ki imajo večino hčera še v fazi prireje.
Ključne besede: govedo / pasme / slovenska črno-bela pasma / dolgoživost / Weibullov model sorazmernih ogroženosti
dine et al., 1996; Strandberg and Soelkner, 1996). With the increase of longevity, the proportion of mature cows that produce more milk increases. For example, Strandberg (1996) estimated that an increase in longevity from three to four lactations increases average milk yield per
1	Univ. of Ljubljana, Biotechnical Fac. Dept. of Animal Science, Groblje 3, SI-1230 Domžale, Slovenia
2	J.J. Strossmayer Univ. in Osijek, Fac. of Agriculture, Trg Svetog Trojstva 3, 31000 Osijek, Croatia
3	Agricultural institute of Slovenia, Hacquetova 17, SI-1000 Ljubljana, Slovenia
Acta argiculturae Slovenica, 98/2, 93-109, Ljubljana 2011
K. POTOČNIK et al.
lactation and profit per year between 11 and 13%. In addition, improvement in longevity decreases replacement costs and somewhat increases selection intensity.
There are several ways to implement selection on longevity in the breeding goal, directly or indirectly. Direct longevity can be represented as the length of (productive) life (LPL) or stayability. In cattle breeding LPL is usually defined as the elapsed time between the first calving and culling, while stayability is defined as a binary trait that measures cow survival (live or culled) at a certain point in time. The use of LPL is preferred since stayability as a discrete trait provides less information. Unfortunately, LPL, as well as stayability, can be quantified only after the cows are culled, though both approaches provide partial information when cow survives to the next "period" in life. Therefore, the information on the longevity of daughters of a sire becomes available with the increasing age of a sire. This inherently leads to the prolonged generation interval. Low heritability for longevity (Short and Lawlor, 1992; Vollema and Groen, 1996) induces unreliable estimation of breeding values (BV) based only on the information of parents or grandparents.
Due to long generation interval, breeding programmes also include indirect measures of longevity via correlated traits such as fertility, health, and conformation traits (Burnside et al., 1984). Additional gain is due to the fact that the data on these indirect traits can be collected relatively early in the life of a cow. Nonetheless, both representations of longevity (direct and indirect) have a merit in a modern breeding goal (Essl, 1998).
Analysis of indirect representations of longevity is to a large extent done with a standard linear model based on the Gaussian (normal) distribution. Specific approach is needed for a proper analysis of the LPL, due to the presence of live animals at the time of analysis (censored records) and changes in culling criteria over the productive life of cows (time varying covariates) (e.g. Ducrocq et al., 1988a). Exclusion of censored records from the analysis, or treating them as uncensored leads to biased results (Ducrocq, 1994). Additionally, relationship between longevity and its effects is rather multiplicative than additive (e.g. Ducrocq et al., 1988a). Survival analysis can handle this kind of data. In the last years several countries introduced direct longevity in the routine genetic evaluation of cattle and most of them use the Weibull proportional hazard model (INTERBULL, 2009), which represents a class of models in the field of survival analysis. Other statistical approaches (models) can also be used, but proportional hazard model have better properties (e.g. Caraviello et al., 2004; Jamrozik et al., 2008; Potocnik et al., 2008).
The aim of this study was to present the results of
genetic evaluation for the length of productive life in Slovenian Holstein population using a Weibull proportional hazard model.
2 MATERIAL AND METHODS
2.1 DATA
Raw data for 126,716 Slovenian Holstein cows born from 1982 to 2008 were provided by the Agricultural Institute of Slovenia. In order to use old data but to avoid modelling the data up to the year 1991, the truncation date was set at January 1st 1991. On the other side, the date of last data collection was January 29th 2010. For cows alive at that time longevity was treated as right censored. Longevity was defined as the length of productive life (LPL) and was calculated as the number of days from the first calving to culling (uncensored/complete records) or to the moment of data collection (incomplete/censored records). The LPL of cows surviving beyond the sixth lactation was also censored in order to avoid the effect of preferential treatment and to focus on early culling in the life of a cow. Cows with missing or inconsistent data within the defined limits were removed (29,252 cows): culling before the date of truncation, calving date after the date of culling, no information for 600 days after calving, missing data for the first three lactations, daughters of sires with less than 20 daughters, and missing covariate or factor data.
The structure of used data and descriptive statistics are given in Table 1. Altogether LPL for 116,200 cows from 3891 herds were used in the analysis. Cows in the analysis were daughters of 707 sires, while the whole sire-maternal grandsire pedigree consisted of 2,284 sires. Cows were on average culled in the third lactation, which
Table 1: Structure of data and descriptive statistics (± standard deviation)
Preglednica 1: Struktura podatkov in opisna statistika (± standardni odklon)
Cows, no.	116,200
Sires, no.	707
Pedigree, no.	2,284
Censored records, %	41.0
Number of lactations in life	
uncensored records	3.0
censored records	3.0
Length of productive life, days	
uncensored records	1,095 ± 660
censored records	1,129 ± 754
94 
Acta agriculturae Slovenica, 98/2 - 2011
ANALYSIS OF LONGEVITY IN SLOVENIAN HOLSTEIN CATTLE
amounted to 1,095 days of productive life. Percentage of censored records was 41.0%. Censored records had about the same means, but larger variability.
2.2 SURVIVAL ANALYSIS
Weibull proportional hazards model was used for the analysis of LPL. This model is built upon the Weibull distribution, whose density (1) and hazard (2) function for the ¿-th record t are:
f(t. | A, p) = Ap(At)p-1 exp(-(At)p), h(t. | A, p) = Ap(At)p,
(1) (2)
where X (scale) and p (shape) are strictly positive parameters. In proportional hazard model it is assumed that the baseline hazard function changes proportionally with change in covariate(s) or factor levels. For the analysis of LPL the hazard function was modelled as:
h(t..,, | X, p, else) =
v ijklmnop 1 '	'
h„(t..,, | X, p) exp(c. + l. + y + h + d + s + 1/2s ), (3)
0 v ¿jklmnop' r '	j ' k l m n	o/ v/
where:
h(t ¡jktmnop I X, P, else) = hazard of culling p-th cow given other parameters,
h (t,,, 0 s ijklmnop
yk h
s + 1/2s
| A, p) = baseline Weibull hazard function (2),
= ¿-th age at first calving: 0 (unknown) and
from 19 to 50 months, = j-th lactation stage (1-60 days, 61-150 days, 151-270 days, 271-days till drying, and dry period) within parity - altogether 30 levels; time varying factor, = k-th season defined as year (1990-2010);
time varying factor, = l-th herd (3891 levels); time varying factor, = m-th herd size deviation in comparison to previous year (< -70%, (-70%, -40%], (-40%, -10%], (-10%, 10%], (10%, 40%], (40%, 70%], and > 70%); time varying factor, = n-th sire and the o-th maternal grandsire (onwards both effects are termed sire effect) of the p-th cow.
Levels of time varying factors (lactation stage within parity, year, herd, and herd size deviation) changed with cow "status" changes in time creating subsequent elementary records, while levels for others effects were constant over whole lifetime of a cow. Altogether, there were 1,839,307 elementary records. Herd and sire effects were modelled hierarchically: log-gamma distribution for herd effect and multivariate normal for sire effect with additive genetic covariance matrix build from the
pedigree. The used Weibull proportional hazards model and the corresponding assumptions can be sketched in matrix form as:
y | b, h, s, p ~ Weibull (Xb + Zh + Ws, p), h | y ~ Log - Gamma (y, y), s | G ~ Normal (0, G), where:
(4)
(5)
(6)
b = a vector with intercept p ln (y)) and parameters for the following effects: age at first calving, stage of lactation within parity, year, and the deviation of herd size from year to year, h = the vector of parameters for herd effect, s = the vector of parameters for sire effect, y = Log-Gamma distribution parameter,
G = additive genetic covariance matrix among sires - a product of numerator relationship matrix between sires A and additive genetic variance between sires (ff^).
Heritability according to the model (3, 4-6) was calculated following Meszaros et al. (2010):
h2 =
4<t:
oî+ytoî+</'1(r)+1 Xtoî+oî+ï
(7)
where:
ff2 + /ff2 is variance due to sire and maternal grandsire effects (3), S d2
vm(x) = — log r(x)	it.
ox	is a trigamma function used to evaluate the vari-
ance of log-gamma process (5) giving between herd variance (ff2), while the value of 1 is the underlying residual variance.
Data processing was done with SAS software package (SAS Institute, 2000), while Survival Kit version 3.10 (Ducrocq and Soelkner, 1998) was used for modelling and parameter estimation. In the first step a series of log-likelihood ratio tests were performed for effects that were not modelled hierarchically - the importance of each effect was tested as a comparison between the full model and the model where the effect under testing was excluded. In the next step herd and sire effects were added to the model to obtain estimates for all model parameters. In results relative risk is equal to the value of solutions for model parameters on exponential scale (3) proportional to some specified baseline level that has a solution equal to 1 (e.g. 26 months for the age at first calving). Each plot of relative risks is also augmented with the number of censored and uncensored records per level of a factor. In the case of time varying factors only the last elementary record of a cow was considered for computing the number of records per level of a factor.
c
d
m
Acta agriculturae Slovenica, 98/2 - 2011 95
K. POTOČNIK et al.
Figure 1: Relative risk of culling and number of records by age at first calving Slika 1: Relativno tveganje za izločitev in število zapisov glede na starost ob prvi telitvi
3 RESULTS AND DISCUSSION
All effects included in the model were highly significant (P < 0.001), which is not surprising given the size of data set and the previous knowledge of effect importance for LPL. Distribution of age at first calving was expected with the majority of cows in the range between two and three years of age (Figure 1). Relative risk
of culling increased almost linearly with the increasing age at first calving. Similar results were obtained also by Vollema and Groen (1998) and Rogers et al. (1991), while others did not found significance (Ducrocq et al., 1988a; Ducrocq, 1994) or concluded that this effect was not important (Vukasinovic et al., 1997). This can be at least partially attributed to the fact that our results do not directly imply causal relationship between LPL and the
Figure 2: Hazard of culling and number of records by stage of lactation within parity
Slika 2: Ogroženost za izločitev in število zapisov glede na stadij laktacije in zaporedno laktacijo
96 Acta agriculturae Slovenica, 98/2 - 2011
ANALYSIS OF LONGEVITY IN SLOVENIAN HOLSTEIN CATTLE
Table 2: Estimates of model hyperparameters and derived quantities
Preglednica 2: Ocene parametrov modela in izpeljanih količin
Hyperparameter / Quantity	Estimate
Shape, p	2.00
Log-gamma parameter, y	8.50
Between herd variance, o2 = f(1)(y)	0.12
Between sire variance, o2 s	0.04
Additive genetic variance, o2 = 4o2	0.16
	0.14
age at first calving. Results only imply that there is association between LPL and the age at first calving in our population, which indicates that cows that had late first calving had also some other problems (likely related to reproductive success) that increased risk of being culled early. Estimates at the start and the end of considered age interval were very variable due to the smaller number of records. Given almost linear relationship, variable results at margins, and that age is time independent effect a possible approach would be to model this effect with linear regression. Regression is not appropriate for time dependent effects due to the explosion of number of elementary records.
The stage of production has a significant effect on risk of a cow being culled due to biological (increased probability of mastitis at the start of lactation) or tech-
nological factors (owners' decisions in the dry period). Since the stage of lactation within parity changes with increasing age (dependent variable) we represented this effect using baseline hazard function (3) multiplied with the corresponding relative risk for each stage within parity (Figure 2) for a fixed calving interval of 400 days. Number of culled cows was highest in the second parity and dropped in later parities. Hazard increased over time with considerable changes at the end of lactation - that is in the period between 271 days after lactation and dry-off and in the dry period. Virtually the same results were obtained also in other studies (e.g. Ducrocq, 1994; Vuka-sinovic et al., 1997; Potocnik et al., 2010). The first parity showed unique pattern with increased hazard in the first two periods of lactation (1-60 and 61-150 days), which might be due to the higher incidence of health disorders during early lactation. Similar estimates were obtained also by Ducrocq (1994) and Vukasinovic et al. (1997).
Herd size dynamics has also influence on culling criteria. In general herd expansion lowers risk, while risk is higher in shrinking herds. Herds in Slovenia are in general small, so there is substantial variability in herd size changes from year to year. Majority of records (censored or not) were in the range of -40 to 40% of herd size change. Relative risk for culling was rather constant for herd size change levels above -40%, while it increased for the two levels bellow this threshold as expected - cows from herds with decreasing size have larger probability of being culled. Weigel et al. (2003) calculated the relative culling risk of high producing (top 20%) and low
Figure 3: Relative risk of culling and number of records by levels of variation in herd size between years Slika 3: Relativno tveganje za izločitev in število zapisov glede na letno spremembo velikosti črede
Acta agriculturae Slovenica, 98/2 - 2011 97
K. POTOČNIK et al.
producing (bottom 20%) cows relative to average cows in the same herd with regard to herd size changes. They determined that, before herd expansion, low producing cows were 4.2 times more likely to be culled than average cows, while high producing cows were only 0.5 times as likely to be culled as average cows. After herd expansion, the relative risk for low producing cows dropped to 2.6 times that of average cows, and the risk for high producing cows increased slightly to 0.7 times that of average cows, which clearly shows the effect of herd expansion on the reduced risk of culling.
Culling criteria also change with time. Possible reasons are disease outbreak in population, change of prices, milk quota, etc. In order to capture such variations we included in the model effect to season defined as year. Our data spanned period between 1991 and 2010 (Figure 4). Relative risk of culling was very low in the first years due to lack of records in that period, but reached overall level after the year 1995. After this year we observe overall increase in relative risk of culling over years. A sharp decrease in risk was estimated for the year 2002, which can be attributed to farmers' decision to keep more animals on farm in order to get higher milk quota and subsidies per farm with the forthcoming new quota and subsidy system in Slovenia at that period. Risk was logically very low in the last year (2010) due to the large number of live (censored) animals in the analysis.
Based on the used statistical model, between sire variance was estimated to 0.04, while between herd variance was estimated to 0.12 (Table 2). These values are on exponential scale of the Weibull model (3) and do not
o —:
u
O
tu
ai
LTi ©
Figure 4: Relative risk of culling and number of records by year Slika 4: Relativno tveganje za izločitev in število zapisov glede na leto
have meaningful units related to the analysed variable (LPL). Estimated heritability using the approach of Mes-zaros et al. (2010) was 0.14, which is similar to the values reported in Austria (0.18) and Germany (0.16) and relatively high comparing with other countries that are members of INTERBULL (INTERBULL, 2009).
Breeding values for LPL were presented on scale with mean 100 and standard deviation 12 with favourable values (longer LPL) above 100. Genetic trend by year of birth for sires in the period 1984-2005 shows that there was no selection on longevity (Figure 5). Overall trend was negative (-0.12 ± 0.14), but not significant (p = 0.385). Differences between years were minimal except for the last four years. This could be attributed to the smaller number of evaluated sires and to the fact that these sires have a lot of daughters with censored records.
The accuracy of genetic evaluation highly depends on the ratio of censored and uncensored records. As the proportion of censored records decreases, the evaluation accuracy increases. Also, it is necessary to have sufficient number of daughters per sire. Vukasinovic et al. (1997) stated that more than 30 to 40% of censored records would lead to inaccurate results. Same authors stated that small number of daughters per sire without any or only few uncensored records biases sires ranking. Egger-Danner et al. (1993) performed retrospective study where they compared the ranking of sires from a full data file without censored records and from a truncated data with a different proportion of censored records. They observed that rank correlations between breeding values dropped as the proportion of censoring increased. Fur-
98 Acta agriculturae Slovenica, 98/2 - 2011
ANALYSIS OF LONGEVITY IN SLOVENIAN HOLSTEIN CATTLE
Figure 5: Genetic trend for the length of productive life in Slovenian Holstein sires Slika 5: Genetski trend za dolgoživost pri plemenjakih slovenske črno-bele pasme
ther research is needed in our population to determine the impact of censored records on the accuracy of genetic evaluation as well as to determine how many daughters per sire are needed.
4	CONCLUSION
Survival analysis methodology was applied for the genetic evaluation of longevity (defined as the length of productive life) in Slovenian Holstein cattle. Statistical model included the effect of age at first calving, lactation stage within parity, yearly herd size deviation, year, herd, and additive genetic effect as captured by sire and maternal grandsire effects. Parameter estimates were similar to studies in other countries. Genetic trend was slightly negative (unfavourable), yet not significant. Relatively high differences between average breeding values were observed in last years. As accuracy of genetic evaluation highly depends on the number of daughters per evaluated sire and on the ratio of censored and uncensored records further investigations are needed.
5	REFERENCES
Burnside E.B., McClintock A.E., Hammond K. 1984. Type, production and longevity in dairy cattle: a review. Animal Breeding Abstracts, 52, 711-719 Caraviello D.Z., Weigel K.A., Gianola D. 2004. Comparison between a Weibull proportional hazards model and a lin-
ear model for predicting the genetic merit of US Jersey sires for daughter longevity. Journal of Dairy Science, 87, 1469-1476. http://jds.fass.org/cgi/content/abstract/87/5/ 1469 (26. Mar. 2010) Charffeddine N., Alenda R., Carabano M.J., Bejar F. 1996. Selection for total merit in the Spanish Holstein-Friesian. INTERBULL Bulletin, 12, 142-146 Ducrocq V. 1994. Statistical analysis of length of productive life for dairy cows of the Normande breed. Journal of Dairy Science, 77, 855-866. http://jds.fass.org/cgi/content/ab-stract/77/3/855 (26. Mar. 2010) Ducrocq V., Quaas L.R., Pollak E.J., Casella G. 1988a. Length of productive life of dairy cows. 1. Justification of a Weibull model. Journal of Dairy Science, 71, 3061-3070. http://jds.fass.org/cgi/content/abstract/71/11/3061 (26. Mar. 2010)
Ducrocq V., Quaas L.R., Pollak E.J., Casella G. 1988b. Length of productive life of dairy cows. 2. Variance component estimation and sire evaluation. Journal of Dairy Science 71, 3071-3079. http://jds.fass.org/cgi/content/abstract/71/11/ 3071 (26. Mar. 2010) Ducrocq V., Soelkner J. 1998. The Survival Kit - V3.0: A Package for Large Analysis of Survival Data. In: Proceedings of the 6th World Congress on Genetics Applied to Livestock Production, Armidale, Australia, 22, 51-52 Egger-Danner C., Soelkner J., Essl A. 1993. Zuchtwertschätzung für Merkmale der Langlebigkeit: Vergleich der Proportional Hazards Analyse für das Ausfallsrisiko mit BLUP für die Verbleiberate. Vortragstagung der DGfZ/GfT, Göttingen, Germany Essl A. 1998. Longevity in dairy cattle breeding: a review. Livestock Production Science, 57, 79-89. http://dx.doi.org/10.1016/S0301-6226(98)00160-2 (26. Mar. 2010)
Acta agriculturae Slovenica, 98/2 - 2011 99
K. POTOČNIK et al.
INTERBULL 2009. Interbull Routine Genetic Evaluation for Direct Longevity January 2009 - Table 3. http://www-interbull.slu.se/longevity/l-table3-091.html (20. Jan. 2009)
Jamrozik J., Fatehi J., Schaeffer L.R. 2008. Comparison of models for genetic evaluation of survival traits in dairy cattle: a simulation study. Journal of Animal Breeding and Genetics, 125, 75-83.
http://dx.doi.org/10.1111/j.1439-0388.2007.00712.x (26. Mar. 2010)
Meszaros G., Palos J., Ducrocq V., Soelkner J. 2002. Heritability of longevity in Large White and Landrace sows using continuous time and grouped data models. Genetics Selection and Evolution, 42, 13.
http://dx.doi.org/10.1186/1297-9686-42-13 (26. May 2010)
Potočnik K., Krsnik J., Štepec M., Gorjanc G. 2008. Comparison between methods for estimation of breeding values for longevity in Slovenian Holstein population. INTERBULL Bulletin, 38, 61-65
Potočnik K., Gorjanc, G., Štepec, M., Krsnik J. 2010) Genetic evaluation for longevity in Slovenian Simmental cattle. Acta Agraria Kaposvariensis, 14, 153-160
Rogers G.W., Hargrove G.L., Cooper J.B., Barton E.1991. Relationships among survival and linear type traits in Jerseys. Journal of Dairy Science, 74, 286-291. http://jds.fass.org/cgi/content/abstract/74/1/286 (26. Mar. 2010)
SAS Institute. 2000. The SAS System. Version 8. Cary, SAS Institute
Short T.H., Lawlor T.L. 1992. Genetic parameters of conformation traits, milk yield, and herd life in Holsteins. Journal of Dairy Science, 75, 1987-1998.
http://ids.fass.Org/cgi/content/abstract/75/7/1987 (26. Mar. 2010)
Strandberg E. 1996. Breeding for longevity in dairy cows. In: Progress in Dairy Science. Phillips, C.J.C. (ed.). Oxon, CAB International, Wallingford, 125-144
Strandberg E., Soelkner J. 1996. Breeding for longevity and surival in dairy cattle. In: Proceedings of the International Workshop on Genetic Improvement of Functional Traits in Cattle. Gembloux, Belgium, 12, 111-119
Vollema A.R., Groen A.F. 1996. Genetic parameters of longevity traits of an upgrading population of dairy cattle. Journal of Dairy Science, 79, 2261-2267. http://ids.fass.org/cgi/content/abstract/79/12/2261 (26. Mar. 2010)
Vollema A.R., Groen A.F. 1998. A Comparison of Breeding Value Predictors for Longevity Using a Linear Model and Survival Analysis. Journal of Dairy Science, 81, 3315-3320. http://ids.fass.org/cgi/content/abstract/81/12/3315 (26. Mar. 2010)
Vukasinovic N., Moll J., Kuenzi N. 1997. Analysis of Productive Life in Swiss Brown Cattle. Journal of Dairy Science 80, 2372-2579.
http://ids.fass.org/cgi/content/abstract/80/10/2572 (26. Mar. 2010)
Weigel K.A., Palmer R.W., Caraviello D.Z. 2002. Investigation of factors affecting voluntary and involuntary culling in expanding dairy herds in Wisconsin using survival analysis. Journal of Dairy Science, 86, 1482-1486. http://ids.fass.org/cgi/content/abstract/86/4/1482 (26. Mar. 2010)
100 Acta agriculturae Slovenica, 98/2 - 2011