library(tidyverse)
heroes <- read_csv("https://raw.githubusercontent.com/Pozdniakov/tidy_stats/master/data/heroes_information.csv",
na = c("-", "-99", "NA"))
heroes$good <- heroes$Alignment == "good"23 Общая линейная модель и ее расширения
23.1 Общая линейная модель
Обобщением множественной линейной регрессии можно считать общую линейную модель (general linear model). Она может предсказывать не одну, а сразу несколько объясняемых переменных в отличие от множественной линейной регрессии.
\[Y = XB + E,\] где \(Y\) — матрица объясняемых переменных, \(X\) — матрица предикторов, \(B\) — матрица параметров, \(E\) — матрица ошибок.
Почти все пройденные нами методы можно рассматривать как частный случай общей линейной модели: \(t\)-тесты, коэффициент корреляции Пирсона, линейная регрессия, ANOVA.
23.2 Обобщенная линейная модель
Обобщенная линейная модель (generalized linear model) была придумана как обобщение линейной регрессии и ее сородичей: логистической регрессии и пуассоновской регрессии.
Общая линейная модель задается формулой \[Y = XB + E\]
Обобщенная оборачивает предиктор \(XB\) связывающей функцией (link function), которая различается для разных типов регрессионных моделей.
Давайте попробуем построить модель, в которой объясняемой переменной будет то, является ли супергерой хорошим или плохим.
Обычная линейная модель в случае с бинарной объясняемой переменной работает плохо: ее предсказания не ограничены диапазоном от 0 до 1 (а именно так — как вероятность — хочется трактовать предсказание для бинарного исхода), да и ошибки такой модели заведомо далеки от нормальных. Эту проблему решает логистическая регрессия, которая является частным случаем обобщенной линейной модели.
Но прежде уберем строки с пропусками в интересующих нас переменных: скоро мы будем сравнивать несколько моделей, а сравнивать их можно, только если они построены на одних и тех же данных.
heroes_good <- heroes %>%
drop_na(good, Weight, Gender)Для этого нам понадобится функция glm(), а не lm() как раньше. Ее синтаксис очень похож, но нам теперь нужно задать еще один важный параметр family = для выбора связывающей функции (в данном случае это логит-функция, которая является связывающей функцией по умолчанию для биномиального семейства функций в glm()).
heroes_good_glm <- glm(good ~ Weight + Gender, heroes_good, family = binomial())
summary(heroes_good_glm)
Call:
glm(formula = good ~ Weight + Gender, family = binomial(), data = heroes_good)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.763917 0.235410 7.493 6.73e-14 ***
Weight -0.004253 0.001124 -3.783 0.000155 ***
GenderMale -0.760310 0.245851 -3.093 0.001984 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 605.22 on 477 degrees of freedom
Residual deviance: 570.88 on 475 degrees of freedom
AIC: 576.88
Number of Fisher Scoring iterations: 4Результат очень похож по структуре на вывод lm(), однако вместо \(R^2\) перед нами AIC. AIC расшифровывается как информационный критерий Акаике (Akaike information criterion) — это критерий, использующийся для выбора из нескольких моделей. Чем он меньше, тем лучше модель. Как и adjusted \(R^2\), AIC «наказывает» за большое количество параметров в модели.
Поскольку AIC — это относительный показатель качества модели, нам нужно сравнить его с AIC другой модели, построенной на тех же данных, — например, более простой, без массы супергероев.
heroes_good_glm_noweight <- glm(good ~ Gender, heroes_good, family = binomial())
summary(heroes_good_glm_noweight)
Call:
glm(formula = good ~ Gender, family = binomial(), data = heroes_good)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.4040 0.2109 6.657 2.80e-11 ***
GenderMale -0.9311 0.2389 -3.898 9.72e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 605.22 on 477 degrees of freedom
Residual deviance: 588.52 on 476 degrees of freedom
AIC: 592.52
Number of Fisher Scoring iterations: 4AIC стал больше, следовательно, мы выберем модель с массой супергероев.
23.3 Модель со смешанными эффектами
Модели со смешанными эффектами (mixed-effects models) — это то же самое, что и иерархическая регрессия (hierarchical regression) или многоуровневое моделирование (multilevel modelling). Этому методу повезло иметь много названий — в зависимости от области, в которой он используется. Модели со смешанными эффектами позволяют включать в линейную регрессию не только фиксированные эффекты (fixed effects), но и случайные эффекты (random effects).
Для экспериментальных дисциплин это интересно тем, что в таких моделях можно не усреднять показатели по испытуемым или образцам, а учитывать влияние соответствующей группирующей переменной как случайный эффект. В отличие от обычного фактора как в линейной регрессии или дисперсионном анализе (здесь он называется фиксированным), случайный эффект не интересует нас сам по себе, а его значения считаются случайной переменной.
Смешанные модели используются в самых разных областях. Они позволяют решить проблему зависимости наблюдений без усреднения значений по испытуемым или группам, что повышает статистическую мощность.
Для работы со смешанными моделями в R есть два известных пакета: {nlme} и более современный {lme4}.
install.packages("lme4")library(lme4)Для примера возьмем данные исследования влияния депривации сна на время реакции.
data("sleepstudy")Данные представлены в длинном формате: каждая строчка — это усредненное время реакции для одного испытуемого в соответствующий день эксперимента.
sleepstudy %>%
ggplot(aes(x = Days, y = Reaction)) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = 0:9) +
facet_wrap(~Subject) +
theme_minimal()
Можно заметить, что, в среднем, время реакции у испытуемых повышается от первого к последнему дню. С помощью смешанных моделей мы можем проверить, различается ли скорость возрастания времени реакции от дня к дню у разных испытуемых.
Для этого мы сравниваем две вложенные модели: более простая получается из более сложной, если убрать часть параметров. Более простая модель предполагает, что время реакции увеличивается у всех испытуемых одинаково, а испытуемые различаются только средним временем реакции. Более сложная модель дополнительно разрешает, чтобы скорость увеличения (наклон по дням) различалась между испытуемыми.
sleep_lme0 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
sleep_lme1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)Визуализируем предсказания двух моделей:
sleepstudy$predicted_by_sleep_lme0 <- predict(sleep_lme0)
sleepstudy$predicted_by_sleep_lme1 <- predict(sleep_lme1)sleepstudy %>%
rename(observed_reaction_time = Reaction) %>%
pivot_longer(cols = c(observed_reaction_time, predicted_by_sleep_lme0, predicted_by_sleep_lme1), names_to = "model", values_to = "Reaction") %>%
ggplot(aes(x = Days, y = Reaction)) +
geom_line(aes(colour = model)) +
geom_point(data = sleepstudy, alpha = 0.4) +
scale_x_continuous(breaks = 0:9) +
facet_wrap(~Subject) +
theme_minimal()
Линия более простой модели (sleep_lme0) имеет везде один и тот же наклон, а линия более сложной (sleep_lme1) — разный наклон у каждого испытуемого.
Есть несколько способов сравнивать модели, например, уже знакомый нам AIC. Кроме того, можно сравнить две модели с помощью теста хи-квадрат, воспользовавшись функцией anova().
anova(sleep_lme0, sleep_lme1)refitting model(s) with ML (instead of REML)Data: sleepstudy
Models:
sleep_lme0: Reaction ~ Days + (1 | Subject)
sleep_lme1: Reaction ~ Days + (Days | Subject)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
sleep_lme0 4 1802.1 1814.8 -897.04 1794.1
sleep_lme1 6 1763.9 1783.1 -875.97 1751.9 42.139 2 7.072e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Модель со случайным наклоном прямой оказалась лучше, о чем нам говорят как более низкие AIC и BIC1, так и тестирование с помощью хи-квадрат.
BIC — это байесовский информационный критерий (Bayesian information criterion), родственник AIC, но с немного отличающейся формулой.↩︎