# In multiple linear regression does heteroskedasticity need to be evaluated for each predictor or just for the overall model?

If we want to make valid inferences about a simple linear regression estimator with one term, then homogeneity of variance of the residuals is a common assumption for making valid inferences.

However, I'm confused on what the assumption actually is in the multiple regression context. Do we assume homogeneity of variance for each predictor relative to the response in partial regression form or do we simply assume it for the entire model?

I'm not clear on how we would actually make valid predictor inferences just based on the entire model's residuals being normally distributed. When it comes to estimation uncertainty, I'd think we'd care about the specific regressor's unmeasurable error against the response.

• What inferences do you want to make?
– Dave
Oct 16, 2023 at 3:45
• @Dave On individual parameters in MLR ($H_0: \beta \neq 0$). So for example squared parameters or interactions. Oct 16, 2023 at 3:53
• Or just use robust standard errors, like many people do, and do not worry about the issue any further. Oct 16, 2023 at 10:13

First, notice that heteroskedasticity relates to the variances of the latent error terms, and not the computed residuals. Also, normality does not really play a role in the discussion. Further, heteroskedasticity or the lack thereof relates to if so, and how, the error variances relate to the regressors, as in $$E(u_i^2|x_i)=\sigma^2h(x_i)$$ for $$x_i=(x_{i1},\ldots,x_{iK})'$$ and for some nonnegative function $$h(x_i)\geq0$$. We would get homoeskedasticity if $$h(x_i)=c$$.

Now, I interpret your question as considering a case where $$h(x_i)=h(\tilde{x}_{i})$$ with $$\tilde{x}_{i}\subset x_i$$ in the sense that only the elements of $$\tilde x_i$$ "drive" the heteroskedasticity and then asking if inference on the coefficients on $$x_i\backslash \tilde x_i$$ is affected when those variables do not generate heteroskedasticity.

Here is an illustration that this indeed need not be the case (but is just an example without a claim to generality!).

It in particular shows that the conventional and robust standard errors (so the ones that are correct even in the presence of heteroskedasticity) on the "non-driving" X[,2] variable are the same, while they differ for the "driving" one X[,1], and that even though the two regressors are correlated so that the effect could maybe be expected to propagate:

library(mvtnorm)
library(lmtest)
library(sandwich)

n <- 20000

V1 <- 3
C12 <- 0.5
V2 <- 1

CovMat <- matrix(c(V1, C12, C12, V2), ncol=2)
X <- rmvnorm(n, sigma = CovMat)
e <- rnorm(n)
u <- X[,1]*e
y <- u

linreg <- lm(y~X-1)

> summary(linreg)\$coef
Estimate  Std. Error     t value  Pr(>|t|)
X1  0.0006061996 0.007442212  0.08145422 0.9350815
X2 -0.0174413564 0.012886659 -1.35344286 0.1759295

> coeftest(linreg, vcov = vcovHC(linreg, "HC0"))

t test of coefficients:

Estimate Std. Error t value Pr(>|t|)
X1  0.0006062  0.0125739  0.0482   0.9615
X2 -0.0174414  0.0132955 -1.3118   0.1896

> S <- matrix(c(3*V1^2, 3*V1*C12, 3*V1*C12, V1*V2+2*C12^2), ncol=2)

> # the plims
> sqrt(diag(solve(CovMat)%*%S%*%solve(CovMat))/n)
[1] 0.01243163 0.01279204

> sqrt(V1*diag(solve(CovMat))/n)
[1] 0.007385489 0.012792043


Somewhat more formally, we have that the asymptotic variance of OLS under heteroskedasticity is (a proof of this result may be found in, e.g., Hayashi, Econometrics, Proposition 2.1). $$Avar(\hat\beta)=\Sigma^{-1}S\Sigma^{-1},$$ with $$\Sigma=E(x_ix_i')$$ and $$S=E(u_i^2x_ix_i')$$. Note that, in our bivariate case and given that we generated heteroskedasticity via $$u_i=x_{i1}e_i$$, $$S=E\begin{pmatrix} e_i^2x_{i1}^4&e_i^2x_{i1}^3x_{i2}\\ e_i^2x_{i1}^3x_{i2}&e_i^2x_{i1}^2x_{i2}^2 \end{pmatrix}$$ By the law of iterated expectations, $$S=E\left[E\begin{pmatrix} e_i^2x_{i1}^4&e_i^2x_{i1}^3x_{i2}\\ e_i^2x_{i1}^3x_{i2}&e_i^2x_{i1}^2x_{i2}^2 \end{pmatrix}\biggl|x_i\right]=E(e_i^2)\begin{pmatrix} E(x_{i1}^4)&E(x_{i1}^3x_{i2})\\ E(x_{i1}^3x_{i2})&E(x_{i1}^2x_{i2}^2) \end{pmatrix}$$ From Isserlis' Theorem (see also here or here), letting $$c=E(x_{i1}x_{i2})$$ and $$v_j=V(x_{ij})$$ (notice our regressors are multivariate normal, have mean zero for simplicity; we also took $$E(e_i^2)=1$$), so that $$S=\begin{pmatrix} 3v_1&3v_1c\\ 3v_1c&v_1v_2+2c^2 \end{pmatrix}$$ and hence $$Avar(\hat\beta)=\frac{1}{v_1 v_2-c^2}\begin{pmatrix}3 v_1 v_2-2 c^2 & -c v_1\\ -c v_1 & v_1^2\end{pmatrix}$$

This is what will be consistently estimated by robust standard errors (more precisely, after squaring these and multiplying by $$n$$ to undo the effect of taking the square root and dividing by $$n$$ to get standard errors).

The standard asymptotic variance estimate in turn will tend to $$Var(u_i)\Sigma^{-1}=E(e_i^2)E(x_{i1}^2)\Sigma^{-1}=\frac{v_1}{v_1 v_2 - c^2}\begin{pmatrix}v_2 & -c\\ -c & v_1\end{pmatrix}$$ Hence, in the present example, the nonrobust and robust variance estimates will indeed agree for the regressor that does not "produce" heteroskedasticity (the $$(2,2)$$-elements), while they do not for the "producing" one (the $$(1,1)$$-elements).

The example however relies on, e.g., specific features of the multivariate normal. I for example replaced the above regressors with multivariate t ones as in X <- rmvt(n, sigma = CovMat, df=3) and got different standard errors for both coefficients.