What is the meaning of the residual standard error in linear ordinary least squares output? In the linear model where $Y=X\beta+e$ and $e\sim N(0,\sigma^2 I_n)$, I understand that the expected value of $RSS$ is $\sigma^2 (n-p)$, so we can estimate $\sigma$ with $\sqrt(RSS/(n-p))$, and this is what R calls "residual standard error" in the summary of a linear model. 
My first question is what the phrase residual standard error actually refers to. I think this was asked before but the answer did not seem specific enough. I think residual standard error should refer to an estimate of the standard deviation of the estimates $\hat{e}_i$, but it doesn't seem obvious to me that this (sd($\hat{e_i}$)) would be $\sigma$. Or does it mean an estimate of the standard deviation of $e_i$ and is just misnamed. 
Also, even though I understand (on paper) why $RSS/(n-p)$ is an unbiased estimate of $\sigma^2$, why exactly does the formula for the sample standard deviation $\sum (\hat{e_i}-0)^2/(n-1)$ not apply? It seems like $\hat{e_i}$ is a "sample from a population with variance $\sigma^2$", but there are some conditions that prevent the usual formula from applying. 
Sorry if these have been asked before, but none of the posts I checked exactly answered my questions.
 A: The error is adjusted for the number of parameters in the model. This is just like standard deviation. For standard deviation one is adjusting for the fact that the mean value is an estimate by using $n-1$ in the denominator as opposed to $n$. If one has other parameters as well then one has to subtract $p$ from $n$ to adjust for that. Suppose that $p=n$, then there is no "play" in the solution, and no regression. What you would have in that case is an exact solution, not a regression per se. Then if you have $n=p+1$, there is only one degree of freedom, so all of the variability is only applied to one surfeit measurement, which has so little variability that it accounts the entire root RSS. As $n$ becomes progressively larger compared to $p$, then in general (for non-linear and linear modelling) we are adjusting by $n-p$.
BTW, residual sum squared is defined as $RSS = \sum_{i=1}^n (y_i - f(x_i))^2$ for any $f(x_i)$, whether linear or not.
A: We have the following linear model
$$\mathrm y = \mathrm X \beta + \mathrm e$$
where $\mathrm y \in \mathbb R^n$ and $\mathrm X \in \mathbb R^{n \times p}$ are given, $\beta \in \mathbb R^p$ is to be estimated, and $\mathrm e \sim N (0_n, \sigma^2 \mathrm I_n)$.
Left-multiplying both sides by $\mathrm X^{\top}$,
$$\mathrm X^{\top} \mathrm y = \mathrm X^{\top} \mathrm X \beta + \mathrm X^{\top} \mathrm e$$
If $\mathrm X$ has full column rank, then $\mathrm X^{\top} \mathrm X \in \mathbb R^{p \times p}$ is invertible. Left-multiplying both sides by $(\mathrm X^{\top} \mathrm X)^{-1}$, we obtain an unbiased estimate of $\beta$
$$\hat{\beta} := (\mathrm X^{\top} \mathrm X)^{-1} \mathrm X^{\top} \mathrm y = \beta + (\mathrm X^{\top} \mathrm X)^{-1} \mathrm X^{\top} \mathrm e$$
Hence,
$$\mathrm y - \mathrm X \hat{\beta} = \mathrm X \beta + \mathrm e - \mathrm X \beta - \mathrm X (\mathrm X^{\top} \mathrm X)^{-1} \mathrm X^{\top} \mathrm e = \left( \mathrm I_n - \mathrm X (\mathrm X^{\top} \mathrm X)^{-1} \mathrm X^{\top} \right) \mathrm e$$
where


*

*$\mathrm P_{\mathrm X} := \mathrm X (\mathrm X^{\top} \mathrm X)^{-1} \mathrm X^{\top}$ is the projection matrix that projects onto the column space of $\mathrm X$, which is $p$-dimensional (as we assumed that $\mathrm X$ has full column rank).

*$\mathrm I_n - \mathrm X (\mathrm X^{\top} \mathrm X)^{-1} \mathrm X^{\top} = \mathrm I_n - \mathrm P_{\mathrm X}$ is the projection matrix that projects onto the orthogonal complement of the column space of $\mathrm X$, which is the $(n-p)$-dimensional left null space of $\mathrm X$.


Thus,
$$\boxed{\mathrm y - \mathrm X \hat{\beta} = \left( \mathrm I_n - \mathrm P_{\mathrm X} \right) \mathrm e}$$
gives us the $(n-p)$-dimensional component of the noise that is orthogonal to the column space of $\mathrm X$. Taking the expected value of the $2$-norm of this projected component of the noise vector,
$$\begin{array}{rl} \mathbb E \left( \| \mathrm y - \mathrm X \hat{\beta} \|_2^2 \right) &= \mathbb E \left( \mbox{tr} \left( (\mathrm y - \mathrm X \hat{\beta}) (\mathrm y - \mathrm X \hat{\beta})^{\top} \right) \right)\\ &= \mathbb E \left( \mbox{tr} \left( \left( \mathrm I_n - \mathrm P_{\mathrm X} \right) \mathrm e \mathrm e^{\top} \left( \mathrm I_n - \mathrm P_{\mathrm X} \right) \right) \right)\\ &=  \mbox{tr} \left( \left( \mathrm I_n - \mathrm P_{\mathrm X} \right) \mathbb E (\mathrm e \mathrm e^{\top}) \left( \mathrm I_n - \mathrm P_{\mathrm X} \right) \right)\\ &=  \mbox{tr} \left( \left( \mathrm I_n - \mathrm P_{\mathrm X} \right) \sigma^2 \mathrm I_n \left( \mathrm I_n - \mathrm P_{\mathrm X} \right) \right)\\ &=  \sigma^2 \, \mbox{tr} \left( \left( \mathrm I_n - \mathrm P_{\mathrm X} \right)^2 \right)\\ &= \sigma^2 \, \mbox{tr} \left( \mathrm I_n - \mathrm P_{\mathrm X} \right)\\ &= \sigma^2 \, ( n - \mbox{tr} (\mathrm P_{\mathrm X}) )\\ &= \sigma^2 \, ( n - \mbox{rank} (\mathrm P_{\mathrm X}) )\\ &= ( n - p ) \, \sigma^2\end{array}$$
where we used the symmetry and idempotence of the projection matrix $\mathrm P_{\mathrm X}$. We also used the fact that the trace of a projection matrix is equal to its rank. 
Note that the expected value of the $2$-norm of the unprojected $n$-dimensional noise is
$$\mathbb E \left( \| \mathrm e \|_2^2 \right) = \mathbb E \left( \mbox{tr} \left( \mathrm e \mathrm e^{\top} \right) \right) = \mbox{tr} \left( \mathbb E \left( \mathrm e \mathrm e^{\top} \right) \right) = \mbox{tr} \left( \sigma^2 \mathrm I_n \right) = n \, \sigma^2$$
