Загружаем нужные пакеты:

library("knitr") # создание отчётов
opts_chunk$set(fig.align = 'center') # выравнивание картинок по центру

library("ggplot2") # графики
library("sandwich") # оценка Var для гетероскедастичности
library("lmtest") # тест Бройша-Пагана
library("dplyr") # манипуляции с данными

# library("Greg") # пакет Greg Макса Гордона с готовой функцией для построения доверительных интервалов при гетероскедастичности; 
# ставится с github командой devtools::install_github("gforge/Greg")

theme_set(theme_bw()) # чёрно-белая тема оформления графиков

Гетероскедастичность — \(Var(\epsilon_i | x_i)\neq const\).

Последствия

Что делать?

\[ \frac{y_i}{\hat{\sigma}_i}=\beta_1 \frac{1}{\hat{\sigma}_i}+\beta_2 \frac{x_i}{\hat{\sigma}_i}+\frac{\epsilon_i}{\hat{\sigma}_i} \]

Файл с данными по стоимости квартир в Москве доступен по ссылке goo.gl/zL5JQ

filename <- "../datasets/flats_moscow.txt"
flats <- read.table(filename, header = TRUE)


ggplot(flats) + geom_point(aes(x = totsp, y = price)) + 
  labs(x = "Общая площадь квартиры, кв.м", y = "Цена квартиры, $1000", 
  title = "Стоимость квартир в Москве")

Строим регрессию стоимости на общую площадь.

m1 <- lm(price ~ totsp, data = flats)
summary(m1)
## 
## Call:
## lm(formula = price ~ totsp, data = flats)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -120.58  -17.44   -3.56   10.99  444.52 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -62.04484    3.71178  -16.72   <2e-16 ***
## totsp         2.59346    0.04973   52.15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.96 on 2038 degrees of freedom
## Multiple R-squared:  0.5716, Adjusted R-squared:  0.5714 
## F-statistic:  2719 on 1 and 2038 DF,  p-value: < 2.2e-16

Покажем доверительные интервалы для коэффициентов из регрессии.

confint(m1)
##                  2.5 %     97.5 %
## (Intercept) -69.324117 -54.765570
## totsp         2.495927   2.690998

Посмотрим на \(\widehat{Var}(\hat{\beta}|X)\).

vcov(m1)
##             (Intercept)        totsp
## (Intercept)  13.7772925 -0.180775186
## totsp        -0.1807752  0.002473516

Если не хотим смотреть на всё summary, а только на то, что касается бет с крышкой.

coeftest(m1)
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -62.044844   3.711778 -16.716 < 2.2e-16 ***
## totsp         2.593462   0.049734  52.146 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Как было сказано ранее: если в данных присутствует гетероскедастичность, но нужно проверить гипотезы и строить доверительные интервалы, то нужно использовать правильную оценку ковариационной матрицы, она выглядит так: \[ \widehat{Var}_{HC}(\hat{\beta}|X)=(X'X)^{-1}X'\hat{\Omega}X(X'X)^{-1} \]

\[ \hat{\Omega}=diag(\hat{\sigma}^2_1,\ldots,\hat{\sigma}^2_n) \]

Есть разные способы состоятельно оценить отдельные дисперсии \(\hat{\sigma}^2_i\), подробно можно почитать в описании пакета sandwich. Наиболее популярный на сегодня — это так называемый “HC3”: \[ \hat{\sigma}^2_i=\frac{\hat{\epsilon}_i^2}{(1-h_{ii})^2} \]

Здесь \(h_{ii}\) — диагональный элемент матрицы-шляпницы (hat matrix), \(\hat{y}=Hy\), \(H=X(X'X)^{-1}X'\). Матрицу-шляпницу также называют проектором, т.к. она проецирует вектор \(y\) на линейную оболочку регрессоров.

Посмотрим на матрицу \(\widehat{Var}_{HC}(\hat{\beta}|X)\) в R.

vcovHC(m1)
##             (Intercept)       totsp
## (Intercept)  61.7662920 -0.89503034
## totsp        -0.8950303  0.01303563

Применяя правильную \(\widehat{Var}_{HC}(\hat{\beta}|X)\), проверим гипотезы.

coeftest(m1, vcov = vcovHC(m1))
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -62.04484    7.85915 -7.8946 4.716e-15 ***
## totsp         2.59346    0.11417 22.7151 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Доверительный интервал надо строить руками. За основу мы возьмем табличку с правильными стандартными ошибками и прибавим и вычтем их от оценок коэффициентов.

conftable <- coeftest(m1, vcov = vcovHC(m1))
ci <- data_frame(estimate = conftable[, 1], 
                 se_hc = conftable[, 2],
                 left_95 = estimate - qnorm(0.975)*se_hc,
                 right_95 = estimate + qnorm(0.975)*se_hc)
ci
## # A tibble: 2 × 4
##     estimate     se_hc    left_95   right_95
##        <dbl>     <dbl>      <dbl>      <dbl>
## 1 -62.044844 7.8591534 -77.448501 -46.641186
## 2   2.593462 0.1141737   2.369686   2.817239

Ещё есть малоизвестный пакет Greg Макса Гордона и в нём есть готовая функция:

Greg::confint_robust(m1)

Пакет Greg не доступен на официальном репозитории CRAN, ставится с github командой devtools::install_github("gforge/Greg").

В условиях гетероскедастичности обычная \(t\)-статистика не имеет ни \(t\)-распределения при нормальных \(\varepsilon_i\), ни даже асимптотически нормального. И такая же проблема возникает с \(F\)-статистикой для сравнения вложенных моделей. Отныне привычная \(F\)-статистика

\[ F=\frac{(RSS_R-RSS_{UR})/r}{RSS_{UR}/(n-k_{UR})} \] не имеет ни \(F\)-распределения при нормальных остатках, ни \(\chi^2\)-распределения асимптотически. В частности, для проверки незначимости регрессии в целом некорректно использовать \(F\)-статистику равную \(F=\frac{ESS/(k-1)}{RSS/(n-k)}\).

Если для теста Вальда использовать корректную оценку ковариационной матрицы, то его можно применять в условиях гетероскедастичности. Допустим, мы хотим проверить гипотезу \(H_0\): \(A\beta=b\), состоящую из \(q\) уравнений, против альтернативной о том, что хотя бы одно равенство нарушено. Тогда статистика Вальда имеет вид:

\[ W = (A\beta - b)'(A\cdot \widehat{Var}_{HC}(\hat\beta)A')^{-1}(A\beta - b) \] и имеет асимптотически \(\chi^2\)-распределение с \(q\) степенями свободы.

В R реализуется с помощью опции для команды waldtest. Покажем на примере гипотезы о незначимости регрессии в целом:

m0 <- lm(price ~ 1, data = flats) # оцениваем ограниченную модель (регрессия на константу)
waldtest(m0, m1, vcov = vcovHC(m1))
## Wald test
## 
## Model 1: price ~ 1
## Model 2: price ~ totsp
##   Res.Df Df      F    Pr(>F)    
## 1   2039                        
## 2   2038  1 515.97 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Переходим ко взвешенному МНК. Взвешенный МНК требует знания структуры гетероскедастичности, что редко встречается на практике. Предположим, что в нашем случае \(Var(\epsilon_i | totsp_i) = const \cdot totsp_i\).

В R вектор весов weights должен быть обратно пропорционален дисперсиям, т.е. \(w_i=1/\sigma^2_i\).

m2 <- lm(price ~ totsp, weights = I(1 / totsp), data =flats)
summary(m2)
## 
## Call:
## lm(formula = price ~ totsp, data = flats, weights = I(1/totsp))
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.571  -1.994  -0.591   1.235  39.088 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -52.50302    3.53967  -14.83   <2e-16 ***
## totsp         2.46290    0.04931   49.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.58 on 2038 degrees of freedom
## Multiple R-squared:  0.5504, Adjusted R-squared:  0.5501 
## F-statistic:  2495 on 1 and 2038 DF,  p-value: < 2.2e-16

Тесты на гетероскедастичность

Графическое обнаружение гетероскедастичности:

В R:

m1.st.resid <- rstandard(m1) # получаем стандартизированные остатки
ggplot(aes(x = totsp, y = abs(m1.st.resid)), data =flats) + geom_point(alpha = 0.2) + 
  labs(x = "Общая площадь квартиры, кв.м", y = expression(paste("Стандартизированные остатки, ", s[i])),
  title = "Графическое обнаружение гетероскедастичности")