What does having "constant variance" in a linear regression model mean?

Question

What does having "constant variance" in the error term mean? As I see it, we have a data with one dependent variable and one independent variable. Constant variance is one of the assumptions of linear regression. I am wondering what homoscedasticity means. Since even if I have 500 rows, I would have a single variance value which is obviously constant. With what variable should I compare the variance?

Answer 1

It means that when you plot the individual error against the predicted value, the variance of the error predicted value should be constant. See the red arrows in the picture below, the length of the red lines (a proxy of its variance) are the same.

Answer 2

This is a place where I've found looking at some formulas helps, even for people with some math anxiety (I'm not suggesting that you do, necessarily). The simple linear regression model is this:
$$ Y=\beta_0+\beta_1X+\varepsilon \\ \text{where } \varepsilon\sim\mathcal N(0, \sigma^2_\varepsilon) $$ What's important to note here is that this model explicitly states once you've estimated the meaningful information in the data (that's the "$\beta_0+\beta_1X$") there is nothing left over but white noise. Moreover, the errors are distributed as a Normal with a variance of $\sigma^2_\varepsilon$.

It's important to realize that $\sigma^2_\varepsilon$ is not a variable (although in junior high school level algebra, we would call it that). It doesn't vary. $X$ varies. $Y$ varies. The error term, $\varepsilon$, varies randomly; that is, it is a random variable. However, the parameters ($\beta_0,~\beta_1,~\sigma^2_\varepsilon)$ are placeholders for values we don't know--they don't vary. Instead, they are unknown constants. The upshot of this fact for this discussion is that no matter what $X$ is (i.e., what value is plugged in there), $\sigma^2_\varepsilon$ remains the same. In other words, the variance of the errors / residuals is constant. For the sake of contrast (and perhaps greater clarity), consider this model:
$$ Y=\beta_0+\beta_1X+\varepsilon \\ \text{where } \varepsilon\sim\mathcal N(0, f(X)) \\ ~ \\ \text{where } f(X)=\exp(\gamma_0+\gamma_1 X) \\ \text{and }\gamma_1\ne 0 $$ In this case, we plug in a value for $X$ (starting on the third line), pass it through the function $f(X)$ and get the error variance that obtains at that exact value of $X$. Then we move through the rest of the equation as usual.

---

The above discussion should help with understanding the nature of the assumption; the question also asks about how to assess the it. There are basically two approaches: formal hypothesis tests and examining plots. Tests for heteroscedasticity can be used if you have experimental-ish data (i.e., that only occur at fixed values of $X$) or an ANOVA. I discuss some such tests here: Why Levene test of equality of variances rather than F-ratio. However, I tend to think looking at plots is best. @Penquin_Knight has done a good job of showing what constant variance looks like by plotting the residuals of a model where homoscedasticity obtains against the fitted values. Heteroscedasticity can also possibly be detected in a plot of the raw data, or in a scale-location (also called spread-level) plot. R conveniently plots the latter for you with a call to plot(linear.model, which=3); it is the square root of the absolute values of the residuals against the fitted values, with a lowess curve helpfully overlaid. You want the lowess fit to be flat, not sloped.

Consider the plots below, which compare how homoscedastic vs. heteroscedastic data might look in these three different types of figures. Note the funnel shape for the upper two heteroscedastic plots, and the upward sloping lowess line in the last one.

For completeness, here is the code that I used to generate these data:

set.seed(5)

N  = 500
b0 = 3
b1 = 0.4

s2 = 5
g1 = 1.5
g2 = 0.015

x        = runif(N, min=0, max=100)
y_homo   = b0 + b1*x + rnorm(N, mean=0, sd=sqrt(s2            ))
y_hetero = b0 + b1*x + rnorm(N, mean=0, sd=sqrt(exp(g1 + g2*x)))

mod.homo   = lm(y_homo~x)
mod.hetero = lm(y_hetero~x)