A Bayesian generalized random regression model for estimating heritability using overdispersed count data

Background Faecal egg counts are a common indicator of nematode infection and since it is a heritable trait, it provides a marker for selective breeding. However, since resistance to disease changes as the adaptive immune system develops, quantifying temporal changes in heritability could help improve selective breeding programs. Faecal egg counts can be extremely skewed and difficult to handle statistically. Therefore, previous heritability analyses have log transformed faecal egg counts to estimate heritability on a latent scale. However, such transformations may not always be appropriate. In addition, analyses of faecal egg counts have typically used univariate rather than multivariate analyses such as random regression that are appropriate when traits are correlated. We present a method for estimating the heritability of untransformed faecal egg counts over the grazing season using random regression. Results Replicating standard univariate analyses, we showed the dependence of heritability estimates on choice of transformation. Then, using a multitrait model, we exposed temporal correlations, highlighting the need for a random regression approach. Since random regression can sometimes involve the estimation of more parameters than observations or result in computationally intractable problems, we chose to investigate reduced rank random regression. Using standard software (WOMBAT), we discuss the estimation of variance components for log transformed data using both full and reduced rank analyses. Then, we modelled the untransformed data assuming it to be negative binomially distributed and used Metropolis Hastings to fit a generalized reduced rank random regression model with an additive genetic, permanent environmental and maternal effect. These three variance components explained more than 80 % of the total phenotypic variation, whereas the variance components for the log transformed data accounted for considerably less. The heritability, on a link scale, increased from around 0.25 at the beginning of the grazing season to around 0.4 at the end. Conclusions Random regressions are a useful tool for quantifying sources of variation across time. Our MCMC (Markov chain Monte Carlo) algorithm provides a flexible approach to fitting random regression models to non-normal data. Here we applied the algorithm to negative binomially distributed faecal egg count data, but this method is readily applicable to other types of overdispersed data.


Background
Faecal egg count is a commonly used indicator of susceptibility to gastrointestinal nematode infection and provides a marker for selective breeding programs. Infection has traditionally been controlled with anthelmintic drugs, but resistance to these drugs has directed attention towards selective breeding as a sustainable and viable alternative [1]. Selective breeding programs rely on estimates of the animals' breeding values, which represent the sum of the additive effect of the genes received from both parents [2]. Thus, designing an effective selective breeding scheme requires accurate assessment of the heritability.
The analysis of faecal egg count data is challenging for two reasons. First, the data are overdispersed, which has led previous studies to use transformed faecal egg counts [3,4]. Transformations of faecal egg count data (commonly a log transformation) can result in bimodal data [5] and therefore may not be appropriate [6,7]. However, any transformation can be avoided by modelling the raw faecal egg counts as a negative binomial distribution [3]. The second challenge is that, as the adaptive immune response develops, the sources of variation and the heritability of faecal egg counts are expected to change over time, which suggests that a multivariate approach may be appropriate. The goal of this paper was to estimate the change in heritability of faecal egg count over the grazing season.
Random regression models are commonly used to model changes in quantitative traits measured over a continuous scale such as time or age [8]. In particular, they can be used to estimate changes in genetic and environmental variance components as continuous functions over a time frame by specifying time-dependent functions φ 1 , . . . , φ k [9]. However, these models can involve the estimation of a large number of parameters that may exceed the number of observations and become computationally intractable, which prompts the use of reduced rank random regression [8]. By estimating each covariance matrix using relatively few principal components, or eigenfunctions, the number of parameters to be estimated can be significantly reduced [8]. Then, the reduced rank random regression models estimate continuous covariance functions using a small number of the largest eigenvalues [10].
Random regression models have been widely used to estimate genetic parameters of repeated measurements over time [11] and, previously, Bayesian methods have been used to capture the skewed distribution of faecal egg count data [12,13], but because of the methodological challenges, to date, these two approaches have not been combined.
Software is available to fit random regression models, such as WOMBAT [14], ASReml [15] and RRGIBBS [16]. Each program can be used to estimate variance components. WOMBAT fits linear mixed models, of full or reduced rank, through restricted maximum likelihood (REML). ASReml can be used to fit generalised linear mixed models and RRGIBBS fits random regression models using a Gibbs sampler. However, none of these packages can conduct random regression analyses of generalised mixed models that are also reduced rank. For these reasons, we ultimately took a Bayesian approach by fitting a multivariate negative binomial reduced rank random regression model using a Metropolis Hastings algorithm implemented in R [17].
This paper deals with heritability estimates on two scales. Using faecal egg count data collected from five consecutive cohorts of 200 Scottish Blackface lambs, our overall goal was to provide an alternative method to estimate the heritability of faecal egg count on the link scale and identify the short-comings of estimating heritability on a latent scale by giving specific attention to log transformations [7,18,19].
The specific steps in this work were to: (1) demonstrate, using our data, the short-comings of the approaches that are commonly used to handle faecal egg count data, in particular, the use of univariate analyses when the data exhibit correlations that are significantly different from zero between variance components over time, and the use of a log transformation when the data follow a negative binomial distribution; (2) we used standard software (WOMBAT) to demonstrate and compare full and reduced rank analyses of the log transformed data; (3) we applied a reduced rank random regression approach, which assumes that the data are negative binomially distributed. Using residual diagnostics, we showed that the random regression model using the untransformed data provided a better model than the random regression based on the log transformed data. Heritability estimates from the random regression model using the untransformed data were similar to those from the univariate analyses in the later months; however, in the early months, the random regression model explained a much greater proportion of the variance and resulted in substantially higher values of the heritability.

Data source
Data were collected over consecutive years from five cohorts of 200 straight-bred Scottish Blackface lambs between 1992 and 1996 [4]. Within each year, lambs were born over a two-week period and monthly faecal egg counts were estimated using standard procedures between May and October with additional post-mortem samples. As is common practice on this farm, immediately after collection of the faecal samples each month, the lambs were given the recommended dose of anthelmintic.
Since eggs are counted in a 1 50 th of a gram of faeces, each egg counted represents 50 eggs per gram [20]. Following standard practice, we use the term faecal egg count (FEC) for eggs per gram of faeces, and subsequently use the term raw egg count to refer to the distribution of counts per 1 50 th gram.
Seven categories of adult nematodes were detected with variable prevalences across the five years. Teladorsagia circumcincta was identified for nearly all the lambs analysed [21]. Across the five years, 36 rams and 485 dams were examined.
The data are unbalanced with varying degrees of missingness between months and as much as two-thirds of the data missing in May (Table 1). Post-mortem counts were not measured in the last year of the study and most of the data were collected from male lambs. Due to the unbalanced nature of the data and level of missingness, the first and last time points were removed from all multivariate analyses in this paper.

Modelling the data
Quantitative traits are often assumed to be normally distributed [22,23]. Consequently, it has become common practice to fit a normal distribution to log transformed faecal egg counts (FEC) of the form log(FEC + c) for some constant c [4,24] or to Box-Cox transform faecal egg counts [25]. However, such transformations can result in bimodal data [5].
Transformation of faecal egg count data may not be necessary [6,7,26] since such count data often follow a negative binomial distribution [3]. The negative binomial distribution is parameterised by the arithmetic mean μ and a positive exponent r [27] and has also been shown to be a good fit to the distribution of faecal egg counts for a wide variety of parasites [5,28]. In this parameterisation, the variance of the distribution is μ(1 + μ r ) and thus approaches a Poisson distribution as the dispersion parameter r increases and so smaller values reflect more dispersion.
We began our analyses of these data by fitting a series of univariate nested half-sib design mixed models to log transformed faecal egg counts and untransformed raw faecal egg counts. We then used a multitrait animal model to quantify temporal correlations, which prompted us to use a random regression model. We used standard software (WOMBAT) to implement full and reduced rank analyses of the log transformed data and finally applied our negative binomial reduced rank random regression model.

Univariate analysis
In previous studies, heritabilities of log(FEC + 1) and log(FEC + 25) transformations of these data were estimated by treating each month independently [4,24]. We examined the use of log transformed faecal egg count data by comparing univariate heritability estimates of log(FEC + c) for a range of values of c over the sevenmonth period. We assumed that the log transformed data were normally distributed. For each model, we fitted year of birth and sex as fixed effects. The nlme package in R [29] was used to estimate heritabilities on the latent scale using nested half-sib design mixed models. This model takes the form: where Y i j k is an observation from the i-th lamb, S k is the random effect of the k-th sire and D j k is the random effect of the j-th dam within the k-th sire. Coefficients of fixed effect are denoted by θ θ θ, X i is the corresponding design matrix and i j k is the residual variance. We used similar univariate nested half-sib design mixed models to estimate month by month heritabilities for the raw egg count data. The residual variation in these models was estimated using the method described by Tempelman et al. [30]. The R package glmmadmb [31] was used to fit negative binomial mixed effects models.
The MCMCglmm (MCMC generalised linear mixed models) package [32] was used to fit transformed data (log(FEC + 1)) to a multitrait animal model to estimate correlations between June and October. Using this method, we estimated pairwise phenotypic, maternal and genetic correlations between the five months simultaneously. Due to the large number of parameters estimated, we compared these results to those that were obtained using a series of bivariate animal models as an informal check of convergence [33]. Any major deviations between the two sets of models would indicate that the chain had not fully converged.
Each model was run for 1 000 000 iterations with a burnin period of 50 000 iterations. Inverse-Wishart priors were used for each covariance matrix. We constructed these priors such that faecal egg counts accross months were a priori independent with the total variance being spread evenly across months [32].

Multivariate analysis
A log(FEC + 1) transformation was used for the multivariate analyses of the transformed faecal egg count. This transformation was chosen for consistency with previous multivariate analyses of these data [4].
Random regression models are mixed effect models with individual functions of continuous covariates fitted as random effects and are commonly fitted by specifying time-dependent basis functions φ 1 (t), . . . , φ K (t) [8]. In this paper, we aimed at modelling individual changes in additive genetic, maternal and permanent environmental effects over time. Time refers to month and we assumed that all fixed effects remained constant across months. Sex and year of birth were included as fixed effects in all models. For where T denotes the number of time points), a random regression model takes the form: The matrix contains the set of basis functions evaluated across the time frame. The vector θ θ θ contains fixed effect regression coefficients and X i is the corresponding data matrix for lamb i.
are individual regression coefficients relating to the additive genetic, permanent environmental and maternal random effects. Note that the number of basis functions can differ between random effects. Addition residual variance is denoted by vector Since any covariance matrix can be decomposed as = EλE T EλE T EλE T , where λ λ λ is a diagonal matrix of ordered eigenvalues and E E E an orthonormal matrix of eigenvectors, the matrix can be approximated by setting the smallest eigenvalues (in the matrix λ λ λ) to zero. Using this decomposition, the random regression model can be written as: This transformation forms the basis of a reduced random random regression model, however a full description is in Meyer and Kirkpartick [8].

Random regression using WOMBAT
The implications of a reduced rank analysis of these data were assessed using full and reduced rank random regression models in WOMBAT [14]. This program assumes normally distributed traits and, thus, we used transformed faecal egg count data (log(FEC + 1)).
The number of time-dependent basis functions and the number of eigenfunctions are specified by the user. We used K to denote the number of basis functions and M to denote the number of eigenfunctions used in any reduced rank analyses [8]. Legendre polynomials were used as basis functions. The first three Legendre polynomials are [34]: We began with a second order analysis, using two Legendre polynomials (K = 2) for each variance component. A first order analysis was not considered since including only the first Legendre polynomial would not allow variance components to be functions of time. In total, we considered two, three and four Legendre polynomials (K = 2, 3 or 4) in full rank analyses (M = K) and then considered the same range of Legendre polynomials in reduced rank analyses, each with two eigenfunctions (M = 2). The results with K = 2 and M = 1 are also presented.

Negative binomial random regression model
We describe a method to estimate heritability using raw faecal egg count data multivariately using a negative binomial random regression model. For lamb i at time t, the negative binomial model with mean μ it and dispersion r t can be parameterised as: [35]. The dispersion of this negative binomial distribution is controlled by the value of r t with smaller values producing a higher variance [36]. Using a log link function, we set: The vector θ θ θ contains the fixed effect regression coefficients corresponding to the overall mean, the sex of the lamb and the year the lamb was born. The matrix X i is the corresponding data matrix. The matrix is of dimension 5 × K, where K is the number of Legendre polynomials functions used and 5 is the number of time points in this study. The tk-th entry of the matrix is φ k (t).
Three variance components are included: an additive genetic, permanent environmental and a maternal effect.
to the maternal component. The suffix dam i is used to denote the dam corresponding to lamb i and n d is the total number of dams. M is the number of eigenfunctions included in the reduced rank model. Here, we included the same number of polynomial functions and eigenvalues for each of the three variance components. However, the model can easily be extended to include polynomials of different degrees and values of M.
Estimated maternal effects can depend on the number of progeny per dam, the number of dams with recorded egg counts and the number of generations of recorded data [37]. In this dataset, there were rarely more than two progeny per dam and for each lamb, only the sire and dam are known. The structure of this dataset is a source of difficulty in distinguishing between maternal additive genetic and maternal environmental effects [4,37]. Therefore, the maternal effect modelled here may include both genetic and environmental effects. However, the model is easily adapted to specifically model genetic and environment effects [38].
Coefficients γ γ γ = {γ 1 γ 1 γ 1 , . . . , γ L γ L γ L } are assumed independent between lambs: Matrices I L I L I L and I n d I n d I n d are identity matrices of dimension L × L and n d × n d respectively.
The values λ α 1 , . . . , λ α M are the M largest eigenvalues of the additive genetic covariance matrix, and the matrix E α M E α M E α M contains the corresponding eigenvectors. Similarly, the values λ γ 1 , . . . , λ γ M and λ β 1 , . . . , λ β M are the M largest eigenvalues of the permanent environmental and maternal covariance matrices respectively, with the corresponding eigenvectors stored in Legendre polynomials are orthogonal [40] and subsequently the matrices and corresponding eigenfunctions are constrained to be orthogonal [41].
All estimates presented in this paper are based on posterior median values from the last 20 000 samples based on 1 000 000 iterations.

Specifying prior distributions
Inference was based on Markov chain Monte Carlo simulation and the parameters were updated with a Metropolis-Hastings algorithm. Prior distributions were specified for θ θ θ, r r r = {r 1 , . . . , We chose to use non-informative priors for regression coefficients relating to fixed effects and incorporated shrinkage priors for variance components [42][43][44]. Although other priors could naturally be incorporated into this type of model, results presented in this paper are based on the priors described in this section.
Normal priors were used for each element of θ θ θ, that is: Using a large variance ensures that the prior distribution is flat over a wide range of values and is centered around zero.
The likelihood function can be penalized to improve estimation of covariance matrices and variance components by shrinking sample eigenvalues towards a given value [42], a property that can be enforced through the prior distribution. We used normal priors for log(λ α m ) [43,45]: Similar prior distributions were used for log(λ γ m ) and log(β α m ). In this case, we chose to shrink our eigenvalues towards zero with little variation, setting μ λ = 0 and σ λ =5. Similarly, variance components can be shrunk towards a value by setting flat gamma priors on corresponding precision parameters [44]. We used gamma priors for each of the five dispersion parameters: . The values in the matrix are used to form a matrix of polar coordinates and then transformed to form an orthogonal matrix, E E E, using a Gram-Schmidt transformation. Following Martinez et al. [45], we set

Estimating variance components
To estimate variance components as functions of time, we followed the model of Meyer and Kirkpartick [8] and set: for times t 1 and t 2 . The functionĜ is the additive genetic covariance function,P is the permanent environmental covariance function andM the maternal covariance function. Given these functions, the heritability of raw faecal egg counts can be estimated on a continuous time scale. Specifically, at time t, the heritability of raw faecal egg counts can be estimated [25] as: whereˆ 2 e t is the estimated residual variance at time t on a continuous scale. The residual variance in a negative binomial regression can be estimated using the methods described by Tempelman et al. [30]. The negative binomial model is an extension to a Poisson model and there are two sources of residual variance. There is variation in the Poisson sampling and some additional variation due to the overdispersion that is not captured in the Poisson model [7,19]. Following these methods, at observation time t, we approximated the residual variation from the Poisson sampling byλ t −1 and estimated the additional variance as: We assumed the residual covariance matrix 2 2 2 to be a 5 × 5 diagonal matrix with the t-th diagonal entry equal toσ 2 e t = φ (1) (r t ) +λ −1 t . Given the estimated covariance matrix, a continuous residual variance functionˆ 2 e t was estimated by the methods described by Kirkpatrick et al. [34]. For consistency, we used the same number of Legendre polynomials and eigenfunctions to estimateˆ 2 e t as that used to estimateĜ,P andM.
The R codes that were used to analyse these data are available on request to the corresponding author.

Application to other distributions
It should be noted that this model can be applied to other distributions. In the case of a normally distributed trait, we would assume: The prior distribution specified for σ 2 μ t corresponds to a flat prior that shrinks the residual variances towards zero [44]. Heritability in this case is estimated as:

Model selection
By fixing the value of M (number of eigenfunctions) and increasing the value of K (number of Legendre polynomials), we increased the number of polynomials used but fixed the number of components used to estimate each covariance matrix. Therefore, there was little change in the number of parameters estimated between models. However, since a large number of parameters was estimated compared to the number of observations, we compared these models using a bias corrected Akaike information criterion (AIC c ) [46]: where L(y) is the log likelihood of the data given the estimated parameters, p is the number of parameters estimated and n the total number of observations. We used this model selection criterion because a large number of parameters were estimated with respect to the sample size and because, in such a complex model, AIC c enforces a higher penalty than AIC, BIC (Bayesian information criterion) or DIC (deviance information criterion).

Results
Faecal egg counts (FEC) were quite dispersed (Fig. 1, black  bars). Exposing the lower 95th percentile of each distribution (shaded grey area) shows that as the grazing season progresses, the data are more variable and higher levels of dispersion are observed in May, August, September and in the post-mortem counts (Fig. 1, black bars). In particular, raw faecal egg counts are not normally distributed (Fig. 2a) and performing a log transformation i.e. log(FEC + 1) produced a mass around zero in our data (Fig. 2b). Figure 2a shows that the raw faecal egg counts are consistent with a negative binomial distribution (χ 2 goodness of fit test p-value = 0.7276). The log transformed FEC data (Fig. 2b) are significantly different from a normal distribution (Shapiro-Wilk test p-value < 0.001). In addition, log transformed FEC were not normally distributed for c = 25 or 0.2 (Shapiro-Wilk test p-value < 0.001, in both cases).
Treating each month as a separate trait, we estimated the additive genetic ( G G G corr ), maternal ( M M M corr ) and phenotypic ( Total Total Total corr ) correlation matrices as: An asterisk highlights correlations that differ significantly from zero. We also computed a series of bivariate animal models to estimate pairwise phenotypic, maternal and genetic correlations and found no significant differences between correlations computed in the bivariate models and in the full multitrait model. In addition to running multiple chains, this provided confidence that the chain had converged. Maternal correlations were not significant. Genetic correlations between July, August, September and October were strong and positive and between June and the remaining four months were strong and negative. Genetic correlations between June and October and July and October were not significant. Phenotypic correlations increased as the season progressed with a strong phenotypic correlation between September and October.

Fixed effects
In all analyses, sex and year were significant fixed effects. Considering each month separately, female lambs consistently had lower faecal egg counts than male lambs with the exception of the post-mortem counts. These differences were significant for the last three months (Fig. 3a).
Differences in faecal egg counts across the five years during June to October and in the post-mortem counts were significant (Fig. 3b). Generally, lower faecal egg counts were recorded for years 1992 and 1994 and higher egg counts were recorded for years 1995 and 1996.

Univariate analysis
Heritabilities estimated on the latent scale depended on the additional constant used in the log transformation (Fig. 4a). Therefore, heritability on this transformed scale depended on the particular transformation used and thus, did not provide robust heritability estimates. Heritability on the link scale, using a log link function in the negative binomial model of the untransformed raw count data, showed that the heritability of raw faecal egg counts increased univariately from 0 to 0.4 with a dip in August (Fig. 4b).

Random regression using WOMBAT
We compared the full rank analyses with K = M = 2, K = M = 3 and K = M = 4, and also compared each reduced rank analysis (with M = 1 and K = 2 or M = 2 and K = 3 and 4) with the corresponding full rank analysis. The reduced rank analysis model using four Legendre polynomials (K = 4) and two eigenfunctions (M = 2) minimised the AICc marginally (Table 2).
In general, including more Legendre polynomials gave less smooth curves but smaller values of AIC c ( Table 2). We found little difference in variance components between setting M = 2 and M = 1 with K = 2 (Fig. 5a) and small discrepancies between setting M = 2 and M = 3 with K = 3 (Fig. 5b). However, differences were greater between setting M = 2 and M = 4 with K = 4 (Fig. 5c). Therefore, a reduced rank analysis with higher degree polynomials had a greater effect on variance component estimates.

Negative binomial random regression model
The full rank scenarios were computationally infeasible due to the length of time required for convergence. Therefore, we considered reduced rank scenarios that were similar to those presented in the previous section. We used two eigenfunctions (M = 2) with two, three and four Legendre polynomials (K = 2, 3 and 4). The model using K = 3 minimised AIC c and thus, gave a better fit compared to K = 2 or K = 4 (Table 3).
We expected our variance curves to behave similarly to the WOMBAT curves with possible inflated variance estimates at the boundaries of the dataset (Fig. 5). The phenotypic variance increased between July and October (Fig. 6a, black lines). Likewise, the additive genetic a b  (Fig. 6a, brown and pink lines, respectively). The permanent environmental variance also increased slightly between July and October (Fig. 6, green lines). A heritability of 0.25 was estimated for June (Fig. 6b), but it may be influenced by inflated boundary estimates. It rose from 0.25 between June and July to 0.4 in October with a dip in August. This dip coincided with a similar pattern found in the univariate analysis (Fig. 4b).

Residual diagnostics
By plotting residuals against the fitted values (Fig. 7), we visually assessed the validity of two key assumptions of the normal log(FEC+1) and the negative binomial linear regression models: linearity and homoscedasticity (constant residual variance). In the negative binomial model, the assumption of homoscedasticity applies only to Pearson residuals, which adjust for the expected meanvariance relationship.
The residuals from the best fitting model of log(FEC+1) (setting K = 4 and M = 2 in WOMBAT) showed a strong nonlinear relationship with the fitted values (Fig. 7a, red line) whereas only a very weak trend was observed in the Pearson residuals from the best fitting model of raw faecal egg counts (setting K = 3 and M = 2 in the negative binomial random regression model; Fig. 7b, red line). Homoscedasticity was difficult to assess in the normal model due to the distinct stripe of points on the lower left corner of the plot. Pearson residuals from the negative binomial model appeared to vary more for low fitted values, apparently violating the assumption of homoscedasticity. However, interpretation was also difficult due to the much higher density of points at low fitted values. To informally assess homoscedasticity, we separated the ordered fitted values into ten equally sized groups and calculated the standard deviation of the corresponding residuals within each group. Variation in the standard deviations along the fitted values was substantial for the normal model and was twice as large for mid-level fitted values compared to low and high fitted values (Fig. 7c). By contrast, Pearson residuals from the negative binomial model showed little variation in standard deviation along the fitted values (Fig. 7d).
These observations indicate that the negative binomial random regression model provides a considerably better fit to raw faecal egg count data than the normal random regression of log transformed data.

Discussion
The overall goal of this work was to present some of the challenges that are inherent in the estimation of heritabilities from faecal egg count data and to demonstrate the use of random regression models to estimate heritability over time on a link scale using untransformed data. Faecal egg count data are typically overdispersed. Consequently, these data have previously been modelled both univariately and multivariately by transforming faecal egg count data and assuming that they are normally distributed. Our univariate analysis showed that heritability estimates depend on the additional constant used in the log transformation, as was demonstrated in previous studies [47]. The effectiveness of log and Box-Cox transformations decreases as the number of animals with zero egg counts increases [47]. We found that by increasing the additional constant in the log transformation, the heritabilities on the latent scale approached heritabilities on the link scale (Fig. 4). In a separate analysis, we estimated significant phenotypic and genetic correlations between faecal egg counts taken at different times during the season. This suggested that a multivariate analysis such as random regression was more appropriate for these data.
Our paper presents a method to fit a random regression that is based on MCMC. This may be viewed as computationally inefficient compared to REML estimates (for example, as implemented in WOMBAT [14]); however, Bayesian MCMC methods are appealing since they allow for more flexible fitting to complex models which is advantageous in the genetic analysis of non-normal data [18]. Here, we modelled untransformed raw faecal egg counts assuming that they were negative binomially distributed although this method can be easily adapted for fitting longitudinal data following other distributions. Various studies have compared heritability estimates from sire or animal models based on REML and MCMC methods using real and simulated data [44,[48][49][50]. Studies on real data have reported that REML estimates are more conservative than the corresponding MCMC estimates [48,49]. Previous studies showed the importance of prior specification of variance components and adequate mixing in Bayesian analysis [44]. Incorporating more informative priors, for example shrinkage priors, can produce more precise and accurate variance component estimates. However, imposing strict penalties in shrinkage priors inevitably affects parameter estimates. Generally, using more informative priors is beneficial in this type of analysis given the large number of parameters to be estimated. We chose non-informative priors for fixed effect regression coefficients and shrinkage priors for variance components. It is possible to impose stricter penalties [43] and a range of other prior distributions, however, this type of sensitivity analysis was outside the scope of this paper. Random regressions require fitting polynomial functions, but larger amounts of well-distributed data are required to fit high order polynomials [51]. Consequently, these methods are typically used in systems in which it is straightforward to generate much larger datasets than the one used here [52][53][54]. To address this issue, and reduce the number of parameters to be estimated, we conducted reduced rank analyses that estimated variance functions by eigen-decomposition of the corresponding covariance matrices of the random regression coefficients.
The reduced rank approach fits a number of polynomial functions (K) over time and estimates covariance functions using the largest few eigenfunctions. Analyses of the log transformed faecal egg counts, using the software WOMBAT, showed that decreasing the number of eigenfunctions (M) produced more differences compared to the full range analysis (Fig. 5) [55]. For all WOMBAT models considered here, we found possible inflated estimates of variance components at the boundaries of our dataset and consequently, we believe that there are inaccuracies in our reduced rank curves at the boundaries of the data (Figs. 5 and 6). However, this could also be the result of insufficient data to adequately estimate each of the covariance matrices.
Our analyses confirm findings that were obtained with other systems i.e., that heritabilities are not necessarily constant [23] and our results showed that the heritability of faecal egg counts increased as lambs got older, which is expected since exposure of grazing lambs to nematodes helps build immunity to infection. Heritabilities estimated with the MCMC random regression developed here ranged from 0.2 to 0.4 (Fig. 6) and were consistent with those on the latent scale from previous studies [4,24].
Faecal egg count is a commonly used marker in selective breeding programs for resistance to gastrointestinal nematodes and the heritability of such resistance markers predicts the effectiveness of breeding programs. Our negative binomial model provides estimates of heritability of raw faecal egg counts on a different scale to that used in previous analyses of log transformed data, however, breeding values on different scales tend to be highly correlated [56] and therefore can be interpreted in a similar way to those obtained using the normal animal model.
There is no correct scale on which to measure heritability. Statistical convenience and correspondence with the infinitesimal model in quantitative genetics has often led to calculate heritabilities for faecal egg counts on the latent scale (i.e. following log transformation). However, overdispersed count data such as that obtained here may also arise from multiplicative processes, providing a biological justification for the log link function [18]. We found that the negative binomial random regression model to provided a good fit to raw faecal egg counts whereas residual assumptions of the normal random regression model using log transformed faecal egg counts were not satisfied, which provided evidence in favour of the negative binomial model. Overall, our results agree with previous analyses and demonstrate that raw faecal egg count is a highly heritable trait.
In summary, random regression is a useful tool to analyse variance components from multivariate traits or, in this instance, multiple records of a trait per animal spread over a trajectory. To our knowledge, this is the first time random regression analyses have been used to handle non-normal data on the link scale. We have demonstrated the use of a Bayesian MCMC approach to apply random regression models to negative binomially distributed data. However, the approach taken here can be easily adapted to model data that are consistent with other non-normal distributions.