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.

• How is the RSS defined? Oct 23, 2016 at 22:22
• @Rodrigo It would be defined as $(Y-X\hat{\beta})'(Y-X\hat{\beta})$
– user135912
Oct 24, 2016 at 0:59
• @Rocky Residual sum squared (RSS) is not limited to the linear case.
– Carl
Nov 10, 2016 at 16:36

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.

• And terminology-wise, does residual standard error refer to this "error associated with the linear model"?
– user135912
Oct 24, 2016 at 0:50
• So basically, more samples used to estimate other parameters means less samples contribute to estimating $\sigma$, and we need to divide by a smaller amount? And can you expand on this idea of "play"?
– user135912
Oct 24, 2016 at 0:51
• Play. Well not a typically statistical expression. When there are as many parameters as samples being fit, the solution is exact, there are no residuals, or if you wish to even define residuals, they are all zero. However, the solution may be complex valued, if a real solution doesn't exist. When it is complex, there are a transfinite number of them. However, there is always at least one exact solution, and if there are more of them, they are all exact, ergo, no residuals.
– Carl
Oct 24, 2016 at 2:33
• Terminology-wise, the models can be anything, exponentials, gamma variates, linear logarithms-whatever. They can be multiple linear regression or just one linear model, anything. Whatever the model is, count up the number of parameters, $p$. Subtract that from the number of samples, $n$.
– Carl
Oct 24, 2016 at 2:37
• When there is one more sample than the number of parameters, there is no exact solution, and, there will be non-zero residuals, but those residuals will be smaller in magnitude and closer to the samples than if there were either more samples or fewer parameters.
– Carl
Oct 24, 2016 at 2:40

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$$

• This only applies to linear models, which is an unnecessarily restrictive condition. It is nice to provide a theorem structure, but one should be careful that this is not just a long winded half-truth.
– Carl
Nov 10, 2016 at 16:29
• And isn't the question about linear models? Nov 10, 2016 at 17:25
• Yes, but the answer is more general than the question. For example, $RSS = \sum_{i=1}^n (y_i - f(x_i))^2$ in general. To answer otherwise is not as general.
– Carl
Nov 10, 2016 at 17:59
• @RodrigodeAzevendo More general contexts include non-linear OLS, multiple linear OLS regression, weighted OLS regression, Non-OLS regression; Theil-Sen, Deming, Tikhonov regularization and ridge regression.
– Carl
Nov 10, 2016 at 18:20
• The ordinary least squares in $y$ regression minimizes error in predicting $y$ given $x$ on the closed interval $\{Min(X),Max(X)\}$ of overdefined equations of the form $y=mx+b$, where $\{x_i,y_i\}$ are the data, $1\leq i\leq n$ and $n\geq3$.
– Carl
Nov 10, 2016 at 21:53