Семинар 5. Множественная регрессия
Июнь, 14, 2018
Шаманское заклинание для настройки глобальных опций отчёта:
1 Регрессия при гомоскедастичных ошибках
library(tidyverse) # обработка данных, графики...
library(skimr) # описательные статистики
library(rio) # импорт фантастического количества форматов данных
library(broom) # метла превращает результаты оценивания моделей в таблички
library(GGally) # больше готовых графиков
library(sjPlot) # ещё больше графиков
library(lmtest) # диагностика линейных моделей
library(Ecdat) # много-много разных наборов данных
library(sjstats) # удобные мелкие функции для работы с моделями
library(sandwich) # оценка Var для гетероскедастичности
library(AER) # работа с инструментальными переменнымиВозьмём набор данных об изменении пульса студентов до и после упражнений. Его описание доступно по ссылке.
Загрузим данные и посмотрим на описательные статистики.
Установите в качестве текущей папки ту папку, где лежит данный .Rmd-файл: Session — Set working directory — To source file location.
Правильно укажите путь к данным относительно .Rmd-файла:
pulse <- import('data/pulse.txt')
skim(pulse)## Skim summary statistics
## n obs: 110
## n variables: 11
##
## Variable type: integer
## variable missing complete n mean sd p0 p25 p50 p75 p100
## Age 0 110 110 20.56 3.92 18 19 20 21 45
## Alcohol 0 110 110 1.38 0.49 1 1 1 2 2
## Exercise 0 110 110 2.21 0.65 1 2 2 3 3
## Gender 0 110 110 1.46 0.5 1 1 1 2 2
## Height 0 110 110 171.58 16.08 68 165.25 172.5 180 195
## Pulse1 1 109 110 75.69 13.3 47 68 76 82 145
## Pulse2 1 109 110 96.8 31.57 56 72 84 125 176
## Ran 0 110 110 1.58 0.5 1 1 2 2 2
## Smokes 0 110 110 1.9 0.3 1 2 2 2 2
## Year 0 110 110 95.63 1.75 93 95 96 97 98
## hist
## ▇▁▁▁▁▁▁▁
## ▇▁▁▁▁▁▁▅
## ▂▁▁▇▁▁▁▅
## ▇▁▁▁▁▁▁▇
## ▁▁▁▁▁▂▇▅
## ▁▇▇▃▁▁▁▁
## ▅▇▃▁▃▂▂▁
## ▆▁▁▁▁▁▁▇
## ▁▁▁▁▁▁▁▇
## ▇▁▁▇▆▁▇▆
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25 p50 p75 p100 hist
## Weight 0 110 110 66.33 15.16 27 56.25 63 75 110 ▁▁▇▇▅▃▁▁- Сколько в данных наблюдений? Сколько переменных?
- Какие переменные на самом деле являются факторными?
Вернём факторным переменным положенный статус!
pulse_fct <- pulse %>%
mutate_at(vars(-Weight, -Height, -Age, -Pulse1, -Pulse2), factor)И теперь можно приступать к регрессиям. Начнём с парной и построим регрессию пульса после упражнений Pulse2 на константу и пульс до упражнений Pulse1. Для этого воспользуемся командой lm() и передадим ей данные и формулу.
model_r <- lm(data = pulse_fct, Pulse2 ~ Pulse1)С помощью команды summary() посмотрим на описание модели.
summary(model_r)##
## Call:
## lm(formula = Pulse2 ~ Pulse1, data = pulse_fct)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.09 -23.07 -17.92 28.08 72.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.8476 16.4415 1.937 0.055379 .
## Pulse1 0.8581 0.2140 4.010 0.000113 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29.57 on 107 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1307, Adjusted R-squared: 0.1225
## F-statistic: 16.08 on 1 and 107 DF, p-value: 0.0001125- На сколько в среднем отличается пульс после упражнений у двух студентов, если пульс до упражнений у них отличается на единицу?
- Какие коэффициенты значимы?
Какая доля разброса пульса после упражнений объясняется пульсом до упражнений?
Упражнение 1.
Упражнения будем выполнять на встроенном наборе данных Housing. Посмотрите описание переменных в справке :) Затем «бросьте взгляд» на данные, есть ли в них факторные переменные, которые объявлены иначе?
house <- Housing
glimpse(house)## Observations: 546
## Variables: 12
## $ price <dbl> 42000, 38500, 49500, 60500, 61000, 66000, 66000, 6900...
## $ lotsize <dbl> 5850, 4000, 3060, 6650, 6360, 4160, 3880, 4160, 4800,...
## $ bedrooms <dbl> 3, 2, 3, 3, 2, 3, 3, 3, 3, 3, 3, 2, 3, 3, 2, 2, 3, 4,...
## $ bathrms <dbl> 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1,...
## $ stories <dbl> 2, 1, 1, 2, 1, 1, 2, 3, 1, 4, 1, 1, 2, 1, 1, 1, 2, 3,...
## $ driveway <fct> yes, yes, yes, yes, yes, yes, yes, yes, yes, yes, yes...
## $ recroom <fct> no, no, no, yes, no, yes, no, no, yes, yes, no, no, n...
## $ fullbase <fct> yes, no, no, no, no, yes, yes, no, yes, no, yes, no, ...
## $ gashw <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, n...
## $ airco <fct> no, no, no, no, no, yes, no, no, no, yes, yes, no, no...
## $ garagepl <dbl> 1, 0, 0, 0, 0, 0, 2, 0, 0, 1, 3, 0, 0, 0, 0, 0, 1, 0,...
## $ prefarea <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, n...Постройте парную регрессию цены дома price на константу и площадь дома lotsize.
house_r <- lm(data = house, price ~ lotsize)
summary(house_r)##
## Call:
## lm(formula = price ~ lotsize, data = house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -69551 -14626 -2858 9752 106901
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.414e+04 2.491e+03 13.7 <2e-16 ***
## lotsize 6.599e+00 4.458e-01 14.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22570 on 544 degrees of freedom
## Multiple R-squared: 0.2871, Adjusted R-squared: 0.2858
## F-statistic: 219.1 on 1 and 544 DF, p-value: < 2.2e-16- Совпадает ли знак коэффициента при
lotsizeс ожидаемым? - Можно ли считать коэффициента при
lotsizeзначимым на уровне 1%?
Другой способ посмотреть на описание модели — воспользоваться функциями пакета broom. Оценки коэффциентов, стандартные ошибки, p-значения выводит команда tidy():
tidy(model_r)## term estimate std.error statistic p.value
## 1 (Intercept) 31.8476001 16.4414717 1.937029 0.0553785865
## 2 Pulse1 0.8581347 0.2139792 4.010365 0.0001125123А функция glance() покажет общие характеристики модели, для которых достаточно одной строчки, — коэффциент детерминации, значение лог-функции правдоподобия, AIC, BIC…
glance(model_r)## r.squared adj.r.squared sigma statistic p.value df logLik
## 1 0.1306681 0.1225435 29.5705 16.08303 0.0001125123 2 -522.8137
## AIC BIC deviance df.residual
## 1 1051.627 1059.701 93562.33 107Если из общей таблицы описания модели нужна только информация о коэффициентах, то поможет команда coeftest() из пакета lmtest().
coeftest(model_r)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.84760 16.44147 1.9370 0.0553786 .
## Pulse1 0.85813 0.21398 4.0104 0.0001125 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Нарисуем линию регрессии для модели, в которой зависимая переменная — это пульс после упражнений, а объясняющая — пульс до упражнений. Для этого к диаграмме рассеяния добавим дополнительный слой geom_smooth(). Внутри него за линейную регрессию отвечает method = 'lm'.
ggplot(data = pulse_fct, aes(x = Pulse1, y = Pulse2)) +
geom_point() +
geom_smooth(method = 'lm')
- Упражнение 2.
Визуализируйте линию регрессии для модели house_r. В ней зависимой переменной была цена дома price, а объясняющей — площадь lotsize.
ggplot(data = house, aes(x = lotsize, y = price)) +
geom_point() +
geom_smooth(method = 'lm')
Теперь оценим более сложную модель. К объясняющим переменным добавим вес студента Weight и бинарные переменные Ran, где 2 соответствует студенту, который выполнял упражнения , и Smokes, где значение 2 — это студент-курильщик.
model_ur <- lm(data = pulse_fct, Pulse2 ~ Weight + Pulse1 + Ran + Smokes)
summary(model_ur)##
## Call:
## lm(formula = Pulse2 ~ Weight + Pulse1 + Ran + Smokes, data = pulse_fct)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.631 -3.717 0.170 4.360 42.354
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.12572 12.14842 4.620 1.10e-05 ***
## Weight 0.02187 0.09113 0.240 0.811
## Pulse1 0.89417 0.10421 8.580 9.77e-14 ***
## Ran2 -52.20287 2.75215 -18.968 < 2e-16 ***
## Smokes2 1.90805 4.52173 0.422 0.674
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.1 on 104 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8078, Adjusted R-squared: 0.8004
## F-statistic: 109.3 on 4 and 104 DF, p-value: < 2.2e-16- Упражнение 3.
Оцените более сложную модель для набора данных о домах. Среди объясняющих переменных укажите площадь дома lotsize, количсевто спален bedrooms и бинарную переменную prefarea, которая равна 1 для домов, расположенных в хороших районах.
house_ur <- lm(data = house, price ~ lotsize + bedrooms + prefarea)
summary(house_ur)##
## Call:
## lm(formula = price ~ lotsize + bedrooms + prefarea, data = house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -58243 -13182 -2681 9315 89082
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.515e+03 3.982e+03 1.636 0.102
## lotsize 5.485e+00 4.225e-01 12.984 < 2e-16 ***
## bedrooms 1.024e+04 1.211e+03 8.457 2.57e-16 ***
## prefareayes 1.273e+04 2.142e+03 5.944 4.97e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20590 on 542 degrees of freedom
## Multiple R-squared: 0.4088, Adjusted R-squared: 0.4055
## F-statistic: 124.9 on 3 and 542 DF, p-value: < 2.2e-16Есть ли проблемы у модели? Диагностический график в студию!
ggnostic(model = model_ur)
На этом графике четыре линии! * Первая: зависимость величины остатка \(\hat u_i = y_i - \hat y_i\) от каждого регрессора. * Вторая: зависимость величины \(\hat \sigma_i\) от каждого регрессора.
Величина \(\hat \sigma_i\) показывает, какой была бы оценка стандартного отклонения случайной составляющей, если бы \(i\)-ое наблюдение выкинули из модели.
Слишком маленькие значения \(\hat \sigma_i\) говорят о том, что при выкидывании наблюдения резко падает оценка стандартного отклонения. Значит, наблюдение могло быть выбросом.
- Третья: зависимость величины \(h_{ii}\) от каждого регрессора.
Величина \(h_{ii}\) имеет две интерпретации. Во-первых, она показывает, насколько изменится прогноз \(\hat y_i\) при увеличении фактического наблюдения \(y_i\) на единицу. Во-вторых, эта величина показывает отношение дисперсии прогноза \(\hat y_i\) к дисперсии случайной составляющей: \[ h_{ii} = \frac{Var(\hat y_i)}{Var(u_i)} \]
- Четвертая: зависимость расстояния Кука от каждого из регрессоров.
Расстояние Кука показывает, насколько сильно изменятся прогнозы по всему набору данных, при выкидывании \(i\)-го наблюдения,
\[ D_i = \frac{\sum_{j=1}^n(\hat y_j - \hat y_j^{(-i)})^2}{k\hat\sigma^2}, \] где \(k\) — число регрессоров, считая единичный столбец.
Подробнее про расстояние Кука можно прочитать в маленьких заметках.
- Какие наблюдения являются потенциальными выбросами?
Построим 5%-ые доверительные интервалы для всех коэффициентов модели model_ur. Для этого будем использовать команду coefci(). Поменять уровень значимости можно аргументом level. Затем визуилизируем получившиеся доверительные интервалы командой sjp.lm() из пакета sjPlot.
coefci(model_ur)## 2.5 % 97.5 %
## (Intercept) 32.0349494 80.2164943
## Weight -0.1588415 0.2025909
## Pulse1 0.6875121 1.1008223
## Ran2 -57.6604911 -46.7452493
## Smokes2 -7.0587150 10.8748163plot_model(model_ur, ci.lvl = 0.9)
Загляните в справку :) Какая опция функции
coefciотвечает за уровень доверия?Упражнение 4.
Постройте 95%-ные доверительные интервалы для всех коэффициентов модели house_ur и визуализируйте их.
coefci(house_ur)## 2.5 % 97.5 %
## (Intercept) -1307.098806 14336.341933
## lotsize 4.655317 6.315047
## bedrooms 7863.559266 12622.144408
## prefareayes 8524.476076 16938.849973plot_model(house_ur, ci.lvl = 0.95)
Если есть две вложенных модели, то есть одна является частным случаем другой, то можно провести тест Вальда, чтобы выбрать одну из них.
\(H_0\): верна ограниченная (короткая) модель;
\(H_a\): верная неограниченная (длинная) модель;
waldtest(model_r, model_ur)## Wald test
##
## Model 1: Pulse2 ~ Pulse1
## Model 2: Pulse2 ~ Weight + Pulse1 + Ran + Smokes
## Res.Df Df F Pr(>F)
## 1 107
## 2 104 3 122.16 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Сравниваем p-значение с уровнем значимости.
- Упражнение 5.
Проведите тест Вальда для моделей house_r и huose_ur.
waldtest(house_r, house_ur)## Wald test
##
## Model 1: price ~ lotsize
## Model 2: price ~ lotsize + bedrooms + prefarea
## Res.Df Df F Pr(>F)
## 1 544
## 2 542 2 55.804 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Какую модель следует выбрать по результатам теста на 5%-ом уровне значимости?
2 О бедном t-тесте замолвите слово!
Иногда душа хочет проверки самых простых гипотез!
Если хочется построить доверительный интервал для математического ожидания пульса до упражнений, то достаточно построить регрессию на константу:
model_mu <- lm(data = pulse_fct, Pulse1 ~ 1)
coefci(model_mu)## 2.5 % 97.5 %
## (Intercept) 73.16341 78.21274Если хочется построить доверительный интервал для разности математического ожидания пульса у курящих и некурящих, то можно построить регрессию на константу и дамми-переменную.
model_diff <- lm(data = pulse_fct, Pulse1 ~ Smokes)
coefci(model_diff)## 2.5 % 97.5 %
## (Intercept) 69.56907 85.521838
## Smokes2 -10.47800 6.346272- Упражнение 6.
Постройте доверительный интервал для математического ожидания цены и доверительный интервал для разности математических ожиданий цены дома в хорошем и обычном районе.
house_mu <- lm(data = house, price ~ 1)
summary(house_mu)##
## Call:
## lm(formula = price ~ 1, data = house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43122 -18997 -6122 13878 121878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68122 1143 59.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26700 on 545 degrees of freedomhouse_diff <- lm(data = house, price ~ prefarea)
summary(house_diff)##
## Call:
## lm(formula = price ~ prefarea, data = house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52086 -17263 -4486 10664 111737
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63264 1234 51.248 < 2e-16 ***
## prefareayes 20723 2550 8.128 2.96e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25240 on 544 degrees of freedom
## Multiple R-squared: 0.1083, Adjusted R-squared: 0.1067
## F-statistic: 66.06 on 1 and 544 DF, p-value: 2.959e-15Дисперсионный анализ легким движением руки также заменяется на построение единственной регрессии.
Проверим, что две качественные переменные, Smokes и Gender, не оказывают влияния на частоту пульса.
model_2anova <- lm(data = pulse_fct, Pulse1 ~ Smokes * Gender)
summary(model_2anova)##
## Call:
## lm(formula = Pulse1 ~ Smokes * Gender, data = pulse_fct)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.298 -8.375 0.196 6.196 71.196
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.3750 4.7203 16.180 <2e-16 ***
## Smokes2 -2.5711 5.0771 -0.506 0.614
## Gender2 4.2917 9.0387 0.475 0.636
## Smokes2:Gender2 -0.7977 9.4333 -0.085 0.933
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.35 on 105 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.01995, Adjusted R-squared: -0.008052
## F-statistic: 0.7124 on 3 and 105 DF, p-value: 0.5467anova_stats(model_2anova)## # A tibble: 4 x 12
## term df sumsq meansq statistic p.value etasq partial.etasq omegasq
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Smok… 1 4.22e1 42.2 0.237 0.628 0.002 0.002 -0.007
## 2 Gend… 1 3.38e2 338. 1.89 0.172 0.018 0.018 0.008
## 3 Smok… 1 1.27e0 1.27 0.007 0.933 0 0 -0.009
## 4 Resi… 105 1.87e4 178. NA NA NA NA NA
## # ... with 3 more variables: partial.omegasq <dbl>, cohens.f <dbl>,
## # power <dbl>- Упражнение 7.
Проверьте, что наличие подъезда к дому driveway и расположение в хорошем районе prefarea, не влияют на цену дома.
house_2anova <- lm(data = house, price ~ driveway * prefarea)
summary(house_2anova)##
## Call:
## lm(formula = price ~ driveway * prefarea, data = house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -41441 -16441 -4441 11163 108559
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48732 2822 17.271 < 2e-16 ***
## drivewayyes 17709 3115 5.685 2.14e-08 ***
## prefareayes -6782 17508 -0.387 0.699
## drivewayyes:prefareayes 24995 17692 1.413 0.158
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24440 on 542 degrees of freedom
## Multiple R-squared: 0.1672, Adjusted R-squared: 0.1626
## F-statistic: 36.27 on 3 and 542 DF, p-value: < 2.2e-163 Регрессия при гетероскедастичных ошибках
Если не предполагать, что дисперсии ошибок \(Var(u_i)\) одинаковы для всех наблюдений, то построенные нами доверительные интервалы и выполненная проверка гипотез — полный отстой :)
Повторяем табличку с тестами на значимость отдельных коэффициентов с учётом поправки на гетероскедастичность.
coeftest(model_ur, vcov. = vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.125722 21.191573 2.6485 0.009345 **
## Weight 0.021875 0.095912 0.2281 0.820039
## Pulse1 0.894167 0.217629 4.1087 7.954e-05 ***
## Ran2 -52.202870 3.252667 -16.0492 < 2.2e-16 ***
## Smokes2 1.908051 2.516461 0.7582 0.450029
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Также легко получить доверительные интервалы:
coefci(model_ur, vcov. = vcovHC)## 2.5 % 97.5 %
## (Intercept) 14.1020380 98.1494057
## Weight -0.1683218 0.2120713
## Pulse1 0.4626001 1.3257342
## Ran2 -58.6530314 -45.7527089
## Smokes2 -3.0821863 6.8982876Результаты в R чуть отличаются от результатов в stata. По-умолчанию, R использует корректировку HC3 для подсчёта стандартных ошибок коэффициентов, а stata — менее удачную, устаревшую HC1. Сравнение есть у Achim Zeileis. Если хочется воспроизвести именно корректировку HC1, то запросто :)
coeftest(model_ur, vcov. = vcovHC(model_ur, type = "HC1"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.125722 16.388352 3.4247 0.0008819 ***
## Weight 0.021875 0.084090 0.2601 0.7952746
## Pulse1 0.894167 0.165376 5.4069 4.111e-07 ***
## Ran2 -52.202870 3.121163 -16.7255 < 2.2e-16 ***
## Smokes2 1.908051 2.250136 0.8480 0.3984014
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Сравним две вложенных модели с учётом поправки на гетероскедастичность:
waldtest(model_r, model_ur, vcov = vcovHC)## Wald test
##
## Model 1: Pulse2 ~ Pulse1
## Model 2: Pulse2 ~ Weight + Pulse1 + Ran + Smokes
## Res.Df Df F Pr(>F)
## 1 107
## 2 104 3 91.751 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Упражнение 8.
Получите таблицу с тестами на значимость коэффициентов в модели house_ur, используя корректировку HC3. А затем сравните две модели house_r и house_ur с учётом поправки на гетероскедастичность.
coeftest(house_ur, vcov. = vcovHC, type = 'HC3')##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.5146e+03 3.8045e+03 1.7123 0.08741 .
## lotsize 5.4852e+00 5.2896e-01 10.3697 < 2.2e-16 ***
## bedrooms 1.0243e+04 1.2472e+03 8.2128 1.593e-15 ***
## prefareayes 1.2732e+04 2.2766e+03 5.5924 3.555e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1waldtest(house_r, house_ur, vcov = vcovHC)## Wald test
##
## Model 1: price ~ lotsize
## Model 2: price ~ lotsize + bedrooms + prefarea
## Res.Df Df F Pr(>F)
## 1 544
## 2 542 2 52.281 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14 Инструментальные переменные
Исследователи всегда мечтают обнаружить не просто статистическую связь, а причинно-следственную. Проще всего было бы обнаружить связь с помощью рандомизированного эксперимента: поделить случайно индивидов случайно на две группы, одну попросить стирать обычным порошком, а вторую — новым чудопорошком. Важно упаковать порошки совершенно одинаково и говорить одни и те же слова :)
Мы часто имеем с данными наблюдений. Рандомизированного эксперимента при этом не было. В редких случаях можно попытаться обнаружить причинно-следственную связь на данных наблюдений.
Рассмотрим набор данных по предложению труда женщин Mroz из пакета Ecdat.
mroz <- Ecdat::Mroz
glimpse(mroz)## Observations: 753
## Variables: 18
## $ work <fct> no, no, no, no, no, no, no, no, no, no, no, no, no,...
## $ hoursw <int> 1610, 1656, 1980, 456, 1568, 2032, 1440, 1020, 1458...
## $ child6 <int> 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, ...
## $ child618 <int> 0, 2, 3, 3, 2, 0, 2, 0, 2, 2, 1, 1, 2, 2, 1, 3, 2, ...
## $ agew <int> 32, 30, 35, 34, 31, 54, 37, 54, 48, 39, 33, 42, 30,...
## $ educw <int> 12, 12, 12, 12, 14, 12, 16, 12, 12, 12, 12, 11, 12,...
## $ hearnw <dbl> 3.3540, 1.3889, 4.5455, 1.0965, 4.5918, 4.7421, 8.3...
## $ wagew <dbl> 2.65, 2.65, 4.04, 3.25, 3.60, 4.70, 5.95, 9.98, 0.0...
## $ hoursh <int> 2708, 2310, 3072, 1920, 2000, 1040, 2670, 4120, 199...
## $ ageh <int> 34, 30, 40, 53, 32, 57, 37, 53, 52, 43, 34, 47, 33,...
## $ educh <int> 12, 9, 12, 10, 12, 11, 12, 8, 4, 12, 12, 14, 16, 12...
## $ wageh <dbl> 4.0288, 8.4416, 3.5807, 3.5417, 10.0000, 6.7106, 3....
## $ income <int> 16310, 21800, 21040, 7300, 27300, 19495, 21152, 189...
## $ educwm <int> 12, 7, 12, 7, 12, 14, 14, 3, 7, 7, 12, 14, 16, 10, ...
## $ educwf <int> 7, 7, 7, 7, 14, 7, 7, 3, 7, 7, 3, 7, 16, 10, 7, 10,...
## $ unemprate <dbl> 5.0, 11.0, 5.0, 5.0, 9.5, 7.5, 5.0, 5.0, 3.0, 5.0, ...
## $ city <fct> no, yes, no, no, yes, yes, no, no, no, no, no, no, ...
## $ experience <int> 14, 5, 15, 6, 7, 33, 11, 35, 24, 21, 15, 14, 0, 14,...Отберём работавших женщин.
labor <- filter(mroz, wagew > 0)Мы хотим оценить причинно-следственный эффект от дополнительного года обучения на заработную плату женщины.
На первый взгляд разумно попробовать обычный МНК:
model_lm <- lm(data = labor, log(wagew) ~ educw + experience + I(experience^2))
summary(model_lm)##
## Call:
## lm(formula = log(wagew) ~ educw + experience + I(experience^2),
## data = labor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.22168 -0.22193 0.00529 0.24018 1.17754
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0310084 0.1295490 0.239 0.81098
## educw 0.0834844 0.0090775 9.197 < 2e-16 ***
## experience 0.0238019 0.0087690 2.714 0.00699 **
## I(experience^2) -0.0003113 0.0002532 -1.229 0.21983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3784 on 332 degrees of freedom
## Multiple R-squared: 0.2524, Adjusted R-squared: 0.2457
## F-statistic: 37.37 on 3 and 332 DF, p-value: < 2.2e-16Эта регрессия прекрасно оценивает силу статистической взаимосвязи и показывает, насколько различается зарплата двух женщин, у которых разное образование, и одинаковый опыт работы.
Эта регрессия вовсе не говорит, насколько увеличится зарплата данной женщины, если она надумает ещё год поучиться.
Проблема состоит в ненаблюдаемых характеристиках, например, в способностях жещнины. Логично предположить, что ненаблюдаемые способности связаны положительно и с зарплатой и с продолжительностью обучения.
При всей спорности предположения, будем считать, что \[ \log(wagew_i) = \beta_1 + \beta_2 educw_i + \beta_3 experience_i + \beta_4 experience_i^2 + \beta_5 ability_i + u_i, \] где ошибка \(u_i\) некоррелирована с остальными переменными в правой части.
Предположим, что образование матери и отца, educwm и educwf не коррелированы с ошибкой \(u_i\). Тогда их можно использовать как инструменты для переменной educw_i.
В качестве инструментальных переменных для educw возьмём образование матери и отца, educwm и educwf:
model_iv <- ivreg(data = labor,
log(wagew) ~ educw + experience + I(experience^2) |
experience + I(experience^2) + educwm + educwf)
summary(model_iv, diagnostics = TRUE)##
## Call:
## ivreg(formula = log(wagew) ~ educw + experience + I(experience^2) |
## experience + I(experience^2) + educwm + educwf, data = labor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2646 -0.2287 0.0176 0.2501 1.1828
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0924839 0.2642086 -0.350 0.72653
## educw 0.0934865 0.0207397 4.508 9.1e-06 ***
## experience 0.0232873 0.0088372 2.635 0.00881 **
## I(experience^2) -0.0002929 0.0002560 -1.144 0.25335
##
## Diagnostic tests:
## df1 df2 statistic p-value
## Weak instruments 2 331 39.395 4.49e-16 ***
## Wu-Hausman 1 331 0.288 0.592
## Sargan 1 NA 0.136 0.713
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3791 on 332 degrees of freedom
## Multiple R-Squared: 0.2497, Adjusted R-squared: 0.2429
## Wald test: 15.91 on 3 and 332 DF, p-value: 1.085e-09После вертикальной палочки | в команде ivreg указывают экзогенные переменные, используемые в качестве инструментов.
Если верить в то, что образование матери и отца некоррелировано с ненаблюдаемыми способностями женщины, тогда коэффициент при образовании женщины educw показывает, насколько изменится зарплата женщины при увеличении образования на один год, постоянном опыте и постоянных способностях.
Автоматически команда summary с опцией diagonostics = TRUE проводит три теста:
- Тест на слабые инструменты
Нулевая гипотеза: \(H_0\): инструментальные переменные не коррелированы с потенциально эндогенным регрессором
- Тест Хаусмана на равенство оценок IV и обычного МНК.
Нулевая гипотеза: \(H_0\): ковариация потенциально эндогенного регрессора с ошибкой равна нулю.
- Тест Саргана
Нулевая гипотеза: \(H_0\): инструментальные переменные правильно подобраны, то есть они не коррелированы с ошибкой, а невключённые инструменты и не должны быть включены в модель.
Какие гипотезы отвергаются, какие не отвергаются в нашем случае?
Упражнение 9.
Для упражнения возьмём данные о потреблении сигарет в США CigaretteeSW из пакета AER. Его описание можно прочесть в справке!
Возьмём наблюдения только за 1995 год и создадим несколько новых переменных: логарифм реальных цен lrprice, логарифм дохода на душу населения lrincome, логарифм количества пачек сигарет на человека lquant и реальный налог на сигареты tdiff.
data("CigarettesSW")
cig <- subset(CigarettesSW, year == 1995)
cig$lrprice <- log(cig$price / cig$cpi)
cig$lrincome <- log(cig$income / cig$population / cig$cpi)
cig$lquant <- log(cig$packs)
cig$tdiff <- (cig$taxs - cig$tax) / cig$cpi
glimpse(cig)## Observations: 48
## Variables: 13
## $ state <fct> AL, AR, AZ, CA, CO, CT, DE, FL, GA, IA, ID, IL, IN,...
## $ year <fct> 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 199...
## $ cpi <dbl> 1.524, 1.524, 1.524, 1.524, 1.524, 1.524, 1.524, 1....
## $ population <dbl> 4262731, 2480121, 4306908, 31493524, 3738061, 32652...
## $ packs <dbl> 101.08543, 111.04297, 71.95417, 56.85931, 82.58292,...
## $ income <dbl> 83903280, 45995496, 88870496, 771470144, 92946544, ...
## $ tax <dbl> 40.50000, 55.50000, 65.33333, 61.00000, 44.00000, 7...
## $ price <dbl> 158.3713, 175.5425, 198.6075, 210.5047, 167.3500, 2...
## $ taxs <dbl> 41.90467, 63.85917, 74.79082, 74.77133, 44.00000, 8...
## $ lrprice <dbl> 4.643604, 4.746543, 4.869992, 4.928169, 4.698749, 4...
## $ lrincome <dbl> 2.558416, 2.498898, 2.605622, 2.777178, 2.792119, 3...
## $ lquant <dbl> 4.615966, 4.709917, 4.276029, 4.040580, 4.413803, 4...
## $ tdiff <dbl> 0.9216975, 5.4850193, 6.2057067, 9.0363074, 0.00000...skim(cig)## Skim summary statistics
## n obs: 48
## n variables: 13
##
## Variable type: factor
## variable missing complete n n_unique top_counts ordered
## state 0 48 48 48 AL: 1, AR: 1, AZ: 1, CA: 1 FALSE
## year 0 48 48 1 199: 48, 198: 0, NA: 0 FALSE
##
## Variable type: numeric
## variable missing complete n mean sd p0
## cpi 0 48 48 1.52 0 1.52
## income 0 48 48 1.3e+08 1.4e+08 1e+07
## lquant 0 48 48 4.54 0.24 3.9
## lrincome 0 48 48 2.68 0.14 2.42
## lrprice 0 48 48 4.78 0.13 4.56
## packs 0 48 48 96.33 23.78 49.27
## population 0 48 48 5426461.81 5800282.87 478447
## price 0 48 48 183.25 24.05 145.98
## tax 0 48 48 53.69 15.68 26.5
## taxs 0 48 48 61.87 18.48 34.44
## tdiff 0 48 48 5.36 2.85 0
## p25 p50 p75 p100 hist
## 1.52 1.52 1.52 1.52 ▁▁▁▇▁▁▁▁
## 3.5e+07 8.4e+07 1.6e+08 7.7e+08 ▇▃▁▁▁▁▁▁
## 4.38 4.53 4.69 5.15 ▁▁▃▇▇▇▁▁
## 2.56 2.67 2.78 3.04 ▂▇▅▇▇▃▂▁
## 4.69 4.75 4.87 5.06 ▂▃▇▃▅▂▃▁
## 80.23 92.84 109.27 172.65 ▂▅▇▆▅▁▁▁
## 1670598.25 3796654.5 6197099.25 3.1e+07 ▇▃▁▁▁▁▁▁
## 166.61 176.15 198.51 240.85 ▅▇▅▃▃▂▂▁
## 42 51.25 62.5 99 ▂▇▆▇▃▃▁▁
## 48.75 59.84 74.78 112.63 ▇▇▅▇▅▃▂▁
## 4.64 5.63 7.28 10.26 ▅▁▁▂▇▅▂▂Постройте регрессию логарифма количества пачек сигарет на человека lquant на логарифм реальной цены одной пачки lrprice, используя в качестве инстурмента реальный налог tdiff.
cig_iv1 <- ivreg(data = cig, lquant ~ lrprice | tdiff)
summary(cig_iv1, diagnostics = TRUE)##
## Call:
## ivreg(formula = lquant ~ lrprice | tdiff, data = cig)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64619 -0.07732 0.02981 0.11283 0.41904
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.7199 1.5141 6.420 6.79e-08 ***
## lrprice -1.0836 0.3166 -3.422 0.00131 **
##
## Diagnostic tests:
## df1 df2 statistic p-value
## Weak instruments 1 46 40.956 7.27e-08 ***
## Wu-Hausman 1 45 0.314 0.578
## Sargan 0 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1904 on 46 degrees of freedom
## Multiple R-Squared: 0.4011, Adjusted R-squared: 0.3881
## Wald test: 11.71 on 1 and 46 DF, p-value: 0.001313Добавьте к объясняющим переменным логарифм реального дохода на душу населения lrincome, а к инструментам — реальный налог на сигареты tax / cpi.
cig_iv2 <- ivreg(data = cig, lquant ~ lrprice + lrincome | tdiff + lrincome + I(tax/cpi))
summary(cig_iv2, diagnostics = TRUE)##
## Call:
## ivreg(formula = lquant ~ lrprice + lrincome | tdiff + lrincome +
## I(tax/cpi), data = cig)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6006931 -0.0862222 -0.0009999 0.1164699 0.3734227
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.8950 1.0586 9.348 4.12e-12 ***
## lrprice -1.2774 0.2632 -4.853 1.50e-05 ***
## lrincome 0.2804 0.2386 1.175 0.246
##
## Diagnostic tests:
## df1 df2 statistic p-value
## Weak instruments 2 44 244.734 <2e-16 ***
## Wu-Hausman 1 44 3.068 0.0868 .
## Sargan 1 NA 0.333 0.5641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1879 on 45 degrees of freedom
## Multiple R-Squared: 0.4294, Adjusted R-squared: 0.4041
## Wald test: 13.28 on 2 and 45 DF, p-value: 2.931e-05