360-glm.qmd

# Общая линейная модель и ее расширения {#sec-glm_general}

## Общая линейная модель {#sec-general_linear_model}

Обобщением множественной линейной регрессии можно считать **общую линейную модель (general linear model)**. Общая линейная модель может предсказывать не одну, а сразу несколько объясняемых переменных в отличие от множественной линейной регрессии.

$$Y = XB$$ где $Y$ --- матрица объясняемых переменных, $X$ --- матрица предикторов, $B$ --- матрица параметров.

Почти все пройденные нами методы можно рассматривать как частный случай общей линейной модели: t-тесты, коэффициент корреляции Пирсона, линейная регрессия, ANOVA.

![](https://lindeloev.github.io/tests-as-linear/linear_tests_cheat_sheet.png)

## Обобщенная линейная модель {#sec-glm_model}

**Обобщенная линейная модель (generalized linear model)** была придумана как обобщение линейной регрессии и ее сородичей: логистической регрессии и пуассоновской регрессии.

Общая линейная модель задается формулой $$Y = XB$$

Обобщенная оборачивает предиктор $XB$ **связывающей функцией (link function),** которая различается для разных типов регрессионных моделей.

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

```{r}
library(tidyverse)
heroes <- read_csv("https://raw.githubusercontent.com/Pozdniakov/tidy_stats/master/data/heroes_information.csv",
                   na = c("-", "-99"))
heroes$good <- heroes$Alignment == "good"
```

Обычная линейная модель нам не подходит, если распределение наших ошибок далеко от нормального. А это значит, что мы не можем использовать общую линейную модель с бинарной объясняемой переменной. Эту проблему решает логистическая регрессия, которая является частным случаем обобщенной линейной модели.

Для этого нам понадобится функция `glm()`, а не `lm()` как раньше. Ее синтаксис очень похож, но нам теперь нужно задать еще один важный параметр `family =` для выбора связывающей функции (в данном случае это логит-функция, которая является связующей функцией по умолчанию для биномиального семейства функций в `glm()`).

```{r}
heroes_good_glm <- glm(good ~ Weight + Gender, heroes, family = binomial()) 
summary(heroes_good_glm)
```

Результат очень похож по своей структуре на `glm()`, однако вместо $R^2$ перед нами AIC. AIC расшифровывается как информационный критерий Акаике (Akaike information criterion) --- это критерий использующийся для выбора из нескольких моделей. Чем он меньше, тем лучше модель. Как и Adjusted R^2^, AIC "наказывает" за большое количество параметров в модели.

Поскольку AIC --- это относительный показатель качества модели, нам нужно сравнить его с AIC другой, более общей модели. Можно сравнить с моделью без веса супергероев.

```{r}
heroes_good_glm_noweight <- glm(good ~ Gender, heroes, family = binomial()) 
summary(heroes_good_glm_noweight)
```

AIC стал больше, следовательно, мы выберем модель с весом супергероев.

## Модель со смешанными эффектами {#sec-lme_model}

Модели со смешанными эффектами (mixed-effects models) --- это то же самое, что и иерархическая регрессия (hierarchical regression) или многоуровневое моделирование (multilevel modelling). Этому методу повезло иметь много названий - в зависимости от области, в которой он используется. Модели со смешанными эффектами позволяет включать в линейную регрессию не только фиксированные эффекты (fixed effects), но и случайные эффекты (random effects).

Для экспериментальных дисциплин это интересно тем, что в таких моделях можно не усреднять показатели по испытуемым или образцам, а учитывать влияние соотвествующей группирующей переменной как случайный эффект. В отличие от обычного фактора как в линейной регрессии или дисперсионном анализе (здесь он называется фиксированным), случайный эффект не интересует нас сам по себе, а его значения считаются случайной переменной.

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

Для работы со смешанными моделями в R есть два известных пакета: `nlme` и более современный `lme4`.

```{r, eval = FALSE}
install.packages("lme4")
```

```{r}
library(lme4)
```

Для примера возьмем данные исследования влияния депривации сна на время реакции.

```{r}
data("sleepstudy")
```

Данные представлены в длинном формате: каждая строчка --- это усредненное время реакции для одного испытуемого в соответствующий день эксперимента.

```{r}
sleepstudy %>%
  ggplot(aes(x = Days, y = Reaction)) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = 0:9) +
  facet_wrap(~Subject) +
  theme_minimal()
```

Можно заметить, что, в среднем, время реакции у испытуемых повышается от первого к последнему дню. С помощью смешанных моделей мы можем проверить, различается ли скорость возрастания времени реакции от дня к дню у разных испытуемых.

Для этого мы сравниваем две модели, одна из которых является "вложенной" в другую, то есть усложненной версией более общей модели. В данном случае, более общая модель предполагает, что время реакции увеличивается у всех испытуемых одинаково, а испытуемые различаются только средним временем реакции.

```{r}
sleep_lme0 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
sleep_lme1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
```

Визуализируем предсказания двух моделей:

```{r}
sleepstudy$predicted_by_sleep_lme0 <- predict(sleep_lme0)
sleepstudy$predicted_by_sleep_lme1 <- predict(sleep_lme1)
```

```{r}
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_line(aes(y = predicted_by_M1), colour = "orange") + 
  #geom_line(aes(y = predicted_by_M2), colour = "purple") +   
  geom_point(data = sleepstudy, alpha = 0.4) +
  scale_x_continuous(breaks = 0:9) +
  facet_wrap(~Subject) +
  theme_minimal()
```

Зеленая линия (нулевая модель) имеет везде один и тот же наклон, а синяя (альтернативная модель) имеет разный наклон у всех испытуемых.

Есть несколько способов сравнивать модели, например, уже знакомый нам AIC. Кроме того, можно сравнить две модели с помощью теста хи-квадрат, восполльзовавшись функцией `anova()`.

```{r}
anova(sleep_lme0, sleep_lme1)
```

Модель со случайным наклоном прямой оказалась лучше, о чем нам говорят как более низкие AIC и BIC, так и тестирование с помощью хи-квадрат.