Краткое summary, что мы уже знаем

Знаем базовый R и умеем работать в пакете tidyverse
Умеем тестировать статистические гипотезы и проверять значимость рассчитанных статистик
Умеем подбирать модель простой и множественной линейной регрессии и тестировать ее качество
Шарим в дисперсионном анализе и понимаем, что это та же саммая регрессия, только без количественного предиктора.

Ограничения общих линейных моделей

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

Вспомним, как выглядит уравнение такой модели: \[y = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p + \varepsilon\]

Предикторы в такой модели могут быть как дискретными, так и непрерывными.

Важнейшим допущением / требованием этой модели является распределение ошибки: \[\varepsilon \thicksim \mathcal{N}(0, \, \sigma)\]

Поскольку ошибка модели должна быть распределена нормально, а моделируется среднее значение, то можно сформулировать более общее допущение / требование: \[y_i \thicksim \mathcal{N}(\mu_i, \, \sigma)\]

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

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

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

Биномиальное распределение

Биномальное распределение описывает вероятность выпадения определенного числа «успехов» в серии испытаний Бернулли.

Здесь можно поиграться.

Бинарные переменные

Эти переменные достаточно широко распространены как повседневности, так и в науке. Блюдо вкусное или невкусное, команда выиграла или проиграла, пациент в результате медицинских манипуляций выжил или умер, в ходе эксперимента была выбрана какая-то опция или нет, сдал студент экзамен или не сдал — и т.д.

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

Морские звезды и мидии

Данные Khaitov et al., 2018

Тут есть датасет про мидий. Это экспериментальные данные, которые были собраны, чтобы понять, различают ли морские звезды два вида мидий.

Ценные промысловые моллюски, атлантические мидии, обитают в Белом море, но недавно туда вселились мидии другого вида — тихоокеанские. Вселенцы имеют меньшую промысловую значимость, и важно понять, что регулирует их численность. Наиболее значимый фактор — морские звезды, которые ими питаются.

Грузим данные.

library(tidyverse)

## ── Attaching packages ──────────────────────────────────────────────── tidyverse 1.2.1 ──

## ✔ ggplot2 3.2.0     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   1.0.2     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0

## ── Conflicts ─────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

theme_set(theme_bw())

mussels <- read_csv('https://raw.githubusercontent.com/angelgardt/hseuxlab-andan/master/aster_mussel.csv')

## Parsed with column specification:
## cols(
##   Box = col_double(),
##   L = col_double(),
##   Sp = col_character(),
##   Outcome = col_character()
## )

Смотрим их структуру.

str(mussels)

## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 286 obs. of  4 variables:
##  $ Box    : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ L      : num  33.1 24.8 33 18.8 34 22.6 34.5 17.8 29.9 21.4 ...
##  $ Sp     : chr  "Tr" "Ed" "Ed" "Ed" ...
##  $ Outcome: chr  "not_eaten" "not_eaten" "not_eaten" "not_eaten" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   Box = col_double(),
##   ..   L = col_double(),
##   ..   Sp = col_character(),
##   ..   Outcome = col_character()
##   .. )

У нас есть следуюшие переменные:

Outcome — съели мидию или нет
Sp — вид мидий (Ed — корреной вид, Tr — вселенец)
L — длина мидий (мм)
Box — номер контейнера, в котором прододился эксперимент

Поизучаем и подправим данные.

colSums(is.na(mussels))

##     Box       L      Sp Outcome 
##       0       0       0       0

table(mussels$Box)

## 
##  1  2  3  4 
## 66 68 74 78

Подправим данные. Сделаем факторными переменные, которые должны быть таковыми, и перекодируем Outcome.

mussels %>% mutate(Box = as_factor(Box),
                   Sp = as_factor(Sp),
                   Out = ifelse(mussels$Outcome == 'eaten', 1, 0)) -> mussels

Перекодирование нужно для того, чтобы алгоритм подборки модели корректно отработал.

Задача построения линейной модели

Чтобы ответить на вопрос, различают ли звезды разные виды мидий, нам нужно построить модель, описывающую связь Outcome (=Out) с предикторами Sp и L.

Попробуем обычную линейную регрессию.

fit0 <- lm(Out ~ Sp * L, mussels)
summary(fit0)

## 
## Call:
## lm(formula = Out ~ Sp * L, data = mussels)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41803 -0.11835 -0.08124 -0.04458  0.96140 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.725105   0.135773   5.341 1.91e-07 ***
## SpEd        -0.479496   0.178580  -2.685 0.007681 ** 
## L           -0.019811   0.005119  -3.870 0.000135 ***
## SpEd:L       0.013519   0.006715   2.013 0.045024 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3146 on 282 degrees of freedom
## Multiple R-squared:  0.09123,    Adjusted R-squared:  0.08156 
## F-statistic: 9.436 on 3 and 282 DF,  p-value: 5.82e-06

mussels %>% ggplot(aes(L, Out, color = Sp)) +
  geom_point(size=2, alpha=.5) +
  geom_line(aes(y = fit0$fitted.values))

Но получается что-то странное…

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

Логистическая кривая

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

Сама зависимая переменная в зависимости от предиктора распределена примерно так:

Введем новые обозначения:

\(p_i\) — вероятность события \(y_i = 1\) при данных значениях предиктора,
\(1 - p_i\) — вероятность альтернативного события \(y_i = 0\).

Получается непрерывная величина \(0 \leq p_i \leq 1\).

И вроде бы как ее можно уже моделировать. Но нужно помнить, что вероятность изменяется в пределах от нуля до единицы, а прямая ничем не ограничена. Поэтому прямая — не лучший вариант.

Такая закономерность моделируется логистической кривой.

Она лежит как раз в пределах от 0 до 1. Наша логистическая кривая задается уравнением \[p_i = \frac{e^{\beta_0 + \beta_1 x_i}}{1 + e^{\beta_0 + \beta_1 x_i}}\]

Логистическую кривую мы больше никогда не увидим [печаль :(], но она используется внутри функций, которыми мы строим модель.

Шансы и логиты

Теперь нам надо побороться с ограниченностью логистической кривой. Для этого можно перейти от вероятностей к шансам.

Шанс (отношение шансов, odds, odds ratio) — это отношение вероятности успеха к вероятности неудачи. Их величина варьируется от \(0\) до \(+\infty\).

Уже лучше, но все еще не самый лучший вариант… Последний шаг, необходимый нам, чтобы все было хорошо, юзануть логарифм, который преобразуем шансы в логиты.

\[ \mathrm{logit}(p) = \ln\left(\frac{p_i}{1-p_i}\right) \]

Значения логитов — трансформированные оценки вероятностей события. Они варьируют от \(-\infty\) до \(+\infty\), симметричны относительно нуля, и их удобно брать в качестве зависимой переменной для построения модели. Кроме того, logit-преобразование еще и линеаризует логистическую кривую. Доказательство.

Функция, используемая для линеаризации связи между предиктором и зависимой переменной, называется связывающей функцией (linked function). Функция logit-преобразоввания — одна из нескольких связывающих функций, применяемых для анализа бинарных переменных отклика.

Итак, summary от всего, что было выше: 1. От дискретной оценки событий (0 и 1) переходим к оценке вероятностей. 2. Связь вероятностей с предиктором описывается логистической кривой. 3. Если при помощи функции связи перейти от вероятностей к логитам, то связь будет описываться прямой линией. 4. Параметры линейной модели для такой прямой можно оценить с помощью регрессионного анализа.

Формулирование модели в математическом виде

GLM с биномиальным распределением отклика

\[y_i \thicksim \mathrm{Bin}(n = 1, p_i) \\ \mathbb{E}(y_i) = p_i = \frac{e^{\beta_0 + \beta_1 x_1 + \dots + \beta_p x_p}}{1 + e^{\beta_0 + \beta_1 x_1 + \dots + \beta_p x_p}}\]

Функция связи (linked function): \[ \ln \left(\frac{p_i}{1 - p_i}\right) = \eta_i \\ \eta_i = \beta_0 + \beta_1 x_{1i} + \dots + \beta_p x_{pi} \]

Для перехода от логитов к вероятностям используется функция, обратная функции связи. В данном случае, логистическое преобразование: \[ p_i = \frac{e^{\eta_i}}{1 + e^{\eta_i}} \]

Подбор линейной модели

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

fit <- glm(Out ~ Sp*L, family = binomial(link = 'logit'), mussels)
summary(fit)

## 
## Call:
## glm(formula = Out ~ Sp * L, family = binomial(link = "logit"), 
##     data = mussels)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1053  -0.4917  -0.4066  -0.3128   2.5156  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  1.82408    1.08147   1.687  0.09167 . 
## SpEd        -1.98930    1.73332  -1.148  0.25110   
## L           -0.12879    0.04601  -2.799  0.00513 **
## SpEd:L       0.03895    0.07146   0.545  0.58576   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 212.58  on 285  degrees of freedom
## Residual deviance: 190.94  on 282  degrees of freedom
## AIC: 198.94
## 
## Number of Fisher Scoring iterations: 5

Модель сошлась, все ровно. Вызываем функцию summary(), чтобы посмотреть статистики модели. Вывод функции нам уже знаком, он аналогичен тому, с которым мы встречались, однако все-таки отличается некоторыми деталями.

Во-первых, нет F-статистики. Во-вторых, нет R². Это связано с тем, что алгоритм GLM не работает с дисперсией и суммой квадратов — модель подбирается методом максимального правдоподобия. Данный метод позволяет учесть разную форму распределений, которые мы можем моделировать с использованием GLM.

Метод максимального правдоподобия

Правдоподобие (likelihood) — способ измерить соответствие имеющихся данных тому, что можно получить при определенных значениях параметров модели. Оно представляет собой произведение вероятностей получения каждой из точек данных:

\[ L(\theta | \mathrm{data}) = \prod_{i=1}^n f(\mathrm{data}|\theta), \] где \(f(\mathrm{data}|\theta)\) — функция распределения с параметрами \(\theta\).

Параметры модели должны максимизировать значения [логарифма] правдоподобия, т.е. \[ L(\theta | \mathrm{data}) \rightarrow \max \] Однако для упрощения вычислений используют логарифмы правдоподобий (loglikelihood) и максимизируют их \[ \ln L(\theta | \mathrm{data}) \rightarrow \max \] Аналитически такие задачи решаются редко, чаще используются методы численной оптимизации.

Анализируем модель

Итак, продолжим изучать аутпут функции summary(). Обратимся к ней еще раз.

summary(fit)

## 
## Call:
## glm(formula = Out ~ Sp * L, family = binomial(link = "logit"), 
##     data = mussels)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1053  -0.4917  -0.4066  -0.3128   2.5156  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  1.82408    1.08147   1.687  0.09167 . 
## SpEd        -1.98930    1.73332  -1.148  0.25110   
## L           -0.12879    0.04601  -2.799  0.00513 **
## SpEd:L       0.03895    0.07146   0.545  0.58576   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 212.58  on 285  degrees of freedom
## Residual deviance: 190.94  on 282  degrees of freedom
## AIC: 198.94
## 
## Number of Fisher Scoring iterations: 5

Мы видим знакомую нам табличку с оценками коэффициентов. В нашем случае есть intercept, которые содержит один из уровеней категориального предиктора (SpTr), второй уровень категориального предиктора (SpEd) и количественный предиктор на каждом из двух уровней (L и SpEd:L). У каждого предиктора есть информация о значении коэффициента при нем (Estimate), значение стандартной ошибки, z-value и p-value, рассчитанный для последнего.

Судя по аутпуту, значимым является предиктор длины только для вида-вселенца. Но так ли это?

Внимательно посмотрим на эту табличку. Раньше у нас были значения t-статистики — теперь z-теста. К чему это может приводить?

Чтобы ответить на этот вопрос, надо разобраться, что такое этот z-тест.

Значение z-value является результатом подсчета теста Вальда и позволяет оценить значимость коэффициента модели. Расчет статистики похож на подсчет t-теста, только распределение данной статистики ассимптотически стремиться к нормальному (отсюда и z). Ассимптотика приводит к тому, что опеределение значимости коэффициентов на маленьких выборках будет неточным.

\[ H_0: \beta_k = 0 \\ H_1: \beta_k \neq 0 \\ \\ z = \frac{b_k - \beta_k}{SE_{b_k}} = \frac{b_k}{SE_{b_k}} \thicksim \mathcal{N}(0,\,1) \]

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

Анализ девиансы и упрощение модели

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

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

Дла понимания механики анализа нам потребуется несколько полезных абстракций:

Насыщенная модель (saturated model) — каждое уникальное наблюдение (сочетание предикторов) описывается одним из \(n\) параметров.

\[ \ln L_{\mathrm{sat}} = 0 \\ \mathrm{df}_{\mathrm{sat}} = n - p_{\mathrm{sat}} = n - n = 0 \]

Нулевая модель (null model) — все наблюдения описываются одним параметром (средним значением)

\[ \hat y_i = \beta_0 \\ \ln L_{\mathrm{null}} \neq 0 \rightarrow -\infty \\ \mathrm{df}_{\mathrm{null}} = n - p_{\mathrm{null}} = n - 1 \]

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

\[ \hat y_i = \beta_0 + \beta_1 x_{1i} + \dots + \beta_p x_{pi} \\ \ln L_{\mathrm{model}} \neq 0 \\ \mathrm{df}_\mathrm{model} = n - p_\mathrm{model} \]

Девианса является мерой различия правдоподобий двух моделей (оценка разницы логарифмов правдоподобий). [см. рисунок выше]

Остаточная девианса \[ d_\mathrm{resid} = 2\big(\ln L_\mathrm{sat} - \ln L_\mathrm{model} \big) = -2 \ln(L_\mathrm{model}) \]

Нулевая девианса \[ d_\mathrm{null} = 2 \big(\ln L_\mathrm{sat} - \ln L_\mathrm{null}\big) = -2 \ln(L_\mathrm{null}) \]

Сравнение нулевой и остаточной девианс позволяет судить о статистической значимости модели в целом. Такое сравнение проводится с помощью теста отношения правдоподобий (likelihood ratio test, LRT).

\[ d_\mathrm{null} - d_\mathrm{resid} = -2 ( \ln L_\mathrm{null} - \ln L_\mathrm{model} ) = 2 (\ln L_\mathrm{model} - \ln L_\mathrm{null}) = 2 \ln \left( \frac{L_\mathrm{model}}{L_\mathrm{null}} \right) \\ \mathrm{LRT} = 2 \ln \left( \frac{L_\mathrm{M_1}}{L_\mathrm{M_2}} \right) = 2 (\ln L_\mathrm{M_1} - L_\mathrm{M_2}), \] где \(M_1,\, M_2\) — вложенные модели (\(M_1\) — более полная, \(M_2\) — уменьшенная), \(L_\mathrm{M_1}, \, L_\mathrm{M_2}\) — правдоподобия.

Распределенеи разницы логарифмов правдоподобий аппроксиммируется распределением \(\chi^2\) с числом степеней свободы \(\mathrm{df} = \mathrm{df}_\mathrm{M_2} - \mathrm{df}_\mathrm{M_1}\).

LRT для тестирования значимости модели в целом:

\[ \mathrm{LRT} = 2 \ln \left( \frac{L_\mathrm{model}}{L_\mathrm{null}} \right) = 2 (\ln L_\mathrm{model} - \ln L_{\mathrm{null}}) = d_\mathrm{null} - d_\mathrm{resid} \\ \mathrm{df} = p_\mathrm{model} - 1 \]

LRT для тестирования значимости отдельных предикторов:

\[ \mathrm{LRT} = 2 \ln \left( \frac{L_\mathrm{model}}{L_\mathrm{reduced}} \right) = 2 (\ln L_\mathrm{model} - \ln L_\mathrm{reduced}) \\ \mathrm{df} = p_\mathrm{model} - p_\mathrm{reduced} \]

Реализуем изученное в R.

Протестируем значимость модели в целом. Для этого неоходимо создать нулевую модель с одним предиктором — среднее — и передать эту модель и тестируемую в функцию anova() с указанием параметра test = 'Chi'. Получим нужным нам аутпут:

fit_null <- glm(Out ~ 1, family = binomial(link = 'logit'), mussels)
anova(fit_null, fit, test = 'Chi')

## Analysis of Deviance Table
## 
## Model 1: Out ~ 1
## Model 2: Out ~ Sp * L
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       285     212.58                          
## 2       282     190.94  3   21.631 7.783e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Как можно увидеть, наша модель получилась статистически значима.

Для тестирования значимости отдельных предикторов можно (и лучше) использовать функцию Anova() из пакета car.

Anova(fit)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Out
##      LR Chisq Df Pr(>Chisq)    
## Sp     7.8938  1  0.0049604 ** 
## L     11.7480  1  0.0006091 ***
## Sp:L   0.2960  1  0.5863920    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Мы получаем табличку, похожую на результат работы summary(), но с другими статистиками.

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

drop1(fit, test = 'Chi')

## Single term deletions
## 
## Model:
## Out ~ Sp * L
##        Df Deviance    AIC     LRT Pr(>Chi)
## <none>      190.94 198.94                 
## Sp:L    1   191.24 197.24 0.29601   0.5864

Например, в нашем случае взаимодействие факторов оказалось статистически незначимым. Следовательно, его можно удалить из модели. Обновим имеющуюся модель в помощью функции update() и запишем в новый объект.

fit1 <- update(fit, .~. -Sp:L)
drop1(fit1, test = 'Chi')

## Single term deletions
## 
## Model:
## Out ~ Sp + L
##        Df Deviance    AIC     LRT  Pr(>Chi)    
## <none>      191.24 197.24                      
## Sp      1   199.13 203.13  7.8938 0.0049604 ** 
## L       1   202.99 206.99 11.7480 0.0006091 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Повторное тестирование модели показывает, что оба оставшихся предиктора статистически значимы — больше модель упростить нельзя. Таким образом, мы получили нашу финальную модель.

На прошлом занятии мы упоминали информационный критерий Акаике (Akaike information criterion), позволяющий сделать выбор из нескольких моделей. Его можно использовать и в нашем случае как дополнительную проверку.

AIC(fit, fit1)

##      df      AIC
## fit   4 198.9443
## fit1  3 197.2404

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

Интерпретация коэффициентов модели

Еще раз позовем summary модели:

summary(fit1)

## 
## Call:
## glm(formula = Out ~ Sp + L, family = binomial(link = "logit"), 
##     data = mussels)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0583  -0.5099  -0.4096  -0.2894   2.5934  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  1.46896    0.85072   1.727  0.08421 . 
## SpEd        -1.06973    0.37879  -2.824  0.00474 **
## L           -0.11328    0.03497  -3.239  0.00120 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 212.58  on 285  degrees of freedom
## Residual deviance: 191.24  on 283  degrees of freedom
## AIC: 197.24
## 
## Number of Fisher Scoring iterations: 5

Что значат эти непонятные коэффициенты при предикторах? В общих линейных моделях все ясно как белый день, а тут какая-то дичь…

Посмотрим на математическую запись нашей модели:

\[ \mathrm{Out}_i \thicksim \mathrm{Bin}(n = 1, \, p_i) \\ \mathbb{E}(\mathrm{Out}_i) = p_i \\ \ln \left(\frac{p_i}{1 - p_i}\right) = \eta_i \\ \eta_i = 1.469 - 1.070 \, \mathrm{Sp}_{\mathrm{Ed} \, i} - 0.113 \, \mathrm L_{i} \]

Наша зависимая переменная — логарифм отношения шансов. От этого и будем толкаться.

Тогда,

\(b_0\), интерсепт, показывает логарифм отношения шансов для базового уровня (Tr) дискретного предиктора (Sp);
\(b_1\) показывает, на сколько единиц изменяется логарифм отношения шансов для данного уровня (Ed) дискретного предиктора (Sp) по сравнению с базовым уровнем (Tr);
\(b_2\) показывает, на сколько единиц изменится логарифм отношения шансов при изменении значения предиктора (L) на единицу.

Корректно, но непонятно… Чтобы разобраться лучше, немного потупим в алгебру.

Немного алгебры для понимания сути коэффициентов

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

\[ \eta = b_0 + b_1 x \]

В терминах логитов её можно записать следующим образом:

\[ \eta = \ln \left( \frac{p}{1 - p} \right) = \ln (\mathrm{odds}) \]

Как изменится предсказание модели при изменении предиктора на единицу?

\[ \eta_{x+1} - \eta_x = \ln (\mathrm{odds}_{x+1}) - \ln (\mathrm{odds}_{x}) = \ln \left( \frac{\mathrm{odds}_{x+1}}{\mathrm{odds}_x} \right) \] \[ \eta_{x+1} - \eta_x = b_0 + b_1 (x + 1) - b_0 - b_1x = b_1 \] \[ \ln \left( \frac{\mathrm{odds}_{x+1}}{\mathrm{odds}_{x}} \right) = b_1 \] \[ \frac{\mathrm{odds}_{x+1}}{\mathrm{odds}_{x}} = e^{b_1} \]

Таким образом, \(e^{b_1}\) показывает, во сколько раз изменится отношение шансов при увеличении предиктора на единицу. Для дискретных предикторов \(e^{b_1}\) покажет, во сколько раз различается отношение шансов для данного уровня предиктора по сравнению с базовым.

Теперь мы можем тракровать полученные коэффициенты нормальным языком.

\[ \eta_i = 1.469 - 1.070 \, \mathrm{Sp}_{\mathrm{Ed} \, i} - 0.113 \, \mathrm L_{i} \]

При увеличении длины тела мидии на 1 мм отношения шансов быть съеденной увеличатся в \(e^{−0.113} = 0.893\) раза. То есть мидия, имеющая больший размер, имеет меньше шансов быть съеденной (больше шансов на выживание).

Отношение шансов быть съеденной у мидии, относящейся к группе Ed дискретного фактора Sp, отличается в \(e^{-1.07} = 0.343\) раза от аналогичной характеристики мидии, относящейся к базовому уровню (Tr). То есть вероятность выжить у мидии из группы Ed больше, чем у мидии из группы Tr.

Предсказательная сила модели

Проверим, насколько эффективна наша модель для предсказаний. Воспользуется функцией predict() и создадим новую колонку с предсказанными вероятностями.

mussels %>% mutate(predicted = predict(fit1, type = 'response')) -> mussels

Далее нам необходимо выбрать критерий — пороговое значение вероятность —- после которого мы будем решать, что данную мидию съели. Выберем очень либеральный критерий. Результаты заносим в новый столбец.

mussels %>% mutate(out_pred = ifelse(mussels$predicted > .6, 'eaten_p', 'not_eaten_p')) -> mussels

Теперь построим confusion matrix и посмотрим на качество предсказаний:

table(mussels$Outcome, mussels$out_pred)

##            
##             not_eaten_p
##   eaten              35
##   not_eaten         251

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

GLM и биномиальная регрессия

A. Angelgardt