HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Erratum An erratum has been published Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Cantet, R. J. C. Articles by Fioretti, C. Search for Related Content PubMed PubMed Citation Articles by Cantet, R. J. C. Articles by Fioretti, C.

Semiparametric animal models via penalized splines as alternatives to models with contemporary groups¹

R. J. C. Cantet^*,,2, A. N. Birchmeier^*, A. W. Canaza Cayo and C. Fioretti

^* Departamento de Producción Animal, Universidad de Buenos Aires, C1417DSE Buenos Aires, Argentina; and Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Argentina; and Universidad del Altiplano, Puno, Perú; and and Estancias y cabaña, Las Lilas, S.A., Argentina

	Abstract

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

 
Contemporary groups (CG) are used in genetic evaluation to account for systematic environmental effects of management, nutritional level, or any other differentially expressed group effect; however, because the functional form of the distribution of those effects is unknown, CG serve as an approximation to a time-varying mean. Conversely, in semiparametric models, there is no need to assume any functional form for the time-varying effects. In this research, we present a semiparametric animal model (AMS) using the covariate day of birth (DOB) by means of penalized splines (P-splines), as an alternative to fitting CG. In the AMS, the functionality of the data on DOB is expressed by means of a Basic segmented polynomial line (B-spline) basis, and proper covariance matrices are used to reflect the serial correlation among the points of support (or knots) at different times. Three different covariance matrices that reflect either short- or long-range dependences among knots are discussed. Different models were fitted to birth weight data from Polled Hereford calves. Models compared were an animal model with CG, an animal model with CG and the covariate DOB nested within CG (CG + DOB), and P-splines with the first difference penalty matrix and three different AMS with 20, 40, 60, 80, or 120 knots. Models were compared using a modified Akaike information criterion (AIC_C), which was calculated as a byproduct of the estimation of variance components by REML using the expectation maximization algorithm. All three AMS had smaller (better) values of AIC_C than the regular model with CG, while producing almost the same ranking of predicted breeding values and similar average predicted error variance. In all AMS, the inference and all measures of comparison were similar when the number of knots was equal 40. The model CG + DOB had analogous performance to the AMS, but at the expense of using more parameters. It is concluded that the use of penalized regression splines using a B-spline basis with proper covariance matrices is a competitive method to the fitting of CG into animal models for genetic evaluation, without having to assume any functional form for the covariate DOB.

Key Words: Contemporary Groups • Penalized Splines • Semiparametric Models

	Introduction

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

 
In genetic evaluation, contemporary group (CG) is a discrete explanatory variable that accounts for a high percentage of the total phenotypic variation in most traits. An important issue is the definition of CG by setting cut-off dates to form new CG (Crump et al., 1997; Carabaño et al., 2004); bias may be introduced for records close to a cut-off date between two adjacent CG (Sivajarasingam, 1993). The reason is that the functional form of the relationship between the evaluated trait and the effects of herd environment in time is unknown. Thus, CG is an approximation to a time-varying mean in an attempt to describe the unknown distribution of those effects.

Semiparametric models are employed when the functional form of a covariate is unknown. In this situation, the underlying smooth function is usually a nuisance parameter, and the interest lies in accounting for the effects of the regressor variable (Altman, 2000). Thus, an alternative to CG is to fit a time covariate (e.g., day of birth [DOB]) without imposing any particular functional form on its effect on the trait, as suggested by Cantet (2002). Under this approach, one or more functions are fitted to the data that account for the irregular behavior of DOB. To fit such a function, Eilers and Marx (1996) proposed penalized splines (P-splines), a methodology that is closely connected to mixed models (Ruppert et al., 2003; Wand, 2003). Cantet (2002) suggested that regression splines may be used to substitute CG in genetic evaluation models. His presentation was expository and attempted to use truncated rather than Basic segmented polynomial line (B-spline) basis functions without details on how to perform the numerical computations. Thus, the goals of the present research were 1) to describe an animal model with DOB using P-splines and proper covariance matrices (AMS), and 2) to fit the AMS to birth weight records of beef cattle and to compare them with the fit obtained by the regular model with CG for genetic evaluation.

	Materials and Methods

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

 
P-Splines, B-Splines, and Semiparametric Animal Models
The process we attempt to model is a particular case of functional data in which the time (t) varying covariate DOB follows an underlying continuous time process f(t). This process is described by some random trajectory or trend within a herd, which results from environmental effects: management, grass availability, or weather characteristics such as temperature or rainfall, or a combination of all effects. The idea is to model this trajectory with a smooth function that is fitted using P-splines (Eilers and Marx, 1996). In contrast to the prediction of random regression models for functional breeding values (White et al., 1999; Huisman et al., 2002; Druet et al., 2003; Iwaisaki et al., 2005), in which the unit with repeated observations in time is the individual animal, in the approach used here, the herd is the unit with repeated observations over time. Other time-varying covariates such as age of dam may be modeled using the approach presented here, although a low-order polynomial can be used to describe such data in a more parsimonious fashion. With suitable modifications, the fit of trajectories in two or more dimensions (for example, in space), can be dealt with by using a two-dimensional B-splines basis, but we do not pursue the topic any further.

The P-splines are a combination of regression on a spline basis with a discrete roughness penalty to smooth data series or scatterplots (Eilers and Marx, 1996). The penalty controls the degree of smoothness while fitting the function. For P-splines, Eilers and Marx (1996) used equally spaced B-splines basis functions (De Boor, 1993). When fitted to time-dependent data, a covariate expressed on a B-spline basis results in the union of k-degree polynomial segments that have k – 1 continuous derivatives at the joining points, or knots. The resulting fit is smooth, with better numerical properties than a polynomial fit of high degree (e.g., k > 5; Green and Silverman, 1994). The B-splines of degree k have local support, which means that they are defined only in a small part of the real line between k + 1 knots, being equal to 0 outside this small segment of the real line. Compared with least squares, this prevents any observation that may affect the entire shape of the function. Eilers and Marx (1996) observed that, in P-splines, the number of knots is not a critical parameter as long as "there are enough of them." This lack of sensitivity to the number of knots is due to the control exercised by the penalization. Thus, problems of multi-collinearity can be fixed by changing the penalty parameter (P. Eilers, Leiden Univ., The Netherlands and B. D. Marx, Louisiana State Univ.; personal communication). Moreover, Ruppert and Carroll (2000) found that increasing the number of knots to >40 resulted in slight differences in fit after extensive simulation studies. Using a number of knots in the order of 40 has the effect of dramatically decreasing the dimensionality of the parameter space compared with other smoothing spline methods. Notwithstanding this, increases in the number of records may require enlarging the parameter space. We dealt with this issue by comparing models with different numbers of knots.

Covariate DOB Expressed Using a B-Spline Basis
Consider a random vector y (n x 1) of records, which is functionally dependent on time (t) through a general vector valued function f(t). We will consider t to be DOB in Julian days. It is desired that f(t) be a smooth function of t that can be fitted using mixed model theory. To do so, let

[1]

In [1], f(t) is a linear combination of parameters (b_i) that are to be estimated, and elements of B-spline basis functions B^(k)_i i = 1, 2, ..., nx (De Boor, 1993). The number of basis functions B^(k)_i needed to express DOB is nx + 6 (Eilers and Marx, 1996); three additional knots are added to each extreme. Usually the degree of a B-spline is 3 (k 3). In the present research, cubic splines (k = 3) were used because they are sufficiently flexible to retain most features of the data.

Calculations of the B⁽³⁾_i coefficients were performed with the recursive algorithm of De Boor (1993). A brief explanation follows. Let t₁, t₂, ..., t_nx+6 be the set of knots that expand the range of DOB. Then, B⁽⁰⁾_i = 1 if DOB is in the interval between t_i and t_i+1; otherwise, B⁽⁰⁾_i = 0. For example, suppose t₁₀ = 3303.5 and t₁₁= 3390.8, then for DOB = 3353, B⁽⁰⁾₁₀ = 1, whereas B⁽⁰⁾₁ , ..., B⁽⁰⁾₉ , B⁽⁰⁾₁₁, ..., and B⁽⁰⁾_nx+6 are all 0. The next step in the algorithm is to calculate the B⁽¹⁾_is, then the B⁽²⁾_is and finally the B⁽³⁾_is, using the scheme displayed in Figure 1, which is called the blossom. Four B⁽³⁾_i elements (B⁽³⁾_i–3, B⁽³⁾_i–2, B⁽³⁾_i–1, and B⁽³⁾_i) are needed to express one value of DOB in terms of a cubic B-splines basis. These four values add up to 1, and the recursion formula (De Boor, 1993) to calculate them is equal to

View larger version (12K): [in this window] [in a new window]   Figure 1. Order of calculation of the Basic segmented polynomial line coefficients.

[2]

For the sake of completeness, truncated bases are an alternative to B-splines bases. Using truncated bases, the function in [1] is expressed as

where t_ijis the DOB measure for herd i at time j, m is the order of fit (usually m = 1, 2, or 3), and b₀, b_j, and the u_kl are the parameters to estimate. The values of the knots are_k (k = 1, 2 ..., nx), and the quantity (t_ij– _k)₊ is taken to be equal to t_ij – _k whenever t_ij > _k, or 0 otherwise. If l = 1, the fit is linear; it is quadratic for l = 2 and cubic for l = 3.

Expression [1] can be written in matrix form as Bb, where B is the n x nx matrix that contains the B⁽³⁾_i, and b is the parametric vector (nx x 1) containing the b_i for f(t). Each row of B has all elements equal to zero except for basis coefficients B⁽³⁾_i–3, B⁽³⁾_i–2, B⁽³⁾_i–1, and B⁽³⁾_i in columns i – 3, i – 2, i – 1, and i, respectively. Thus, each value of the covariable DOB is transformed into four B-spline coefficients in the interval (0, 1) for each animal with a record in y.

P-Splines
To obtain the estimators of b in the model

[3]

Eilers and Marx (1996) proposed maximizing the likelihood while penalizing the sum of squares of the differences between adjacent b_i, either for first or for second squared differences . The resulting function is then proportional to

[4]

The scalar controls the amount of smoothing. For first differences, matrices D (nx – 1 x nx) and D'D are respectively equal to

[5]

Notice that D' D (nx x nx) is singular. If second differences are considered, then matrices D and D'D are

[6]

In [6], D is of order (nx – 2 x nx), and D'D is again singular. Eilers and Marx (1996) obtained the penalized estimator of b as the solution of the following system of equations:

[7]

The effect of the penalty is to shrink b in an amount proportional to . The connection between P-splines and mixed models (e.g., Ruppert et al., 2003; Wand, 2003) is now apparent. In a mixed-model setting, b can be viewed as random effects, to the ratio of the error variance to the variance of the bs (Ruppert et al., 2003), whereas D'D may be interpreted as some g-inverse of the variance-covariance matrix of the linear spline parameters. From a Bayesian viewpoint, a singular D'D induces the prior distribution of the linear spline parameters to be improper, which, in turn, causes the posterior distribution to be improper (Hobert and Casella, 1996; Lang and Brezger, 2004). A way around this problem is to formulate an equivalent mixed model (Henderson, 1984) such as the one discussed by Currie and Durban (2002); however, the fitting of such a model is computationally involved, as it makes the mixed model equations extremely dense and does not behave in a numerically stable manner. Alternatively, proper covariance matrices of the bs at the knots provide a better fit to the model than singular penalty matrices as shown subsequently.

Animal Models with a Functional Covariate Using a B-Spline Fit
To specify an animal model suitable for genetic evaluation, we now incorporate fixed effects (e.g., age of dam and sex) and breeding values to the specification of DOB in [1] and [3]. The resulting animal model is equal to

[8]

The incidence matrix X (n x p) associates data to the p x 1 parametric vector ß of fixed effects; Z (n x q) relates elements of y to the random vector a (q x 1) of breeding values, whereas e is the n x 1 vector of random errors. We assume a full-rank parametrization in ß so that rank [X] = p. The Gaussian vectors a and e are stochastically independent; both have zero expectations and covariance matrices equal to A²_A and I²_e, respectively. Matrix A contains the additive relationships. The random variables that relate to the knots are assumed to be distributed as b ~ N(0, G_S²_b), where ²_b is the variance component between knots and G_S is a positive definite matrix that portrays the covariance structure in time for the knots. Three alternative specifications for G_S are described in the next section. The expectations and the variances and covariances of all random vectors in [8] are equal to

[9]

where V = ZAZ' ²_A + BG_SB' ²_b + I²_e. The variance components are the additive genetic variance where ²_A, the variance between knots where ²_b, and the error variance ²_e. Mixed-model equations for [8] are:

[10]

where = ²_e/²_b and = ²_e/²_A

Covariance Matrices of the Elements in b
There are several considerations in the search for a suitable covariance matrix G_S. First, it should be able to account for "serial correlation" (Diggle et al., 1994) in such a way that pairs of knots that are closer in time are likely to be more strongly correlated than pairs farther distant in time. A second consideration is that correlation structures that, after inversion, are similar to the banded matrices of differences [5] and [6] would behave similarly to the original formulation of P-splines. It also was desirable to choose a linear covariance structure (i.e., G_S = P²_b) from a regular multivariate normal density rather than from a "curved" normal density (Lehmann, 1983). In this latter case, the covariance structure depends on dispersion parameters in a nonlinear fashion, such as in an autoregressive process. Linear dispersion structures allow the use of REML-expectation maximization or Gibbs sampling, avoiding the need for more involved algorithms such as Metropolis-Hastings. Finally, it was considered to be important that inversion of G_S not be computationally expensive. As a result, we used the covariance structure (P) that accounts for a linear decay of the correlation with time originally proposed by Cantet (2002). Off-diagonal elements of P (P_ij) are functions of time lag t_j – t_i, between the knots i and j. Then, for any pair of knots at times t and t + w, . This covariance structure is stationary, as it depends on time lag w but not on the time moments t or t + w. If the time measures (expressed in days, weeks, or months) are t₁ < t₂ < ... < t_nx–1 < t_nx, diagonal elements of P are equal to 1 and off-diagonals are for t_j > t_i. This has the effect of mapping the difference t_j – t_i into the interval [0,1) such that, as t_j – t_i decreases (i.e., knots are positioned closer), the correlation between the spline effects increases linearly. This formulation of P also allows dealing with irregular timings, whereas the inverse of P can be computed proportional to O(nx) calculations such as for A^–1, a development shown in Appendix A. To model P with equally spaced knots, let j be a nx x 1 vector with all elements equal to 1; then t_nx = (t_j – t_i) nx, and P is equal to

[11]

To exemplify, let nx = 5, then P is

which is a matrix with Toeplitz structure (Marin and Dhorne, 2002).

The second formulation considered is related to the stochastic interpretation of f(t) given by Wahba (1990), who showed that a polynomial smoothing spline is the solution to the stochastic differential equation of a Wiener process. Wecker and Ansley (1983) and De Jong and Mazzi (2001) pursued this idea and obtained an expression for the variance-covariance matrix of the spline function and its derivative. Guo (2002) used this formulation to write down a mixed model for functional data. Finally, Hyndman et al. (2005) used the expression for the variance of the stochastic process involving splines that was obtained by De Jong and Mazzi (2001) to derive the following expression for the covariance between knots i and j: . In terms of a covariance matrix for the cubic spline functions, we have

[12]

Matrix also shows a decay of correlation in time, but this is not linear as in P. We were not able to find an algorithm that allows inverting as simply as P.

In the correlation structures induced by matrices P and , all off-diagonal elements are non-zero, indicating dependency among all random variables in b, although long-range covariances may be quite small. Alternatively, one may want to model a correlation structure in which there is covariance only with the next neighbor knot. For this formulation, Durban et al. (2001) presented a covariance matrix for the spline function using a decomposition of the penalty matrix discussed by Green and Silverman (1994). The resulting matrix is tridiagonal, with main diagonal elements equal to U_i,i = 2/3 i = 1, 2, . . ., nx, and off-diagonals U_i+1,i = U_i,i+1 = 1/6, for i = 1, 2, . . ., nx – 1, being 0; otherwise,

[13]

The matrix U can be viewed as the covariance structure of a first-order moving average process with correlation between adjacent knots equal to ¹/3 and 0 thereafter. All AMS discussed in this section, plus the regular animal model with CG and the original P-spline formulation with matrix D'D in [5], were fit to a data set containing birth weights of beef cattle.

Number of Knots
Although knot selection methods exist, P-splines rely on using a large number of knots (Eilers and Marx, 1996), and on limiting the effect of the knots, while constraining the size of the spline coefficients (the b_i). However, for genetic evaluation purposes, the variable DOB expands several years and calving seasons, and its range increases as new data arrive. Therefore, the dimension of the parameter space (the number of knots) may have to increase when the number of records increases. Brumback and Rice (1998) observed that the curve is a finite nx-dimensional approximation of the function f(t), which lies in an infinite dimensional parameter space. An appropriate reduction of the parameter space can be produced using the method of sieves (Chen, 2005). By sieves, it is meant a sequence of approximating, and significantly less complex parameter spaces, which optimize some criterion function. The consistency of the method is ensured by increasing the dimension of the parameter space with the increase in sample size. Whereas the asymptotics of the method of sieves are complicated because of the dual approach to infinity of the dimension of the parameter space and the number of data (nx and n ), the implementation is easily achieved by fitting models with increased numbers of knots and by comparing them using, for example, the Akaike information criterion (AIC_C) as a criterion function. In general, if too many knots are used, the fit tend to be unstable and highly variable (the curve is too "wiggly"). Conversely, fitting too few knots leads to bias. To quantify the effect of the increase in the number of knots, models with 20, 40, 80, or 120 equally spaced knots may be fitted and compared according to convergence, AIC_C, and the way they handle periods in which DOB are not observed.

Differential Management Effects
For multiple herd evaluations or to model interruptions of measures, different herds, or types of management, G_S can be taken to be block-diagonal. Whenever some records in a herd cannot be considered as part of the stochastic process related to herd-time effects, the fitting of DOB requires a modification of model [8]–[9]. For example, suppose the trait of interest is weaning weight and, throughout the years, a few animals received a different nutritional management than the majority of them. One possible solution is to modify [8] to incorporate unrelated time-varying parameters to the function in [1]. Let the vector of these CG effects be b_G, which is related to the records by the incidence matrix B_G. Rows in B_G will be equal to zero, except for those belonging to animals from the differential management, which will have a 1 in the column related to the CG in b_G. This vector may be treated as either fixed or random. In the latter case, b_G ~ N(0, I²_c), with ²_c either being equal to or different from ²_b, in which case this variance has to be estimated too.

Data
Records were 5,175 birth weights from a purebred Polled Hereford herd, belonging to Estancias y Cabaña Las Lilas S.A. and located in Pasteur, western Buenos Aires province, Argentina. Data were collected from 1972 to 2000. A total of 9,742 animals were included in the pedigree file. The records were from 2,739 males and 2,436 females calves, respectively, sired by 177 bulls and from 1,825 dams. The cow herd was kept in cultivated pastures without supplemental feeding throughout the entire year. During most years, there were calving seasons in spring and fall. Therefore, a CG was defined for each calving season within year, creating 53 CG. The average CG day span was 54 d (range = 12 to 107 d), and the average number of calves per CG was 98 (range = 9 to 204).

Models of Analysis.
Six animal models were fitted to the records. All models included fixed effects of sex of calf and age of dam and random breeding values and error terms. Four models had DOB as a covariate expressed with a cubic B-spline basis, such as in [1] and [8]. The fifth model included CG as a fixed classification variable, and the sixth also included the covariate DOB nested within CG, as suggested by a reviewer. For X to have full column rank, the sex effect of females was set to 0. In the two models with CG, the last age of dam also was set to 0.

B-Spline Fit of DOB.
Covariate DOB was calculated using Julian days. Day 1 was the day the first registered calf in the herd was born (July 2, 1972). The B-spline coefficients for DOB were calculated using the FORTRAN subroutine BSLPVN (De Boor, 1977), with 46 (= nx + 6, or nx = 40) equally spaced knots; however, the order of the vector b, as well as of the penalty matrix D'D or the covariance matrices P, , and U, was equal to nx = 40.

Estimation of Variance Components.
The procedure used to estimate the dispersion parameters was REML (Patterson and Thompson, 1971) by means of the expectation maximization algorithm (Dempster et al., 1977). A program was written in FORTRAN to perform all calculations. The estimator of the additive variance ²_A at the r-iteration was equal to

where â is BLUP(a), and C^aa is the portion of the inverse of the mixed model equations associated with a in the inverse matrix:

[14]

The REML-EM estimator of ²_b is

where = BLUP(b), and C^bb is as in [14]. Finally, the error variances in all AMS were estimated using the formula

whereas in the two models with GC, the estimated error variance was equal to

In all cases n = 5,175, and ê = BLUP(e). The algorithm was performed until the difference in the estimates from two successive iterates for any variance component was <10^–2 kg². To make the estimates from all models comparable, regardless of the definition of either fixed CG or random DOB effects, heritability (h²) was estimated as ² = ²_A/(²_A + ²_e).

To compare the fit obtained with the different models, AIC_C, as adapted by Hurvich et al. (1998) to local smoothing spline estimators, was calculated as

[15]

the estimated error variance and H is a symmetric matrix such that = Hy. The statistic AIC_C was calculated as a byproduct of the EM algorithm when estimating the variance components using the expression (see Appendix B for its derivation):

The model with the smallest value of AIC_C is to be selected (Hurvich et al., 1998).

	Results

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

 
To improve convergence for the model with G_S = P, P^–1 was divided by its (1, 1) or (nx, nx) element, which is equal to 2/nx (see A1 in Appendix A). This is equivalent to Var(b) = P²_b', where ²_b' = (nx/2)²_b. Notice that the resulting matrix is similar to the first difference penalty matrix D'D in [5], except for the non-zero elements in (1, nx) or (nx, 1), which makes P be nonsingular. In addition, the matrix was multiplied by the cube of the constant interval between equally spaced knots (t_i+1 – t_i)³ to improve convergence. Having done that, the four models with DOB and 20, 40, and 60 knots and the two models with CG met the criterion of convergence in <200 rounds of iteration. The convergence criterion was not met after 1,000 rounds of iteration for the models that included matrix with 80 or 120 knots. The REML-EM estimates of ²_A, ²_b, ²_e, h and the value of AIC_C from all models are displayed in Table 1. Within any AMS, estimates of ²_A and ²_e were similar across models with different numbers of knots, and the model with matrix D'D had the largest range of estimates. As a result, ² values in all AMS were similar to the values observed for both models with CG, whereas the model with D'D showed the largest value of ². The estimates of ²_b were different in the four AMS and tended to increase as the number of knots increased. This was not the case for the model with P, where ²_b showed a somewhat erratic downward trend.

Figure 2 displays the solutions of CG effects (the thickest line with irregularities), and the BLUP(b) of all AMS from the models with 40 knots. The shapes of the fitted curves were affected by the model. The curve from the AMS with matrix U was more wiggly, whereas the one from the model with penalty matrix D'D was smoother than the curves from the other models. The diagonal lines in the curve for the model with CG are caused by a lack of observations at that time. Of particular interest is the period expanding from 5,012 to 5,219 d of DOB, in which no birth weights were recorded. Notice that the solution of the CG on the left of 5,012 d is lower than the one on the right of 5,219 d. The curve from the model with penalty matrix D'D passes above the solution of the CG on the left and below the CG on the right. The path from the AMS with P went up in a more or less straight fashion; however, the curves from the models with and U went down after 5,012 d to a local minimum at 5,127 d and then increased to 5,219 d. Thus, this effect seems to be an artifact of both models, which is not supported by the data. During most time intervals with observed birth weights, the curves from the AMS with matrices P and tended to have a similar shape.

View larger version (17K): [in this window] [in a new window]   Figure 2. Differences in birth weight plotted against Julian day of birth (DOB) for different models that differ in the covariance matrices for B-spline coefficients. Covariance matrices are D'D, first-order penalty matrix; P, as proposed by Cantet (2002); , as proposed by Hyndman et al. (2005); U, as proposed by Durban et al. (2001), and the usual model with fixed contemporary group (CG).

View larger version (18K): [in this window] [in a new window]   Figure 3. Differences in birth weight plotted against Julian day of birth (DOB) for the models with covariance matrix for B-spline coefficients P 2b, as proposed by Cantet (2002), for different numbers of knots.

	Discussion

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

Most other applications of splines in animal breeding are based on modeling functional breeding values (White et al., 1999; Huisman et al., 2002; Druet et al., 2003; Iwaisaki et al., 2005). The literature on smoothing methods with splines has become abundant, and the fit of semiparametric methods using mixed linear model software has become popular (Ruppert et al., 2003; Wand, 2003). This is especially true for P-splines (Eilers and Marx, 1996), either viewed from the frequentist (Wand, 2003) or the Bayesian (Lang and Brezger, 2004) camps. There are many reasons for the appeal of this methodology that generalizes ordinary smoothing splines to knot sequences, which are usually much smaller than the response variable. The P-splines have the stable numerical properties of B-splines compared with other basis function approaches, such as truncated basis (Eilers and Marx, 2005). We have shown the flexibility of P-splines to accommodate different specifications of covariance matrices as penalties. Moreover, the fit of the function does not require any assumption on the parametrization of the curve.

In this research, we compared curves fitted with P-splines in models with different numbers of knots but at equal intervals or spacings. Within any type of penalty matrix or covariance structure, the values of all statistics used for comparison (²_A, ²_e, ², AIC_C, average PEV(a), average accuracy, and ranking of predicted breeding values) were similar for models with different numbers of knots. The only parameter that showed sizeable changes with the increase in the number of knots was ²_b, which is the denominator of the smoothing parameter . These results are consistent with the observation of Eilers and Marx (1996) that the number of knots is not a critical parameter in P-splines but the smoothing parameter is. After extensive simulation, Ruppert and Carroll (2000) found that increasing the number of knots to >40 resulted in marginal decreases in mean square error for all models and examples they worked out. Furthermore, setting the knots at equal spacings performed slightly better than knot sequential procedures. In the words of Altman (2000), "penalized regression splines appear to be robust to the choice of knots (as long as there are sufficiently many)." The curves fitted with matrix P (Figure 3) suggest that 20 knots are not adequate to take care of all features of the trend, whereas 80 or 120 are less parsimonious than 40 knots, but give some inconsistencies in intervals of DOB with no birth weight recorded. In cases where more than one curve has to be fit (for example, a curve per herd or per system of management within a herd), the B-spline coefficients and the knots have to be calculated on a herd-by-herd basis, and Var(b) will be a block diagonal matrix.

The choice of covariance matrix in the AMS attempted to take into account the situations of either complete dependence among the random variables related to the knots (covariance structures from matrices P and ) or to reflect covariance only with the next neighbor knot (matrix U) in the other extreme. The performance of all three AMS, in terms of the estimates of the variance components and h², PEV(a), and AIC_C, were similar. However, there were differences in the shape of the fitted curves (model with matrix U produced more wiggly curves) and in the approach to convergence when the number of knots increased (model with matrix not converging when 80 or 120 knots were fitted). All things considered, the AMS with matrix P and 40 knots seems to be the model of choice for this data.

In this research, the estimation of dispersion parameters and the model comparison procedure (AIC_C) followed a likelihood approach; however, to estimate the variance components, a Bayesian method and Gibbs sampling (Sorensen and Gianola, 2002) also are employed in any of the models considered. In fact, models [8]–[9] with covariance matrices [11], [12], or [13] are Bayesian P-splines models (Lang and Brezger, 2004), which results from an nx-variate normal prior for b such that a priori b ~ N(0, G_S²_b) with G_S = P, or U.

	Implications

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

 
Semiparametric animal models with a penalized-spline fit of the covariate day of birth of the animal and a proper covariance matrix of the spline coefficients had a better fit than the regular model with contemporary groups. Moreover, the penalized-splines models produced similar ranking of predicted breeding values and average accuracy for birth weights of beef cattle. Thus, the use of penalized regression splines using a B-spline basis with proper covariance matrices is a competitive method to avoid fitting contemporary groups (or herd-year-seasons) into animal models for genetic evaluation without having to assume any parametric form for the environmental effects of herd in time.

	Appendix A: Elements of P–1

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

 
Let the times of measure (expressed in days, weeks, or months) be t₁ < t₂ < ... < t_nx–1 < t_nx. All diagonal elements of P are 1, whereas off-diagonals are equal to 1 – [(t_j – t_i)/t_nx] whenever t_j > t_i. Let the matrix of knot differences be

So we can write P = jj' – (t_nx)^–1, where j is a nx x 1 vector with all elements equal to 1. To invert P, we use the formula for the inverse of a sum of matrices (Harville, 1997):

Rybicki and Press (1992) observed that ^–1 is equal to

Although P may have all elements different from 0, its inverse shares the structure of ^–11. Again, on using the inverse of a sum formula and after algebra, the elements of P^–1 are

With equally spaced knots, all differences t_j– t_i are equal, and P^–1 is equal to

[A1]

	Appendix B: Corrected Akaike Information Criterion (AICC)

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

 
Hurvich et al. (1998) wrote the AIC_C statistics as follows:

[B1]

where ²_e is the estimated error variance, n is the number of data in y, and H is a symmetric matrix such that the predicted data vector is = Hy. Under model [8], is equal to

[B2]

The left-hand side of [B2] can be written as

where and as in [14].

Then, take and rotate it so as to obtain

[B3]

Now, let S be equal to

and add and subtract S from [B3] in tr(H). We then have

where MS is equal to

Finally,

or

[B4]

Expression [B4] is the operational formula to calculate tr(H).

	Footnotes

 
¹ Funding for this research was provided by grants from Secretaría de Ciencia y Técnica, UBA (UBACyT G018—2004–2007) and from Agencia Nacional de Ciencia y Tecnología (FONCyT PICT 08-9502) of Argentina.

² Correspondence: Agronomía UBA, Av. San Martin 4453 (phone: 54-011-4514-8000, ext. 8784; fax: 54-011-4514-8000, ext. 8785; e-mail: rcantet{at}agro.uba.ar).

Received for publication April 29, 2005. Accepted for publication July 12, 2005.

	Literature Cited

Top Abstract Introduction Materials and Methods Results Discussion Implications Appendix A: Elements of... Appendix B: Corrected Akaike... Literature Cited

 


Altman, N. 2000. Krige, smooth, both or neither? Aust. N. Z. J. Stat. 42:441–461.

Brumback, B. A., and J. A. Rice. 1998. Smoothing spline models for the analysis of nested and crossed samples of curves. J. Am. Stat. Assoc. 93:961–976.

Cantet, R. J. C. 2002. B-spline fitting of time as an alternative to contemporary group effect. Communication No. 20-11 in Proc. 7th World Cong. Genet. Appl. Livest. Prod., Montpellier, France.

Carabaño, M. J., A. Moreno, P. López-Romero, and C. Díaz. 2004. Comparing alternative definitions of the contemporary group effect in Avileña Negra Ibérica beef cattle using classical and Bayesian criteria. J. Anim. Sci. 82:3447–3457.[Abstract/Free Full Text]

Chauhan, V. P. S., and R. Thompson. 1986. Dairy sire evaluation using a ‘rolling months’ model. J. Anim. Breed. Genet. 103:321–333.

Chen, X. 2005. Large sample sieve estimation of semi-nonparametric models. Handbook of Econometrics volume 6. Available: http://athens.src.uchicago.edu/jenni/handbookv6/Chen/. Accessed June 28, 2005.

Crump, R. E., N. R. Wray, R. Thompson, and G. Simm. 1997. Assigning pedigree beef performance records to contemporary groups taking account of within-herd calving patterns. Anim. Sci. 65:193–198.

Currie, I., and M. Durban. 2002. Flexible smoothing with P-splines: A unified approach. Stat. Model. 4:333–349.

De Boor, C. 1977. Package for calculating with B-splines. SIAM J. Numerical Anal. 14:441–472.

De Boor, C. 1993. B(asic)-spline basics. Pages 27–49 in Fundamental Developments of Computer-Aided Geometric Modeling. L. Piegl, ed. Academic Press, San Diego, CA.

De Jong, P., and S. Mazzi. 2001. Modeling and smoothing unequally spaced sequence data. Stat. Inference Stochastic Processes 4:53–71.

Dempster, A. P., N. M. Laird, and D. B. Rubin. 1977. Maximum likelihood from incomplete data via the EM algorithm. J. Royal Stat. Soc., Series B 39:1–38.

Diggle, P. J., K. Y. Liang, and S. L. Zeger. 1994. Analysis of Longitudinal Data. Clarendon Press, Oxford, UK.

Druet, T., F. Jaffrézic, D. Boichard, and V. Ducrocq. 2003. Modeling lactation curves and estimation of genetic parameters for first lactation test-day records of French Holstein cows. J. Dairy Sci. 86:2480–2490.[Abstract/Free Full Text]

Durban, M., I. Currie, and R. Kempton. 2001. Adjusting for fertility and competition in variety trials. J. Agric. Sci. (Camb.) 136:129–140.

Eilers, P. H. C., and B. D. Marx. 1996. Flexible smoothing with B-splines and penalties (with comments and rejoinder). Stat. Sci. 11:89–121.

Eilers, P. H. C., and B. D. Marx. 2005. Splines, knots and penalties. Technical report. Available: http://www.stat.lsu.edu/faculty/marx/. Accessed June 23, 2005.

Green, P. J., and B. W. Silverman. 1994. Nonparametric Regression and Generalized Linear Model. Chapman & Hall, London, UK.

Guo, W. 2002. Functional mixed models. Biometrics 58:121–128.[Medline]

Harville, D. A. 1997. Matrix Algebra from a Statistician’s Perspective. Springer-Verlag, New York, NY.

Henderson, C. R. 1984. Applications of Linear Models in Animal Breeding. Univ. of Guelph, Ontario, Canada.

Hobert, J. P., and G. Casella. 1996. The effect of improper priors on Gibbs sampling in hierarchical linear mixed models. J. Am. Stat. Assoc. 91:1461–1473.

Huisman, A. E., R. F. Veerkamp, and J. A. M. van Arendonk. 2002. Genetic parameters for various random regression models to describe the weight data of pigs. J. Anim. Sci. 80:575–582.[Abstract/Free Full Text]

Hurvich, C. M., J. S. Simonoff, and C. Tsai. 1998. Smoothing parameter selection in nonparametric regression using an improved Akaike Information Criterion. J. Royal Stat. Soc., Series B 60:271–293.

Hyndman, R. J., M. L. King, I. Pitrun, and B. Billah. 2005. Local lineal forecasts using cubic smoothing splines. Aust. N. Z. J. Stat. 47:87–99.

Iwaisaki, H., S. Tsuruta, I. Misztal, and J. K. Bertrand. 2005. Genetic parameters estimated with multitrait and linear spline-random regression models using Gelbvieh early growth data. J. Anim. Sci. 83:757–763.[Abstract/Free Full Text]

Lang, S., and A. Brezger. 2004. Bayesian P-splines. J. Comput. Graphic. Stat. 13:183–212.

Lehmann, E. L. 1983. Theory of Point Estimation. J. Wiley & Sons, New York, NY.

Marin, J. M., and T. Dhorne. 2002. Linear Toeplitz covariance structure models with optimal estimators of variance components. Linear Algebra Appl. 354:195–212.

Patterson, H. D., and R. T. Thompson. 1971. Recovery of inter-block information when block sizes are unequal. Biometrika 58:545–554.[Abstract/Free Full Text]

Ruppert, D., and R. J. Carroll. 2000. Spatially-adaptive penalties for spline fitting. Aust. N. Z. J. Stat. 42:205–253.

Ruppert, D., M. P. Wand, and R. J. Carroll. 2003. Semiparametric Regression. Cambridge Univ. Press, Cambridge, UK.

Rybicki, G. B., and W. H. Press. 1992. Interpolation, realization, and reconstruction of noisy, irregularly sampled data. Astrophys. J. 398:169–176.

Sivarajansingam, S. 1993. Comparison of alternative methods of handling contemporary group effects in animal model prediction. J. Anim. Breed. Genet. 110:401–411.

Sorensen, D., and D. Gianola. 2002. Likelihood, Bayesian and MCMC Methods in Genetics. Springer Verlag, Berlin, Germany.

Taylor, J. M. G., and M. Yu. 2002. Bias and efficiency loss due to categorizing an explanatory variable. J. Multivariate Anal. 83:248–263.

Wand, M. P. 2003. Smoothing and mixed models. Comput. Stat. 18:223–249.

Wahba, G. 1990. Spline Models for Observational Data. Soc. Ind. Appl. Math., Philadelphia, PA.

Wecker, W. E., and C. F. Ansley. 1983. The signal extraction approach to nonlinear regression and spline smoothing. J. Am. Stat. Assoc. 78:81–89.

White, I. M. S., R. Thompson, and S. Brotherstone. 1999. Genetic and environmental smoothing of lactation curves with cubic splines. J. Dairy Sci. 82:632–638.[Abstract]