# Why does heteroskedasticity not affect $R^2$ and why does it make estimated regression more statistically significant?

I am studying what the consequences of heteroskedasticity are. And i found that assuming that the model is linear in the parameters (i.e $$Y=X\beta+\epsilon$$), is identifiable, has no perfect collinearity (i.e $$X$$ is full column rank), is exogenous (i.e. $$E[\epsilon|X]=0$$) and has no serial correlations (i.e. $$Var[\epsilon|X]=\Sigma$$, where $$\Sigma$$ is diagonal) we have :

1. $$E[\hat{\beta}] = \beta$$, because $$\hat{\beta} = \beta + (X^TX)^{-1}X^T\epsilon$$ and $$E[\epsilon|X]=0$$

2. OLS is no longer BLUE because the homoskedasticity assumption is absent (and hence Gauss-Markov theorem cannot be applied).

3. Estimated variances of regression coefficients will be biased (because $$E\bigg[\hat{Var}[\hat{\beta}_i|X]\bigg|X\bigg] = (X^TX)^{-1}_{ii}E\bigg[MS_{res}\bigg|X\bigg]\neq (X^TX)^{-1}_{ii}\Sigma_{ii}=Var[\hat{\beta}_i|X]$$, where $$MS_{res} = \dfrac{e^Te}{n-p}$$, where $$e = Y-\hat{Y}=$$ the residuals)

4. The corresponding t-statistics of the estimated regression coefficients will appear to be more significant than they really are (page 4). This implies that $$MS_{res}<\Sigma_{ii}$$. I do not understand why this should be the case.

5. Since $$R^2=1-\dfrac{\epsilon^T(I-H)\epsilon}{Y^T(I-J_n)Y}$$ is based on overall sums of squares, it is unaffected by heteroskedasticity (page 5). Again I do not understand why this should be the case. Specifically, in what sense is $$R^2$$ unaffected? Like if I somehow made a different linear model to deal with the heteroskedasticity, then would this new linear model have the same $$R^2$$? If so how?

Could you please verify points 1,2 and 3, and give non-heuristic explanations for points 4 and 5?

## Edit1:

@ChristophHanck,
I dont understand why $$R^2$$ would converge at all under heteroskedasticity. I mean, under homoskedasticity, assuming $$e_i\stackrel{iid}{\sim}(\text{mean}=0,\text{variance}= \sigma^2)$$ and $$n>>p$$, we have $$\dfrac{e^Te}{n-p}\approx\dfrac{e^Te}{n}=\dfrac{e_1^2+\cdots+e_n^2}{n}\rightarrow E[\epsilon_i^2]=\sigma^2$$, by weak law of large numbers. But under heteroskedasticity, assuming $$e_i\stackrel{independent}{\sim}(\text{mean}=0,\text{variance}= \sigma^2_i)$$ and $$n>>p$$, why should $$\dfrac{e_1^2+\cdots+e_n^2}{n}$$ be expected to converge at all?
And in terms of this answer to your cited question, why should $$\dfrac{\sum_{i=1}^N{m_{ii}\sigma^2_i}}{N - K}$$ be expected to converge at all?
also, under heterskedasticity, the $$y_i$$'s are not iid. So, why should $$\sum(y_i-\bar{y})^2$$ converge at all?

• Do you think the definition of $R^2$ is correct as you have specified it? I don't think so (and you might want to confirm this) because you have used error terms rather than the residuals in defining it. The OLS estimate doesn't make any assumption about the error distribution, and thus, should be agnostic to distributional assumptions about errors. Thus, residuals are also agnostic to assumptions about errors as long as you are using OLS estimates, and thus, $R^2$ shouldn't be affected by heteroscedasticity in residuals. Sep 22 at 9:49
• @medium-dimensional $e=Y-\hat{Y}=(I-H)Y=(I-H)\epsilon \implies e^Te=\epsilon^T(I-H)\epsilon$ Sep 22 at 10:18
• could some one explain the 4th point also? Sep 22 at 10:22
• I understand the point as follows: if your errors have different variances, so will your Y's, and since R^2 considers, as your formula shows, the ratio of the two, the effect cancels out in R^2 Sep 22 at 11:04

As to 5, $$R^2$$ will tend to $$1-\frac{Var(u)}{Var(y)}$$ under homoskedasticity. Under heteroskedasticity (no detailed derivation here), the numerator will be related to some function of the average variance (by suitably modifying the derivation here: How do I use the standard regression assumptions to prove that $\hat{\sigma}^2$ is an unbiased estimator of $\sigma^2$?)

In a linear model, we have, under suitable assumptions, $$Var(y)=\beta^2Var(x)+Var(u),$$ so that $$R^2\to_p1-\frac{Var(u)}{\beta^2Var(x)+Var(u)}$$ So provided the average variance under heteroskedasticity is similar to the variance under homosekdasticity, we may expect similar results for $$R^2$$.

Illustration:

n <- 10000

u.hom <- rnorm(n)
u.het <- rnorm(n, sd=runif(n, 0.5, 1.5)) #  i.e., on average, the s.d. will also be one
x <- rnorm(n)

y.hom <- 2*x + u.hom
y.het <- 2*x + u.het

> summary(lm(y.hom~x))

Call:
lm(formula = y.hom ~ x)

Residuals:
Min      1Q  Median      3Q     Max
-3.5527 -0.6754 -0.0023  0.6854  3.8040

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.008325   0.010040  -0.829    0.407
x            2.009258   0.009825 204.494   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.004 on 9998 degrees of freedom
Multiple R-squared:  0.807, Adjusted R-squared:  0.807
F-statistic: 4.182e+04 on 1 and 9998 DF,  p-value: < 2.2e-16

> summary(lm(y.het~x))

Call:
lm(formula = y.het ~ x)

Residuals:
Min      1Q  Median      3Q     Max
-5.0801 -0.6242 -0.0082  0.6446  4.9125

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.003788   0.010510    0.36    0.719
x           2.011002   0.010286  195.51   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.051 on 9998 degrees of freedom
Multiple R-squared:  0.7927,    Adjusted R-squared:  0.7927
F-statistic: 3.823e+04 on 1 and 9998 DF,  p-value: < 2.2e-16

• please refer to edit 1 of my post. Sep 22 at 13:37
• also under heterskedasticity, the $y_i$'s are not iid. So, why should $\sum(y_i-\bar{y})^2$ converge at all? Sep 22 at 13:41
• As I said, my argument surely is heuristic, as I did not have time to think about the issue rigorously. I expect convergence to some average variance. The sample code appears to illustrate that there seems to be something to the heuristic. Sep 22 at 15:20
• Yes, of course. Also, your answer provides some much needed intuition. +1 for that. Could you also give some explanation for the 4th point. Thanks. Sep 22 at 16:16
• Please see my separate answer for four, I did not want to clutter this one even more Sep 23 at 14:20

As to 4, that is indeed, while empirically often the case, not necessarily true.

Consider as a tractable example the expressions for the standard and robust variance estimator (the correct one under heteroskedasticity) for a regression on a constant and a dummy investigated in these answers: Eicker-Huber-White Robust Variance Estimator and How to prove equality of standard errors for two-sample t-test and linear regression?

In particular, for $$n_1$$ "treated" observations for which the dummy equals one and the rest of $$n$$ being the other $$n_0$$ observations, we have for the robust variance estimator of the slope coefficient that $$\widehat{\mathbb{V}}_{\text{EHW}}=\frac{n_1^2s_0^2+n_0^2s_1^2}{n_0^2n_1^2},$$ with $$s_t^2$$, $$t=0,1$$, the sum of squared residuals belonging to the two groups of observations, and $$s^2=s_0^2+s_1^2$$.

Also (here, we work with the ML estimator for the error variance for simplicity which divides the sum of squared residuals with $$n$$ instead of $$n-2$$ to ease the algebra; that does not matter for large $$n$$; also note the definition of $$s^2$$ in the second link is different from the one employed here, just scale by $$n$$), the standard OLS slope variance estimator is $$\widehat{\mathbb{V}}_{\text{standard}}=\frac{s_0^2+s_1^2}{n_0n_1},$$

We have $$\widehat{\mathbb{V}}_{\text{EHW}}<\widehat{\mathbb{V}}_{\text{standard}}$$ when (unless $$n_0= n_1$$, in which case we always have $$\widehat{\mathbb{V}}_{\text{EHW}}=\widehat{\mathbb{V}}_{\text{standard}}$$) $$s_1^2>\frac{n_1}{n_0}s_0^2$$ or $$\frac{s_1^2}{n_1}>\frac{s_0^2}{n_0},$$ i.e., when the variance of the "treated" observations is higher than that of the untreated.

Take, e.g., $$n_0=n/4$$ and hence $$n_1=3n/4$$. Then, $$\widehat{\mathbb{V}}_{\text{EHW}}<\widehat{\mathbb{V}}_{\text{standard}}$$ will result when $$s_1^2>3s_0^2,$$ i.e., when the variance associated with the dummy's observations equal to 1 is sufficiently larger that the others, unless there are relatively too many such observations (i.e., when $$n_1/n_0$$ is large). Illustration:

library(sandwich)
library(lmtest)

n <- 2000
n0 <- 1/4*n
n1 <- 3/4*n

y <- c(rnorm(n0, sd=1), rnorm(n1, sd=1.05))  # i.e. more variance when dummy=1
x1 <- c(rep(0, n0), rep(1,n1)) # a dummy regressor

reg <- lm(y~x1)

> summary(reg)

Call:
lm(formula = y ~ x1)

Residuals:
Min      1Q  Median      3Q     Max
-4.3243 -0.6771 -0.0006  0.6766  3.0950

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03402    0.04677  -0.727    0.467
x1           0.02879    0.05400   0.533    0.594

Residual standard error: 1.046 on 1998 degrees of freedom
Multiple R-squared:  0.0001422, Adjusted R-squared:  -0.0003582
F-statistic: 0.2842 on 1 and 1998 DF,  p-value: 0.594

> robust.se <- coeftest(reg, vcov = vcovHC(reg, "HC0"))

> robust.se

t test of coefficients:

Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.034021   0.045556 -0.7468   0.4553
x1           0.028788   0.053065  0.5425   0.5875

> residuals <- resid(reg)

> 3*sum(residuals[1:n0]^2) # n_1/n_0*s_0^2
 1556.493

> sum(residuals[(n0+1):n]^2) # s_1^2
 1666.235


Interesting question! @christoph hanck has treated a lot of stuff here. In your edit, you ask a question: why does the average of $$e_i^2$$, with iid $$e_i\sim(0,\sigma_i^2)$$ converge? Here I provide an answer to this question. First we have \begin{align} \frac 1n \sum_i e_i^2 &= \frac 1n \sum_i (e_i^2-\sigma_i^2) + \frac 1n \sum_i \sigma_i^2 \end{align} Let's look at the first term. Given the indepence assumed on $$e_i$$, $$z_i=e_i^2-\sigma_i^2$$ is also indepndent and has a zero mean. Therefore as long as you have a moment of $$2+\delta$$ for $$e_i$$, the result (that $$\frac 1n \sum_i e_i^2$$ converges) would follow. For simplicity, however, let's assume finite fourth moments for $$e_i$$ such that we could invoke elementary results here; for some real $$a>0$$, we have \begin{align} \Pr\left(|\bar z|>a\right)&\leq \frac{var(\bar z)}{a^2}\\ &= \frac{1}{n^2a^2}\sum_ivar(z_i)\\ &\leq \frac{1}{n^2a^2}nc\\ &=\frac c{a^2}n^{-1}\end{align} for some constant c. This in turn implies that $$\frac 1n \sum_i (e_i^2-\sigma_i^2) \overset{p}{\to}0$$. Therefore, $$\frac 1n \sum_i e_i^2\overset{p}{\to} \lim_{n\to \infty}\frac 1n \sum_i \sigma_i^2.$$

Note, however, that in a regression framework and when the focus is the OLS estimators, one focuses often on the standard errors of the OLS estimator and that would be slightly different.