In the random regression model (RRM) for milk yield, by replacing empirical lactation curves with the five-order Legendre polynomial to fit fixed groups, the RRM can be transformed to a hierarchical model that consisted of a RRM in the first hierarchy with Legendre polynomials as individuals’ lactation curves resolved by restricted maximum likelihood (REML) software, and a multivariate animal model for phenotypic regression coefficients in the second hierarchy resolved by DMU software. Some empirical lactation functions can be embedded into the RRM at the first hierarchy to well fit phenotypic lactation curve of the average observations across all animals. The functional relationship between each parameter and time can be described by a Legendre polynomial or an empirical curve usually called submodel, and according to three commonly used criteria, the optimal submodels were picked from linear and nonlinear submodels except for polynomials. The so-called hierarchical estimation for the RRMs in dairy cattle indicated that more biologically meaningful models were available to fit the lactation curves; moreover, with the same number of parameters, the empirical lactation curves (MIL1, MIL5, and MK1 for 3, 4, and 5 parameters, respectively) performed higher goodness of fit than Legendre polynomial when modelling individuals’ phenotypic lactation curves.