Joint modelling of repeated measurements and event time: Application to performance traits and survival of lambs bred in sub-humid tropics

We considered the analysis of a study for Dorper, Red Maasai and crossbred lambs born over a period of 6 years at the Diani Estate, Kenya. The study was designed to compare survival and performance traits of genotypes with differing susceptibilities to helminthiasis. The available data include information on time to death and repeated measurements of body weight, packed cell volume (PCV) and faecal egg count (FEC) of the animals. In the paper, we consider joint modelling of the survival time and the repeated measurements. Such an approach allows to account for the possible association between the survival and repeated measurement processes. The advantages and limitations of the joint modelling are discussed and illustrated using the Diani Estate study data.


INTRODUCTION
The data used in this study came from an animal breeding experiment carried out by the International Livestock Research Institute (ILRI) from 1991 to 1996 [1] at the Diani Estate, Kenya Coast. The objective of the experiment was to study genetic resistance to naturally acquired gastro-intestinal nematodes in different breeds of sheep, namely the Red Maasai (found in East Africa and perceived to be resistant to helminthiasis), Dorper (originating from South Africa and presumed to be susceptible) and their cross breeds. Three ewe genotypes (Dorper, Red Maasai, and Red Maasai × Dorper) were mated to Dorper and Red Maasai rams in single-sire mating groups in a diallel design to generate two pure-bred and four cross-bred genotypes. The data collected during the study consisted of repeated measurements of the traits: body weight (BWT), packed cell volume (PCV) and faecal egg count (FEC), which were measured periodically over a lamb's first year of life. Overall, 1745 lambs were born alive during the six-year period of the study. The repeated measurements, however, were highly unbalanced since 655 (38%) of the lambs died and 92 (5%) were stolen before they reached one year of age.
An assessment of the level of genetic resistance has been carried out by Baker et al. [1,2], by applying linear mixed models to individual measurements of the traits that were recorded from birth to one year of age. In the analysis undertaken by these authors, only the information on the animals that survived to each of the analysed time points was utilised. Nguti et al. [16] used these same traits in a frailty model for survival, where a correlation was induced among the ages at the time of death for lambs from the same sire. In that analysis, PCV and FEC were separately included in the model as time-varying covariates, with either a baseline or time-varying BWT. One handicap of this latter analysis was that the measured traits were only recorded at monthly time points, and were thus unknown at the time of death. To overcome this problem, the nearest preceding measured value of the trait was used to impute for the value at the time of death, resulting in a piecewise covariate process. In addition, the analysis did not correct for the possible measurement error in the recorded values of BWT, PCV and FEC. Such an error can result from both the analytical error and short-term, biological variability. Failure to account for the measurement error and for any missing time-varying covariate observations has been shown to cause the estimated regression parameters in the Cox proportional hazard (PH) model to be biased towards the null [4,18].
In the study described, as in many longitudinal studies where individuals are followed over time, the data can be grouped into three categories: (1) the elapsed time to an event such as death; (2) repeated measurements of timevarying variables (like PCV, FEC, BWT); (3) time-varying (e.g., rainfall) or constant baseline (e.g., sire, breed) covariates that may affect both the repeated measurement and the time-to-event processes. When modelling of the repeated measurements is of interest, one may focus, for instance, on how the measurements change with time, on how the parameter estimates are influenced by the drop-out of individuals during the course of the study, or on how the measurements may be affected by the additional covariates. Looking at things from a time-to-event process point of view, the interest may focus on, e.g., how the time to the event is affected by both the repeated measurement process and the additional covariates. A vast amount of literature exists on the methods suitable for either approach. To analyse event-time data, parametric or non-parametric models can be used, but the Cox PH model is often the method of choice. Repeated measurements are commonly analysed using linear mixed effects models [13]. These models are attractive for several reasons, one of them being the ability to easily accommodate unbalanced designs, especially regarding the timing and frequency of the observations. The models also allow for an explicit partitioning of variability and estimation of individual effects. In particular, at least two sources of variability are readily identified: between-and within-individual variation. The between-individual variability is often modelled by a vector of correlated, individual random effects.
In the last ten years, many methods, which simultaneously use the information available in both the time-to-event and the repeated measurement processes, have been proposed in medical research. In particular, several models have been developed in the area of acquired immunodeficiency syndrome (AIDS) research [5,7,20,22,24] and in schizophrenia studies [9,25]. A detailed review of research work in joint modelling of times to an event and repeated measurements is given in the reference [23]. The need for joint models to model survival and performance traits in animal studies is discussed in [6].
Several advantages of joint modelling of the repeated measurement and the time-to-event processes have been highlighted in the literature: (1) the repeated measurements can be extrapolated from the observed measurement times to the specific event time in a way that utilises the entire measurement history; (2) the time-to-event is allowed to depend on the 'true' but unknown value of the repeated measurement, thus making an adjustment for the measurement error, which in turn leads to reduced bias of the parameter estimates of the Cox model; and (3) the repeated measurement process is adjusted for any loss of information arising from death or loss of individuals.
The objective of the current paper was to use the joint modelling approach to model the time to death and the repeated measurements of PCV, BWT and FEC. For this, the methodology proposed by Henderson et al. [9] was used. The paper is organised as follows. In Section 2 more details on the motivating dataset are given. Section 3 provides a brief background on the linear mixed effects and the Cox PH models, as well as the joint model of Henderson et al. [9]. In Section 4 we adopt the joint model to the analysis of the motivating dataset. The results are presented in Section 5 and discussed in Section 6.

MOTIVATING DATA
Measurements of PCV, FEC and BWT were taken from lambs from birth to one year of age in batches of lambs born in each of the years 1991 to 1996. Packed cell volume and FEC were measured according to the methods reported by Baker et al. [1]. All lambs were weighed at birth and their BWT, PCV and FEC were subsequently recorded at one and two months of age. On either of these latter occasions, when individual lambs had an FEC greater than or equal to 2000 eggs per gram (epg) and/or a PCV less than or equal to 20%, they were treated (drenched) with an anthelmintic drug. Lambs with low PCV were also checked for trypanosome infections. At about three months of age, when weaned, the lambs were again weighed, and blood and faecal samples were collected for PCV and FEC, respectively. All lambs were then drenched. The lambs were then left to graze on pasture, separately from the ewes and rams. Every week, a monitor group of about 50 lambs, made up of approximately equal numbers of lambs of each genotype and sex, was sampled and their mean FEC was recorded. If the mean FEC was over 2000 epg then, during two consecutive days, all lambs were weighed, faeces and blood samples were taken and the lambs were drenched. This procedure was followed until the lambs reached on average one year of age. It resulted in five drenchings per year except for 1994 and 1996. In 1994, the year with the highest rainfall, the lambs were drenched eight times post-weaning, while in 1996 six drenchings occurred.

METHODOLOGY
In this section, we first briefly discuss the methods used in (separate) modelling of repeated measurements and time-to-event data. Then we shortly review joint modelling.

Linear mixed effect models for repeated measures
Data sets resulting from follow-up studies are often highly unbalanced, with subjects having an unequal number of measurements. Moreover, the data have complex correlation structure due to repeated measurements for each individual. As a result, such data are not ideally suited to analysis by classical least squares techniques and linear mixed effects models [13] are now standard tools for analysing such complex hierarchical data.
Let Y i j denote the observed jth measurement for the ith individual recorded at time t i j (i = 1, . . . , N; j = 1, . . . , n i ) and let Then a linear mixed effects model is written as where X 1i and Z i are n i × p and n i × q design matrices, respectively, β L is a p × 1 vector containing the fixed effects, and s i is a q × 1 vector of the random effects. It was assumed that s i is N(0, G), i.e., it is normally distributed with mean zero and variance-covariance matrix G = [g kl ], where g kl = Cov(s ik , s il ). Furthermore, it was assumed that s i is independent from the vector of residual random errors ε i . The residual errors were assumed to be N(0, R i ), with variance-covariance matrix R i depending on i only via its size (n i × n i ). It then follows that, marginally, Y i is normally distributed with mean X 1i β L and variance-covariance matrix V i = Z i GZ T i + R i . In model (1), R i captures the within-individual variability while the between-individual variability is modelled through the random effects s i . In particular, if then model (1) is known as a random intercept and slope model (see reference [21], p. 25). The underlying assumption of this model is that the measurements increase linearly with time, but for each individual the linear trend has its own intercept and slope. Furthermore, if Var(ε i j ) = σ 2 e , the result that the assumed covariance function of the response for this model is where 1(A) is the indicator function of event A. Note that function (2) is quadratic over time.
Model (1) has been used extensively to analyse repeated measurements arising from animal breeding programmes [8,11,14,15]. In these applications more emphasis has been placed on the covariance structure of the random effects (s i , ε i ) in order to capture different sources of variability, such as those due to maternal, paternal and environmental effects.
To estimate the parameters of model (1), various approaches can be applied. One approach is the classical method of maximum likelihood (ML), which results in generalised least squares (GLS) estimates for β L . This method of estimation, however, leads to underestimation of the variance parameters involved in G and R i . As an alternative, the restricted maximum likelihood estimation (REML) is more commonly used, which remedies this problem [21].

The Cox proportional hazard model
Let the ith individual (i = 1, . . . , N) be observed from a time zero to a failure time T i or to a potential right censoring time C i . Let T o i = min(T i , C i ), be the observed time and δ i be the censoring indicator which is equal to 1 if T o i = T i and 0 otherwise. Hence the observed data for the ith subject are (T o i , δ i ). The basic analytical quantities for time-to-event data are the survivor function which is the probability of surviving beyond time t, and the hazard function which is the instantaneous failure rate after surviving up to time t. In most time-to-event studies, interest focuses on how the hazard function is affected by independent covariates. To assess this effect, Cox [3] has proposed the following hazard model: where λ 0 (t) is the baseline hazard function common to all individuals, x 2i is a p × 1 vector of the observed covariates for the ith individual, and β S is the corresponding vector of regression parameters to be estimated. Model (3) is known as the Cox proportional hazard model and has been used extensively over the last three decades in the analysis of failure-time data. The Cox model assumes that the analysed time-to-event observations are independent. Models which allow to account for correlation of the failure times have been studied extensively (see [10,12,19]). With clustered time-to-event data, a common approach is to add a random multiplicative factor to the Cox model, shared by the observations from the same cluster. This factor is called the frailty and accounts for the correlation of the times within a cluster (e.g., the correlation of survival times of lambs from the same sire). The resulting model is called a shared frailty model.

A joint model for repeated measurements and time-to-event
.., Y in i ) be the vector of the repeated measurements for the ith individual measured at times t T i = (t i1 , . . . , t in i ). Let T o i = min(T i , C i ) and δ i denote, respectively, the observed time and censoring indicator for the ith individual. The observed data available for the ith individual are thus , where X 1i denotes the matrix of the observed values of covariates believed to influence the repeated measurements Y i , while x 2i is a vector of the observed values of covariates believed to affect the time-to-event process.
Henderson et al. [9] have proposed a model for the joint analysis of both the time-to-event and repeated measurements. They postulate a latent (unobserved) bivariate Gaussian process W i (t) = {W 1i (t), W 2i (t)} such that the repeated measurements process depends on W 1i (t) and the event time process depends on W 2i (t). In particular, for the repeated measurements process, consider a model of the general form is the systematic component, which can be described by a linear model, e.g., µ i (t i ) = X 1i β L . As a basic example for the latent process W 1i (t), Henderson et al. [9] consider where (U 1i , U 2i ) is a bivariate normal random vector with zero mean and variance-covariance On the contrary, the time-to-event is modelled through a Cox proportional hazard model It is assumed that the repeated measurement and time-to-event processes are conditionally independent given W i (t). However, in order to induce association between the two processes, W 2i (t) is taken to be related to particular components of W 1i (t). This is achieved via the general equation For example, a joint model with W 2i (t) = γ 1 U 1i + γ 2 U 2i , would allow both the random intercept U 1i and slope U 2i , involved in (5), to affect the risk of the event.
The parameters of the models for the repeated measurement process and the time-to-event process are estimated jointly by maximising the observed joint likelihood of the data, as described in the references [23,24].

APPLICATION
We now describe the application of the joint model described in the previous section to the data introduced in Section 2. Separate analyses of the repeated measurements of PCV, FEC and BWT were performed. Survival times of lambs that survived beyond one year, or those of lambs that were stolen, were censored at one year and at the last recorded observation, respectively.
Nguti et al. [16] reported an average lamb mortality of 19% in the preweaning period and 31% in the post-weaning period. The age at death during the post-weaning period ranged from 3 to 12 months (median 6.4 months). The number of repeated measurements recorded from weaning ranged from 1 to 8 (median 6) per lamb with 1994 having the most post-weaning measurements. Variable patterns of each trait were observed over the one year as well as across the six years. For example, Figure 1 shows the scatter plot of the measurements recorded from one to 12 months for PCV. In the plot, individual profiles for a randomly selected sample of 15 lambs are highlighted. Although all animals were weighed and sampled on the same day, the ages varied as a result of the lambs being born within a period of about 20−40 days.
In the joint models with either PCV, BWT or FEC as repeated measurements, the fixed effects of genotype (6 levels), year of birth (6 levels) and sex (2 levels) were included in the repeated measurements component of the joint model. As suggested by Figure 1, PCV was assumed to be curvilinear over time. This curvilinear trend over time was also assumed for FEC and BWT (figures not shown). The age of dam (5 levels) was considered as a baseline covariate for the time-to-event component only, but not for the repeated measurements, where it was not found to be significant.
Consequently, we can define β T . . , 5) are the binary indicators capturing the breed effects, ς k (k = 1, . . . , 5) are the indicators for year of birth effects and l is the binary indicator for males. As a result, the repeated measurements model can be written as where X 1i is the n i × 12 design matrix corresponding to β L , (η 1 , η 2 ) are the parameters associated with the time trend, t * T i = (t 2 i1 , . . . , t 2 in i ) is the vector of the quadratic times, and W 1i ( (8) can then be re-written as To specify the survival component of the joint model, , where a r (r = 1, . . . , 4) are the binary indicators coding the dam age groups (with levels: ≤ 2 years, 3, 4, 5 and ≥ 6 years). The model for survival time is then given by where x 2i is the n i × 15 design matrix associated with β S . For all three traits the following settings for W 1i and W 2i were considered: Settings (S1) and (S2) assume independence between the repeated measurement and survival processes. Settings (S3) and (S4) correspond to (S1) and (S2), respectively, with respect to the structure of W 1i (t), but allow for dependence between the processes (joint models).
To obtain parameter estimates for the fixed effects, variance components and the association parameters of the joint models (S3) and (S4) specified above, a programme was written in SAS  . Estimates from either setting (S1) and (S2) were computed using PROC MIXED (for repeated measurements) and PROC PHREG (for survival time) in SAS  .
Estimates of the standard errors for all parameter estimates in the joint models were obtained by using the jackknife method. This was achieved by leaving out the observations for lambs from the same sire and then re-fitting the model to the remaining observations. Classically, jackknife provides reliable estimates of the standard errors if the observations omitted are independent from those that are left in. When observations for lambs from the same sire are left out, so is the genetic component of these lambs. This genetic component is assumed here to provide more individual contribution to the lamb characteristics (e.g. survival, BWT) than the common environmental components, which are shared by lambs born in the same year. This is supported by the findings in reference [2]. These authors show that for the analysed data set, the differences observed in heritability estimates of PCV and FEC for Dorper-compared with Red Maasai-sired lambs were more likely due to the differences in genetic variance rather than in environmental variability.

RESULTS
In this section the results of the fitted models for PCV, BWT and FEC for settings (S1)-(S2) are reported. For each trait the results of the joint model are compared with those of the corresponding independence model for both the repeated measurements and survival estimates.

Packed cell volume from one month
Initially, the analysis of the repeated measurements for PCV from one month until one year of age was considered. When fitting the model corresponding to setting (S2), a non-positive definite estimate of the variance-covariance matrix G 1 (see Sect. 3.3) was obtained. On further investigation, it was discovered that the PCV repeated measurements were negatively correlated, with the correlation increasing in absolute value over time. This negative serial correlation cannot be captured by a model with a random intercept and random slope, as specified under setting (S2). In particular, by taking into account the magnitude of the correlation and low variability of the slopes (g 22 in Eq. (2)) of individual profiles, the variance of the random slope is estimated to be less than or equal to 0, which is obviously an inadmissible value. Therefore, for the PCV measurements collected from one month to one year of age, the models for settings (S2) and (S4) could not be fitted. On the contrary, the models for settings (S1) and (S3), which do not account for the negative correlation, could be used. The results for these models are given in Table I. Strictly speaking, one should treat the results with caution since they are based on models with a possibly misspecified variance-covariance structure.
Setting (S1): Under this setting, which assumes independence between PCV measurements and survival time, the Dorper (D×D) breed had the lowest mean PCV from one month to one year of age, which was between 0.1 to 1.9% units lower than for other genotypes (see Tab. I). This difference increased as the Red Maasai genotype in the lambs increased, with the Red Maasai having the highest mean PCV. The linear and quadratic time effects were both significant (P < 0.001) implying an average non-linear trend in PCV. The trend is as indicated in Figure 1, which shows a general sharp decline in PCV after one month in all years except 1996 followed by a slight rise. The lambs born in 1992-1995 had on average lower PCV (0.3 to 3.2%) units than those born in 1991. The mean PCV was the highest in 1996. On average, male lambs had lower PCV than female lambs.
By exponentiating the estimates given in the 'Survival model -S1' column in Table I, one can see that, as compared to the Dorper, the relative Table I. Estimates from settings S1 (W 1i (t) = U 1i , W 2i (t) = 0) and S3 (W 1i (t) = U 1i , W 2i (t) = γW 1i (t)) for repeated measurements of PCV(%) from one month to 12 months and survival.

Repeated measurements model
Survival model S1 S3 S1 S3 est s.e. est s.e. est s.e. est s.e. mortality hazard of the other genotypes ranged from exp(−0.476) = 0.62 to exp(−1.34) = 0.26. The R × (R × D) and R × R breeds had the lowest, and similar, mortality. The hazard of the lambs born in the years 1993-1996 was statistically significantly higher than that of the lambs born in 1991 and ranged from 2.2 to 4.0. Male lambs had a higher mortality hazard than females while the hazard ratio decreased with increasing age of the dam. Setting (S3): This setting corresponds to (S1), but assumes dependence between PCV measurements and survival time. As compared with (S1) the differences in the mean PCV, relative to the Dorper breed, increased slightly for all other genotypes. For instance, the estimated mean PCV for the non-Dorper genotypes was 0.2−2.2% units higher than for the Dorper breed (Tab. I). This increase might be the result of the adjustment of the analysis of the repeated measurements for the variation in death rates. The estimated time trend parameters for the repeated measurements model for the (S1) setting were similar to those obtained for (S3).
Relative to the mortality hazard for the Dorper breed, the hazard ratio for the non-Dorper genotypes now ranged between 0.60 to 0.24, as compared to the (S1) setting (Tab. I). Significant negative estimates (P < 0.001) were obtained for the association parameters (γ in Tab. I) for the survival model under (S3). This indicates that the mortality hazard decreased with increasing PCV.

Packed cell volume from weaning
Since the critical period for assessing genetic resistance to endoparasites in lambs is between weaning and 12 months of age [1], the analysis of the PCV repeated measurements for the period from weaning onwards was also considered. In this analysis, the survival time was redefined by using weaning as the time of origin. Consequently in this analysis only the animals alive at the time of weaning were considered. Models under settings (S1) to (S4) were considered. The results for settings (S1) and (S3) (not shown) were similar to those reported for the analysis of data from one month of age (Tab. I). However, unlike the period from one month of age, the time trend had a more moderate negative slope estimate. This corresponds to Figure 1, which shows a more gradual decline in PCV after weaning than before. Furthermore, the relative mortality hazard in the post-weaning period exhibited similar patterns as in the analysis of data from one month of age, but now the R × (R × D) had the lowest mortality (exp(−1.641) = 0.19) when compared to the Dorper breed. A detailed description of the results from (S2) and (S4) (Tab. II) is given below.

Settings (S2):
In the repeated measurements model with both random intercepts and slopes, (Tab. II), similar trends for the fixed effects parameter estimates were observed as in the simpler, random-intercept-only model. However the ranges of the estimates were reduced. For example for setting (S1), the genotype parameter ± s.e. estimates were the following: 0.361 ± 0.287 (D × (D × R) The random intercept and slope were negatively correlated (σ 12 = −1.82, see Tab. II). This implies that lambs with a high PCV at weaning had a more rapid decline in PCV than those with a low PCV. The estimated variance component for the random intercept (σ 2 1 ) was 15.89 (s.e. 0.99) compared to 6.47 (s.e. 0.43) in the simpler model (S1) (results not shown). Note, however, that, according to equation (2), the negative correlation between random intercept and slope will reduce the total variability under the (S2) setting, making it similar to the sum of σ 2 1 and σ 2 e in the simpler (S1) model. It is also worth noting that, according to equation (2), under (S2) the covariance between two PCV measurements depends on time, and can be negative if the correlation between random intercept and slope is negative. Since it is possible that the correlation between PCV measurements may change with time and may be negative (see remarks in Sect. 5.1), this gives more credibility to the results obtained for (S2) because the model used in an (S1) settting fimposes a constant positive correlation between PCV measurements. Setting (S4): As compared with (S1) and (S2) settings, the differences in the mean PCV for all other genotypes, as compared with the Dorper breed, increased slightly in the joint models constructed under both (S3) (results not shown) and (S4) (Tab. II) settings. For instance, for the (S3) setting, the estimated mean PCV from weaning for the non-Dorper genotypes was 0.6−3.4% units higher than for the Dorper breed. For the (S4) setting (Tab. II), the difference was between 0.3% and 3.0%. This increase might be attributed to the adjustment of the analysis of the repeated measurements for the variation in death rates.
Relative to the mortality hazard for the Dorper breed, the hazard ratio for the non-Dorper genotypes now ranged from 0.55 to 0.16 for the (S4) setting (Tab. II). These estimates were lower than the estimates obtained for the corresponding (S1) and (S2) survival models. For (S3) setting the mortality hazard for lambs born in 1993-1995 was about five times higher than in 1991   (with parameter ± s.e. estimates as: 1.637±0.209 (1993), 1.640±0.252 (1994), 1.677±0.250 (1995)). For (S4) the ratio was about five for 1993, while for 1994 and 1995 it was about four times, a lower ratio than for (S3). This is supported by Figure 1. The lambs born in 1994 and 1995 had much lower PCV measurements at weaning than those born in 1991, and PCV has been shown to correlate with survival. However for the former, PCV increased over time while those for the latter decreased. Thus, the adjustment for the PCV evolution over time results in a slight decrease in the mortality hazard for 1994-1995. On the contrary, the lambs born in 1993 and 1991 had almost similar PCV measurements at weaning. The decrease in 1993 over time was, however, much sharper (larger negative slope) than in 1991. Adjusting for this decrease translates into a higher mortality hazard for 1993. Finally, higher hazard ratios are observed for 1994 and 1995 in (S3) than (S4), since the latter model adjusts the risk only for the overall level of PCV over time (without adjusting for the rate of change). Significant negative estimates (P < 0.001) were obtained for all the association parameters:γ = −0.303 ± 0.0203 for setting (S3), and γ 1 − γ 3 as in Table II for setting (S4). This indicates that the mortality hazard decreased with increasing PCV. Thus, after weaning in the (S3) setting, lambs with PCV measurements higher than the average had a lower mortality hazard than those with lower PCV measurements. The standard deviation of the distribution of the random intercepts in the repeated measurements part of the joint model for the (S3) setting was estimated to be equal to 2.51. Thus, the model predicts that for a (random) increase of PCV due to an increase of the random intercept W 1i = U 1i by one standard deviation, the risk of death decreases by exp(−γ * W 1i ) = exp(−0.303 * 2.51) = 0.47. For the (S4) model, a large negative estimate was obtained for γ 2 (= −1.92), which corresponds to the random individual slope. Thus, the model indicates that at any time lambs that had a large decrease in PCV had an increased risk of death.

Body weight
Measurements of body weight from birth to one year of age were used in this analysis. The models were fitted using all four settings. The parameter estimates of the fixed effects were similar for settings (S1)-(S3) and (S2)-(S4). We thus report only the estimates of the models constructed under settings (S2) and (S4) which are shown in Table III.  Table III. Estimates from settings S2 (W 1i (t) = U 1i +U 2i t, W 2i (t) = 0) and S4 (W 1i (t) = U 1i + U 2i t, W 2i (t) = γ 1 U 1i + γ 2 U 2i + γ 3 W 1i (t)) for repeated measurements of BWT(kg) from birth to 12 months and survival.

Setting (S2):
The Dorper (D × D) breed had the highest mean BWT, which was between 0.04 and 0.72 kg higher than for the other genotypes (Tab. III). There was a non-linear trend in the change of body weight over time. On average, the lambs born in 1994-1996 were lighter than those born in 1991. Male lambs were on average 0.11 kg heavier than the females.
Compared with the Dorper, the relative mortality hazard of the other genotypes ranged from 0.61 to 0.27, with the Red Maasai (R × R) and R × (R × D) having the lowest hazard. As for the previous PCV analyses, an increased mortality hazard was noted for the years 1993-1996. Lambs born to ewes ≥ 3 years of age had lower hazard than those born to younger ewes.
Setting (S4): Adjusting the repeated measurement process for the variation in death rates only had a slight effect on the parameter estimates of the (S4) models when compared to setting (S2).
The relative mortality of the genotypes now ranged from 0.61 to 0.24 with the R×R genotype having the lowest hazard mortality. This result indicates that despite being lighter in body weight when compared to the other genotypes, the Red Maasai demonstrated better performance in terms of survival. The age of the dam effect was non-significant. This could be due to the fact that in this analysis we accounted for the low body weight of lambs born to young dams, which biologically may be due to low milk production of the dam in her first parity.
Negative estimates of the parameters relating the random components of the repeated measurements model to the survival model were observed (see the estimates for γ 1 − γ 3 in Tab. III). In particular, there was a significant negative association only between random growth rate (γ 2 (= −3.28), P < 0.001) and risk of death. This shows that the animals that had weight profiles with increasing slope had a reduced risk of death. In a reduced model with γ 1 and γ 3 constrained to zero (results not shown), the estimate obtained for the association parameter γ 2 wasγ 2 = −3.262 (s.e. = 0.389). The standard deviation of the random slope in this reduced model was equal to 0.22. Thus for every increase by one standard deviation, the mortality hazard associated with the rate of change in BWT was reduced by exp(−3.262 × 0.22) = 0.49.

Faecal egg count
The repeated measurements for FEC were log-transformed into LFEC = ln(FEC + 25) since the data were highly skewed. Initially, an analysis of the repeated measurements from one month to one year of age was considered. As in the case with PCV, a negative correlation between repeated measurements, that increased over time, was observed. This problem precluded the estimation of any of the settings (S1)-(S4). Consideration of only the LFEC measurements from weaning onwards did not resolve the problem, since the negative correlation continued to persist among measurements collected towards the end of the one year period. The problem may be due to strong oscillations of the individual patterns that were observed in the LFEC measurements from one month onward, as well as from the time of weaning, that could be linked with the treatment protocol instigated in the study. A possible solution might be to extend the random structure (W 1i ) of the repeated measurements model by including a serial correlation component in the model for the repeated measurements process. Owing to the limitation of software availability, this modelling option was not attempted.

DISCUSSION
The data discussed in this paper were previously analysed by Baker et al. [2] and Nguti et al. [16]. Baker et al. [2] analysed the repeated measurements of BWT, PCV and FEC without taking into account the survival pattern of the animals, and chose to analyse the data for each time point separately. They confirmed the higher resistance (lower FEC) and higher resilience (higher PCV) of Red Maasai than the Dorpers. Nguti et al. [16], studied the survival of each genotype and introduced the effects of BWT, PCV and FEC as time-varying covariates in a shared frailty model, with the frailty defined as a random effect of sire. Introduction of PCV and FEC as time-varying covariates in that analysis in models with BWT (time-invariant or time-varying) reduced the magnitude of the sire variance, confirming the moderate levels of heritability reported by Baker et al. [2].
In the analysis presented in the current paper, individual repeated measurements were analysed jointly with the survival process. By doing so, parameter estimates in both components of the joint model in general increased in absolute order of magnitude, as compared with the models assuming independence between the two processes. For instance, in the joint model with PCV as the repeated measurement, the range of the mortality hazard ratios for different genotypes relative to the Dorper changed from 0.61−0.19 for (S1) and (S2) models to 0.55−0.16 for (S3) and (S4) models. Thus, the adjustment for the evolution of PCV resulted in a clearer separation between the Dorpers and the Red Maasai. On the contrary, adjustment for the death rates widened the ranges of the parameters reflecting the differences among the various genotypes in PCV measurements. As shown by Baker et al. [2] the PCV values among the genotypes were similar at the first sampling at one month, prior to any sign of the disease. The subsequent reductions in PCV were inversely correlated with the levels of FEC, suggesting that the differences in overall mean PCV values among the genotypes were more likely the result of the different levels of resistance between the Dorpers and the Red Maasai, rather than any genetic differences between breeds in their normal PCV values.
In general, repeated measurements such as PCV, BWT or FEC are only recorded at specific time points. When such variables are used in a proportional hazards model as time-varying covariates, the standard method is to impute the missing observations by using the last observed value. This results in a piecewise constant profile. This approach has been shown to lead to biased model parameters [18], with the presence of measurement error in the covariate attenuating the estimates towards zero. On the contrary, in the joint analysis, the repeated measurements are imputed by the 'true' values predicted under the model for the repeated measurement process. This reduces the attenuation and can explain the increase in the magnitude of parameter estimates observed for the joint models. In addition, the joint modelling allows other characteristics of the repeated measurements pattern, such as the rate of change (slope), to influence the risk of death. The effects of these characteristics is well demonstrated in this study when repeated measurements of PCV recorded from weaning were considered. Including a random slope had a major effect on the postweaning risks of death for the years 1993 to 1995 when compared with 1991.
Nguti [17] compares the results from the survival component of the joint model under setting (S4) for the repeated measurements of PCV (from weaning) and BWT (from birth) together with those for a Cox PH model with a time-varying covariate for these measurements. Differences in hazard estimates were observed by the two approaches. However, it would be erroneous to make direct comparisons of the two sets of results since the methodological aspects of the two models are different. The survival component of the joint model can be viewed as a conditional model (conditioned on the random effects which may be partially confounded with genotype), while a Cox PH model with a time-dependent covariate can be viewed as a population average model.
The aforementioned potential of the joint models for a more accurate analysis of the repeated measurements and survival processes is an attractive feature of the methodology. However, it might be that the models do not have the capacity to fit the actual data as well as with the FEC measurements. In such a case a careful reflection on the source of this lack-of-fit should result in proposing a change in the structure of the model. For instance, it should be noted that the post-weaning measurements were obtained at the time lambs were diagnosed for treatment based on the observations of the sentinel group of 50 lambs (see Sect. 2). Treatment will have resulted in decreases in FEC especially, and also increases in PCV. Thus data collected in between treatments, when measurements were more likely to have been within normal limits, were missing. These regular treatments are likely to have influenced the covariance structure resulting in the observed negative correlations. A similar problem occurred pre-weaning when individual animals were treated on the basis of PCV and FEC. It should also be noted here that not all lambs that died did so due to helminthiasis. Although Nguti [16] found a difference in survival rate between Dorper and Red Maasai regardless of whether the deaths due to causes other than helminthiasis were omitted or not, the decreases in PCV and increases in FEC, that are typical of the onset of helminthiasis, may not have been reflected in deaths due to other causes. More frequent sampling and use of the data collected mid-way between treatment interventions may have both helped in better understanding the patterns of changes in PCV and FEC, thus overcoming the problem of a negative correlation.
In the application of any statistical method one needs to be aware of the assumptions made. For example, in the joint models considered in this paper, the models used a bivariate Gaussian process to induce the association between the repeated measurements and the survival process. This, in turn, imposes a particular, hierarchical structure (e.g., random intercepts and slopes) of the model for the repeated measurements process. Such a structure may not necessarily be adequate for a particular dataset. As already noted, the models with random intercepts and random slopes were not able to capture the intrinsic patterns of FEC induced by the treatment interventions. On the contrary, they coped better with the BWT profiles from the time of birth. In general, random intercepts and slopes provide a representation of the dominant part of the evolution of the profiles, but do not capture a more subtle behaviour (e.g. short-term oscillations around an average pattern). Such a behaviour might be captured by using an autocorrelated stochastic process. Henderson et al. [9] has proposed the use of a non-stationary Gaussian process in their approach. Unfortunately, due to the lack of appropriate software, this solution is not yet available in practice.
The joint models applied in this paper did not allow for adjusting of correlations among survival times of lambs from the same sire. For this purpose, Nguti et al. [16] used frailty models (see Sect. 3.2) with a random sire effect. In theory, the models formulated by Henderson et al. [9] allow for the inclusion of a frailty term in the time-to-event component of the joint model, but implementation needs to await the development of appropriate software.

CONCLUSION
Classically, the effect of repeated measurements on the risk of death is assessed by treating the repeated measurement as a time-dependent covariate in a survival model. Estimates from such an analysis have been shown to be biased towards zero, thus showing over-estimated hazard ratios. The bias can be removed by using joint models. This paper illustrates the application of such models. Accounting for the survival (drop-out) mechanism was shown to affect the parameter estimates related to the repeated measurement process.