Why does the sum of residuals equal 0 from a graphical perspective? I've seen the proof for why in least squares regression the sum of residuals is always equal to 0, and I kind of understand why from that algebraic perspective. Basically, you're finding the minimum of the least squares equation, so you take the partial derivative and set it equal to 0, right?
To make things more confusing, I've heard from some sources that the sum doesn't necessarily equal 0, but it's expected value is 0, so there's some variation. (That makes more sense to me, because you're dealing with random variables that naturally have random variation in the points.)
What I can't wrap my head around is that from a graphical perspective, if I imagine a set of points with a regression line through them, how can you guarantee that the sum of residuals ends up being 0? Can someone help me understand from a point of view that isn't necessarily the algebraic proof?
 A: To fit $\hat y = b_1 + b_2 x$ by OLS we minimise the residual sum of squares, $RSS = \sum_{i=1}^n e_i^2$. As the question states, we can do this analytically by setting
$$\frac{\partial RSS}{\partial b_1}=0 , \qquad \frac{\partial RSS}{\partial b_2}=0 $$
and solving the resulting normal equations. Note that from the normal equations we can deduce $\sum_{i=1}^n e_i = 0$ and $\sum_{i=1}^n x_i e_i = 0$; the latter is equivalent (think about the inner or "dot" product) to stating that the vector of observations $\mathbf{x} = (x_1, x_2, \dots, x_n)^t$ in $\mathbb{R}^n$ is orthogonal (perpendicular) to the vector of residuals $\mathbf{e} = (e_1, e_2, \dots, e_n)^t$, and analogously your requirement that the sum of residuals is zero is equivalent to the statement that $\mathbf{e}$ is orthogonal to the vector of ones, $\mathbf{1}_n=(1,1,\dots,1)^t$.

Both these results can be seen geometrically from knowing the design matrix $\mathbf{X}$ includes a column of ones to represent the intercept term and another column for the $x_i$ data, and that the vector of residuals is orthogonal to each column of the design matrix because the hat matrix $\mathbf{H}$ is an orthogonal projection onto the column space of $\mathbf{X}$. For more on how to interpret the diagram, see Geometric interpretation of multiple correlation coefficient $R$ and coefficient of determination $R^2$. But that takes place in the $n$-dimensional subject space; it would be nice to develop some intuition from the scatter plot itself.

To illustrate OLS geometrically, draw a vertical line segment from each point $(x_i, y_i)$ to its fitted value on the regression line $(x_i, \hat y_i)$, then draw a square with this side. The length of this segment is the magnitude of the residual $|e_i$| (note $e_i$ is positive for points above the line, e.g. with the blue square, negative for points below, e.g. with the red square) and the area of the square is $e_i^2$. I have illustrated two examples points, and for convenience placed the squares on whichever side avoids overlapping the regression line. The RSS is the sum of the areas of all the squares (both red and blue squares counting as positive area), and to find the OLS solution we seek $b_1$ and $b_2$ that minimise this area. 
Now take any regression line (not necessarily the OLS one) and consider translating its intercept up by $\delta b_1$. I have drawn the new regression line as dotted, and overlaid the residual squares based on both new and original fitted values.

For the point on the left with a positive residual, this shift in the regression line has reduced the area of the residual square — this is logical since the observed point lay above the line, and the regression line has moved up closer to it, so the fit is better. Two $\color{blue}{\text{blue}}$ rectangular strips (which I consider to extend the whole length of the larger, original $e_1$ by $e_1$ square) have been cut off the left and bottom sides. But these overlap (the grey square at bottom left) so subtracting both rectangles has double-counted the grey square, and we need to add this area back on once again. Overall the new (reduced) residual square is found from the original square, minus the two blue rectangular strips, plus the grey square.
For points below the line with a negative residual, raising the line has worsened the fit and made the residual more negative; the residual square has grown by two $\color{red}{\text{red}}$ rectangular strips (running across the sides of the smaller, original $e_2$ by $e_2$ square) plus the grey square in the upper right.

If we dissect the residual squares for all data points in this way, then total the results,
$$\color{Gold}{\text{New RSS}} = \text{Old RSS} + \sum \color{red}{\text{red rectangles}} - \sum \color{blue}{\text{blue rectangles}} + \sum \color{grey}{\text{grey squares}}$$
From the diagram it's clear $\delta e_i = - \delta b_1$ for all points — raising the intercept increases each fitted $\hat y_i$ by $\delta b_1$ so $e_i = y_i - \hat y_i$ falls by $\delta b_1$. Each $\color{grey}{\text{grey}}$square has area $(\delta b_1)^2$. Each horizontally-aligned rectangles is as tall as the intercept was shifted, and as wide as the original residual for that point. The corresponding vertically-aligned rectangles are congruent, just rotated by a right angle. So I can line up all the $\color{blue}{\text{blue}}$ rectangles to form a single rectangle as high as the change in intercept and twice as wide as the sum of the positive residuals. The $\color{red}{\text{red}}$ rectangles form a single rectangle just as high, and twice the width of the sum of the (absolute values of the) negative residuals. 
Now suppose the original $\hat y = b_1 + b_2 x$ satisfied $\sum_{i=1}^n e_i=0$, which occurs (prove it!) if the line passes through the centroid $(\bar x, \bar y)$. Since the positive and negative residuals cancel out, the blue and red rectangles must too:
$$\color{Gold}{\text{RSS after change in intercept}} = \text{RSS of line through centroid} + \color{grey}{n(\delta b_1)^2}$$
Whatever the gradient of our original line, adjusting the intercept so it avoids the centroid will make the RSS worse, by the area of the grey squares. The least-squares line must therefore pass through the centroid and have  $\sum_{i=1}^n e_i=0$.
This does not tell us anything about which gradient minimises the RSS, but we can adapt our approach to consider a fixed intercept $b_1$ and a change in slope of $\delta b_2$. This time, the fitted values $\hat y_i$ rise by $x_i \delta b_2$, so fitted values rise (and residuals fall) by more the further $x_i$ lies to the right: the rectangles do not all have the same height, nor the grey squares the same area. But otherwise the dissection is much as before.

We find the area of a $\color{grey}{\text{grey}}$ square is $(x_i \delta b_2)^2$ and the area of a rectangle is proportional (by $\pm \delta b_2$) to $x_i e_i$. For the same reasons as before we want the $\color{blue}{\text{blue}}$ and $\color{red}{\text{red}}$ rectangles to balance out, which would require the positive and negative values of $x_i e_i$ to cancel, i.e. $\sum_{i=1}^n x_i e_i =0$. Totalling,
$$\color{Gold}{\text{RSS after change in slope}} = \text{RSS of line for which }\sum_{i=1}^n x_i e_i \text{ is zero} + \color{grey}{(\delta b_2)^2\sum_{i=1}^n x_i^2}$$
Regardless of the intercept, if we draw a line with a slope such that $\sum_{i=1}^n x_i e_i = 0$, then any changes to the slope will result in an RSS which is worse (higher) by the area of the grey squares. The least-squares line 
must therefore satisfy $\sum_{i=1}^n x_i e_i = 0$.

The diagrams oversimplify things: I didn't consider cases like a point with positive residual that develops a negative residual after the line sweeps above it, or where $x_i$ was negative. But we can verify the intuition algebraically:
$$\sum_{i=1}^n (e_i+\delta e_i)^2 = \sum_{i=1}^n e_i^2 + 2 \sum_{i=1}^n e_i \delta e_i + \sum_{i=1}^n (\delta e_i)^2 \tag{1}$$
The middle term represents the rectangular areas; note that the red and blue rectangles will take opposite signs since $e_i$ was positive for the blue case and negative for the red. Writing RSS as a function of $b_1$ and $b_2$,
$$RSS(b_1, b_2) = \sum_{i=1}^n e_i^2 = \sum_{i=1}^n (y_i - b_1 - b_2 x_i)^2 \tag{2}$$
Translating the regression line from $\hat y = b_1 + b_2 x$ to $\hat y = (b_1 + \delta b_1) + b_2 x$ reduces the residuals by the change in the intercept, $\delta e_i = -b_1$, so $(1)$ yields
$$
\begin{align}
RSS(b_1 + \delta b_1, b_2) &= RSS(b_1, b_2) + 2 \sum_{i=1}^n e_i (-\delta b_1) + \sum_{i=1}^n (-\delta b_1)^2 \\
RSS(b_1 + \delta b_1, b_2) &= RSS(b_1, b_2) - 2 \delta b_1 \sum_{i=1}^n e_i + n(\delta b_1)^2 \tag{3}
\end{align}
$$
Switching from $\hat y = b_1 + b_2 x$ to $\hat y = b_1 + (b_2 + \delta b_2) x$ gives $\delta e_i = -x_i \delta b_2$, and $(1)$ yields
$$
\begin{align}
RSS(b_1, b_2 + \delta b_2) &= RSS(b_1, b_2) + 2 \sum_{i=1}^n e_i (-x_i \delta b_2) + \sum_{i=1}^n (-x_i \delta b_2)^2 \\
RSS(b_1, b_2 + \delta b_2) &= RSS(b_1, b_2) - 2 \delta b_2 \sum_{i=1}^n x_i e_i + (\delta b_2)^2 \sum_{i=1}^n x_i^2 \tag{4}
\end{align}
$$
The argument can then proceed as before. Note that $(2)$ shows RSS is a quadratic function of both slope and intercept, so the expansions $(3)$ and $(4)$ could alternatively be derived by using the Taylor expansion
$$f(a+h) = f(a) + f'(a) h + \frac{1}{2} f''(a) h^2$$
with the partial derivatives
$$\frac{\partial RSS}{\partial b_1}= -2 \sum_{i=1}^n e_i, \quad \frac{\partial^2 RSS}{\partial b_1^2}=2n, \quad \frac{\partial RSS}{\partial b_2}=-2\sum_{i=1}^n x_i e_i, \quad \frac{\partial^2 RSS}{\partial b_2^2}= 2\sum_{i=1}^n x_i^2 $$
The $h$ represents our change $\delta b_1$ or $\delta b_2$. Note that setting the coefficient of the linear term in the change to zero (which we did by getting the red and blue rectangles to cancel out) is equivalent to putting $f'(a)=0$, i.e. ensuring the point we are expanding about satisfies the first order conditions to be a turning point. Checking that the coefficient of the quadratic term in the change is positive is equivalent to verifying $f''(a)>0$, which is the second order condition for this turning point to be a minimum: this is equivalent to our argument that, because the grey squares were being added on, the RSS was lower before the change.

Note that we were considering the residuals, $e_i = y_i - \hat y_i$, and not the errors, $\varepsilon_i = y_i - (\beta_1 + \beta_2 x_i)$ (where the betas are the "true" population regression parameters). Both residuals and errors are stochastic, in that if we re-sampled we would get a whole new bunch of random errors, a new OLS regression estimate, and a whole new set of residuals (we'd be fitting a different line to different points). The total of the residuals (as measured from the new OLS line) would still be zero. There's no restriction on the sum of the error terms, though. How could there be, remembering that the errors are generally assumed to be independent (or at least uncorrelated)? However, since the expected value of each error term is zero, the expected value of their sum is zero.
The arguments above should indicate that the sum of residuals would no longer be guaranteed to be zero if the intercept is removed from the model. Nor need the residuals sum to zero if the line was fitted by a method other than OLS. A vivid example of that is provided by Least Absolute Deviations: when moving an outlier up and down, you'll find (try it) that the LAD regression line remains "latched" to other points, and won't budge to take account of this change. The fact you can vary one residual while the others stays constant illustrates very dramatically that the sum of residuals is not invariant.     
A: One geometrical answer that might or might not resonate with you is that 

two legs of any right triangle are perpendicular.

It applies here because the "constant" or "intercept" in the regression results in projecting the response $y$ onto the vector $\hat{y}=\bar{y}(1,1,\ldots,1)$, which is one leg, and the residuals are the other leg $y-\hat{y} = (e_1, e_2, \ldots, e_n)$. 

This picture shows a response $y$ regressed against $x_1 = (1,1,\ldots,1)$. The coefficient is $\alpha = \bar{y}$ and the residual is $y_{\cdot 1}$.  It is geometrically obvious that $y_{\cdot 1}$ and $x$ are perpendicular.
Their perpendicularity, expressed in terms of the dot product, says the sum of the residuals is zero, because
$$\eqalign{
0 &= y_{\cdot 1} \cdot x_1 = (e_1, e_2, \ldots, e_n)\cdot (1, 1, \ldots, 1) = e_1(1) + e_2(1) + \cdots + e_n(1) \\ &= e_1+e_2+\cdots+e_n.}$$
Notice


*

*This is exact, not an expectation.  After all, it's just geometry!  We did not need to make any probabilistic assumptions about $y$.

*The result requires that the regression include a constant.

*The result is true for all multiple regressions that span a constant, because they can be carried out by first regressing $y$ against the constant and then performing a multiple regression of its residuals.

The image is taken from a longer document I have posted on how to control for variables.  It elaborates on point (3).
A: You must keep in mind the difference between the true values of the parameters and your estimates of them. Graphically, you can think of the situation as the true line and the estimated line.  
The usual way of returning the estimated line from the data just happens to correspond to the sum of residuals being zero. The residuals are deviations from the estimated line and the errors are deviations from the true line, the sum of errors doesn't necessarily equal zero but that is its expected value. 
So if you knew what the true line and you generated observations from that you would see that the expected value of the distribution of sum of errors would be zero. 
I.e. you would see the sum of errors distributed around zero if you ran multiple simulations. 
All of linear regression theory is predicated on the true line existing. 
