Лекция 4 - практическая часть

R
Автор

Е. Тымченко

Дата публикации

2 июля 2023 г.

Библиотеки

library(tidyverse)
library(hrbrthemes)
library(stargazer)
library(reshape2)

Данные

Сегодня мы построим модель, которая позволит диагностировать сахарный диабет, на основании различных данных о пациенте.

dta <- read_csv('https://raw.githubusercontent.com/ETymch/Econometrics_2023/main/Datasets/framingham.csv')

test <- sample_n(dta, 800) # тестовая выборка. На ней мы будем проверять модели.
train <- setdiff(dta, test) # на этой выборке мы будем оценивать модели

Почему в данном случае МНК - не лучший вариант для оценки параметров?

Наиболее значимым фактором для постановки такого диагноза является уровень глюкозы в крови. Построим график зависимости диагноза от кровня глюкозы:

dta %>%
  ggplot(aes(x = glucose, y = diabetes)) +
  geom_point() +
  geom_smooth(method = 'lm') + # добавим линию, оценённую с помощью МНК
  geom_vline(xintercept = 99, color = 'red') + # Верхняя граница нормы
  theme_ipsum() +
  labs(x = 'Уровень глюкозы в крови, мг./Дл.',
       y = 'Сахарный диабет')

Линия плохо описывает взаимосвязь. Возможно, такая форма лучше?

dta %>%
  ggplot(aes(x = glucose, y = diabetes)) +
  geom_point() +
  geom_smooth(method = "glm", 
              method.args = list(family = "binomial"), 
              se = FALSE) +
  geom_vline(xintercept = 99, color = 'red') +
  theme_ipsum() +
  labs(x = 'Уровень глюкозы в крови, мг./Дл.',
       y = 'Сахарный диабет')

Оценим модель с помощью МНК.

mod_ols <- lm(diabetes ~ glucose, train)
stargazer(mod_ols, type = 'html', header = F)
Dependent variable:
diabetes
glucose 0.004***
(0.0001)
Constant -0.323***
(0.008)
Observations 3,116
R2 0.375
Adjusted R2 0.375
Residual Std. Error 0.130 (df = 3114)
F Statistic 1,870.258*** (df = 1; 3114)
Note: p<0.1; p<0.05; p<0.01

Но как её интерпретировать? Что значит отрицательный intercept? Если человек не сдавал анализ на глюкозу, его заболеваемость отрицательная? Также, вспомните Теорему Гауccа-Маркова, в которой сказано, что если выполнены несколько условий, одним из которых является линейность связи между зависимой и независимой переменной, то оценка МНК - лучшая несмещённая оценка. В нашем случае связь, очевидно, нелинейная, а оценки МНК - не лучшие. А значит нам нужен иной метод - Метод Максимального Правдоподобия.

Для оценки таких моделей у нас есть более удачная функция, чем обыкновенная прямая. Эта функция - сигмоид.

Логистическая функция

Напишем уравнение для линии.

y <- function(x){
  intercept + b*x
}

Уравнение логистической функции.

S <- function(x){
  1 / ( 1 + exp(-y(x)))
}

Параметры и интервал

intercept = 1
b = 1
x <- seq(-5, 5, by = 0.1)

Теперь нарисуем их:

tibble(x,
       line = y(x),
       sigmoid = S(x)) %>%
  melt(id = 'x') %>%
  ggplot(aes(x = x, y = value, color = variable)) +
  geom_line(lwd = 1, alpha = 0.6) +
  ylim(0, 1) +
  theme_ipsum() +
  labs(title = 'Сигмоид и линия', x = 'X', y = 'Диабет - 1, Здоров - 0')

Чтобы получить больше интеиции о том, как устроена логистическая функция, Сравним две модели:

  • МНК: \(y_i = intercept + \beta_1 x_i + \epsilon_i\)
  • Логистическая регрессия: \(y_i = \frac{1}{1 + e^{-(intercept + \beta_1 x_i)}} + \epsilon_i\)

Зафиксируем \(\beta_1\) и посмотрим, что произойдёт при изменении intercept.

Теперь зафиксируем intercept и посмотрим, что будет при изменении \(\beta_1\):

Оценим две модели:

mod_ml <- glm(diabetes ~ glucose, family = 'binomial', train) # МНК
mod_ols <- lm(diabetes ~ glucose, train) # ММП

stargazer(mod_ols, mod_ml, type = 'html', header = F) # сравним модели
Dependent variable:
diabetes
OLS logistic
(1) (2)
glucose 0.004*** 0.083***
(0.0001) (0.007)
Constant -0.323*** -11.732***
(0.008) (0.720)
Observations 3,116 3,116
R2 0.375
Adjusted R2 0.375
Log Likelihood -185.367
Akaike Inf. Crit. 374.735
Residual Std. Error 0.130 (df = 3114)
F Statistic 1,870.258*** (df = 1; 3114)
Note: p<0.1; p<0.05; p<0.01

Ставим диагнозы при помощи модели

S <- function(x){
  1 / (1 + exp(- (coefficients(mod_ml)[1] + coefficients(mod_ml)[2] * x)))
}

y <- function(x){
  coefficients(mod_ols)[1] + coefficients(mod_ols)[2]*x
}

gl_test <- test$glucose # Столбец с значением анализа для тестовой выборки.
dia_test <- test$diabetes # Столбец с истинным значением диагноза.

# Создаём табличку

result <- tibble(gl_test, 
       dia_test,
       pred_ols = y(gl_test),
       pred_ml = S(gl_test)) %>%
  mutate(pred_ml_01 = ifelse(pred_ml < 0.5, 0, 1)) # создайм ещё одну переменную со значением прогноза на основе предсказанной вероятности.

result %>%
  head(10) # посмотрим на вервые 10 строк в таблице.
# A tibble: 10 × 5
   gl_test dia_test pred_ols pred_ml pred_ml_01
     <dbl>    <dbl>    <dbl>   <dbl>      <dbl>
 1     166        1  0.387   0.892            1
 2      79        0  0.0152  0.00581          0
 3     120        1  0.191   0.151            0
 4      75        0 -0.00195 0.00417          0
 5      72        0 -0.0148  0.00325          0
 6      78        0  0.0109  0.00534          0
 7      68        0 -0.0319  0.00233          0
 8      67        0 -0.0362  0.00214          0
 9      80        0  0.0194  0.00631          0
10      74        0 -0.00622 0.00383          0

В отличие от оценок, сделанных при помощи МНК, оценки ММП смещённые. Это значит, что средняя ошибка в модели не равна 0.

residuals(mod_ols) %>% mean() # Ошибки в модели МНК.
[1] -3.123725e-17
residuals(mod_ml) %>% mean() # Ошибки в модели ММП.
[1] -0.09042032

Для большей наглядности, нарисуем распределение ошибок в моделях:

tibble(id = 1:length(residuals(mod_ols)),
                redid_ols = residuals(mod_ols),
                resid_ml = residuals(mod_ml)) %>%
  melt(id = 'id') %>%
  ggplot(aes(x = value, color = variable, fill = variable)) +
  geom_density(alpha = 0.5) +
  theme_minimal()+
  xlim(-1, 0.5)

Добавим в модель больше параметров и сравним с первой:

mod_ml_1 <- glm(diabetes ~ glucose + age + education + cigsPerDay + totChol, family = 'binomial', train) # Диагноз теперь зависит от уровня глюкозы, возраста, уровня образования, интенсивности курения и уровня холестерина.

stargazer(mod_ml, mod_ml_1, type = 'html', header = F)
Dependent variable:
diabetes
(1) (2)
glucose 0.083*** 0.082***
(0.007) (0.007)
age 0.048**
(0.021)
education -0.055
(0.163)
cigsPerDay -0.008
(0.016)
totChol 0.003
(0.003)
Constant -11.732*** -14.682***
(0.720) (1.523)
Observations 3,116 3,018
Log Likelihood -185.367 -174.440
Akaike Inf. Crit. 374.735 360.880
Note: p<0.1; p<0.05; p<0.01

Сравним, какая модель стваит диагнозы: с однйо объясняющей переменной или с несколькими:

S <- function(x){
  1 / (1 + exp(-x))
}

predict_ml <- predict(mod_ml, test %>% select(-diabetes)) %>%
  S()
predict_ml_1 <- predict(mod_ml_1, test %>% select(-diabetes)) %>%
  S()

predictions <- tibble(true_diagnoses = test$diabetes,
       prob_ml = predict_ml,
       prob_ml_1 = predict_ml_1) %>%
  mutate(dia_ml = ifelse(prob_ml < 0.5, 0, 1),
         dia_ml_1 = ifelse(prob_ml_1 < 0.5, 0, 1)
         ) %>%
  mutate(error_ml = ifelse(dia_ml == true_diagnoses, 0, 1),
         error_ml_1 = ifelse(dia_ml == true_diagnoses, 0, 1)
         )

predictions %>%
  tail(10) # последние 10 наблюдений в табличке
# A tibble: 10 × 7
   true_diagnoses  prob_ml prob_ml_1 dia_ml dia_ml_1 error_ml error_ml_1
            <dbl>    <dbl>     <dbl>  <dbl>    <dbl>    <dbl>      <dbl>
 1              0 NA        NA           NA       NA       NA         NA
 2              0  0.0156    0.0249       0        0        0          0
 3              0  0.00353   0.00360      0        0        0          0
 4              0  0.00275   0.00159      0        0        0          0
 5              0 NA        NA           NA       NA       NA         NA
 6              0  0.00453   0.00195      0        0        0          0
 7              0  0.00181   0.00105      0        0        0          0
 8              0  0.0170    0.00856      0        0        0          0
 9              0  0.00214   0.00105      0        0        0          0
10              0  0.00534   0.00526      0        0        0          0
predictions$error_ml %>% na.omit() %>% sum
[1] 9
predictions$error_ml_1 %>% na.omit() %>% sum
[1] 9

Удивительно, но обе модели ошиблись в диагнозах менее 20 раз из 800 - менее 5%! Модель с одной объясняющей переменной показала отличную способность прогнозировать диагнозы, потому что мы выбрали наиболее значимую объясняющую переменную - уровень инсулина.