Гетероскедастичность в R

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

library(knitr)
opts_chunk$set(cache=FALSE)

library(ggplot2) # графики
library(sandwich) # оценка Var для гетероскедастичности
library(bstats) # тест Уайта, тест Бройша-Пагана
library(lmtest) # тест Бройша-Пагана


theme_set(theme_bw())

Гетероскедастичность — \(Var(\epsilon_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"
h <- read.table(filename,header=TRUE)

qplot(totsp,price,data=h)

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

m1 <- lm(price~totsp,data=h)
summary(m1)
## 
## Call:
## lm(formula = price ~ totsp, data = h)
## 
## 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})\).

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)^{-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})\) в R.

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

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

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],sd_hc=conftable[,2])
ci$left_95 <- ci$estimate-qnorm(0.975)*ci$sd_hc
ci$right_95 <- ci$estimate+qnorm(0.975)*ci$sd_hc
ci
##               estimate     sd_hc    left_95   right_95
## (Intercept) -62.044844 7.8591534 -77.448501 -46.641186
## totsp         2.593462 0.1141737   2.369686   2.817239

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

model<- lm(price~totsp, weights=I(1/totsp), data=h)
summary(model)
## 
## Call:
## lm(formula = price ~ totsp, data = h, 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

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

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

  • Оценивается исходная регрессия и из неё получаются остатки \(\hat{\varepsilon}_i\).
  • Если верить в гомоскедастичность, то дисперсия вектора остатков определяется по формуле \(Var(\hat{\varepsilon})=\sigma^2 (I-X(X'X)^{-1}X')\). Т.е. дисперсия у остатков разная (!) даже если \(Var(\varepsilon_i)=\sigma^2\). Поэтому остатки стандартизуют по формуле: \[ s_i = \frac{\varepsilon_i}{\sqrt{\hat{\sigma}^2(1-h_{ii})}} \] Остатки полученные после такого преобразования называют стандартизированными или иногда стьюдентизированными.
  • Можно построить график зависимости величины \(s_i^2\) или \(|s_i|\) от переменной, предположительно повинной в гетроскедастичности.

В R:

m1.st.resid <- rstandard(m1) # получаем стандартизированные остатки
qplot(totsp,abs(m1.st.resid),data=h,alpha=0.05)

Иногда для простоты пропускают второй шаг со стандартизацией остатков и в результате могут ошибочно принять принять разную \(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-статистику можно получить с помощью вспомогательной регрессии. Авторская процедура:

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

Современная модификация выглядит (неизвестный рецензент Коэнкера) так: * Оценивается регрессия \(y_i=\beta_1+\beta_2 x_i +\varepsilon_i\) * Оценивается регрессия \(s_i^2=\alpha_1+\alpha_2 x_i +\varepsilon_i\) * Достоинством модификации является верный асимпотитический уровень значимости при нарушении нормальности (? Коэнкер в статье критикует что-то про мощность ?)

При верной \(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.9538, df = 1, p-value < 2.2e-16
bptest(m1,studentize=FALSE) # классика Бройша-Пагана
## 
##  Breusch-Pagan test
## 
## data:  m1
## BP = 2365.978, df = 1, p-value < 2.2e-16

Тест Goldfeld-Quandt.

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

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

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

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

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

Тест Уайта. White

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

Процедура: * Оценивается регрессия \(y_i=\beta_1+\beta_2 x_i +\varepsilon_i\) * Строим регрессию \(\hat{\varepsilon_i^2}\) на … * асимптотически \(nR^2\) имеет хи-квадрат распределение

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

white.test(m1)
## 
##  White test for constant variance
## 
## data:  
## White = 247.3515, df = 2, p-value < 2.2e-16

Преимущество теста Уайта над тестами Голдфельда-Квандта и Бройша-Пагана состоит в том, что не требуется нормальность остатков. Вероятно, это самый популярный тест на гетероскедастичность на сегодня.

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

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

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

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

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

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

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

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

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