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

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 = "Графическое обнаружение гетероскедастичности")

Иногда для простоты пропускают второй шаг со стандартизацией остатков и в результате могут ошибочно принять принять разную \(Var(\hat{\varepsilon}_i)\) за гетероскедастичность, т.е. за разную \(Var(\varepsilon_i)\). Впрочем, если гетероскедастичность сильная, то её будет видно и на нестандартизированных остатках.

Тест Бройша-Пагана, Breusch-Pagan:

Есть две версии теста Бройша-Пагана: авторская и современная модификация.

Предпосылки авторской версии:

  • нормальность остатков, \(\varepsilon_i \sim N(0,\sigma^2_i)\)
  • \(\sigma^2_i=h(\alpha_1+\alpha_2 z_{i2}+\ldots+\alpha_{p}z_{ip})\)
  • у функции \(h\) существуют первая и вторая производные
  • тест асимптотический

Суть теста: Используя метод максимального правдоподобия посчитаем LM-статистику. При верной \(H_0\) она имеет хи-квадрат распределение с \(p-1\) степенью свободы.

Оказывается, что LM-статистику можно получить с помощью вспомогательной регрессии. Авторская процедура:

  1. Оценивается регрессия \(y_i=\beta_1+\beta_2 x_i +\varepsilon_i\)
  2. Переходим к \(g_i=\frac{n}{SSR}\hat{\varepsilon_i}\)
  3. Строим регрессию \(g_i\) на \(\alpha_1+\alpha_2 z_{i2}+\ldots+\alpha_{p}z_{ip}\)
  4. \(LM=\frac{ESS}{2}\)

Современная модификация выглядит (неизвестный рецензент Коэнкера) так:

  1. Оценивается регрессия \(y_i=\beta_1+\beta_2 x_i +\varepsilon_i\)
  2. Оценивается регрессия \(\hat{\varepsilon_i^2}\) на переменные вызывающие гетероскедастичность
  3. При верной \(H_0\) асимптотически: \[ nR^2 \sim \chi^2_{p-1}, \] где \(p\) — число оцениваемых коэффициентов во вспомогательной регрессии. По смыслу \((p-1)\) — это количество факторов, от которых потенциально может зависеть дисперсия \(Var(\varepsilon_i)\).

Тест Бройша-Пагана в R:

bptest(m1) # версия Коэнкера
## 
##  studentized Breusch-Pagan test
## 
## data:  m1
## BP = 201.95, df = 1, p-value < 2.2e-16
bptest(m1, studentize = FALSE) # классика Бройша-Пагана
## 
##  Breusch-Pagan test
## 
## data:  m1
## BP = 2366, df = 1, p-value < 2.2e-16

У современной модификации Бройша-Пагана более слабые предпосылки:

  • \(H_0\): \(Var(\varepsilon_i)=\sigma^2\), нормальность не требуется
  • \(E(\varepsilon_i^4)=const\)
  • тест асимптотический
  • (? Коэнкер в статье критикует что-то про мощность ?)

Тест Уайта. White test.

Предпосылки:

  • \(H_0\): \(Var(\varepsilon_i)=\sigma^2\), нормальность не требуется
  • \(E(\varepsilon_i^4)=const\)
  • тест асимптотический

Процедура:

  1. Оценивается регрессия \(y_i=\beta_1+\beta_2 x_i +\varepsilon_i\)
  2. Строим регрессию \(\hat{\varepsilon_i^2}\) на исходные регрессоры и их квадраты
  3. Асимптотически \(nR^2\) имеет хи-квадрат распределение

Тест Уайта в R:

bptest(m1, varformula = ~ totsp + I(totsp^2), data = flats)
## 
##  studentized Breusch-Pagan test
## 
## data:  m1
## BP = 247.35, df = 2, p-value < 2.2e-16

Тест Уайта является частным случаем современной модификации теста Бройша-Пагана. В тесте Бройша-Пагана во вспомогательной регресии можно брать любые объясняющие переменные. В тесте Уайта берутся исходные регрессоры, их квадраты и попарные произведения. Если в R попросить сделать тест Бройша-Пагана и не указать спецификацию вспомогательной регрессии, то он возьмет только все регрессоры исходной.

Тест Goldfeld-Quandt.

Предпосылки:

  • нормальность остатков, \(\varepsilon_i \sim N(0,\sigma^2_i)\)
  • Наблюдения упорядочены по возростанию гетероскедастичности
  • тест точный (неасимптотический)

Процедура:

  1. Упорядочить наблюдения в том порядке, в котором подозревается рост гетероскедастичности
  2. Выкинуть некий процент наблюдений по середине, чтобы подчеркнуть разницу в дисперсии между начальными и конечными наблюдениями.
  3. Оценить исходную модель по наблюдениям из начала выборки и по наблюдениям конца выборки
  4. Получить, соответственно, \(RSS_1\) и \(RSS_2\)

При верной \(H_0\) \[ \frac{RSS_2}{RSS_1}\sim F_{r-k,r-k}, \] где \(r\) — размер подвыборки в начале и в конце.

Тест Голдфельда-Квандта в R:

flats2 <- flats[order(flats$totsp), ] # сменим порядок строк в табличке h
m3 <- lm(price ~ totsp, data = flats2) 
gqtest(m3, fraction = 0.2) # проведем GQ тест выкинув посередине 20% наблюдений
## 
##  Goldfeld-Quandt test
## 
## data:  m3
## GQ = 8.2121, df1 = 814, df2 = 814, p-value < 2.2e-16

Тестирование и борьба с гетероскедастичностью

Чайники часто используют такой ошибочный подход:

Протестировать гетероскедастичность. Если тесты выявили гетероскедастичность, то бороться с ней (использовать устойчивые стандартные ошибки или применять взвешенный мнк). Если тесты не выявили гетероскедастичность, то не бороться с ней.

Этот подход неверен!

Правильно делать так:

  • Если есть теоретические основания ожидать гетероскедастичность и наблюдений достаточно, то использовать устойчивые стандартные ошибки вне зависимости от результатов тестирования.
  • Если есть теоретические основания ожидать гетероскедастичность и наблюдений мало, то отказаться от девичьих грёз об эффективности оценок и проверке гипотез. Довольствоваться возможной несмещенностью оценок.

Тут обычно спрашивают: “А зачем же тогда тестировать на гетероскедастичность, если это решение ни на что не влияет?”. Ответ прост — чтобы узнать, есть ли гетероскедастичность. :) Впрочем, если есть теоретические основания её ожидать, то можно и не тестировать.

Дополнительно: