Загружаем библиотеки

library("mfx") # подсчет предельных эффектов
library("ggplot2") # графики
library("lmtest") # линейные регрессии 
library("foreign") # чтение файлов в некоторых форматах
library("dplyr") # манипуляции с таблицами
library("broom") # описание модели в виде таблички

Загружаем данные в R:

role <- read.dta("../datasets/mesto_jenshini.dta")

Это данные американского социологического опроса National Longitudinal Survey of Youth. Основная переменная, agree, показывает согласна ли опрашиваемая женщина с утверждением “Место женщины у плиты!” Мы оценим logit и probit модели для переменной agree.

Tsvetkov, Mesto jenshini

Андрей Цветков, “Место женщины”

До оценки моделей

Проверяем, что всё корректно загрузилось — .dta это не родной для R формат.

glimpse(role)
## Observations: 3,705
## Variables: 16
## $ obs      (int) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16...
## $ id       (int) 1, 2, 3, 4, 8, 10, 12, 14, 16, 19, 20, 21, 22, 25, 27...
## $ age      (dbl) 20.33333, 20.00000, 17.41667, 16.41667, 20.50000, 18....
## $ age2     (dbl) 413.4445, 400.0000, 303.3403, 269.5070, 420.2500, 333...
## $ agree    (int) 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,...
## $ disagree (int) 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,...
## $ nonint   (int) 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ mw14     (int) 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,...
## $ meduc    (int) 8, 5, 10, 11, 9, 12, 15, 12, 12, 10, 12, 12, 12, 14, ...
## $ adjinc   (int) 13416, 8944, 10013, 10013, 4173, 2048, 14434, 6577, 2...
## $ nsibs    (int) 1, 8, 3, 3, 7, 3, 3, 2, 3, 3, 2, 2, 1, 1, 1, 1, 1, 1,...
## $ fpro     (int) 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ cath     (int) 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,...
## $ so       (int) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ urb      (int) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ tradrole (int) 2, 4, 2, 1, 2, 2, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1, 3, 1,...

Смотрим на описательные статистики:

tidy(role)[c("column", "n", "mean", "median", "min", "max", "sd")]
##      column    n         mean  median       min   max           sd
## 1       obs 3705 1853.0000000 1853.00   1.00000  3705 1069.6857015
## 2        id 3705 5603.1284750 4682.00   1.00000 12676 3881.4924983
## 3       age 3705   18.3615385   18.50  14.08333    22    2.2399992
## 4      age2 3705  342.1623384  342.25 198.34027   484   81.5823005
## 5     agree 3705    0.1479082    0.00   0.00000     1    0.3550569
## 6  disagree 3705    0.8520918    1.00   0.00000     1    0.3550569
## 7    nonint 3705    0.2504723    0.00   0.00000     1    0.4333435
## 8      mw14 3705    0.5033738    1.00   0.00000     1    0.5000561
## 9     meduc 3705   11.6129555   12.00   0.00000    20    2.5282944
## 10   adjinc 3705 8646.2933873 7354.00   0.00000 49497 5910.1002907
## 11    nsibs 3705    3.3284750    3.00   0.00000    17    2.1649853
## 12     fpro 3705    0.2801619    0.00   0.00000     1    0.4491388
## 13     cath 3705    0.3122807    0.00   0.00000     1    0.4634862
## 14       so 3705    0.2898785    0.00   0.00000     1    0.4537671
## 15      urb 3705    0.7581646    1.00   0.00000     1    0.4282529
## 16 tradrole 3705    1.7832659    2.00   1.00000     4    0.7792120

Несколько графиков до построения моделей:

ggplot(role) + geom_bar(aes(factor(agree))) + scale_x_discrete(breaks = c("1", "0"), 
  labels = c("Да", "Нет")) + labs(x = "Согласны ли с утверждением?", y = "Количество женщин", 
  title = "Результаты опроса")

Собственно оценка моделей

Логит, пробит и мнк. Предпосылки применения мнк нарушены, но мы его всё равно применим и посмотрим, что выйдет.

m_logit <- glm(agree ~ age + adjinc + nsibs,
          data = role,
          family = binomial(link = "logit"))
m_probit <- glm(agree ~ age + adjinc + nsibs,
          data = role,
          family = binomial(link = "probit"))
m_ols <- glm(agree ~ age + adjinc + nsibs,
          data = role)

Кратко о логит-модели:

glance(m_logit)
##   null.deviance df.null    logLik      AIC      BIC deviance df.residual
## 1       3105.26    3704 -1534.526 3077.052 3101.922 3069.052        3701

Оценки коэффициентов:

tidy(m_logit)
##          term      estimate    std.error statistic      p.value
## 1 (Intercept) -7.241760e-01 3.926483e-01 -1.844338 6.513396e-02
## 2         age -5.025358e-02 2.066656e-02 -2.431638 1.503071e-02
## 3      adjinc -3.775219e-05 9.283632e-06 -4.066532 4.771795e-05
## 4       nsibs  5.804622e-02 2.039539e-02  2.846047 4.426568e-03

Ковариационную матрицу оценок коэффициентов также раздобыть легко

vcov(m_logit)
##               (Intercept)           age        adjinc         nsibs
## (Intercept)  1.541727e-01 -7.731805e-03 -7.622419e-07 -1.532192e-03
## age         -7.731805e-03  4.271065e-04 -4.348763e-10 -9.594252e-06
## adjinc      -7.622419e-07 -4.348763e-10  8.618583e-11  2.889028e-08
## nsibs       -1.532192e-03 -9.594252e-06  2.889028e-08  4.159717e-04

Коэффициенты в логит и пробит моделях плохо интерпретируемы, т.к. скрытая переменная, в данной задаче склонность к ответу “да”, измеряется в непонятных единицах. В дополнение к коэффициентам рассчитывают предельные эффекты. Предельные эффекты отвечают на вопрос, как изменится примерно вероятность y = 1 с ростом регрессора на единицу. Обычно предельные эффекты считаются для среднего значения каждого регрессора:

logitmfx(agree ~ age + adjinc + nsibs, data = role)
## Call:
## logitmfx(formula = agree ~ age + adjinc + nsibs, data = role)
## 
## Marginal Effects:
##              dF/dx   Std. Err.       z     P>|z|    
## age    -6.2048e-03  2.5443e-03 -2.4387   0.01474 *  
## adjinc -4.6612e-06  1.1316e-06 -4.1193 3.801e-05 ***
## nsibs   7.1669e-03  2.5113e-03  2.8538   0.00432 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Сообщения об ошибках

  • fitted probabilities numerically 0 or 1 occurred. Означает, что при итерациях при поиске максимума правдоподобия для какого-то наблюдения оцененные вероятности были настолько близки к нулю или единице, что у компьютера не хватает ячейки памяти, чтобы их отличить от нуля или единицы. Само по себе это нестрашно.
  • algorithm did not converge. Можно попробовать увеличить количество итераций, иногда это помогает. В сочетании с предыдущим сообщением может быть показателем того, что 0 и 1 идеально разделяются какой-то переменной и возможно «идеальное» прогнозирование.

Прогнозы

Создаем набор данных, для которого мы будем строить прогнозы:

new <- data_frame(age = c(20, 24),
                  adjinc = c(16000, 4000),
                  nsibs = c(2, 5))
new
## Source: local data frame [2 x 3]
## 
##     age adjinc nsibs
##   (dbl)  (dbl) (dbl)
## 1    20  16000     2
## 2    24   4000     5

В логит и пробит-моделях можно прогнозировать значение скрытой переменной (склонности ответить “да”) или вероятность \(y_i=1\) (вероятность ответить “да”). Прогноз скрытой переменной:

augment(m_logit, newdata = new)
##   age adjinc nsibs  .fitted    .se.fit
## 1  20  16000     2 -2.21719 0.09875815
## 2  24   4000     5 -1.79104 0.13486036

Немножко другой командой получаются прогнозы вероятности \(y_i=1\):

augment(m_logit, newdata = new, type.predict = "response")
##   age adjinc nsibs    .fitted     .se.fit
## 1  20  16000     2 0.09821738 0.008747081
## 2  24   4000     5 0.14294530 0.016522004

Заметим приятный факт. Для средних значений регрессоров обычный МНК, который формально нельзя применять из-за нарушения предпосылок теоремы Гаусса-Маркова, даёт вполне корректные прогнозы:

augment(m_ols, newdata = new)
##   age adjinc nsibs    .fitted    .se.fit
## 1  20  16000     2 0.09714633 0.01053247
## 2  24   4000     5 0.14418222 0.01675517

Проблемы МНК видны на краевых значениях регрессоров. Тут запросто могут быть отрицательные вероятности :)

augment(m_ols, newdata = data_frame(age = 100, adjinc = 10^5, nsibs = 0))
##   age adjinc nsibs    .fitted   .se.fit
## 1 100  1e+05     0 -0.7668123 0.2299096

Поскольку в методе максимального правдоподобия оценка любого параметра является асимпотически нормальной, то доверительный интервал для любого параметра \(a\) строится по принципу \([\hat{a}-z_{cr}se(\hat{a});\hat{a}+z_{cr}se(\hat{a})]\).

Поскольку чаще всего нужен доверительный интервал именно для вероятностей, мы построим только его. Если нужен доверительный интервал для ожидаемой склонности ответить да, то достаточно просто убрать опцию prediction.type:

new2 <- augment(m_logit, newdata = new, type.predict = "response")
z_cr <- qnorm(0.975)
new2 <- mutate(new2, 
               ci_left = .fitted - z_cr * .se.fit,
               ci_right = .fitted + z_cr * .se.fit)
new2
##   age adjinc nsibs    .fitted     .se.fit    ci_left  ci_right
## 1  20  16000     2 0.09821738 0.008747081 0.08107342 0.1153613
## 2  24   4000     5 0.14294530 0.016522004 0.11056277 0.1753278

Binary prediction. Best binary prediction?

Выбор между вложенными моделями

Likelihood ratio, LR-test, Тест отношения правдоподобия

m2_logit <- glm(agree ~ age + age2 + adjinc + nsibs,
          data = role,
          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 -1534.5                     
## 2   5 -1534.5  1 0.0225     0.8808

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.0225     0.8808

Кривые ROC

ROC ручками…

ROC с помощью пакета ROCR