Concerning bivariate least squares linear regression, the classical approach pursued for functional models in earlier attempts is reviewed using a new formalism in terms of deviation (matrix) traces. Within the framework of classical error models, the dependent variable relates to the independent variable according to the usual additive model. Linear models of regression lines are considered in the general case of correlated errors in X and in Y for heteroscedastic data. The special case of (C) generalized orthogonal regression is considered in detail together with well known subcases. In the limit of homoscedastic data, the results determined for functional models are compared with their counterparts related to extreme structural models. While regression line slope and intercept estimators for functional and structural models necessarily coincide, the contrary holds for related variance estimators even if the residuals obey a Gaussian distribution, with a single exception. An example of astronomical application is considered, concerning the [O/H]-[Fe/H] empirical relations deduced from five samples related to different stars and/or different methods of oxygen abundance determination. For selected samples and assigned methods, different regression models yield consistent results within the errors for both heteroscedastic and homoscedastic data. Conversely, samples related to different methods produce discrepant results, due to the presence of (still undetected) systematic errors, which implies no definitive statement can be made at present. A comparison is also made between different expressions of regression line slope and intercept variance estimators, where fractional discrepancies are found to be not exceeding a few percent, which grows up to about 20% in presence of large dispersion data.