regression - Why does a linear least squares fit appear to have a bias when applied to simple test data? - Cross Validated
43
$\begingroup$
I used python to generate a correlated data set for testing, and then plotted a basic linear least-squares fit. The result looked a bit strange to me, because the line doesn't really seem to pass "centrally" through the data. It looks a little "tilted":
So, instead, I then diagonalized the covariance matrix to obtain the eigenvector that gives the direction of maximum variance. This is shown by the black arrow in the figure below. This one does indeed point in the direction that I would expect:
So, I am looking for a bit of an intuitive explanation for what is going on here. I know that they are not measuring the same thing - the fit is minimizing the sum of the squared "errors", the vertical distances between the measured values and the fitted model. Whereas the eigenvector is chosen to maximize variance.
But I guess I am surprised by the result - I would've expected the fit to go through the center of the data cluster, and not have this bias. Is it something to do with the fact that minimizing the vertical distances somehow breaks the symmetry in the system? Is it inappropriate to apply a basic linear fit to such data?
Here is the code:
import numpy as np import matplotlib.pyplot as plt
def get_correlated_dataset(n, dependency, mu, scale): latent = np.random.randn(n, 2) dependent = latent.dot(dependency) scaled = dependent * scale scaled_with_offset = scaled + mu return scaled_with_offset[:, 0], scaled_with_offset[:, 1]
Generate Data
dependency = [[30, 30], [30, 2]] mu = [150,-100] scale = 1 x, y = get_correlated_dataset(20000, dependency, mu, scale)
Calculate direction of maximum variance (eigenvector of covariance matrix)
covariance_matrix = np.cov( np.stack((x, y), axis=0) ) EIGVALS,EIGVECS = np.linalg.eig( covariance_matrix ) vec = EIGVECS[:, np.argmax(EIGVALS)] # Get the eigenvector corresponding to largest eigenvalue
Magnitude of vector, 3 standard deviations, and mean (for plotting arrow)
vec_norm = np.sqrt(vec[0]2 + vec[1]2) sigma3 = np.sqrt(EIGVALS.max())*3 mu_x, mu_y = np.mean(x), np.mean(y)
Linear fit to the data
slope, intercept = np.polyfit(x, y, 1)
fig1,ax1 = plt.subplots(figsize=(8, 5)) ax1.plot( x , y ,'o',label="Generated Data") ax1.plot( x , xslope+intercept ,label="Linear Fit") ax1.quiver(mu_x, mu_y, sigma3vec[0]/vec_norm, sigma3*vec[1]/vec_norm, zorder=10,angles='xy', scale_units='xy',scale=1,label="Eigenvector of Covariance Matrix")
Cite
Follow
12.7k1414 gold badges2828 silver badges4747 bronze badges
asked Jan 2 at 11:37
59911 gold badge55 silver badges1313 bronze badges
$\endgroup$
18
- 4
$\begingroup$ I don't have any idea what is going on here, but there is a lot of overplotting in your graphs. You might want to plot either smaller dots or unfilled circles. I don't think this is the real issue here, but it seems like good practice $\endgroup$
Peter Flom
2026-01-02 12:21:00 +00:00
Commented Jan 2 at 12:21
- 7
$\begingroup$ This seems like the difference between total-least-squares and ols (i.e. ordinary least squares). We have a few relevant threads on the topic, including this recent one: stats.stackexchange.com/questions/672804/… $\endgroup$
mkt
– mkt
2026-01-02 12:27:59 +00:00
Commented Jan 2 at 12:27
- 5
$\begingroup$ This paper is also very relevant here: humans tend to overestimate the slope of a regression line when making a "visual fit"; we intuit more towards orthogonal least squares than ordinary (as mkt also mentioned). Importantly, the latter assumes no error exists along the X-axis. Eigen decomposition is effectively total least squares, so this confirms the bias described in the paper. $\endgroup$
PBulls
– PBulls
2026-01-02 12:44:00 +00:00
Commented Jan 2 at 12:44
- 11
$\begingroup$ Its important to know that the ordinary regression line doesnt try to fit the major axis of the bivariate ellipse-shape; instead it estimates the conditional expectation of Y|X=x, which is somewhat nearer to the horizontal. See the second plot here for an illustration: stats.stackexchange.com/a/226054/805 (indeed, the question and answers there may help you) $\endgroup$
Glen_b
– Glen_b
2026-01-02 13:32:54 +00:00
Commented Jan 2 at 13:32
- 3
$\begingroup$ @Ben It duplicates a lot of material (I linked to a search in my earlier comment), but I have been reluctant to close the thread due to the interest it has been generating. $\endgroup$
whuber
– whuber ♦
2026-01-02 19:28:30 +00:00
Commented Jan 2 at 19:28
Show 13 more comments
7 Answers 7
Sorted by: Reset to default
Highest score (default) Date modified (newest first) Date created (oldest first)
47
$\begingroup$
The regression line does, in fact, go through the "center," but since the errors are measured vertically, so is the "center."
This image shows it quite clearly (see the black vertical chord): 
Cite
Follow
57655 silver badges1919 bronze badges
answered Jan 2 at 21:50
60444 silver badges44 bronze badges
$\endgroup$
3
- 9
$\begingroup$ This is the best visualization in this question so far. $\endgroup$
justhalf
– justhalf
2026-01-03 03:47:01 +00:00
Commented Jan 3 at 3:47
- 1
$\begingroup$ A property of OLS is that the regression line goes through the point at $(\bar x, \bar y)$, which is the center of gravity of the cloud... $\endgroup$
jginestet
2026-01-05 02:41:22 +00:00
Commented Jan 5 at 2:41
- $\begingroup$ See Fig.1 and discussion in this paper of the same problem: crocodile.org/lord/baroaltitude.pdf $\endgroup$
krokodil
– krokodil
2026-01-05 04:13:50 +00:00
Commented Jan 5 at 4:13
Add a comment |
32
$\begingroup$
There are already loads of great answers here but I feel like it's necessary to add that TLS/the main PCA/ eigenvector, however you want to call it, can be something very different from what we normally understand as a linear regression line. They just converge towards each other the closer you are to the conditions in @Dominic Scoccheras Answer.
Consider:
$ X \sim \mathcal N(0,1) \text{ and } Y = 0.5X+e, e\sim\mathcal N(0,2) $
The resulting picture with n = 20000 looks like this:
If you wanted to use the main PCA, i.e. the red arrow, to predict $E[Y| X=2]$ you'd end up way of the scale at 26, while correct value is of course 1!
Here's the R-code:
n <- 20000 x <- rnorm(n) res_sd <- 2 y <- 0.5*x + rnorm(n, sd = res_sd)
just for checking syntax
eigen(M)[["vectors"]][, 1]
draw_eigen_arrow <- function(M, color){ eigen_vec <- eigen(M)[["vectors"]][, 1]*res_sd arrows(0, 0, eigen_vec[1], eigen_vec[2], col = color, lwd = 2) }
the following commented commands are for ploting with appropiate resolution
setwd("picture_directory")
png(width = 500sd(x), height = 500sd(y))
plot(x, y, col = rgb(0,0,0, alpha = 0.25)) M <- matrix(c(1, 0.25, 0.25, 0.25 + res_sd^2), nrow = 2) abline(a=0, b=0.5, lwd = 2) points(ellipse(M), col = 2, lwd = 2, type = "l") draw_eigen_arrow(M, 2)
dev.off()
eigen_vec <- eigen(M)[["vectors"]][, 1] eigen_vec[2]/eigen_vec[1]*2
Cite
Follow
answered Jan 2 at 19:10
5,69388 silver badges2929 bronze badges
$\endgroup$
Add a comment |
26
$\begingroup$
So, I am looking for a bit of an intuitive explanation for what is going on here.
The "intuitive" TL;DR is that OLS as a model inherently treats the dependent and the independent variable asymmetrically: only $y$ is assumed to be noisy, while $x$ is assumed to be exact.
I would've expected the fit to go through the center of the data cluster.
I'm not sure where that expectation is violated. The OLS fit does go through the "center", as long as you define the center as the point $(\mu_x, \mu_y)$.
Cite
Follow
answered Jan 2 at 13:13
2,27388 silver badges1515 bronze badges
$\endgroup$
3
- $\begingroup$ When I say "go through the center", I don't just mean passes through a center point at one place. What I really mean is "line is placed through the region of largest density", or "passes through the most number of points". This is what I think of as being a "best-fit" line. $\endgroup$
teeeeee
– teeeeee
2026-01-02 13:17:22 +00:00
Commented Jan 2 at 13:17
- 6
$\begingroup$ That's no doubt your expectation, but the least-squares fit pays no attention to local density and it is not guaranteed to go through any data points. The criterion is just in terms of a sum of squared errors to be minimised. $\endgroup$
Nick Cox
– Nick Cox
2026-01-02 13:25:22 +00:00
Commented Jan 2 at 13:25
- 1
$\begingroup$ @teeeeee Ah, I see. Yes, that can indeed be another definition of center, although I'd say that the more common interpretation is a centroid point. $\endgroup$
arpad
– arpad
2026-01-02 13:27:28 +00:00
Commented Jan 2 at 13:27
Add a comment |
23
$\begingroup$
If you want an intuitive understanding of what you observed, you first need to realize, as mentioned already by several contributors, that OLS treats Y and X asymmetrically.
To convince you of this, just plot, on your existing graph, in addition to the OLS line of Y against X, the OLS line of X against Y. You should see a different regression line. The two lines do in fact intersect at the point $(\bar x, \bar y)$ (that is a property of OLS; and that point is the center of gravity of your 2D cloud of points); and you will find that your "intuitive" best fit is somewhere between these 2 lines.
This is because OLS regression of Y against X minimizes the sum of vertical errors, while OLS regression of X against Y minimizes the sum of horizontal errors, but our visual intuition (our pattern matching brain) goes for minimizing the Euclidean (orthogonal) distance: hence the cognitive disconnect.
The plot below illustrates the scenario (it is based on some data I had, looking at the 6-minute walk distance (6MWD) as a function of age), when I was making that very point.
You can clearly see that both OLS regressions are slightly "biased" wrt what one might call the "main axis".
If you want a regression which matches your intuition, then look (as others have mentioned) at total least squares regression, or at Deming regression, or at orthogonal regression.
So the short answer to the question posed in your title is, "because it is supposed to have a bias". But the OLS does go through the center of gravity; just not at the angle our brains would have picked, and that is intentional, given what the OLS does.
Cite
Follow
answered Jan 2 at 18:08
17.4k22 gold badges1212 silver badges4747 bronze badges
$\endgroup$
2
- $\begingroup$ Great answer. When does it make sense, from a perspective of prediction, to "want" OLS vs TLS? $\endgroup$
Jonah
– Jonah
2026-01-04 22:54:46 +00:00
Commented Jan 4 at 22:54
- 2
$\begingroup$ @Jonah, OLS should be used when you can treat X or Y as "noise free" (even when it may not truly be so); e.g. you have a noisier but cheaper measurement system, and you want to calibrate it against a "gold standard" (but more $$) system, which has its own (but tolerated) issues. TLS is when you just want the "best" model which minimizes total noise (and X & Y are both treated as noisy). E.g. the 6MWD is used to predict a "physical age" (not real age), so the distance is considered "noise free" (even though, if the same subject was tested k times, there would be some noise). $\endgroup$
jginestet
2026-01-04 23:12:58 +00:00
Commented Jan 4 at 23:12
Add a comment |
15
$\begingroup$
As @mkt said this is the difference between OLS and total least squares. The solution to the slope of the line in OLS is $\beta_{ols}= {\sigma_{xy}} / {\sigma_{x}^2}$. What you also did was find the eigenvector with the maximum corresponding eigenvalue (PCA). For TLS we find the line with the smallest eigenvalue, which is the same line as the eigenvectors are orthogonal. The direction line from PCA is $\mathcal{l}(t)=\mu+t v$ ($\mu=(\mu_{x}, \mu_{y})^{T}$, $v$ is the eigenvector $v=(a,b)^{T}$) and the TLS line is $\mathcal{l}(t)=\mu+t (-b,a)^{T}$. These lines differ from OLS as the slope of the PCA line is $\beta_{PCA}=\bigg[\sigma_{y}^{2}-\sigma_{x}^{2}+\sqrt{(\sigma_{y}^{2}-\sigma_{x}^{2})^{2}+4\sigma_{xy}^{2}}\bigg]/{2\sigma_{xy}}$, which is only the same as the OLS slope if $(\sigma_{x}^{2}=\sigma_{y}^{2})$ or $\sigma_{xy}=\pm\sigma_{x}\sigma_{y}$ (perfect linear relationship $|\rho|=1$).
Here is your code edited to show this:
import numpy as np import matplotlib.pyplot as plt
def get_correlated_dataset(n, dependency, mu, scale): latent = np.random.randn(n, 2) dependent = latent.dot(dependency) scaled = dependent * scale scaled_with_offset = scaled + mu return scaled_with_offset[:, 0], scaled_with_offset[:, 1]
----------------------------
Generate data
----------------------------
np.random.seed(0)
dependency = [[30, 30], [30, 2]] mu = [150, -100] scale = 1
x, y = get_correlated_dataset(20000, dependency, mu, scale)
mu_x, mu_y = np.mean(x), np.mean(y)
OLS
ols_slope, ols_intercept = np.polyfit(x, y, 1)
PCA
cov = np.cov(x, y) eigvals, eigvecs = np.linalg.eig(cov) v = eigvecs[:, np.argmax(eigvals)] pca_slope = v[1] / v[0]
TLS
X = np.column_stack((x - mu_x, y - mu_y)) , , Vt = np.linalg.svd(X, full_matrices=False) direction = Vt[0] tls_slope = direction[1] / direction[0]
print("OLS:", ols_slope) print("PCA:", pca_slope) print("TLS:", tls_slope)
Cite
Follow
12.7k1414 gold badges2828 silver badges4747 bronze badges
answered Jan 2 at 13:02
17144 bronze badges
$\endgroup$
6
- $\begingroup$ Thanks, I'm not sure what it is about my generated test data which requires a method that allows errors in the x variable as well (i.e. total least squares)? $\endgroup$
teeeeee
– teeeeee
2026-01-02 13:27:59 +00:00
Commented Jan 2 at 13:27
- 3
$\begingroup$ +1. Welcome to Cross Validated! @teeeeee there's nothing about your data that requires a "method that allows errors in the x variable." But if you choose to use ordinary least squares instead of total least squares (TLS) you will get the "surprising" result that you found. Your intuition agrees with the TLS result, which treats variance inx and y symmetrically. $\endgroup$
EdM
– EdM
2026-01-02 16:38:59 +00:00
Commented Jan 2 at 16:38
- $\begingroup$ Thanks Dominic, I like what you are getting at in this answer. You are giving some explicit expressions for the slopes of the three methods (OLS, TLS, PCA), and making the point that they coincide under certain conditions. However, I am struggling to follow your notation, and what all your symbols mean. Are you parametrizing with the variable $t$? Also what is $\rho$? If you can explain a little more some of these symbols (and give references for your expressions) that would be helpful - please assume I know very little about this! $\endgroup$
teeeeee
– teeeeee
2026-01-05 09:49:09 +00:00
Commented Jan 5 at 9:49
- $\begingroup$ Some people have said in other answers that the difference between TLS and OLS comes if you assume there are errors in the $x$ variable or not - if $x$ has no error then OLS can be used. How can I connect this idea to your condition $\sigma_x^2=\sigma_y^2$ ? $\endgroup$
teeeeee
– teeeeee
2026-01-05 09:52:51 +00:00
Commented Jan 5 at 9:52
- $\begingroup$ sorry for not defining it explicitly in my answer but $\rho$ is the Pearson correlation coefficient en.wikipedia.org/wiki/Pearson_correlation_coefficient#Inference. To quote the Wikipedia article "The correlation coefficient can be derived by considering the cosine of the angle between two points representing the two sets of x and y co-ordinate data.[12] This expression is therefore a number between -1 and 1 and is equal to unity when all the points lie on a straight line." So it just means you have a perfect linear relationship when $\rho=\pm1$. $\endgroup$
Dominic Scocchera
2026-01-06 12:59:43 +00:00
Commented Jan 6 at 12:59
Show 1 more comment
14
$\begingroup$
From a geometric point of view the linear fit line is the shear/skew angle of an ellipse, while the axis that your brains wants to see is the axis of a rotated ellipse. Although shearing an ellipse creates a rotated ellipse, the angle of the resulting ellipse's axis does not equal the shear angle.
The image below you can see angle difference between the the bounds of the sheared/skewed ellipse (blue parallelogram) and the axes of the rotated ellipse (red dots and lines).
Also, an equivalent transformation (and resulting ellipse) can be achieved by shearing along the x-axis. In this case the shear angle is angle of the regression line where the x-variable contains variation, but the y-variable is fixed.
Cite
Follow
answered Jan 5 at 7:21
24155 bronze badges
$\endgroup$
1
- 2
$\begingroup$ +1 for pointing out the difference between shear transform and rotation. If the random distribution is given as a symmetric covariance matrix A, you can perform the symmetric eigen decomposition RDR'=A where R is a rotation matrix and D is diagonal (independent variances for X and Y). Alternatively, you can do the Cholesky decomposition LDL'=A where L is a shear matrix (lower triangular) and D is diagonal. Different decompositions, different results. $\endgroup$
Rainer P.
2026-01-05 14:53:34 +00:00
Commented Jan 5 at 14:53
Add a comment |
8
$\begingroup$
The answer lies in the data generating process. As other answers have pointed out, the simple OLS model assumes
$$ Y = \beta_0 + \beta_1 X + \varepsilon $$
and that $X$ and $\varepsilon$ are independent. It also assumes that we know $X$ with certainty.
Let’s consider the case where we have a perfectly fit regression and all error terms are zero. Because all the data points lie along the regression line, there is no variance in any other direction. So if we were to do PCA on this, the principal component would be in exactly the same direction as our regression line, matching your intuition. Once we start to add noise (i.e. increase the error term), which is strictly in the vertical direction, PCA is going to capture some of that for the first principal component, and leave what’s left over for the second. That is, the direction of maximum variance is going to be more extreme in the $Y$ direction. That’s what you’re observing in your own plots (only talking about the slope discrepancy here, not the centroid, which is covered by other answers).
Now consider the opposite extreme, when error dominates. Here it’s useful to pull apart the variance of $Y$:
$$ \operatorname{Var}(Y) = \operatorname{Var}(\beta_1 X + \varepsilon) = \operatorname{Var}(\beta_1 X) + \operatorname{Var}(\varepsilon) = \beta_1^2 \operatorname{Var}(X) + \operatorname{Var}(\varepsilon) $$
So in the case where
$$ \operatorname{Var}(\varepsilon) \gg \beta_1^2 \operatorname{Var}(X), $$
nearly all variation is vertical. The principal component direction approaches vertical as error variance increases.
Another source of confusion may be that the covariance ellipse only serves as a meaningful data contour when $(X, Y)$ are jointly normal. Here are examples where $X$ has different distributions:
Again, we see the same pattern: when the $X$ variance dominates, the principal component lines up far more with the regression line. But when the ellipsoid doesn’t coincide with such a tight contour of the data, we can more easily see that we’re dealing with different things in PCA and OLS.
Notice that when
$$ \operatorname{Var}(X) \gg \operatorname{Var}(\varepsilon), $$
PCA and OLS align. But the ellipse no longer describes the data's shape. It indicates direction, not density.
Here’s the code used to generate this data. Notice that the make_y_from_x function follows the precise assumptions of OLS. $Y$ is a linear function of $X$ plus i.i.d. Gaussian noise (if you want to have correlated residuals to try to violate the assumption, you can create that in the data generating process). Importantly, $X$ isn’t necessarily a random variable. $X$ could be from an experiment we designed with fixed values, or non-random for any other reason. It can also be random (but still known with certainty), so it’s important to think about this when specifying a model. (See What is a random variable and what isn't in regression models). Notice that the covariance ellipse treats $X$ and $Y$ symmetrically, but the data generating process is quite different.
Quick summary:
- When $\operatorname{Var}(X) \gg \operatorname{Var}(\varepsilon)$: PCA aligns with OLS (the textbook case)
- When $\operatorname{Var}(\varepsilon) \gg \operatorname{Var}(X)$: PCA approaches vertical
- $X$ is not necessarily random, and definitely doesn’t have to be normally distributed.
import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Ellipse
SLOPE = 1.0 INTERCEPT = -100.0 X_MEAN = 150.0 X_SCALE = 50.0
def make_y_from_x(x, noise_std, slope=SLOPE, intercept=INTERCEPT): return intercept + slope * x + np.random.randn(len(x)) * noise_std
def deterministic_binary_x(n): x_low = X_MEAN - X_SCALE x_high = X_MEAN + X_SCALE x = np.repeat([x_low, x_high], n // 2) if len(x) < n: x = np.append(x, x_high) return x
def deterministic_uniform_x(n): return np.linspace(X_MEAN - X_SCALE, X_MEAN + X_SCALE, n)
def bimodal_x(n): mix = np.random.rand(n) < 0.5 x = np.empty(n) x[mix] = np.random.randn(mix.sum()) * (X_SCALE / 3) + (X_MEAN - X_SCALE) x[~mix] = np.random.randn((~mix).sum()) * (X_SCALE / 3) + (X_MEAN + 3 * X_SCALE) return x
def triangular_x(n): return np.random.triangular( X_MEAN - X_SCALE, X_MEAN + X_SCALE, X_MEAN + X_SCALE, size=n )
def plot_from_cov(ax, x, y, title): x = (x - np.mean(x)) / np.std(x) y = (y - np.mean(y)) / np.std(y)
mu_x, mu_y = np.mean(x), np.mean(y)
cov = np.cov(x, y) eigvals, eigvecs = np.linalg.eigh(cov) idx = np.argmax(eigvals)
v = eigvecs[:, idx] angle = np.degrees(np.arctan2(v[1], v[0]))
a = 3 * np.sqrt(eigvals[idx]) b = 3 * np.sqrt(eigvals[1 - idx])
ax.plot(x, y, '.', alpha=0.05)
slope, intercept = np.polyfit(x, y, 1) xs = np.array([x.min(), x.max()]) ax.plot(xs, slope * xs + intercept, lw=2)
ax.quiver( mu_x, mu_y, a * v[0], a * v[1], angles='xy', scale_units='xy', scale=1, width=0.01, color='black', zorder=20 )
ell = Ellipse((mu_x, mu_y), 2 * a, 2 * b, angle=angle, edgecolor='red', facecolor='none', lw=2) ax.add_patch(ell)
xmin = min(x.min(), mu_x - a) xmax = max(x.max(), mu_x + a) ymin = min(y.min(), mu_y - a) ymax = max(y.max(), mu_y + a)
ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) ax.set_aspect('equal', adjustable='box')
ax.text(0.02, 0.98, title, transform=ax.transAxes, ha='left', va='top', fontsize=11, weight='bold')
-------- Joint normal: variance domination --------
fig, axs = plt.subplots(2, 2, figsize=(12, 8)) axs = axs.ravel()
n = 20000 beta = 1.0
settings = [ (0.3, "var(X) >> var(ε)"), (5.0, "var(ε) >> var(X)") ]
for ax, (sigma, label) in zip(axs, settings): X = np.random.randn(n) Y = beta * X + np.random.randn(n) * sigma plot_from_cov(ax, X, Y, label)
for ax in axs[len(settings):]: ax.axis("off")
plt.tight_layout()
-------- Different X distributions --------
fig, axs = plt.subplots(2, 2, figsize=(12, 8)) axs = axs.ravel()
make_xs = [ deterministic_binary_x, deterministic_uniform_x, bimodal_x, triangular_x ]
sigma = 30.0
for ax, make_x in zip(axs, make_xs): X = make_x(n) Y = beta * X + np.random.randn(n) * sigma plot_from_cov(ax, X, Y, make_x.name)
plt.tight_layout() plt.show()
Cite
Follow
answered Jan 5 at 16:25
18133 bronze badges
$\endgroup$
Add a comment |
Your Answer
Thanks for contributing an answer to Cross Validated!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
Draft saved
Draft discarded
Sign up or log in
Sign up using Google
Sign up using Email and Password
Submit
Post as a guest
Name
Required, but never shown
Start asking to get answers
Find the answer to your question by asking.
Explore related questions
See similar questions with these tags.