Why are the residuals in $\mathbb{R}^{n-p}$? I am currently studying Faraway's book "Practical Regression and Anova using R". In the very beginning, he says this:

and then

Now, I understand that $y\in\mathbb{R}^n$ since we would have $n$ obervations, then $\beta\in\mathbb{R}^p$ since we have $p$ variables and hence $p$ weights attached to those variables, but why are the residuals in $\mathbb{R}^{n-p}$? I don't understand, since by my interpretation, the residuals are $$\hat{\epsilon_i}=y_i-\hat{y_i}$$ for $i = 1,...,n$. Hence, they would be in $\mathbb{R}^n$? Am I missing something?
 A: $\newcommand{\rank}{\mathrm{rank}}$
$\newcommand{\tr}{\mathrm{tr}}$
$\newcommand{\real}{\mathbb{R}}$
$\newcommand{\eps}{\epsilon}$
Write the linear model in the matrix form
\begin{align*}
y = X\beta + \epsilon
\end{align*}
where $y \in \mathbb{R}^n, X = \begin{pmatrix} x_1 & \cdots & x_{p}\end{pmatrix} \in \mathbb{R}^{n \times p}, \beta \in \mathbb{R}^p, \epsilon \in \mathbb{R}^n$. By convention, $x_1 \equiv e$, where $e \in \real^n$ is a column vector of all ones, and $\rank(X) = p$.
Recall that $H = X(X'X)^{-1}X'$ is the "hat matrix", and the residual vector can be written as $\hat{\epsilon} = (I_{(n)} - H)y$. First note that $y$ lies in $\mathbb{R}^n$, therefore $\hat{\epsilon}$ lies in the image space of the matrix $I_{(n)} - H$ (you can view the matrix $I_{(n)} - H$ as a linear operator on $\mathbb{R}^n$), say $U$. Linear algebra theory asserts that $\dim(U) = \rank(I_{(n)} - H)$. On the other hand, since $I_{(n)} - H$ is idempotent (i.e., $(I_{(n)} - H)^2 = I_{(n)} - H$), its rank is equal to its trace, i.e.,
\begin{align*}
  & \rank(I_{(n)} - H) = \tr(I_{(n)} - H) = \tr(I_{(n)}) - \tr(H) \\
= & n - \tr(X(X'X)^{-1}X') = n - 
\tr(X'X(X'X)^{-1}) = n - \tr(I_{(p)}) = n - p,
\end{align*}
where we used properties $\tr(A + B) = \tr(A) + \tr(B)$ and $\tr(AB) = \tr(BA)$ of trace.
There is another more geometric flavor argument to derive $\dim(U)$, as your textbook tries to convey. Denote the space spanned by columns of $X$ by $W$, then it is easy to see that
\begin{align*}
\real^n = W \oplus W^\bot, \tag{$*$}
\end{align*}
where $W^\bot$ stands for the orthogonal complement of $W$. For any $y \in \real^n$, the decomposition
$$y = Hy + (I_{(n)} - H)y = \hat{y} + \hat{\eps}. $$
and $\hat{\eps}'Xa = y'(I - H)Xa = 0$ for any $a \in \real^p$ show that $\hat{\eps} \in W^\bot$, hence $U = W^\bot$. By $(*)$, it follows that
\begin{align*}
\dim(U) = \dim(W^\bot) = \dim(\real^n) - \dim(W) = n - p,
\end{align*}
since obviously $\dim(W) = \rank(X) = p$.
A: Orthogonal component
The space for the residuals is the orthogonal complement of the space spanned by $X$.
This is because the residual $\mathbf{y}-\mathbf{\hat{y}}$ is a vector that is perpendicular to the fit $\mathbf{\hat{y}}$. In the geometric representation it is a perpendicular projection of $y$ into the space spanned by the vectors in $X$ namely $\mathbf{\hat y} = X(X^TX)^{-1}X^T \mathbf{y}$.
This orthogonal component has dimension $n-p$

Intuitive illustration
The image in your book is abstract. It might maybe help to draw it for an actual example.
In the example below, you see an illustration for the fitting of $\mathbf{y} = a + b\mathbf{x}$ with only three points.

The error is a vector perpendicular to the surface spanned by $x_1$ and $x_2$. For any observation, the error will point in the same direction and can be seen to be a multiple of a line (a 1D space).

Linear algebra
You can describe any observation $\mathbf{y}$ (a sample of size $n$) in the space of potential observations as a sum of any $n$ orthogonal vectors. If $p$ of those form the vector $\mathbf{\hat{y}}$ then the remainder (the error $\mathbf{\epsilon}$) is a sum of the $n-p$ other ones.
$$\mathbf{y} = \overbrace{\underbrace{x_1 + x_2 + \dots + x_p}_{\substack{\mathbf{\hat{y}}\\\text{These $p$ vectors/regressors summed}\\\text{form the vector $\mathbf{\hat{y}}$}}} + \underbrace{ e_1 + e_2 + \dots + e_{n-p}}_{\substack{\mathbf{\epsilon}\\\text{The error $\mathbf{\epsilon}$ is a sum of}\\\text{the remaining $n-p$}}}}^{\text{$n$ orthogonal vectors of which $p$ are the regressors $x_i$}}$$
The $n-p$ stems from splitting the space into two orthogonal subspaces, one with dimension $p$, the other with dimension $n-p$.
A: I think the constraints on the residuals mean that they can be represented by $n-p$ numbers. Take a trivial case,
$$
y = \bar{y} + \varepsilon
$$
Here, there will be n error terms, but they will sum to zero. So, if you tell me $n-1$ error terms, I can work out the other one.
I think there is something analogous for the full model, but I'm not sure how to explain it. I guess, if you told me $n-p$ error terms, as well as the $X$ matrix and $\beta$ vector, I would be able to work out the rest.
