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

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal 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 Bohmanova, J. Articles by Sullivan, P. G. Search for Related Content PubMed Articles by Bohmanova, J. Articles by Sullivan, P. G.

Comparison of Random Regression Models with Legendre Polynomials and Linear Splines for Production Traits and Somatic Cell Score of Canadian Holstein Cows

J. Bohmanova^*,1, F. Miglior^,, J. Jamrozik^*, I. Misztal and P. G. Sullivan

¹ *Centre for Genetic Improvement of Livestock, Department of Animal and Poultry Science, University of Guelph, Guelph, Ontario, N1G 2W1, Canada
Agriculture and Agri-Food Canada, Dairy and Swine Research and Development Centre, Sherbrooke, Quebec, J1M 1Z3, Canada
Canadian Dairy Network, Guelph, Ontario, N1K 1E5, Canada
Department of Animal and Dairy Science, University of Georgia, Athens 30602

¹ Corresponding author: jbohmano{at}uoguelph.ca

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
A random regression model with both random and fixed regressions fitted by Legendre polynomials of order 4 was compared with 3 alternative models fitting linear splines with 4, 5, or 6 knots. The effects common for all models were a herd-test-date effect, fixed regressions on days in milk (DIM) nested within region-age-season of calving class, and random regressions for additive genetic and permanent environmental effects. Data were test-day milk, fat and protein yields, and SCS recorded from 5 to 365 DIM during the first 3 lactations of Canadian Holstein cows. A random sample of 50 herds consisting of 96,756 test-day records was generated to estimate variance components within a Bayesian framework via Gibbs sampling. Two sets of genetic evaluations were subsequently carried out to investigate performance of the 4 models. Models were compared by graphical inspection of variance functions, goodness of fit, error of prediction of breeding values, and stability of estimated breeding values. Models with splines gave lower estimates of variances at extremes of lactations than the model with Legendre polynomials. Differences among models in goodness of fit measured by percentages of squared bias, correlations between predicted and observed records, and residual variances were small. The deviance information criterion favored the spline model with 6 knots. Smaller error of prediction and higher stability of estimated breeding values were achieved by using spline models with 5 and 6 knots compared with the model with Legendre polynomials. In general, the spline model with 6 knots had the best overall performance based upon the considered model comparison criteria.

Key Words: random regression model • Legendre polynomial • linear spline

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Random regression model (RRM) methodology is commonly used for routine genetic evaluation of production traits in dairy cattle (Interbull, 2007). The basic idea of RRM consists of fitting average lactation curves for a subpopulation and animal specific curves describing deviations from the average curves. Several functions can be used to fit fixed and random regressions. Early applications made use of parametric functions and lactational shape functions such as the Ali and Schaeffer (1987) and Wilmink (1987) functions. However, those functions were subsequently replaced by orthogonal Legendre polynomials (LP) because the former models were unable to properly model the peak of the lactation curve and tended to converge slowly because of high correlations among coefficients (Schaeffer, 2004). Models with LP as regressions are orthogonal and, therefore, have better convergence properties than models with parametric or lactational shape functions (Pool et al., 2000). However, to fit the shape of lactation appropriately, higher order polynomials are required. In this case, implausibly high additive genetic variances at the extremes of lactation and negative correlations between the most distant test-days can be obtained (Jamrozik et al., 2001; López-Romero et al., 2004). The overestimation of variances at the extremes of the lactation is often explained by Legendre polynomials’ lack of asymptotes; poor performance of fitting data at the extremes of lactations (López-Romero et al., 2004); or not accounting for between-herd variation in the shape of the lactation (Gengler and Wiggans, 2001; Jamrozik et al., 2001; de Roos et al., 2004). Recently, splines have been advocated as a good alternative to LP for analyzing test-day yields in RRM (White et al., 1999; Druet et al., 2003; Silvestre et al., 2005), mainly because of their limited sensitivity to the data (individual records influence only a specific part of the function) and higher flexibility of fitting lactation curves (Misztal, 2006).

Splines are piecewise functions consisting of independent segments that are connected in so-called knots (de Boor, 1978). The segments are described by lower order polynomials. The simplest case of a spline function is a linear spline, where the segments are fitted by linear polynomials. Coefficients of linear splines are simple interpolation coefficients between the 2 knots adjacent to the record and 0 between all other knots. Because at most 2 coefficients are nonzero for a given record, the system of equations with splines is sparser than with LP (Misztal, 2006). White et al. (1999) presented an application of cubic splines for analysis of test-day (TD) records of milk production. These functions consist of 2 terms: a) linear interpolation common for all segments and b) second derivatives at internal knots defining the deviation from the overall slope in each interval. However, while the first 2 coefficients (intercept and slope) are allowed to have different variances and are assumed to be correlated, the derivatives are assumed to have the same variances and to be independently distributed. This, in contrast to linear splines, leads to a highly structured covariance matrix. A review of the theory of spline functions and descriptions of different families of splines is given by Meyer (2005).

The Canadian test-day model (CTDM) is a 12-trait random regression model in which TD milk, fat and protein yields, and SCS from the first 3 lactations are analyzed (Schaeffer et al., 2000). Random regressions (additive genetic and permanent environmental) are based on LP of order 4. Only TD records from 5 through 305 DIM are considered in the model. Because less than 40% of Canadian Holstein cows end their lactation before 305 DIM, it is of interest to include into the CTDM the additional TD records recorded after 305 DIM. Higher accuracy of the overall 305-d EBV is expected from such evaluation because of utilization of more information. However, it can also lead to an increase in biasedness of the evaluation because the additional records are more likely to be influenced by the effect of cow’s pregnancy, which is currently not considered in the CTDM.

The objectives of this study were 1) to estimate variance components; 2) to test genetic evaluations from RRM based on LP versus RRM based on linear splines with 4, 5, and 6 knots in terms of a wide range of model comparison criteria; and 3) to identify the best model for a routine genetic evaluation of TD milk, fat, protein yield, and SCS records from 5 through 365 DIM from the first 3 lactation records.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Data
Variance components were estimated using a data set created by random sampling of 50 herds with at least 50 lactating cows in one year from the Canadian Holstein database, provided by Canadian Dairy Network. The data consisted of 96,756 TD milk, fat and protein yields, and SCS records of 6,094 cows from their first 3 lactations initiated from 1988 through 2006. Only TD records with all 4 traits present on a TD from 5 through 365 DIM were included. All cows were required to have a first lactation record. The pedigree file contained 18, 178 animals (Table 1).

Models
The data were analyzed with 4 multiple-trait, multiple-lactation RRM. The general formula common for all models was (Schaeffer et al., 2000):

where y_ijklmo was the oth test-day record of the mth cow for a trait i (TD milk, fat, protein yield, or SCS) in lactation j, HTD_ijk was the kth herd-test-date effect for a trait i and lactation j, _ijln was the nth fixed regression coefficient for a trait i and lactation j specific to the lth region-age-season class, q is the number of covariates, β_ijmn was the nth random regression coefficient for the additive genetic effect of cow m separate for trait i and lactation j, _ijmn is the nth random regression coefficient for the permanent environmental effect of cow m for trait i and lactation j, z(t) was a vector of covariates of size q describing the shape of lactation curve of fixed and random regressions evaluated at t DIM, and e is the residual. Twelve classes of residual variance with 30-d classes were defined for each lactation. The fixed and random regressions were fitted either with LP of order 4 (LEG) or with linear splines with 4 (SPL4), 5 (SPL5), or 6 (SPL6) knots. In the SPL4 model, the 4 knots were located at 5, 125, 245, and 365 DIM (Table 2). The knots were equally spaced (120 DIM) across lactation. The SPL5 model was created by adding one knot at 65 DIM, which created 2 shorter (60 DIM) intervals at the beginning of lactation and gave the spline function flexibility to better fit the peak of the lactation. The SPL6 model had compared with the SPL5 model one additional knot at 305 DIM. The SPL6 model was the most and the SPL4 the least complex model; both the LEG and SPL5 model had the same number of parameters.

For records located at a knot i covariables were defined as z_i = 1 and z_{1....i-1, i+1....q} = 0. Vector z had at most 2 nonzero elements, and sum of all elements of the vector was equal to 1.

For instance, for a record at t = 30, using a model with 4 knots located at t = {5, 65, 245, 365}, the vector of covariables was

In matrix notation each model can be described as

where y was a vector of observations; b was a vector of HTD and regression coefficients for region-age-season class; and u, p, and e were vectors of additive genetic, permanent environmental, and residual values, respectively. X, Z, W were incidence matrices.

The data were assumed to follow

and

where R was a block diagonal matrix with blocks of size 4 x 4 specific to a particular DIM interval and lactation. Residuals were assumed to have the same variance within residual intervals and heterogeneous variance between intervals, and residuals for different DIM were assumed to be uncorrelated. Matrices G and P were the random regression covariance matrices for the genetic and permanent environmental effects, respectively. The size of both matrices was (q x 3 x 4) by (q x 3 x 4), where q is the size of the vector of random regression coefficients, 3 is the number of lactations, and 4 is the number of traits. The matrix A was the additive genetic relationship matrix, and I was an identity matrix. A Bayesian approach via Gibbs sampling was applied to estimate model parameters. A single long chain of 100,000 samples was generated for each model. The first 20,000 samples were discarded as a burn-in, and the remaining samples were used to compute posterior means of model parameters. Convergence of Gibbs chains was monitored by visual inspections of plots of samples. The prior distributions of parameters of the model were

where _g and _p were genetic and permanent environmental covariance matrices, respectively; MVN denoted the multivariate normal distribution, and IW denoted the inverse Wishart distribution with degrees of belief and the scaling factor S. Values adopted for the scaling factors (

and

) were obtained by combining results from previous studies that considered TD with DIM 305 and assumptions made on variances past 305 DIM. Conservative degrees of belief were chosen to represent the vague prior information.

Two genetic evaluations were carried out for each of the 4 models that were compared by using variance components previously estimated from the variance component data. The first evaluation (EVAL01) was based on D01 and the second (EVAL06) used D06 data. Mixed model equations were solved by iteration on data with a preconditioned conjugate gradient algorithm using a block diagonal preconditioner (Lidauer et al., 1999). The convergence criterion was defined as the average relative difference between left and right hand sides of mixed model equations and was required to be <9.9 x 10^–13.

Model Comparison
The overall goodness of fit of all models was investigated by computing percentage of squared bias (PSB), as introduced by Ali and Schaeffer (1987), the correlation (RHO) of observed records in D06 with predicted records by EVAL06, and the residual variance (RV) using D06 data.

The PSB for the jth trait and nth lactation was defined as

where y_jnr was the rth observed record of jth trait and nth lactation,

was the rth predicted record of jth trait and nth lactation and o was the number of records. Values of the PSB, the RHO, and the RV of competing models were contrasted against values of the LEG model. However, none of these criteria takes into account the complexity of the model.

The deviance information criterion (DIC) introduced by Spiegelhalter et al. (2002) is widely used statistic for comparing models in a Bayesian framework. It combines both goodness of fit and degree of parameterization of the model and can be calculated easily from the samples generated by a Markov Chain Monte Carlo (MCMC) simulation. The DIC was defined by Spiegelhalter et al. (2002) as

where

was the posterior expectation of the Bayesian deviance

and

represented a measure of the fit of the model, and p_D was the effective number of parameters representing penalty for increasing model complexity defined as

where

is the vector of parameters of the model,

was the Bayesian deviance evaluated at the posterior mean of the parameters. Models having smaller DIC are favored because it indicates a better fit at lower degree of complexity.

Models were further compared by an error of prediction (ERP) of young bulls’ EBV defined by Sullivan et al. (2005) as

where ebv06 was the EBV calculated from EVAL06, pa01 is the average of animal’s parent EBV calculated by EVAL01 and n was the number of bulls with no daughters in D01 and at least 25 daughters in D06. For milk, fat, and protein yield, ebv06 was defined as 305-d EBV and for SCS as average daily EBV. To account for differences in average of EBV between the 2 genetic evaluations, EBV from EVAL06 were shifted by subtracting the average change in EBV from EVAL01 to EVAL06 for a set of bulls whose average EBV was not expected to change. The adjustment was based on 1,929 bulls with at least 25 daughters in EVAL01, no new daughters, and no more than 10 new granddaughters between D01 and D06.

To compare stability of EBV, average decline of EBV (AD) of 97 "newly proven bulls" [bulls that received their first official EBV at EVAL06 (i.e., had at least 20 daughters with TD records >120 DIM)] was calculated as follows:

Similarly as in ERP, ebv06 were shifted by subtracting the average change in EBV from EVAL01 to EVAL06.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Trends of daily additive genetic variances of milk and SCS are shown in Figures 1 and 2, respectively. Additive genetic variance functions from LEG followed the typical U shape with high variances at the beginning and at the end of lactation compared with the middle of lactation. Additive genetic variance functions of SPL had a parabolic shape between knots. The overall shape of the spline variance function was influenced by number of knots and correlations between random regression coefficients at knots. Functions were smoother in intervals with higher correlations between knots. The parabolic shape of the variance function does not have a biological background and is an artifact specific to linear splines. However, because test-day records are evenly distributed across lactation, the overestimation and underestimation of variances will be cancelled across lactation and, therefore, should not have a significant effect on genetic evaluation. The overall shapes of the variance functions were similar among models for milk, fat, and protein yield. However, LEG had noticeably higher variances at the extremes of lactation than SPL. Variances of LEG for SCS were higher than variances of SPL along the entire lactation. The implausible high variances at the extremes of lactations were reduced by SPL but were not completely removed. This finding suggests that the edge effect is not only caused by the type of function used for regressions but could also be caused by smaller number of TD records at the periphery of the lactation and factors not accounted for in the model such as preferential treatment, between herd variation in the shape of the lactation curve (Gengler and Wiggans, 2001; Jamrozik et al., 2001; de Roos et al., 2004), and the effect of pregnancy.

View larger version (17K): [in this window] [in a new window]   Figure 1. Posterior mean estimates of additive genetic variance of daily milk yield in first, second, and third lactation (LEG = model with Legendre polynomials; SPL4–6 = model with linear splines with 4 to 6 knots).

View larger version (18K): [in this window] [in a new window]   Figure 2. Posterior mean estimates of additive genetic variance of daily somatic cell score in first, second, and third lactation (LEG = model with Legendre polynomials; SPL4–6 = model with linear splines with 4 to 6 knots).

View larger version (16K): [in this window] [in a new window]   Figure 3. Posterior mean estimates of residual variance of daily milk yield in first, second, and third lactation (LEG = model with Legendre polynomials; SPL4–6 = model with linear splines with 4 to 6 knots).

View larger version (17K): [in this window] [in a new window]   Figure 4. Posterior mean estimates of residual variance of daily somatic cell score in first, second, and third lactation (LEG = model with Legendre polynomials; SPL4–6 = model with linear splines with 4 to 6 knots).

View this table: [in this window] [in a new window]   Table 3. Estimates of posterior expectation of the Bayesian deviance (), Bayesian deviance evaluated at the posterior mean of the parameters (), effective number of parameters (pD), deviance information criterion (DIC), and rank of models (Rank) by DIC1

View larger version (12K): [in this window] [in a new window]   Figure 5. Residual variance (RV) of first parity milk yield calculated as variance of differences between observed and predicted records in data set D06 (LEG = model with Legendre polynomials; SPL4–6 = model with linear splines with 4 to 6 knots).

View larger version (11K): [in this window] [in a new window]   Figure 6. Lactation curve of first age-parity class in D06 data set (LEG = model with Legendre polynomials; SPL4–6 = model with linear splines with 4 to 6 knots).

As given in Table 6, differences among models in ERP were not significant in first lactation for milk, fat, and protein. In the second and third lactation, all SPL were better than LEG, and SPL6 had smallest ERP. The SPL4 model, a model with lowest number of parameters, had surprisingly smaller ERP than LEG and also SPL5. For SCS, ERP values of SPL were smaller compared with LEG not only in the second and third but also in the first lactation. The lowest ERP for SCS were found with SPL6, being 6, 6, and 10% smaller than with LEG in first, second, and third lactation, respectively.

This study showed that spline functions with 6 knots are superior to LP of order 4 for fitting TD records with DIM up to 365 d in terms of goodness of fit and stability of 305-d EBV. Models with spline functions and knots located at different DIM than the functions used in this study could provide a better fit, but large changes in overall model performance are not expected with data where test-day records are evenly distributed across lactation because the differences will cancel out on average. Additional research is needed to identify optimal numbers and placement of knots for both fixed and random regression. Fixed regression with denser knots at the beginning and sparser knots in the middle and end of lactation could provide smoother curves than in this study.

Cubic splines or linear splines with a smoothing parameter can create smoother curves than linear splines only (de Boor, 1978). However, these functions were not tested against currently used LP because their covariance matrix is highly structured and implementation of such functions into the CTDM would not be straightforward.

Use of different functions for additive genetic and permanent environmental effects (spline functions for one random component, Legendre polynomials for the other component, or splines for both components but with different knots) were not investigated. López-Romero et al. (2004) suggested that only functions of the same order of fit should be considered to provide equal opportunity of variation for both components. Models with LP of higher order than 4 were not considered. In a preliminary study (unpublished), such models showed extremely high additive genetic variances in early and late lactation and large oscillatory patterns in the middle of lactation. This was in agreement with results by López-Romero et al. (2004).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
The SPL6 model had the best goodness of fit, measured by PSB, RHO, and RV, among the 4 compared models. No significant differences were found in ERP for milk, fat, and protein among models in first lactation, but noticeably smaller ERP were obtained with SPL compared with LEG in second and third lactation. In the case of SCS, SPL performed better in terms of ERP in all 3 lactations. The most stable EBV were obtained with SPL6. The DIC criterion also identified the SPL6 model as the most plausible model. Additive genetic variances at the extremes of lactations were smaller with SPL compared with LEG. The SPL5 and SPL6 were better models than LEG, especially for SCS. In general, SPL6 had the best overall performance among all compared models. Based on results from this study, the best model for genetic evaluation of TD records with = 365 DIM is the RRM using linear splines with 6 knots.

	ACKNOWLEDGEMENTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET: http://www.sharcnet.ca). Funding was provided by DairyGen Council of Canadian Dairy Network and NSERC of Canada. Three anonymous reviewers are thanked for their helpful comments and suggestions.

Received for publication December 12, 2007. Accepted for publication June 3, 2008.

	REFERENCES

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 


Ali, T. E., and L. R. Schaeffer. 1987. Accounting for covariances among test day milk yields in dairy cows. Can. J. Anim. Sci. 67:637–644.

de Boor, C. 1978. A Practical Guide to Splines. Springer-Verlag, NewYork, NY.

de Roos, A. P., A. G. Harbers, and G. de Jong. 2004. Random herdcurves in a test-day model for milk, fat, and protein production of dairy cattle in the Netherlands. J. Dairy Sci. 87:2693–2701.[Abstract/Free Full Text]

Druet, T., F. Jaffrezic, D. Boichard, and V. Ducrocq. 2003. Modelinglactation 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]

Gengler, N., and G. R. Wiggans. 2001. Variance of effects of lactationstage within herd by herd yield. J. Dairy Sci. 84(Suppl. 1):216. (Abstr.)

Interbull. 2007. Interbull routine genetic evaluation for dairy production traits, August 2007. http://www-interbull.slu.se/eval/aug07.html Accessed Nov. 30, 2007.

Jamrozik, J., D. Gianola, and L. R. Schaeffer. 2001. Bayesian estimation of genetic parameters for test day records in dairy cattle using linear hierarchical models. Livest. Prod. Sci. 71:223–240.[CrossRef]

Lidauer, M., I. Stranden, E. A. Mantysaari, J. Poso, and A. Kettunen. 1999. Solving large test-day models by iteration on data and preconditioned conjugate gradient. J. Dairy Sci. 82:2788–2796.[Abstract]

López-Romero, P., R. Rekaya, and M. J. Carabano. 2004. Bayesiancomparison of test-day models under different assumptions of heterogeneity for the residual variance: The change point technique versus arbitrary intervals. J. Anim. Breed. Genet. 121:14–25.[CrossRef]

Meyer, K. 2005. Random regression analyses using B-splines to model growth of Australian Angus cattle. Genet. Sel. Evol. 37:473–500.[CrossRef][Medline]

Misztal, I. 2006. Properties of random regression models using linear splines. J. Anim. Breed. Genet. 123:74–80.[CrossRef][Medline]

Pool, M. H., L. L. Janss, and T. H. Meuwissen. 2000. Genetic parameters of Legendre polynomials for first parity lactation curves. J. Dairy Sci. 83:2640–2649.[Abstract]

Schaeffer, L. R. 2004. Application of random regression models inanimal breeding. Livest. Prod. Sci. 86:35–45.[CrossRef]

Schaeffer, L. R., J. Jamrozik, G. J. Kistemaker, and B. J. Van Doormaal. 2000. Experience with a test-day model. J. Dairy Sci. 83:1135–1144.[Abstract]

Silvestre, A. M., F. Petim-Batista, and J. Colaco. 2005. Genetic parameter estimates of Portuguese dairy cows for milk, fat, and protein using a spline test-day model. J. Dairy Sci. 88:1225–1230.[Abstract/Free Full Text]

Spiegelhalter, D. J., N. G. Best, B. P. Carlin, and A. van der Linde. 2002. Bayesian measures of model complexity and fit. J. R. Stat. Soc. Ser. B Stat. Methodol. 64:583–639.[CrossRef]

Sullivan, P. G., J. W. Wilton, L. R. Schaeffer, G. J. Jansen, J. A. B. Robinson, and O. B. Allen. 2005. Genetic evaluation strategiesfor multiple traits and countries. Livest. Prod. Sci. 92:195–205.[CrossRef]

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]

Wilmink, J. B. M. 1987. Adjustment of lactation yield for age at calving in relation to level of production. Livest. Prod. Sci. 16:335–345.[CrossRef]

This article has been cited by other articles:

				J. Bohmanova, F. Miglior, and J. Jamrozik Use of test-day records beyond three hundred five days for estimation of three hundred five-day breeding values for production traits and somatic cell score of Canadian Holsteins J Dairy Sci, October 1, 2009; 92(10): 5314 - 5325. [Abstract] [Full Text] [PDF]

				J. Bohmanova, J. Jamrozik, and F. Miglior Effect of pregnancy on production traits of Canadian Holstein cows J Dairy Sci, June 1, 2009; 92(6): 2947 - 2959. [Abstract] [Full Text] [PDF]

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal 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 Bohmanova, J. Articles by Sullivan, P. G. Search for Related Content PubMed Articles by Bohmanova, J. Articles by Sullivan, P. G.