Should we always minimize squared deviations if we want to find the dependency of mean on features?

Question

It is known that for any distribution a minimization of squared deviations from a value $c$ gives mean of the distribution. In other words, if we generate many values $y_i$ using our distribution, and then try to find $c$ that minimizes the following expression:

$ \lim_{N \to \infty} \frac{1}{N}\sum_{i}^N (y_i - c)^2 $

then we find that the minimum is achieved when $c = \mu$, where $\mu$ is the mean of the given distribution. And for any $N$ the above minimization gives the sample mean.

---

Now, let's assume tat we have a distribution which depends on a single real-valued feature $x$, and we have a data-set generated by this distribution: ($x_1$, $y_1$), ($x_2$, $y_2$), $\dots$, ($x_n$, $y_n$). Then we try to find a function that gives dependency of the mean of this distribution on the feature:

$ \mu = \mu(x, p_1, p_2, \dots, p_k) $

where $p_i$ are model parameters that we want to find.

Obviously, that we can get the desired function (values of the parameters) if we choose them by minimization of squared deviations:

$ \frac{1}{N} \sum_{i}^{N} [y_i - \mu(x_i, p_1, p_2, \dots, p_k)]^2 $

In other words, we can find the correct answer (correct values of the parameters $p_i$) under condition that $N$ goes to infinity and that the real function, that we are searching, is covered by the space of functions that we consider.

To summarize, the minimization of squared deviations can give not only mean of a given constant distribution, it can also gives correct dependency of the mean on a feature.

---

In reality $N$ is never infinity and can be rather small. As a result minimization of squared deviation will not give exactly correct values of the model parameters $p_i$. For another sample of the same size we will give other estimations of the parameters $p_i$. In general case, we can say that our estimations of the parameters are somehow distributed.

If we, instead of squared deviations minimize some other measure, we will get other distributions of our model parameters $p_i$. So, different measure of accuracy give different distributions in the parameters space for a given fixed sample size.

These distributions are characterized by their "width" and possible a "shift" relative to the correct point (true values of the model parameters).

Can it be the case that a use of something different from squared deviation gives a better estimate of model parameters (for example more accurate (smaller width) and with smaller or no systematic from the correct model parameters?

Answer 1

Minimizing the MSE in the cases you describe indeed produces a consistent estimator for the model parameters. The consistency is related to the fact that the derivative of the MSE, and therefore the first order optimality condition, is linear in the observations $y_i$. So, if you want a consistent estimator without any additional information on the distribution of the samples, MSE is the only option.

On the other hand, if you have a parametric form for the distribution of $y_i$, then the maximum likelihood estimator (MLE) is also consistent, but also asymptotically unbiased and efficient, meaning that it has the smallest variance among all unbiased estimators (in the limit $n \to \infty$). In that sense it is the most "accurate" estimator possible.

Answer 2

NO

It is important to keep in mind that an estimator of a parameter can take on many forms. In fact, constants can be estimators! Consequently, we might find that calculating something other than the empirical mean might have desirable properties for estimating the mean.

In this post, I give an example where estimating the conditional median, by minimizing MAE, gives better estimates of the regression parameters than the OLS estimates, in the sense that the estimator is unbiased (as is the case with the OLS estimator) but has lower variance than the OLS estimator.

Another answer here mentions the Gauss-Markov theorem. As Richard Hardy explains in the answer to my linked question, the MAE minimizer is nonlinear. Thus, Gauss-Markov does not apply. It is fine for minimization of MAE to result in an unbiased estimator that has lower-variance than OLS.

EDIT

Another answer of mine shows when minimizing in-sample MAE results in lower out-of-sample MSE than minimizing in-sample MSE.

EDIT 2

Let's check out a simulation and visualization. In the code, I simulate a heavy-tailed $t_{1.1}$ error term. At each iteration, I calculate the OLS coefficients (with lm) and the MAE-minimizing coefficients (with rq). Then I plot those regression lines, along with the true regression line.

library(quantreg)
set.seed(2022)

N <- 50
B <- 100
beta0 <- 2
beta1 <- -3
x <- seq(0, 1, 1/(N - 1))
yhat <- beta0 + (beta1 * x)
Q0 <- Q1 <- L0 <- L1 <- rep(NA, B)
for (i in 1:B){
  
  y <- yhat + rt(N, 1.1)
  L <- lm(y ~ x)
  L0[i] <- summary(L)$coef[1, 1]
  L1[i] <- summary(L)$coef[2, 1]
  
  Q <- quantreg::rq(y ~ x, tau = 0.5)
  Q0[i] <- summary(Q)$coef[1, 1]
  Q1[i] <- summary(Q)$coef[2, 1]
  
}

par(mfrow = c(2, 1))
plot(
  x, 
  yhat,
  type = 'l', 
  lty = 2, 
  ylim = c(min(L0 + L1*x, Q0 + Q1*x), max(L0 + L1*x, Q0 + Q1*x)),
  main = "OLS"
  )
for (i in 1:B){
  lines(x, L0[i] + L1[i] * x, col = 'red')
}
lines(x, yhat, type = 'l', lty = 2)
#
plot(
  x, 
  yhat,
  type = 'l', 
  lty = 2, 
  ylim = c(min(L0 + L1*x, Q0 + Q1*x), max(L0 + L1*x, Q0 + Q1*x)),
  main = "MAE"
)
for (i in 1:B){
  lines(x, Q0[i] + Q1[i] * x, col = 'red')
}
lines(x, yhat, type = 'l', lty = 2)
par(mfrow = c(1, 1))

The red estimates of the conditional means are much more reasonable when absolute loss is minimized.

EDIT 4

Another example is in "classification" problems with discrete outcomes (say binary for now). The typical loss function minimized is log loss ("crossentropy" in some circles), which corresponds to maximum likelihood estimation in logistic regression. Our Frank Harrell has a strong opinion about minimizing this loss function as opposed to minimizing square loss.

$$ \text{Log Loss}\\ L(y, p) = -\dfrac{1}{N} \sum_{i = 1}^N \bigg[ y_i\log(p_i) + (1 - y_i)\log(1 - p_i) \bigg] $$

EDIT 5

Finally, there is the James-Stein estimator, which shows that the OLS solution to linear regression is inadmissible for any reasonable sample size for doing regression, despite the Gaussian conditional distribution. That is, even the maximum likelihood estimator is inadmissible due to being dominated by James-Stein.

Answer 3

A similar question (if not the same) is: If the predicted value of machine learning method is E(y | x), why bother with different cost functions for y | x?

The theoretical mean of a distribution minimizes the squared error, but that does not mean that the sample mean is always the best estimator with regards the squared error loss function. The sample mean has a statistical variation.

In an answer to the above question an example is givem that shows how the sample median is performing better than the sample mean when the errors are Laplace distributed. Below is a copy of the image:

• Another example question is: Why is the Median Less Sensitive to Extreme Values Compared to the Mean?

  The median can be a better estimator in the case of distributions with outliers.

• Related is also: Could a mismatch between loss functions used for fitting vs. tuning parameter selection be justified?

  The answer to that question explains how, when we wish to have an estimator that optimizes the mean squared error, then it doesn't mean that we need to use the squared error loss function in fitting/training the model.

Asside from using different estimators like the median or a maximum likelihood estimators, there is also the concept of biased estimators that can improve the expectation of the loss. Examples are regularisation (Ridge regression, lasso regression), Bayesian estimators, shrinking.

Answer 4

Your estimator is the OLS estimator in nonlinear regression

Your problem is essentially just the OLS estimation problem in nonlinear regression. To see this, suppose you have nonlinear regression model of the form:

$$Y_i = \mu(x_i, \mathbf{p}) + \varepsilon_i \quad \quad \quad \quad \quad \varepsilon_i \sim \text{IID Dist}(\text{Mean} = 0),$$

with error terms that are independent of the explanatory variables.$^\dagger$ Taking conditional expectations of both sides yields the equation:

$$\mu(x_i, \mathbf{p}) = \mathbb{E}(Y_i|X_i=x_i, \mathbf{p}).$$

This shows that the function $\mu$ is the true regression function which represents the conditional expected value of the response variable conditional on the explanatory variable. Now, if you have $n$ observations from this model, then the OLS estimator for $\mathbf{p}$ is defined by the optimisation requirement:

$$\mu(x_i, \hat{\mathbf{p}}_n) = \min_\mathbf{p} \sum_{i=1}^n (y_i-\mu(x_i, \mathbf{p}))^2,$$

which is the optimisation requirement in your question. So, your problem is about OLS estimation of the parameter $\mathbf{p}$, and you want to know the consistency properties of the OLS estimator $\hat{\mathbf{p}}_n$.

Consisteny properties in regression models is something that has been examined extensively in the statistical literature. Under broad conditions on the sequence of explanatory variables, the estimator $\hat{\mathbf{p}}_n$ is a consistent estimator for $\mathbf{p}$. For OLS estimation in nonlinear regression, these conditions are extensions of the "Grenander conditions" for consistency in linear regression (see e.g., Richardson and Bhattacharyya 1990). The exact conditions are quite technical, but heuristically, they require the sequence of explanatory variables $x_1,x_2,x_3,...$ to be such that the "influence" of any finite set of data points tends to zero as $n \rightarrow \infty$.

If you have the mathematical background to do so, I recommend you read the linked paper to get an understanding of the consistency conditions in nonlinear regression. For a simpler place to get started you can have a look at some answers on this site that look at the Grenander conditions for OLS consistency in linear regression (see e.g., here).

---

$^\dagger$ In this formulation of the nonlinear regression model I allow any error distribution with zero mean. In the special case where the error distribution is normal the OLS estimator will correspond to the MLE. The general case is used here because you have not specified an error distribution in your problem.

Answer 5

Yes, it is possible for estimators obtained by minimizing some different than squared deviation to give a better estimator of model parameters.

The question of whether a given estimator can be beaten by others is studied in statistical decision theory. I'll lay out the basics below then give two examples.

A framework to compare estimators

Suppose we have a data vector $Z \in \mathbb{R}^{p}$, which gives the covariates and response of a single case. Suppose we are estimating a parameter $\mathbf{p}$. Using i.i.d. data $Z_1, \dots, Z_n$, we can estimate the parameter using $\delta(Z_1, \dots, Z_n)$ for some functional $\delta$. We can evaluate the closeness of an estimate to the parameter using the squared-error loss function $\|\delta(Z_1, \dots, Z_n) - \mathbf{p}\|^2$.

Notice the loss function depends on the data as well as the parameter. We can form the risk $R(\delta, \mathbf{p}) = \mathbb{E} \|\delta(Z_1, \dots, Z_n) - \mathbf{p}\|^2$ by averaging over the data. For a given estimation rule $\delta$, the risk $R$ tells us how close we can expect the estimator to be to the truth. Lower is better. Using a simple computation, we can prove that $$R(\delta, \mathbf{p}) = \| \mathrm{Bias}\, \delta(Z_1, \dots, Z_n) \|^2 + \mathrm{trace} \, \mathrm{Var} \, \delta(Z_1, \dots, Z_n),$$ i.e. that the risk trades off the bias (shift) and the variance (width) of the estimator. This decomposition shows that minimizing the variance among unbiased estimators does not lead to the estimator which is expected to be closest to the parameter - instead this tradeoff must be minimized.

We can compare the quality of different estimators $\delta_1$ and $\delta_2$ by comparing the risk curves $\mathbf{p} \mapsto R(\delta_1, \mathbf{p})$ and $\mathbf{p} \mapsto R(\delta_2, \mathbf{p})$. For example, if $R(\delta_1, \mathbf{p}) \leq R(\delta_2, \mathbf{p})$ for all parameters $\mathbf{p}$, this means that the estimation rule $\delta_1$ is will be closer on average to the parameter $\mathbf{p}$ than $\delta_2$ for all possible parameter values. This means that $\delta_1$ dominates $\delta_2$.

A first example

Let's consider the simple example given as the first in OP's question. Here the data $Z=Y \in \mathbb{R}^1$ so that there is only one variable. Let us further assume that $Y = \mu + \epsilon$ for normally distributed $\epsilon$.

The estimator formed by minimizing the empirical squared error is the sample mean $\bar{Y} = \frac{1}{n} \sum_{i=1}^n Y_i$. It is a classical result that the estimator $\delta(Y_1, \dots, Y_n) = \bar{Y}$ is admissible. This means that there does not exist any other estimator $\delta_2(Y_1, \dots, Y_n)$ which dominates the sample mean.

A second example

Now let's consider a linear regression example. Let the data be given by $Z=(y, x)$, where the outcome $y$ is scalar and the covariates $x \in \mathbb{R}^{p-1}$. Assume that $y = x^T \beta + \epsilon$, where $\epsilon$ is normally distributed. Let $\beta$ be the target of inference.

The estimate formed by minimizing the empirical squared error is the OLS estimator $\hat\beta$. When $p > 3$, it turns out this estimator is not admissible: that is, there are other estimators which are always closer on average to the true parameter value $\beta$, regardless of its (unknown) value. A classical example is the James-Stein estimator, which equals $s(Z_1, \dots, Z_n) \hat\beta$ for a suitably chosen data-dependent shrinkage term $s \in (0,1)$.

Conclusion

Basing an estimating equation on the loss function does not necessarily lead to finite sample optimality of the estimator. OP is right to question the basis of the procedure.

Answer 6

Can it be the case that a use of something different from squared deviation gives a better estimate of model parameters (for example more accurate (smaller width) and with smaller or no systematic from the correct model parameters?

The keyword in your question is "better". To rigorously define "better", you need a loss function.

If your loss function is the squared loss and its argument is the prediction error that you commit by using $\mu$ as a prediction of $y$, then estimating $\mu$ by minimizing squared deviations is the right thing to do. In technical terms, you are approximating the expected loss (or statistical risk) with the sample loss (empirical risk). The latter converges to the former (by the law of large numbers), provided that the minimization does not introduce some serious discontinuities at the limit.

If your loss function is not the squared loss, or it is not defined on prediction errors (e.g., it is defined on parameter estimation errors), then the minimization of squared deviations can be sub-optimal. In the first case (different loss function), you just need to be coherent. For example, if you use the absolute loss, then you need to minimize absolute deviations. In the second case (loss defined on something different from prediction errors), things get complicated, and I dare say that there are very few analytical results. One of the few results is the Gauss-Markov theorem for linear regression models: under certain assumptions, even if you minimize a quadratic loss function defined over prediction errors, you achieve optimality also with respect to a quadratic loss function defined over parameter estimation errors.