How to modify a linear regression model when the variance is not constant? I have a very basic question. I have a model that has residuals that don't follow a constant variance pattern, but rather it has a clustered pattern (see the figure attached).
Any idea of how to deal with this challenge? Should I find what those clusters are and add them as variables to my model? Any idea/link/video is highly appreciated.
Thanks

 A: Variance is the spread of the residuals: their tendency to depart from zero, as measured by their typical (root mean square) distance from zero.  The left hand plot uses height to represent the residuals, but there is no discernible variation the amounts by which heights typically deviate from the zero value (near the middle) as one scans across the plot.  This is homoscedasticity, or constant variance, in the response.
I simulated a dataset with essentially the same plots, and then modified the simulation by making the spread in the right hand cluster twice that of the spread in the left hand cluster.

The change in vertical spreads between the clusters at the left is now apparent, even though it has produced only a subtle tinge of non-normality in the Normal QQ plot at right.
A decent way to assess the possibility of non-homoscedasticity ("heteroscedasticity") is to plot the square roots of the absolute values of residuals against the predicted values, as in this next plot of the same data.  We are now looking for a change in average height rather than a change in spread:

On this plot I have superimposed some least-squares fits, one for each (obvious) group, in order to visualize how average height (root residual size) varies with the predicted value. The averages are indeed at (barely noticeably) different levels (2 at the left and 2.5 at the right).  The apparent trend in group 0 at the left is purely accidental: no trend was built into the simulation, apart from the difference in spreads between the clusters.  However, when you do observe such a trend as shown in the left hand cluster, you do begin to suspect some (additional) heteroscedasticity may be apparent.  Experience with simulations like this one can help you decide whether to react to such evidence.
Suppose you decided to identify these clusters and treat the identifier as a new explanatory variable?  It won't help, because the least squares procedure still gets a good answer.  Here is a summary of the enhanced model.
lm(formula = y ~ x + factor(Group), data = X)

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    620.58150    0.78085 794.746   <2e-16 ***
x                1.04758    0.05058  20.712   <2e-16 ***
factor(Group)1  -5.45333    5.14010  -1.061     0.29 

The relatively large p-value of 0.29 shows the cluster factor is not significant.  The diagnostic plots for this model look so much like the original that I won't take up the space on your screen to show them.

This is the R code used to create the simulated data and visualize them.
#
# Create data.
#
g <- rep(0:1, each=100)
set.seed(17)
x <- rnorm(length(g), 100 * g, 10)
y <- rnorm(length(x), x, 10 * (1+g)/2) + 620
X <- data.frame(Group=factor(g), x=x, y=y)
#
# Model them.
#
fit <- lm(y ~ x, X)
#
# Plot the fit.
#
par(mfrow=c(1,2))
plot(fit, which=c(1,2))
par(mfrow=c(1,1))
#
# Study the residuals more closely.
#
library(ggplot2)
Y <- data.frame(factor(g), predict(fit), sqrt(abs(residuals(fit))))
names(Y) <- c("Group", "Predicted", "Root residual")
ggplot(Y, aes(Predicted, `Root residual`, color=Group)) + 
  geom_point() + 
  geom_smooth(method="lm")
#
# Plot an alternative fit.
#
fit.g <- lm(y ~ x + factor(Group), X)
summary(fit.g)
par(mfrow=c(1,2))
plot(fit.g, which=c(1,2))
par(mfrow=c(1,1))

A: *

*The confidence intervals generated by OLS can be biased downwards in the heteroskedasticity case. So, we need to use robust standard errors if we are concerned about the t-stats of coefficient estimates. More information can be found here.

*The other problem is efficiency (ie. do we use the information embedded in our data well?). An observation with higher variance means that it has lower information. We want to give more weight to observations with lower variance and less weight to observations with higher variance. This can be done using a generalized least squares. More information can be found here.
