Обобщением множественной линейной регрессии можно считать общую линейную модель (general linear model). Общая линейная модель может предсказывать не одну, а сразу несколько объясняемых переменных в отличие от множественной линейной регрессии.
Почти все пройденные нами методы можно рассматривать как частный случай общей линейной модели: t-тесты, коэффициент корреляции Пирсона, линейная регрессия, ANOVA.
23.2 Обобщенная линейная модель
Обобщенная линейная модель (generalized linear model) была придумана как обобщение линейной регрессии и ее сородичей: логистической регрессии и пуассоновской регрессии.
Общая линейная модель задается формулой \[Y = XB\]
Обобщенная оборачивает предиктор \(XB\)связывающей функцией (link function), которая различается для разных типов регрессионных моделей.
Давайте попробуем построить модель, в которой объясняемой переменной будет то, является ли супергерой хорошим или плохим.
library(tidyverse)
Warning: package 'dplyr' was built under R version 4.2.3
Warning: package 'stringr' was built under R version 4.2.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.4.4 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.0
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
dat <- vroom(...)
problems(dat)
Rows: 734 Columns: 11
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (8): name, Gender, Eye color, Race, Hair color, Publisher, Skin color, A...
dbl (3): ...1, Height, Weight
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
heroes$good <- heroes$Alignment =="good"
Обычная линейная модель нам не подходит, если распределение наших ошибок далеко от нормального. А это значит, что мы не можем использовать общую линейную модель с бинарной объясняемой переменной. Эту проблему решает логистическая регрессия, которая является частным случаем обобщенной линейной модели.
Для этого нам понадобится функция glm(), а не lm() как раньше. Ее синтаксис очень похож, но нам теперь нужно задать еще один важный параметр family = для выбора связывающей функции (в данном случае это логит-функция, которая является связующей функцией по умолчанию для биномиального семейства функций в glm()).
heroes_good_glm <-glm(good ~ Weight + Gender, heroes, family =binomial()) summary(heroes_good_glm)
Call:
glm(formula = good ~ Weight + Gender, family = binomial(), data = heroes)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8674 -1.3103 0.6334 0.9155 2.4007
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
(256 observations deleted due to missingness)
AIC: 576.88
Number of Fisher Scoring iterations: 4
Результат очень похож по своей структуре на glm(), однако вместо \(R^2\) перед нами AIC. AIC расшифровывается как информационный критерий Акаике (Akaike information criterion) — это критерий использующийся для выбора из нескольких моделей. Чем он меньше, тем лучше модель. Как и Adjusted R2, AIC “наказывает” за большое количество параметров в модели.
Поскольку AIC — это относительный показатель качества модели, нам нужно сравнить его с AIC другой, более общей модели. Можно сравнить с моделью без веса супергероев.
heroes_good_glm_noweight <-glm(good ~ Gender, heroes, family =binomial()) summary(heroes_good_glm_noweight)
Call:
glm(formula = good ~ Gender, family = binomial(), data = heroes)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8082 -1.4164 0.6586 0.9559 0.9559
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.4178 0.1785 7.944 1.95e-15 ***
GenderMale -0.8716 0.2012 -4.332 1.48e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 873.81 on 698 degrees of freedom
Residual deviance: 853.24 on 697 degrees of freedom
(35 observations deleted due to missingness)
AIC: 857.24
Number of Fisher Scoring iterations: 4
AIC стал больше, следовательно, мы выберем модель с весом супергероев.
23.3 Модель со смешанными эффектами
Модели со смешанными эффектами (mixed-effects models) — это то же самое, что и иерархическая регрессия (hierarchical regression) или многоуровневое моделирование (multilevel modelling). Этому методу повезло иметь много названий - в зависимости от области, в которой он используется. Модели со смешанными эффектами позволяет включать в линейную регрессию не только фиксированные эффекты (fixed effects), но и случайные эффекты (random effects).
Для экспериментальных дисциплин это интересно тем, что в таких моделях можно не усреднять показатели по испытуемым или образцам, а учитывать влияние соотвествующей группирующей переменной как случайный эффект. В отличие от обычного фактора как в линейной регрессии или дисперсионном анализе (здесь он называется фиксированным), случайный эффект не интересует нас сам по себе, а его значения считаются случайной переменной.
Смешанные модели используются в самых разных областях. Они позволяют решить проблему зависимости наблюдений без усреднения значений по испытуемым или группам, что повышает статистическую мощность.
Для работы со смешанными моделями в R есть два известных пакета: nlme и более современный lme4.
install.packages("lme4")
library(lme4)
Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Для примера возьмем данные исследования влияния депривации сна на время реакции.
data("sleepstudy")
Данные представлены в длинном формате: каждая строчка — это усредненное время реакции для одного испытуемого в соответствующий день эксперимента.
Можно заметить, что, в среднем, время реакции у испытуемых повышается от первого к последнему дню. С помощью смешанных моделей мы можем проверить, различается ли скорость возрастания времени реакции от дня к дню у разных испытуемых.
Для этого мы сравниваем две модели, одна из которых является “вложенной” в другую, то есть усложненной версией более общей модели. В данном случае, более общая модель предполагает, что время реакции увеличивается у всех испытуемых одинаково, а испытуемые различаются только средним временем реакции.
sleep_lme0 <-lmer(Reaction ~ Days + (1| Subject), sleepstudy)sleep_lme1 <-lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
Зеленая линия (нулевая модель) имеет везде один и тот же наклон, а синяя (альтернативная модель) имеет разный наклон у всех испытуемых.
Есть несколько способов сравнивать модели, например, уже знакомый нам AIC. Кроме того, можно сравнить две модели с помощью теста хи-квадрат, восполльзовавшись функцией anova().
