Загружаем нужные пакеты:
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(ϵi|xi)≠const.
Последствия
Что делать?
yiˆσi=β11ˆσi+β2xiˆσi+ϵiˆσ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
Посмотрим на ^Var(ˆβ|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
Как было сказано ранее: если в данных присутствует гетероскедастичность, но нужно проверить гипотезы и строить доверительные интервалы, то нужно использовать правильную оценку ковариационной матрицы, она выглядит так: ^VarHC(ˆβ|X)=(X′X)−1X′ˆΩX(X′X)−1
ˆΩ=diag(ˆσ21,…,ˆσ2n)
Есть разные способы состоятельно оценить отдельные дисперсии ˆσ2i, подробно можно почитать в описании пакета sandwich
. Наиболее популярный на сегодня — это так называемый “HC3”: ˆσ2i=ˆϵ2i(1−hii)2
Здесь hii — диагональный элемент матрицы-шляпницы (hat matrix), ˆy=Hy, H=X(X′X)−1X′. Матрицу-шляпницу также называют проектором, т.к. она проецирует вектор y на линейную оболочку регрессоров.
Посмотрим на матрицу ^VarHC(ˆβ|X) в R.
vcovHC(m1)
## (Intercept) totsp
## (Intercept) 61.7662920 -0.89503034
## totsp -0.8950303 0.01303563
Применяя правильную ^VarHC(ˆβ|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-распределения при нормальных εi, ни даже асимптотически нормального. И такая же проблема возникает с F-статистикой для сравнения вложенных моделей. Отныне привычная F-статистика
F=(RSSR−RSSUR)/rRSSUR/(n−kUR) не имеет ни F-распределения при нормальных остатках, ни χ2-распределения асимптотически. В частности, для проверки незначимости регрессии в целом некорректно использовать F-статистику равную F=ESS/(k−1)RSS/(n−k).
Если для теста Вальда использовать корректную оценку ковариационной матрицы, то его можно применять в условиях гетероскедастичности. Допустим, мы хотим проверить гипотезу H0: Aβ=b, состоящую из q уравнений, против альтернативной о том, что хотя бы одно равенство нарушено. Тогда статистика Вальда имеет вид:
W=(Aβ−b)′(A⋅^VarHC(ˆβ)A′)−1(Aβ−b) и имеет асимптотически χ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(ϵi|totspi)=const⋅totspi.
В R вектор весов weights
должен быть обратно пропорционален дисперсиям, т.е. wi=1/σ2i.
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(ˆεi) за гетероскедастичность, т.е. за разную Var(εi). Впрочем, если гетероскедастичность сильная, то её будет видно и на нестандартизированных остатках.
Есть две версии теста Бройша-Пагана: авторская и современная модификация.
Предпосылки авторской версии:
Суть теста: Используя метод максимального правдоподобия посчитаем LM-статистику. При верной H0 она имеет хи-квадрат распределение с p−1 степенью свободы.
Оказывается, что LM-статистику можно получить с помощью вспомогательной регрессии. Авторская процедура:
Современная модификация выглядит (неизвестный рецензент Коэнкера) так:
Тест Бройша-Пагана в 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
У современной модификации Бройша-Пагана более слабые предпосылки:
Предпосылки:
Процедура:
Тест Уайта в 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 попросить сделать тест Бройша-Пагана и не указать спецификацию вспомогательной регрессии, то он возьмет только все регрессоры исходной.
Предпосылки:
Процедура:
При верной H0 RSS2RSS1∼Fr−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
Чайники часто используют такой ошибочный подход:
Протестировать гетероскедастичность. Если тесты выявили гетероскедастичность, то бороться с ней (использовать устойчивые стандартные ошибки или применять взвешенный мнк). Если тесты не выявили гетероскедастичность, то не бороться с ней.
Этот подход неверен!
Правильно делать так:
Тут обычно спрашивают: “А зачем же тогда тестировать на гетероскедастичность, если это решение ни на что не влияет?”. Ответ прост — чтобы узнать, есть ли гетероскедастичность. :) Впрочем, если есть теоретические основания её ожидать, то можно и не тестировать.
sandwich
с очень правильным изложением современного взгляда на гетероскедастичность