# Are the errors in this formulation of the simple linear regression model random variables?

On page 21 of Applied Linear Regression, fourth edition, by Sanford Weisberg, the error $$e_i$$ for case $$i$$ under the simple linear regression model is defined to be $$y_i - E(Y | X = x_i)$$, where $$E(Y | X = x_i)$$ is assumed to equal $$\beta_0 + \beta_1 x_i$$ for some unknown $$\beta_0, \beta_1 \in \mathbb{R}$$. The book says that

The errors $$e_i$$ depend on unknown parameters in the mean function and so are not observable quantities. They are random variables and correspond to the vertical distance between the point $$y_i$$ and the mean function $$E(Y | X = x_i)$$.

It doesn't seem to me like $$e_i$$ is a random variable, because it's a function of $$y_i$$ and $$x_i$$, which are non-random, observed values. Why can $$e_i$$ be considered a random variable?

I looked up your citation (4th edition, page 21) because I found it very alarming and was relieved to find is actually given as:

$$\hat{e}_i = y_i − \widehat{E}(Y|X=x_i) = y_i - (\hat{\beta}_0 + \hat{\beta}_1) \tag{2.3}$$

Which is still confusing, I grant you, and the difference isn't actually germane to your question, but at least it isn't patently false. I'll explain why I found it alarming before discussing your (unrelated, I think) question. The "hat" indicates "estimated", usually by MLE in the context of linear regression, and there is a crucial distinction between "true errors" which are denoted $$\epsilon_i$$ and are normally distributed and i.i.d., and "residuals which are denoted $$e_i$$ and are not i.i.d. The formula without the hats would imply the two are exactly equal which is not the case.

On to your real question, which boils down to, "are the given data $$x_i$$ and $$y_i$$ random or not?"

If you believe the pairs $$(x_i, y_i)$$ are known and not-random, e.g. that is, if you believe that $$\forall\; 1 \leq i \leq n,\, (x_i, y_i) \in \mathbb{R} \times \mathbb{R}$$, then the residuals $$e_i$$ are also known and non-random, e.g. $$\forall\; 1 \leq i \leq n,\, e_i \in \mathbb{R}$$. This is because there is a deterministic function for the "best" parameters $$\hat{\beta_0}$$ and $$\hat{\beta_1}$$ from those observations, and then a deterministic function for the residuals in terms of those parameters. This point of view is useful and allows us to derive the MLE estimators of $$\beta$$, for example. It is also the most intuitive view to take when your sitting in front of a concrete, real-world dataset.

However, it kind of puts the cart before the horse and basically shuts down certain kinds of statistical analysis. For example, we cannot talk about the "distribution" of $$\hat{\beta}_1$$ because it is not a random variable and therefore has no distribution! How can we then talk something like the Wald test? Likewise, how do we talk about the "distribution" of residuals so that we can say whether one is an outlier or not?

The way this is done is treating the dataset itself as random. When we want to do statistical inference on a known dataset, we can then treat the known values as a realization of the random dataset. The exact construction is a little bit pedantic but and is often omitted but it helps to go through it at least once. First, we say that $$X$$ and $$Y$$ are two random variables with some joint probability distribution $$F_{X,Y}(\mathbf{\beta}, \sigma^2)$$ with parameters $$\mathbf{\beta} = [\beta_0, \beta_1]^T$$ and $$\sigma$$. $$F_{X,Y}$$ is specified by the model $$Y = X\beta_0 + \beta_1 + \epsilon, \epsilon \sim \mathcal(0, \sigma^2)$$. Now, imagine that we have $$n$$ i.i.d. copies of $$F_{X,Y}$$ that we combine into one big joint probability function $$F_{X_1,Y_1,X_2,Y_2,...,X_n,Y_n}$$.

Now we can imagine the dataset $$(x_i, y_i)$$ for $$i=1,...,n$$ not merely as some known set of numbers, but as a realization sampled from $$F_{X_1,Y_1,X_2,Y_2,...,X_n,Y_n}$$. Each time we sample, we don't just get one pair of numbers, we get $$n$$ pairs of numbers: a brand new dataset. But that means the parameters $$\hat{\beta}$$ get new estimates, and we then calculate new residuals $$e_i$$, right?

Instead of thinking of this as repeated sampling, which is somewhat crude, we can express this entirely in the algebra of random variables. It can be expressed as two $$n$$-dimensional random vectors $$\vec{X}$$ and $$\vec{Y}$$ drawn from $$F_{X_1,Y_1,X_2,Y_2,...,X_n,Y_n}$$. Now $$\hat{\beta}_0$$ and $$\hat{\beta}_1$$ are random variables because they are functions of $$(\vec{X}, \vec{Y})$$. Likewise, all the $$e_i$$ are random variables because they are functions of $$(\vec{X}, \vec{Y})$$.

This state of affairs is much better, because now we can make statements like "The set of residuals $$e_i$$ cannot be independent because they always sum exactly to zero" or "the standard error of $$\hat{\beta}_1$$ follows a t-distribution." without talking literal nonsense. (Both of these statements only make sense if their subjects are random variables.)

In the real world we can't always go and get a brand-new, randomly sampled dataset. We can approximate this with something like the bootstrap, of course, but doing it for real isn't usually practical. But doing it conceptually allows us to think clearly about how randomness during sampling would affect our regression.

You'll note that I did not introduce new notation for $$e_i$$ and $$\hat{\beta}$$ but simply said, "now these things, which we previously thought of a concrete realizations, will now be treated as random variables." As far as I can tell, you just have to be on your toes for this kind of signposting - the same kind you found in your textbook - to indicate whether symbols are referring to random or non-random variables because while there are conventions (such as using uppercase roman letters for random variables) they are not consistently applied. If the author tells you $$e_i$$ is a random variable, he is telling you is also viewing $$x_i$$ and $$y_i$$ as random variables.

In simple linear regression, we assume that the observations are randomly perturbed from the conditional expected value, i.e. $$E[Y|X=x_i]$$; so, each of your observations are assumed to be generated from a model of the form: $$Y=\beta_0+\beta_1X+\epsilon \ \ , \epsilon\sim N(0,\sigma^2)$$

This makes each $$\epsilon_i$$ a RV by definition. Think about a box where you give $$x_i$$ and get $$y_i$$, and you never know what's inside, how much error is introduced by the box etc. Even if we really know that the relation is of the form given above, we don't know the true $$\beta_0,\beta_1$$. If we had known those quantities, we would easily recover $$\epsilon_i$$. Instead, we estimate those, and get residuals.