Abstract Legendre polynomials of orders 3 to 8 in random regression models (RRM) for first-lactation milk production in Canadian Holsteins were compared statistically to determine the best model. Twenty-six RRM were compared using LP of order 5 for the phenotypic age-season groupings. Variance components of RRM were estimated using Bayesian estimation via Gibbs sampling. Several statistical criteria for model comparison were used including the total residual variance, the log likelihood function, Akaike's information criterion, the Bayesian information criterion, Bayes factors, an information-theoretic measure of model complexity, and the percentage relative reduction in complexity. The residual variance always picks the model with the most parameters. The log likelihood and information-theoretic measure picked the model with order 5 for additive genetic effects and order 7 for permanent environmental effects. The currently used model in Canada (order 5 for both additive and permanent environmental effects) was not the best for any single criterion, but was optimal when considering all criteria.