Загружаем библиотеки
library(erer) # предельные эффекты, заодно подтянет пакеты ggplot2, lmtest
library(ggplot2) # графики
library(lmtest)
library(foreign)
library(xtable)
library(texreg)
Загружаем данные в R:
h <- read.dta("../datasets/mesto_jenshini.dta")
Это данные американского социологического опроса National Longitudinal Survey of Youth. Основная переменная, agree, показывает согласна ли опрашиваемая женщина с утверждением “Место женщины у плиты!” Мы оценим logit и probit модели для переменной agree.
Андрей Цветков, “Место женщины”
Проверяем, что всё корректно загрузилось — .dta это не родной для R формат.
head(h)
## obs id age age2 agree disagree nonint mw14 meduc adjinc nsibs fpro
## 1 1 1 20.33 413.4 0 1 0 0 8 13416 1 0
## 2 2 2 20.00 400.0 1 0 0 0 5 8944 8 0
## 3 3 3 17.42 303.3 0 1 1 1 10 10013 3 0
## 4 4 4 16.42 269.5 0 1 1 1 11 10013 3 0
## 5 5 8 20.50 420.2 0 1 0 1 9 4173 7 1
## 6 6 10 18.25 333.1 0 1 1 0 12 2048 3 0
## cath so urb tradrole
## 1 1 0 1 2
## 2 1 0 1 4
## 3 1 0 1 2
## 4 1 0 1 1
## 5 0 0 1 2
## 6 1 0 1 2
summary(h)
## obs id age age2
## Min. : 1 Min. : 1 Min. :14.1 Min. :198
## 1st Qu.: 927 1st Qu.: 2346 1st Qu.:16.4 1st Qu.:270
## Median :1853 Median : 4682 Median :18.5 Median :342
## Mean :1853 Mean : 5603 Mean :18.4 Mean :342
## 3rd Qu.:2779 3rd Qu.: 9243 3rd Qu.:20.3 3rd Qu.:413
## Max. :3705 Max. :12676 Max. :22.0 Max. :484
## agree disagree nonint mw14
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:1.000 1st Qu.:0.000 1st Qu.:0.000
## Median :0.000 Median :1.000 Median :0.000 Median :1.000
## Mean :0.148 Mean :0.852 Mean :0.251 Mean :0.503
## 3rd Qu.:0.000 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
## meduc adjinc nsibs fpro
## Min. : 0.0 Min. : 0 Min. : 0.00 Min. :0.00
## 1st Qu.:11.0 1st Qu.: 4800 1st Qu.: 2.00 1st Qu.:0.00
## Median :12.0 Median : 7354 Median : 3.00 Median :0.00
## Mean :11.6 Mean : 8646 Mean : 3.33 Mean :0.28
## 3rd Qu.:12.0 3rd Qu.:10667 3rd Qu.: 4.00 3rd Qu.:1.00
## Max. :20.0 Max. :49497 Max. :17.00 Max. :1.00
## cath so urb tradrole
## Min. :0.000 Min. :0.00 Min. :0.000 Min. :1.00
## 1st Qu.:0.000 1st Qu.:0.00 1st Qu.:1.000 1st Qu.:1.00
## Median :0.000 Median :0.00 Median :1.000 Median :2.00
## Mean :0.312 Mean :0.29 Mean :0.758 Mean :1.78
## 3rd Qu.:1.000 3rd Qu.:1.00 3rd Qu.:1.000 3rd Qu.:2.00
## Max. :1.000 Max. :1.00 Max. :1.000 Max. :4.00
Let's look on some specific variables. Income. Probably year income divided by family size.
summary(h$adjinc)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 4800 7350 8650 10700 49500
Несколько графиков до построения моделей:
qplot(factor(agree), data = h)
# more with ggplot2
Логит, пробит и мнк. Предпосылки применения мнк нарушены, но мы его всё равно применим и посмотрим, что выйдет.
m.logit <- glm(agree ~ age + adjinc + nsibs, data = h, x = TRUE, family = binomial(link = "logit"))
m.probit <- glm(agree ~ age + adjinc + nsibs, data = h, x = TRUE, family = binomial(link = "probit"))
m.ols <- glm(agree ~ age + adjinc + nsibs, data = h)
Посмотрим на логит модель:
summary(m.logit)
##
## Call:
## glm(formula = agree ~ age + adjinc + nsibs, family = binomial(link = "logit"),
## data = h, x = TRUE)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.894 -0.595 -0.547 -0.480 2.488
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.24e-01 3.93e-01 -1.84 0.0651 .
## age -5.03e-02 2.07e-02 -2.43 0.0150 *
## adjinc -3.78e-05 9.28e-06 -4.07 4.8e-05 ***
## nsibs 5.80e-02 2.04e-02 2.85 0.0044 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3105.3 on 3704 degrees of freedom
## Residual deviance: 3069.1 on 3701 degrees of freedom
## AIC: 3077
##
## Number of Fisher Scoring iterations: 4
Ковариационную матрицу оценок коэффициентов раздобыть легко
vcov(m.logit)
## (Intercept) age adjinc nsibs
## (Intercept) 1.542e-01 -7.732e-03 -7.622e-07 -1.532e-03
## age -7.732e-03 4.271e-04 -4.349e-10 -9.594e-06
## adjinc -7.622e-07 -4.349e-10 8.619e-11 2.889e-08
## nsibs -1.532e-03 -9.594e-06 2.889e-08 4.160e-04
Коэффициенты в логит и пробит моделях плохо интерпретируемы, т.к. скрытая переменная, в данной задаче склонность к ответу “да”, измеряется в непонятных единицах. В дополнение к коэффициентам рассчитывают предельные эффекты. Предельные эффекты отвечают на вопрос, как изменится примерно вероятность y=1 с ростом регрессора на единицу. Обычно предельные эффекты считаются для среднего значения каждого регрессора:
maBina(m.logit)$out
## effect error t.value p.value
## (Intercept) -0.089 0.048 -1.844 0.065
## age -0.006 0.003 -2.439 0.015
## adjinc 0.000 0.000 -4.119 0.000
## nsibs 0.007 0.003 2.854 0.004
Опция x=TRUE в функции glm нужна только для того, чтобы работала функция maBina. Она сохраняет исходные x в полученной регрессии. Если массив данных огромный, то эту опцию можно убрать, но тогда придется считать предельные эффекты ручками.
Создаем набор данных, для которого мы будем строить прогнозы:
new <- data.frame(age = c(20, 24), adjinc = c(16000, 4000), nsibs = c(2, 5))
new
## age adjinc nsibs
## 1 20 16000 2
## 2 24 4000 5
В логит и пробит-моделях можно прогнозировать значение скрытой переменной (склонности ответить “да”) или вероятность \( y_i=1 \) (вероятность ответить “да”). Прогноз скрытой переменной:
predict(m.logit, new)
## 1 2
## -2.217 -1.791
Немножко другой командой получаются прогнозы вероятности \( y_i=1 \):
predict(m.logit, new, type = "response")
## 1 2
## 0.09822 0.14295
Поскольку в методе максимального правдоподобия оценка любого параметра является асимпотически нормальной, то доверительный интервал для любого параметра \( a \) строится по принципу \( [\hat{a}-z_{cr}se(\hat{a});\hat{a}+z_{cr}se(\hat{a})] \). Для того, чтобы получить стандартные ошибки, как для оценки ожидаемой склонности ответить “да”, так и для оценки вероятности ответить “да”, используется опция se.fit=TRUE.
Поскольку чаще всего нужен доверительный интервал именно для вероятностей, мы построим только его. Если нужен доверительный интервал для ожидаемой склонности ответить да, то он строится аналогично:
forecast <- predict(m.logit, new, type = "response", se.fit = TRUE)
new$hat.p <- forecast$fit
new$se.hat.p <- forecast$se.fit
z.cr <- qnorm(0.975)
new$ci.left <- new$hat.p - z.cr * new$se.hat.p
new$ci.right <- new$hat.p + z.cr * new$se.hat.p
new
## age adjinc nsibs hat.p se.hat.p ci.left ci.right
## 1 20 16000 2 0.09822 0.008747 0.08107 0.1154
## 2 24 4000 5 0.14295 0.016522 0.11056 0.1753
Заметим приятный факт. Для средних значений регрессоров обычный МНК, который формально нельзя применять из-за нарушения предпосылок теоремы Гаусса-Маркова, даёт вполне корректные прогнозы:
Проблемы МНК видны на краевых значениях регрессоров. Тут запросто могут быть отрицательные вероятности :)
Binary prediction. Best binary prediction?
Likelihood ratio, LR-test, Тест отношения правдоподобия
m2.logit <- glm(agree ~ age + age2 + adjinc + nsibs, data = h, family = binomial(link = "logit"))
lrtest(m.logit, m2.logit)
## Likelihood ratio test
##
## Model 1: agree ~ age + adjinc + nsibs
## Model 2: agree ~ age + age2 + adjinc + nsibs
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -1535
## 2 5 -1535 1 0.02 0.88
Lagrange multiplier test, LM-test, тест множителей Лагранжа (? by hands)
Wald-test, Тест Вальда
waldtest(m.logit, m2.logit, test = "Chisq")
## Wald test
##
## Model 1: agree ~ age + adjinc + nsibs
## Model 2: agree ~ age + age2 + adjinc + nsibs
## Res.Df Df Chisq Pr(>Chisq)
## 1 3701
## 2 3700 1 0.02 0.88
ROC ручками…
ROC с помощью пакета ROCR
Одна модель:
xtable(m.logit)
Предельные эффекты одной модели:
xtable(maBina(m.logit)$out)
Табличка со сравнением нескольких моделей:
texreg(list(m.logit, m2.logit))
Табличка со сравнением предельных эффектов