Orthogonality of residuals in linear regression In multiple linear regression, I came across the statement that both $e$(residual) and predicted $y$ are projections of actual y and $e$ is orthogonal to predicted $y$.
I was trying to visualize the same for a simple linear regression. Doesn't this mean that the residual term should be perpendicular to the fitted line?
 A: "$e$  is orthogonal to predicted $y$" is the orthogonal in the $n$ dimensional space, where $n$ is number of the observation. It is hard to visualize because of the high dimension. But it follows the general definition of orthogonal. $e'\hat y = 0$.   http://mathworld.wolfram.com/Orthogonal.html
A: Assumptions
The assumption underlying linear regression is that there is a true
underlying $x$ and $y$ relationship that allows you to predict $\hat{y}$, knowing $x$. In the absence of unmodelled variation we have the following relation:
$$\hat{y} = m x+c$$
However error can exist in the actual observations so the equation is rewritten in various forms depending on
how you describe the errors. If we are talking about residuals then we
specifically mean the variation that remains after the model described
variance is taken into account. So $\epsilon$ means that:
$$y_{obs} = mx + c + \epsilon$$
this means that the orthogonality is expected to exist between predicted
$\hat{y}$ and $\epsilon$. The nice thing about simulation is we can create ground
truth then perturb it in a controlled way then probe what happens
This formulation does not seek to untangle the variance arising due to
errors (or to be more precise, unmodelled variation) in $x$ or in $y$.
It is also important to distinguish what we observe from what the true underlying data generating process is:
$$ x_{obs} = x_{true} +\epsilon_x$$
$$ y_{obs} = y_{true} +\epsilon_y$$
Where $\epsilon_x \perp\!\!\!\perp \epsilon_y$ (if unaccounted variance is not independent between $x_{obs}$ and $y_{obs}$ then $y_{obs}$ will be a transformation of $x_{obs}$ rather than $x_{true}$)
Note:
since $mx+c=\hat{y}$ the difference between these two formulas gives
$$\epsilon = y_{obs}-\hat{y}$$
Methods
The following Matlab code was used to create simulated data which started with a ground truth to which was added variation not explained by the true model. This simulated observed data had variations ('errors') in both $x$ and $y$, which are propagated through the model building and prediction.
I.e. we model the situation where $x_{obs}= x_{true}+\epsilon_x$ and $y_{obs}= y_{true}+\epsilon_y$
with $y$ a linear transformation of $x$ and $\epsilon_x \perp\!\!\!\perp \epsilon_y$ allowing us to disentangle the contributions of $\epsilon_x$ and $\epsilon_y$ in a way that can't be done with real data:
    % create the 'ground truth' - no sampling, technical nor measurement errors
modelTrue = [1.5, 0.5];
xTrue = 1:100;
yTrue = polyval(modelTrue , xTrue);


      % create variation not explained by the underlying linear relationship
xUnex = 0.1*(rand(1,100)-0.5);% 10% range of error, uniformly distributed
yUnex = 3*randn(1,100);% random noise with standard deviation of 3
      % Note that many other perturbations could be considered - varing the
      % distribution of noise, using a different transformation model etc.

      % create observed data
xobs = xTrue + xUnex ; %create the data that will actually be observed
yobs = yTrue + yUnex; %create y data that will actually be observed

modelobs =  polyfit(xobs,yobs,1);
ypred = polyval(modelobs,xobs);
EP = yobs -  ypred;
EPvsTrue = yTrue -  ypred;

Results
Here we see scatter plots of the different data sets (ground truth, observed and modelled). Note there is a noticeable change in slope (which will vary depending on randomisation) in the predicted data compared with the ground truth.

Next we look at the correlation between the $\epsilon$ and various parameters.
      Xobs         XTrue       Yobs         Ypred        YTrue 
    ___________    _______    ________    ___________    _______

    -6.3787e-15    0.31728    0.065877    -6.3897e-15    0.31728

There is no correlation between $\epsilon$ and observed $x_{obs}$ nor predicted $\hat{y}$. There is a weak correlation between $\epsilon$ and the $y_{obs}$, which is because $y_{obs}$ is a summation of $\hat{y}+\epsilon$, so the strength of the correlation will depend on the relative magnitudes of $\hat{y}$ and $\epsilon$. There is a stronger correlation between $\epsilon$ and the true values, due to the deviation of the slope in figure 1 from the ideal.
However, you are not interested in such a proof, this is merely a preamble to answering your question about visualisation.
Visualisation of what is orthogonal

Doesn't this mean that the residual term should be perpendicular to
the fitted line?

The 'fitted line' in linear regression typically can be one of two things depending on how you plot it ($x$ vs $y$ with a best fit, or $y_{obs}$ vs $\hat{y}$), but for neither is the orthogonality between the best fit line and the scatter. Linear regression does not attempt to unpackage $\epsilon_x$ and $\epsilon_y$ contributions to unmodeled variation, which would be needed to determine what angle the outcome error is to the true result. In predicted vs observed we have $y_{obs}$ which we saw above contains an explicit contribution from $\epsilon$, so cannot be orthogonal.
Instead we need to look specifically at the $\hat{y}$ and compare it with $y_{obs}$, as shown in the figure below. In addition to the classic $\hat{y}$ vs $y_{obs}$ plot I have included plots of $\hat{y}$ vs $\epsilon$ and  $\epsilon$ vs $y_{obs}$ along with their best fits in each direction:

You can see that for the $y_{obs}$ there is a slope evident in the trend between $y_{obs}$ vs $\epsilon$ (centre of x axis, blue line with slope of 1). Given our discussion above we would expect to see an interaction between these two. In contrast we see that there is no trend linking $\hat{y}$ and $\epsilon$, the two axis are parallel to the axes of the main plot and perpendicular to each other. this is the visual representation of the orthogonality of $\hat{y}$ and $\epsilon$.
Note on extension to higher dimensions
This same principle applies to higher dimension regression models as they all reduce the data to a $\hat{y}$ value. This means the prediction can be composed of multiple orthogonal components, but all $\epsilon$ is concerned with is what is left after your model is applied - so it doesn't matter how complicated your model was.
