• Принципы выбора лучшей линейной модели
• Сравнение линейных моделей разными способами:
  • частный F-критерий
  • тесты отношения правдоподобий
  • информационные критерии AIC и BIC

• Проверка соответствия наблюдаемых данных предполагаемой функциональной связи между зависимой перменной и предикторами:
  • оценки параметров,
  • тестирование гипотез,
  • оценка объясненной изменчивости (\(R^2\)),
  • анализ остатков
• Построение моделей для предсказания значений в новых условиях:
  • Выбор оптимальной модели
  • Оценка предсказательной способности модели

Все ли хорошо с подобранной моделью?

Какая из этих моделей лучше описывает данные?

В недообученной (underfitted) модели слишком мало параметров, ее предсказания неточны.

В переобученной (overfitted) модели слишком много параметров, она предсказывает еще и случайный шум.

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

При увеличении числа предикторов в модели:

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

Хорошее описание существующих данных

Полный набор переменных:

• большая объясненная изменчивость (\(R^2\)),
• маленькая остаточная изменчивость (\(MS_{Residual}\))
• большие стандартные ошибки
• сложная интерпретация

Принцип парсимонии

Entia non sunt multiplicanda praeter necessitatem

Минимальный набор переменных, который может объяснить существующие данные:

• объясненная изменчивость меньше (\(R^2\)),
• остаточная изменчивость больше (\(MS_{Residual}\))
• стандартные ошибки меньше
• интерпретация проще

Объяснение закономерностей

• Нужны точные тесты влияния предикторов: F-тесты или тесты отношения правдоподобий (likelihood-ratio tests)

Описание функциональной зависимости

• Нужна точность оценки параметров

Предсказание значений зависимой переменной

• Парсимония: “информационные” критерии (АIC, BIC, AICc, QAIC, и т.д.)
• Нужна оценка качества модели на данных, которые не использовались для ее первоначальной подгонки: методы ресамплинга (кросс-валидация, бутстреп)

Не позволяйте компьютеру думать за вас!

• Хорошая модель должна соответствовать условиям применимости

• Другие соображения: разумность, целесообразность модели, простота, ценность выводов, важность предикторов.

От каких характеристик зависит концентрация простат-специфического антигена? (Stamey et al. 1989)

Переменных много, мы хотим из них выбрать оптимальный небольшой набор

97 пациентов:

• lcavol — логарифм объема опухоли
• lweight — логарифм веса
• age — возраст пациента
• lbph — логарифм степени доброкачественной гиперплазии
• svi — поражение семенных пузырьков
• lcp — логарифм меры поражения капсулы
• gleason — оценка по шкале Глисона
• pgg45 — доля оценок 4 и 5 по шкале Глисона
• lpsa — логарифм концентрации простат-специфичного антигена

prost <- read_excel("data/Prostate.xlsx")
mod3 <- lm(lpsa ~ lcavol + lweight + age + lbph + svi + gleason, data = prost)

summary(mod3)

## 
## Call:
## lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + gleason, 
##     data = prost)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8084 -0.4056  0.0057  0.4410  1.5954 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.2709     1.0957    0.25   0.8053    
## lcavol        0.5421     0.0786    6.90    7e-10 ***
## lweight       0.4533     0.1698    2.67   0.0090 ** 
## age          -0.0170     0.0110   -1.55   0.1247    
## lbph          0.1069     0.0583    1.83   0.0702 .  
## svi           0.6943     0.2110    3.29   0.0014 ** 
## gleason       0.1106     0.1159    0.95   0.3425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.708 on 90 degrees of freedom
## Multiple R-squared:  0.648,  Adjusted R-squared:  0.624 
## F-statistic: 27.6 on 6 and 90 DF,  p-value: <2e-16

Не все предикторы влияют, возможно эту модель можно оптимизировать…

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

Удаление предиктора - коэффициент при данном предикторе равен нулю.

Полная модель (full model)

М1: \(y _i = \beta _0 + \beta _1 x _1 + \beta _2 x _2 + \epsilon _i\)

Неполные модели (reduced models)

М2: \(y _i = \beta _0 + \beta _1 x _1 + \epsilon _i\)

М3: \(y _i = \beta _0 + \beta _2 x _2 + \epsilon _i\)

M2 вложена в M1
M3 вложена в M1
M2 и M3 не вложены друг в друга

Нулевая модель (null model), вложена в полную (M1) и в неполные (M2, M3)

\(y _i = \beta _0 + \epsilon _i\)

Для тренировки запишем вложенные модели для данной полной модели

(1)\(y _i = \beta _0 + \beta _1 x _1 + \beta _2 x _2 + \beta _3 x _3 + \epsilon _i\)

Для тренировки запишем вложенные модели для данной полной модели

(1)\(y _i = \beta _0 + \beta _1 x _1 + \beta _2 x _2 + \beta _3 x _3 + \epsilon _i\)

Полная модель

\(y _i = \beta _0 + \beta _1 x _{1\ i} + ... + \beta _k x _{k\ i} + \beta _{l} x _{l\ i} + \epsilon _i\)
\(p = l + 1\) – число параметров \(df _{reduced, full} = p - 1\)
\(df _{error, full} = n - p\)

Уменьшеная модель (без \(\beta _l x _{l\ i}\))

\(y _i = \beta _0 + \beta _1 x _{1\ i} + ... + \beta _k x _{k\ i} + \epsilon _i\)
\(p = k + 1 = l\) – число параметров \(df _{reduced, reduced} = (p - 1) - 1\)
\(df _{error, reduced} = n - (p - 1)\)

Как оценить насколько больше изменчивости объясняет полная модель, чем уменьшенная модель?

• Разница объясненной изменчивости — \(SS _{error,reduced} - SS _{error,full}\)
• С чем, по аналогии с обычным F-критерием, можно сравнить эту разницу объясненной изменчивости?
• Можно сравнить с остаточной изменчивостью полной модели — \(SS _{error, full}\)

Полная модель

\(y _i = \beta _0 + \beta _1 x _{1\ i} + ... + \beta _k x _{k\ i} + \beta _{l} x _{l\ i} + \epsilon _i\)
\(p = l + 1\) – число параметров \(df _{reduced, full} = p - 1\)
\(df _{error, full} = n - p\)

Уменьшеная модель (без \(\beta _l x _{l\ i}\))

\(y _i = \beta _0 + \beta _1 x _{1\ i} + ... + \beta _k x _{k\ i} + \epsilon _i\)
\(p = k + 1 = l\) – число параметров \(df _{reduced, reduced} = (p - 1) - 1\)
\(df _{error, reduced} = n - (p - 1)\)

Частный F-критерий - оценивает выигрыш объясненной дисперсии от включения фактора в модель

\[F = \frac {(SS _{error,reduced} - SS _{error,full}) / (df _{reduced, full} - df _{reduced, reduced})} {(SS _{error, full})/ df _{error, full}}\]

Постепенно удаляем предикторы. Модели обязательно должны быть вложенными! *

Обратный пошаговый алгоритм (backward selection)

• • 1.Подбираем полную модель

• Повторяем 2-3 для каждого из предикторов:
  
• 2.Удаляем один предиктор (строим уменьшенную модель)
  
• 3.Тестируем отличие уменьшенной модели от полной

• 4.Выбираем предиктор для окончательного удаления: это предиктор, удаление которого минимально ухудшает модель. Модель без него будет “полной” для следующего раунда выбора оптимальной модели.

• Повторяем 1-4 до тех пор, пока что-то можно удалить.

• — Важно! Начинать упрощать модель нужно со взаимодействий между предикторами. Если взаимодействие из модели удалить нельзя, то нельзя удалять и отдельно стоящие предикторы, из которых оно состоит. Но мы поговорим о взаимодействиях позже

• Незначимо влияние age, lbph, gleason. Можем попробовать оптимизировать модель

summary(mod3)

## 
## Call:
## lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + gleason, 
##     data = prost)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8084 -0.4056  0.0057  0.4410  1.5954 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.2709     1.0957    0.25   0.8053    
## lcavol        0.5421     0.0786    6.90    7e-10 ***
## lweight       0.4533     0.1698    2.67   0.0090 ** 
## age          -0.0170     0.0110   -1.55   0.1247    
## lbph          0.1069     0.0583    1.83   0.0702 .  
## svi           0.6943     0.2110    3.29   0.0014 ** 
## gleason       0.1106     0.1159    0.95   0.3425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.708 on 90 degrees of freedom
## Multiple R-squared:  0.648,  Adjusted R-squared:  0.624 
## F-statistic: 27.6 on 6 and 90 DF,  p-value: <2e-16

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

mod4 <- update(mod3, . ~ . - age)
anova(mod3, mod4)
mod5 <- update(mod4, . ~ . - lbph)
anova(mod3, mod5)
mod6 <- update(mod3, . ~ . - gleason)
anova(mod3, mod6)
# Удаляем gleason, и потом повторяем все снова...

Но мы пойдем другим путем

Вручную тестировать каждый предиктор с помощью anova() слишком долго. Можно протестировать все за один раз при помощи drop1()

drop1(mod3, test = "F")

## Single term deletions
## 
## Model:
## lpsa ~ lcavol + lweight + age + lbph + svi + gleason
##         Df Sum of Sq  RSS   AIC F value Pr(>F)    
## <none>               45.1 -60.4                   
## lcavol   1     23.84 68.9 -21.2   47.61  7e-10 ***
## lweight  1      3.57 48.6 -55.0    7.13 0.0090 ** 
## age      1      1.20 46.3 -59.8    2.40 0.1247    
## lbph     1      1.68 46.8 -58.8    3.36 0.0702 .  
## svi      1      5.42 50.5 -51.3   10.83 0.0014 ** 
## gleason  1      0.46 45.5 -61.4    0.91 0.3425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Нужно убрать gleason

# Убрали gleason
mod4 <- update(mod3, . ~ . - gleason)
drop1(mod4, test = "F")

## Single term deletions
## 
## Model:
## lpsa ~ lcavol + lweight + age + lbph + svi
##         Df Sum of Sq  RSS   AIC F value  Pr(>F)    
## <none>               45.5 -61.4                    
## lcavol   1     28.77 74.3 -15.9   57.50 2.8e-11 ***
## lweight  1      3.23 48.8 -56.7    6.45 0.01282 *  
## age      1      0.96 46.5 -61.4    1.92 0.16953    
## lbph     1      1.86 47.4 -59.5    3.71 0.05716 .  
## svi      1      5.95 51.5 -51.5   11.90 0.00085 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Нужно убрать age 

# Убрали age 
mod5 <- update(mod4, . ~ . - age )
drop1(mod5, test = "F")

## Single term deletions
## 
## Model:
## lpsa ~ lcavol + lweight + lbph + svi
##         Df Sum of Sq  RSS   AIC F value  Pr(>F)    
## <none>               46.5 -61.4                    
## lcavol   1     27.83 74.3 -17.8   55.08 5.6e-11 ***
## lweight  1      2.80 49.3 -57.7    5.54   0.021 *  
## lbph     1      1.30 47.8 -60.7    2.57   0.112    
## svi      1      5.81 52.3 -51.9   11.49   0.001 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Нужно убрать lbph

# Убрали lbph
mod6 <- update(mod5, . ~ . - lbph)
drop1(mod6, test = "F")

## Single term deletions
## 
## Model:
## lpsa ~ lcavol + lweight + svi
##         Df Sum of Sq  RSS   AIC F value  Pr(>F)    
## <none>               47.8 -60.7                    
## lcavol   1     28.04 75.8 -17.9    54.6 6.3e-11 ***
## lweight  1      5.89 53.7 -51.4    11.5   0.001 ** 
## svi      1      5.18 53.0 -52.7    10.1   0.002 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Больше ничего не убрать

summary(mod6)

## 
## Call:
## lm(formula = lpsa ~ lcavol + lweight + svi, data = prost)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7297 -0.4577  0.0281  0.4640  1.5701 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.2681     0.5435   -0.49    0.623    
## lcavol        0.5516     0.0747    7.39  6.3e-11 ***
## lweight       0.5085     0.1502    3.39    0.001 ** 
## svi           0.6662     0.2098    3.18    0.002 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.717 on 93 degrees of freedom
## Multiple R-squared:  0.626,  Adjusted R-squared:  0.614 
## F-statistic:   52 on 3 and 93 DF,  p-value: <2e-16

Проверьте финальную модель на выполнение условий применимости

1) График расстояния Кука

• Выбросов нет

mod6_diag <- fortify(mod6)
mod6_diag_full <- data.frame(
  mod6_diag, 
  prost[, c("lcp", "pgg45", "age", "lbph", "gleason")])

ggplot(mod6_diag_full, aes(x = 1:nrow(mod6_diag_full), y = .cooksd)) + 
  geom_bar(stat = "identity")

2) График остатков от предсказанных значений

• Выбросов нет
• Гетерогенность дисперсии не выявляется
• Есть намек на нелинейность связей

gg_resid <- ggplot(data = mod6_diag_full, aes(x = .fitted, y = .stdresid)) + 
  geom_point() + geom_hline(yintercept = 0)
gg_resid

3) Графики остатков от предикторов в модели и нет

3) Код для графиков остатков от предикторов в модели и нет

res_1 <- gg_resid + aes(x = lcavol)
res_2 <- gg_resid + aes(x = lweight)
res_3 <- gg_resid + aes(x = age)
res_4 <- gg_resid + aes(x = lbph)
res_5 <- gg_resid + aes(x = svi)
res_6 <- gg_resid + aes(x = lcp)
res_7 <- gg_resid + aes(x = gleason)
res_8 <- gg_resid + aes(x = pgg45)
library(gridExtra)
grid.arrange(res_1, res_2, res_3, res_4,
             res_5, res_6, res_7, res_8, nrow = 2)

4) Квантильный график остатков

• Отклонения от нормального распределения остатков не выявляются

library(car)
qqPlot(mod6)

## [1] 39 96

summary(mod6)

## 
## Call:
## lm(formula = lpsa ~ lcavol + lweight + svi, data = prost)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7297 -0.4577  0.0281  0.4640  1.5701 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.2681     0.5435   -0.49    0.623    
## lcavol        0.5516     0.0747    7.39  6.3e-11 ***
## lweight       0.5085     0.1502    3.39    0.001 ** 
## svi           0.6662     0.2098    3.18    0.002 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.717 on 93 degrees of freedom
## Multiple R-squared:  0.626,  Adjusted R-squared:  0.614 
## F-statistic:   52 on 3 and 93 DF,  p-value: <2e-16

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

Мы оцениваем это как произведение вероятностей получения каждой из точек данных

\[L(\theta| data) = \Pi^n _{i = 1}f(x| \theta)\]

где \(f(data| \theta)\) - функция плотности распределения с параметрами \(\theta\)

\(y_i = \beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i} + \epsilon_i\)

\(y_i = \beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i} + \epsilon_i\)

Пусть в нашей модели остатки нормально распределены (\(\epsilon_i \sim N(0, \sigma^2)\)) и их значения независимы друг от друга:

\(N(\epsilon_i; 0, \sigma^2) = \frac {1} { \sqrt {2\pi\sigma^2} } exp (-\frac {1} {2 \sigma^2} \epsilon_i^2)\)

\(y_i = \beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i} + \epsilon_i\)

Пусть в нашей модели остатки нормально распределены (\(\epsilon_i \sim N(0, \sigma^2)\)) и их значения независимы друг от друга:

\(N(\epsilon_i; 0, \sigma^2) = \frac {1} { \sqrt {2\pi\sigma^2} } exp (-\frac {1} {2 \sigma^2} \epsilon_i^2)\)

Функцию правдоподобия (likelihood, вероятность получения нашего набора данных) можно записать как произведение вероятностей:

\(L(\epsilon_i|\mathbf{y}, \mathbf{x}) = \Pi^n _{n = 1} N(\epsilon_i, \sigma^2) = \frac {1} {\sqrt{(2\pi\sigma^2)^n}} exp(- \frac {1} {2\sigma^2} \sum {\epsilon_i}^2)\)

\(y_i = \beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i} + \epsilon_i\)

Пусть в нашей модели остатки нормально распределены (\(\epsilon_i \sim N(0, \sigma^2)\)) и их значения независимы друг от друга:

\(N(\epsilon_i; 0, \sigma^2) = \frac {1} { \sqrt {2\pi\sigma^2} } exp (-\frac {1} {2 \sigma^2} \epsilon_i^2)\)

Функцию правдоподобия (likelihood, вероятность получения нашего набора данных) можно записать как произведение вероятностей:

\(L(\epsilon_i|\mathbf{y}, \mathbf{x}) = \Pi^n _{n = 1} N(\epsilon_i, \sigma^2) = \frac {1} {\sqrt{(2\pi\sigma^2)^n}} exp(- \frac {1} {2\sigma^2} \sum {\epsilon_i}^2)\)

Поскольку \(\epsilon_i = y_i - (\beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i})\)

то функцию правдоподобия можно переписать так:

\(L(\beta_1...\beta_{p - 1}, \sigma^2| \mathbf{y}, \mathbf{x}) = \frac {1} {\sqrt{(2\pi\sigma^2)^n}}exp(- \frac {1} {2\sigma^2} \sum (y_i - (\beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i}))^2)\)

Чтобы найти параметры модели

\(y_i = \beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i} + \epsilon_i\)

нужно найти такое сочетание параметров \(\beta_0\), \(\beta_1\), … \(\beta_{p - 1}\), и \(\sigma^2\), при котором функция правдоподобия будет иметь максимум:

\(\begin{array}{l} L(\beta_1...\beta_{p - 1}, \sigma^2| \mathbf{y}, \mathbf{x}) &= \frac {1} {\sqrt{(2\pi\sigma^2)^n}} exp(- \frac {1} {2\sigma^2} \sum {\epsilon_i}^2) = \\ &= \frac {1} {\sqrt{(2\pi\sigma^2)^n}}exp(- \frac {1} {2\sigma^2} \sum (y_i - (\beta_0 + \beta_1x_{1i} + ... ...\beta_{p - 1}x_{p - 1\ i}))^2) \end{array}\)

Вычислительно проще работать с логарифмами правдоподобий (loglikelihood)

Если функция правдоподобия

\(\begin{array}{l} L(\beta_1...\beta_{p - 1}, \sigma^2| \mathbf{y}, \mathbf{x}) &= \frac {1} {\sqrt{(2\pi\sigma^2)^n}} exp(- \frac {1} {2\sigma^2} \sum {\epsilon_i}^2) = \\ &= \frac {1} {\sqrt{(2\pi\sigma^2)^n}}exp(- \frac {1} {2\sigma^2} \sum (y_i - (\beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i}))^2) \end{array}\)

то логарифм правдоподобия

\(\begin{array}{l} logLik(\beta_1...\beta_{p - 1}, \sigma^2| \mathbf{y}, \mathbf{x}) &= & \\ ln L(\beta_1...\beta_{p - 1}, \sigma^2| \mathbf{y}, \mathbf{x}) &= &- \frac{n}{2} (ln2\pi + ln\sigma^2) - \frac{1}{2\sigma^2}(\sum \epsilon^2_i) = \\ &= &- \frac{n}{2} (ln2\pi + ln\sigma^2) - \\ & &- \frac{1}{2\sigma^2}(\sum (y_i - (\beta_0 + \beta_1x_{1i} + ... + \beta_{p - 1}x_{p - 1\ i}))^2) \end{array}\)

Чем больше логарифм правдоподобия тем лучше модель

Для подбора параметров методом максимального правдоподобия используют функцию glm()

# Симулированные данные
set.seed(9328)
dat <- data.frame(X = runif(n = 50, min = 0, max = 10))
dat$Y <- 3 + 15 * dat$X + rnorm(n = 50, mean = 0, sd = 1)

# Подбор модели двумя способами
Mod     <-  lm(Y ~ X, data = dat) # МНК
Mod_glm <- glm(Y ~ X, data = dat) # МЛ

# Одинаковые оценки коэффициентов
coefficients(Mod)

## (Intercept)           X 
##        2.95       15.02

coefficients(Mod_glm)

## (Intercept)           X 
##        2.95       15.02

\(LogLik\) для модели можно найти с помощью функции logLik()

logLik(Mod_glm)

## 'log Lik.' -65.5 (df=3)

# Предсказанные значения
dat$predicted <- predict(Mod)
# Оценка дисперсии
SD <- summary(Mod)$sigma 
# Вероятности для каждой точки
dat$Prob <- dnorm(dat$Y, mean = dat$predicted, sd = SD)
# Логарифм вероятностей
dat$LogProb <- log(dat$Prob)
# Логарифм произведения, равный сумме логарифмов
sum(dat$LogProb)

## [1] -65.5

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

\[LRT = 2ln\Big(\frac{L_1}{L_2}\Big) = 2(lnL_1 - lnL_2)\]

• \(L_1\), \(L_2\) - правдоподобия полной и уменьшенной модели
• \(lnL_1\), \(lnL_2\) - логарифмы правдоподобий

Разница логарифмов правдоподобий имеет распределение, которое можно апроксимировать \(\chi^2\), с числом степеней свободы \(df = df_2 - df_1\) (Wilks, 1938)

Для этой полной модели

GLM1 <- glm(lpsa ~ lcavol + lweight + age + lbph + svi + gleason, data = prost)

Подберите оптимальную модель при помощи тестов отношения правдоподобий

Тест отношения правдоподобий можно сделать с помощью тех же функций, что и частный F-критерий:

• по-одному anova(mod3, mod2, test = "Chisq")
• все сразу drop1(mod3, test = "Chisq")

drop1(GLM1, test = "Chisq")

## Single term deletions
## 
## Model:
## lpsa ~ lcavol + lweight + age + lbph + svi + gleason
##         Df Deviance AIC scaled dev.      Pr(>Chi)    
## <none>         45.1 217                              
## lcavol   1     68.9 256        41.2 0.00000000014 ***
## lweight  1     48.6 222         7.4        0.0066 ** 
## age      1     46.3 218         2.6        0.1100    
## lbph     1     46.8 218         3.6        0.0594 .  
## svi      1     50.5 226        11.0        0.0009 ***
## gleason  1     45.5 216         1.0        0.3231    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Нужно убрать gleason

# Убираем gleason
GLM2 <- update(GLM1, . ~ . - gleason)
drop1(GLM2, test = "Chisq")

## Single term deletions
## 
## Model:
## lpsa ~ lcavol + lweight + age + lbph + svi
##         Df Deviance AIC scaled dev. Pr(>Chi)    
## <none>         45.5 216                         
## lcavol   1     74.3 261        47.5  5.5e-12 ***
## lweight  1     48.8 220         6.6  0.00998 ** 
## age      1     46.5 216         2.0  0.15497    
## lbph     1     47.4 218         3.9  0.04893 *  
## svi      1     51.5 226        11.9  0.00056 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Нужно убрать age

# Убираем lbph
GLM3 <- update(GLM2, . ~ . - age)
drop1(GLM3, test = "Chisq")

## Single term deletions
## 
## Model:
## lpsa ~ lcavol + lweight + lbph + svi
##         Df Deviance AIC scaled dev. Pr(>Chi)    
## <none>         46.5 216                         
## lcavol   1     74.3 259        45.5  1.5e-11 ***
## lweight  1     49.3 220         5.7  0.01720 *  
## lbph     1     47.8 217         2.7  0.10190    
## svi      1     52.3 225        11.4  0.00073 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Нужно убрать lbph

# Убираем lbph
GLM4 <- update(GLM3, . ~ . - lbph)
drop1(GLM4, test = "Chisq")

## Single term deletions
## 
## Model:
## lpsa ~ lcavol + lweight + svi
##         Df Deviance AIC scaled dev. Pr(>Chi)    
## <none>         47.8 217                         
## lcavol   1     75.8 259        44.8  2.2e-11 ***
## lweight  1     53.7 226        11.3  0.00078 ***
## svi      1     53.0 225        10.0  0.00158 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Больше ничего убрать не получается

summary(GLM4)

## 
## Call:
## glm(formula = lpsa ~ lcavol + lweight + svi, data = prost)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7297  -0.4577   0.0281   0.4640   1.5701  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.2681     0.5435   -0.49    0.623    
## lcavol        0.5516     0.0747    7.39  6.3e-11 ***
## lweight       0.5085     0.1502    3.39    0.001 ** 
## svi           0.6662     0.2098    3.18    0.002 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.514)
## 
##     Null deviance: 127.918  on 96  degrees of freedom
## Residual deviance:  47.785  on 93  degrees of freedom
## AIC: 216.6
## 
## Number of Fisher Scoring iterations: 2

Финальную модель снова нужно проверить на выполнение условий применимости

\(AIC = -2 logLik + 2p\)

• \(logLik\) - логарифм правдоподобия для модели
• \(2p\) - штраф за введение в модель \(p\) параметров

Чем меньше AIC - тем лучше модель

Важно! Информационные критерии можно использовать для сравнения даже для невложенных моделей. Но модели должны быть подобраны с помощью ML и на одинаковых данных!

• \(logLik\) - логарифм правдоподобия для модели
• \(p\) - число параметров
• \(n\) - число наблюдений

AIC(GLM1, GLM2, GLM3, GLM4)

##      df AIC
## GLM1  8 217
## GLM2  7 216
## GLM3  6 216
## GLM4  5 217

• Судя по AIC, лучшая модель GLM2 или GLM3. Но если значения AIC различаются всего на 1-2 единицу — такими различиями можно пренебречь и выбрать более простую модель (GLM3).

Уравнение модели:

\(lpsa = 0.15 + 0.55lcavol + 0.39lweight + 0.09lbph + 0.71svi\)

BIC(GLM1, GLM2, GLM3, GLM4)

##      df BIC
## GLM1  8 238
## GLM2  7 234
## GLM3  6 231
## GLM4  5 229

• Судя по BIC, нужно выбрать модель GLM4

Уравнение модели:

\(lpsa = - 0.27 + 0.55lcavol + 0.51lweight + 0.67svi\)

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

Не важно, какой из способов выбрать, но важно сделать это заранее, до анализа, чтобы не поддаваться соблазну подгонять результаты.

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

• Must read paper! Zuur, A.F. and Ieno, E.N., 2016. A protocol for conducting and presenting results of regression‐type analyses. Methods in Ecology and Evolution, 7(6), pp.636-645.

• Кабаков Р.И. R в действии. Анализ и визуализация данных на языке R. М.: ДМК Пресс, 2014
• Zuur, A., Ieno, E.N. and Smith, G.M., 2007. Analyzing ecological data. Springer Science & Business Media.
• Quinn G.P., Keough M.J. 2002. Experimental design and data analysis for biologists

• Logan M. 2010. Biostatistical Design and Analysis Using R. A Practical Guide