Is it possible to decompose fitted residuals into bias and variance, after fitting a linear model?

Question

I'd like to classify data points as either needing a more complex model, or not needing a more complex model. My current thinking is to fit all the data to a simple linear model, and observe the size of the residuals to make this classification. I then did some reading about the bias and variance contributions to error, and realized that if I could calculate bias directly, it might be a better measure then working with the total error (residual or standardized residual).

Is it possible to estimate bias directly with a linear model? With or without test data? Would cross validation help here?

If not, can one use an averaged bootstrapping ensemble of linear models (I think it's called bagging) to approximate bias?

Answer 1

You generally can't decompose error (residuals) into bias and variance components. The simple reason is that you generally don't know the true function. Recall that $bias(\hat f(x)) = E[\hat f(x) - f(x)],$ and that $f(x)$ is the unknown thing you wish to estimate.

What about bootstrapping? It is possible to estimate the bias of an estimator by bootstrapping, but it's not about bagging models, and I don't believe there is a way to use the bootstrap to assess the bias in $\hat f(x),$ because bootstrapping is still based on some notion of the Truth and can't, in spite of the origins of its name, create something from nothing.

To clarify: the bootstrap estimate of bias in the estimator $\hat \theta$ is $$\widehat{bias}_B = \hat\theta^*(\cdot) - \hat \theta,$$

with $\hat\theta^*(\cdot)$ being the average of your statistic computed on $B$ bootstrap samples. This process emulates that of sampling from some population and computing your quantity of interest. This only works if $\hat\theta$ could in principle be computed directly from the population. The bootstrap estimate of bias assesses whether the plug-in estimate—ie just making the same computation on a sample instead of in the population—is biased.

If you just want to use your residuals to evaluate model fit, that is entirely possible. If you, as you say in the comments, want to compare the nested models $f_1(x) = 3x_1 + 2x_2$ and $f_2(x) = 3x_1 + 2x_2 + x_1x_2$, you can do ANOVA to check whether the larger model significantly reduces sum of squared error.

Answer 2

One situation where you can get an estimate of the decomposition is if you have replicated points (i.e. to have more than one response for various combinations of the predictors).

This is mostly limited to situations where you have control of the independent variables (such as in experiments) or where they're all discrete (when there are not too many x-combinations and you can take a large enough sample that x-value combinations get multiple points).

The replicated points give you a model-free way of estimating the conditional mean. In such situations there's the possibility of decomposition of the residual sum of squares into pure error and lack of fit, but you also have direct (though necessarily noisy) estimates of the bias at each combination of x-values for which you have multiple responses.

Answer 3

In the somewhat more complex Kalman filtering realm, sometimes people test the residuals (observed measurements minus predicted measurements) to look for model changes or fault conditions. In theory, if the model is perfect, and the noise is Gaussian, then the residuals should also be Gaussian with zero mean and also be consistent with a predicted covariance matrix. People can test for nonzero mean with sequential tests like a Sequential Probability Ratio Test (SPRT). Your situation is different because you have a fixed batch of data rather than a steady stream of new data. But the basic idea of looking at the sample distribution of the residuals might still apply.

You indicate that the process you are modeling might change occasionally. Then, to do more with the data you have, you'd probably need to identify other factors causing that change. Consider 2 possibilities: (1) maybe you need local models rather than one global model, e.g., because there are severe nonlinearities only in some operating regions, or (2), maybe the process changes over time.

If this is a physical system, and your samples aren't taken huge time intervals apart, it's possible that these process changes persist over significant time periods. That is, true model parameters may occasionally change, persisting for some time period. If your data is time stamped, you might look at residuals over time. For instance, suppose you have fit y = Ax + b using all your data, finding A and b. Then go back and test the residual sequence r[k] = y[k] - Ax[k] - b , where k is an index corresponding to times in sequential order. Look for patterns over time, e.g., periods where summary statistics like ||r[k] || stays higher than normal for some time. Sequential tests would be the most sensitive to detecting sustained bias sorts of errors, something like SPRT or even CUSUM for individual vector indices. This could point to time periods where you need to consider more complex models.

Answer 4

The answer is no, because bias and variance are attributes of model parameters, rather than the data used to estimate them. There is a partial exception to that statement that pertains to bias and variance varying (ha!) through the predictor space; more on that below. Note that this has absolutely nothing to do with knowing some "true" function relating the predictors and response variables.

Consider the estimate of $β$ in a linear regression, $\hatβ=(X^TX)^{-1}X^TY$, where $X$ is an $N×P$ matrix of predictors, $\hatβ$ is a $P×1$ vector of parameter estimates, and $Y$ is an $N×1$ vector of responses. Let's assume for argument's sake that we have an infinite population of data from which to draw (this is not completely ridiculous, by the way -- if we were actively recording data from some physical process we could record predictor and response data at a rapid rate, thus practically satisfying this assumption). So we draw $N$ observations, each consisting of a single response value and a value for each of the $P$ predictors. We then compute our estimate of $\hatβ$ and record the values. Let us then take this entire process and repeat it $N_{iter}$ times, each time making $N$ independent draws from the population. We will accumulate $N_{iter}$ estimates of $\hatβ$ over which we can compute the variance of each element in the parameter vector. Note that the variance of these parameter estimates is inversely proportional to $N$ and proportional to $P$, assuming orthogonality of the predictors.

The bias of each parameter can be estimated similarly. While we may not have access to the "true" function, let's suppose we can make an arbitrarily large number of draws from the population in order to compute $\hatβ_{best}$, which will serve as a proxy for the "true" parameter value. We'll assume that this is an unbiased estimate (ordinary least squares) and that the number of observations used was sufficiently large such that the variance of this estimate is negligible. For each of the $P$ parameters, we compute $\hatβ_{best_j}-\hatβ_j$, where $j$ ranges from $1$ to $N_{iter}$. We take the average of these differences as an estimate of the bias in the corresponding parameter.

There are corresponding ways of relating bias and variance to the data itself, but they're a little more complicated. As you can see, bias and variance can be estimated for linear models, but you will require quite a bit of hold-out data. A more insidious problem is the fact that once you start working with a fixed dataset, your analyses will be polluted by your personal variance, in that you'll have already begun wandering through the garden of forking paths and there's no way of knowing how that would replicate out-of-sample (unless you just came up with a single model and ran this analysis and committed to leaving it alone after that).

Regarding the matter of the data points themselves, the most correct (and trivial) answer is that if there is any difference between $Y$ and $\hat{Y}$, you need a more complex model (assuming that you could correctly identify all the relevant predictors; you can't). Without going into a boring treatise on the philosophical nature of "error," the bottom line is that there was something going on that caused your model to miss its mark. The problem is that adding complexity increases variance, which will likely cause it to miss the mark on other data points. Therefore, worrying about error attribution at the individual data point level is not likely to be a fruitful endeavor. The exception (mentioned in the first paragraph) stems from the fact that bias and variance are actually functions of the predictors themselves, so you may have large bias in one part of the predictor space and smaller bias in another (same for variance). You could assess this by computing $Y-\hat{Y}$ many times (where $\hat{Y}=X\hatβ$ and $\hatβ$ was not estimated based on $Y$) and plotting its bias (average) and variance as a function of the values of $X$. However, I think that's a pretty specialized concern.