Пример: птицы в лесах Австралии

От каких характеристик лесного участка зависит обилие птиц в лесах юго-западной Виктории, Австралия (Loyn, 1987)

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

Пример из кн. Quinn, Keugh, 2002, данные из Loyn, 1987)

Модель из прошлой лекции

bird <- read.csv("data/loyn.csv")
bird$logAREA <- log(bird$AREA)
bird$logDIST <- log(bird$DIST)
bird$logLDIST <- log(bird$LDIST)
mod2 <- lm(ABUND ~ logAREA + YRISOL + logDIST + logLDIST + ALT, data = bird)

\[\begin{aligned}\widehat{ABUND}_i = &-226.00 + 3.69 \cdot logAREA_i + 0.12 \cdot YRISOL_i \\ &- 0.10 \cdot logDIST_i - 0.33 \cdot logLDIST_i + 0.03 \cdot ALT_i\\ \end{aligned}\]

• Проверим ещё одним способом, влияют ли предикторы.
• Разберёмся, можно ли оптимизировать модель.

Вложенные модели (nested models)

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

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

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

М1: \(y_i = b_0 + b_1 x_{1i} + b_2 x_{2i} + e_i\)

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

М2: \(y_i = b_0 + b_1 x_{1i} + e_i\)

М3: \(y_i = b_0 + b_2 x_{2i} + e_i\)

• М2 вложена в М1
• М3 вложена в М1
• М2 и М3 не вложены друг в друга

Нулевая модель (null model)

\(y_i = b_0 + e_i\) — вложена в полную (М1) и в неполные (М2, М3)

Задание

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

1. \(y_i = b_0 + b_1 x_{1i} + b_2 x_{2i} + b_3 x_{3i} + e_i\)

Сравнение вложенных линейных моделей при помощи F-критерия

Полная модель (со всеми предикторами)

\(y_i = b_0 + b_1 x_{1i} + \ldots + b_k x_{ki} + b_l x_{li} + e_i\)

\(SS_{r, full}\), \(df_{r, full} = p - 1\)

\(SS_{e, full}\), \(df_{e, full} = n - p\)

Уменьшенная модель с \((p - 1)\) параметрами (например, без предиктора \(b_l x_{li}\))

\(y_i = b_0 + b_1 x_{1i} + \ldots + b_k x_{ki} + e_i\)

\(SS_{r, reduced}\), \(df_{r, reduced} = (p - 1) - 1\)

\(SS_{e, reduced}\), \(df_{e, reduced} = n - (p - 1)\)

Полная модель всегда будет лучше

\(SS_{r, full} > SS_{r, reduced}\)

или то же самое:

\(SS_{e, full} < SS_{e, reduced}\)

Можно оценить, насколько именно полная модель лучше уменьшенной (или насколько уменьшенная модель хуже полной).

Частный F-критерий оценивает разницу объяснённой дисперсии между моделями

\[F = \frac{(SS_{e, reduced} - SS_{e, full})/(df_{e, reduced} - df_{e, full})}{SS_{e, full} / df_{e, full}}\]

Если модель значимо ухудшается от удаления предиктора, то влияние этого предиктора значимо.

Краткое обозначение частного F-критерия

Чтобы проверить значимость предиктора ALT, нам нужно сравнить две модели:

ABUND ~ logAREA + YRISOL + logDIST + logLDIST + ALT
ABUND ~ logAREA + YRISOL + logDIST + logLDIST

Это можно обозначить:

\(F_{(\text{ALT | logAREA, YRISOL, logDIST, logLDIST})}\)

А разницу остаточных сумм квадратов \(SS_{e, reduced} - SS_{e, full}\), которую мы тестируем, можно обозначить:

\(SS_{(\text{ALT | logAREA, YRISOL, logDIST, logLDIST})}\)

Частный F-критерий в R: anova(модель_1, модель_2)

Если высота над уровнем моря входит в модель после учёта влияния остальных предикторов, то её влияние не значимо (\(F = 1.84\), \(p = 0.18\)).

\(F_{(\text{ALT | logAREA, YRISOL, logDIST, logLDIST})}\)

mod2_reduced <- update(mod2, . ~ . - ALT)
anova(mod2_reduced, mod2)

Analysis of Variance Table

Model 1: ABUND ~ logAREA + YRISOL + logDIST + logLDIST
Model 2: ABUND ~ logAREA + YRISOL + logDIST + logLDIST + ALT
  Res.Df  RSS Df Sum of Sq    F Pr(>F)
1     51 2206                         
2     50 2128  1      78.2 1.84   0.18

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

Если высота над уровнем моря — единственный предиктор в модели, то его влияние статистически значимо (\(F = 9.45\), \(p < 0.05\)).

\(F_{(\text{ALT})}\)

mod_ALT1 <- lm(ABUND ~ ALT, data = bird)
mod_null <- lm(ABUND ~ 1, data = bird)
anova(mod_null, mod_ALT1)

Analysis of Variance Table

Model 1: ABUND ~ 1
Model 2: ABUND ~ ALT
  Res.Df  RSS Df Sum of Sq    F Pr(>F)   
1     55 6338                            
2     54 5394  1       944 9.45 0.0033 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Порядок проведения сравнений при F-тестах (типы сумм квадратов)

Y ~ A + B + C

Последовательное тестирование anova(модель)

anova(mod2)

Analysis of Variance Table

Response: ABUND
          Df Sum Sq Mean Sq F value  Pr(>F)    
logAREA    1   3471    3471   81.56 4.4e-12 ***
YRISOL     1    607     607   14.27 0.00042 ***
logDIST    1     28      28    0.67 0.41737    
logLDIST   1     25      25    0.59 0.44687    
ALT        1     78      78    1.84 0.18137    
Residuals 50   2128      43                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Осторожно: при последовательном тестировании результат зависит от порядка предикторов

Если мы поставим ALT первым, то изменится SS для других предикторов, а влияние самого ALT станет значимым (т.к. оно оценено без учёта других предикторов).

mod2_reordered <- lm(ABUND ~ ALT + logAREA + YRISOL + logDIST + logLDIST,
                     data = bird)
anova(mod2_reordered)

Analysis of Variance Table

Response: ABUND
          Df Sum Sq Mean Sq F value        Pr(>F)    
ALT        1    944     944   22.17 0.00002017033 ***
logAREA    1   2755    2755   64.74 0.00000000014 ***
YRISOL     1    502     502   11.81        0.0012 ** 
logDIST    1      4       4    0.09        0.7698    
logLDIST   1      5       5    0.12        0.7290    
Residuals 50   2128      43                          
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Поочерёдное тестирование Anova(модель)

library(car)
Anova(mod2)

Anova Table (Type II tests)

Response: ABUND
          Sum Sq Df F value     Pr(>F)    
logAREA     1614  1   37.92 0.00000012 ***
YRISOL       437  1   10.26     0.0024 ** 
logDIST        0  1    0.01     0.9303    
logLDIST       5  1    0.12     0.7290    
ALT           78  1    1.84     0.1814    
Residuals   2128 50                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

При поочерёдном тестировании результат НЕ зависит от порядка предикторов

Anova(mod2_reordered)

Anova Table (Type II tests)

Response: ABUND
          Sum Sq Df F value     Pr(>F)    
ALT           78  1    1.84     0.1814    
logAREA     1614  1   37.92 0.00000012 ***
YRISOL       437  1   10.26     0.0024 ** 
logDIST        0  1    0.01     0.9303    
logLDIST       5  1    0.12     0.7290    
Residuals   2128 50                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F-тест позволяет упростить модель

Anova(mod2)

Anova Table (Type II tests)

Response: ABUND
          Sum Sq Df F value     Pr(>F)    
logAREA     1614  1   37.92 0.00000012 ***
YRISOL       437  1   10.26     0.0024 ** 
logDIST        0  1    0.01     0.9303    
logLDIST       5  1    0.12     0.7290    
ALT           78  1    1.84     0.1814    
Residuals   2128 50                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Если не все предикторы влияют, то, может, можно удалить лишние — то есть упростить модель?

Упрощение линейных моделей при помощи частного F-критерия

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

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

1. Подбираем полную модель.
2. Удаляем один предиктор (строим уменьшенную модель).
3. Тестируем отличие уменьшенной модели от полной.

Повторяем 2–3 для каждого из предикторов:

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

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

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

Что происходит на каждом шаге обратного пошагового алгоритма?

anova(модель_1, модель_2)

Подбираем полную модель, удаляем предикторы по одному, тестируем, ухудшилась ли модель. Для окончательного удаления на этом шаге выбираем предиктор, удаление которого меньше всего её ухудшает.

Но мы пойдём другим путём. Все эти действия может выполнить одна функция drop1().

Частный F-критерий при помощи drop1() (шаг 1)

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

drop1(mod2, test = "F")

Single term deletions

Model:
ABUND ~ logAREA + YRISOL + logDIST + logLDIST + ALT
         Df Sum of Sq  RSS AIC F value     Pr(>F)    
<none>                2128 216                       
logAREA   1      1614 3742 245   37.92 0.00000012 ***
YRISOL    1       437 2565 224   10.26     0.0024 ** 
logDIST   1         0 2128 214    0.01     0.9303    
logLDIST  1         5 2133 214    0.12     0.7290    
ALT       1        78 2206 216    1.84     0.1814    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

Тестируем предикторы (шаг 2)

Удаляем выбранный предиктор и повторяем тесты снова. И т.д.

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

Single term deletions

Model:
ABUND ~ logAREA + YRISOL + logLDIST + ALT
         Df Sum of Sq  RSS AIC F value      Pr(>F)    
<none>                2128 214                        
logAREA   1      1628 3757 244   39.02 0.000000084 ***
YRISOL    1       437 2565 222   10.48      0.0021 ** 
logLDIST  1         9 2137 212    0.20      0.6532    
ALT       1        81 2209 214    1.94      0.1701    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

Тестируем предикторы (шаг 3)

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

Single term deletions

Model:
ABUND ~ logAREA + YRISOL + ALT
        Df Sum of Sq  RSS AIC F value       Pr(>F)    
<none>               2137 212                         
logAREA  1      2112 4248 248   51.39 0.0000000027 ***
YRISOL   1       502 2639 222   12.23      0.00097 ***
ALT      1       123 2260 213    2.99      0.08978 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

Тестируем предикторы (шаг 4)

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

Single term deletions

Model:
ABUND ~ logAREA + YRISOL
        Df Sum of Sq  RSS AIC F value        Pr(>F)    
<none>               2260 213                          
logAREA  1      2473 4732 252    58.0 0.00000000046 ***
YRISOL   1       607 2867 224    14.2       0.00041 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

Итоговая модель

summary(mod5)

Call:
lm(formula = ABUND ~ logAREA + YRISOL, data = bird)

Residuals:
   Min     1Q Median     3Q    Max 
-14.96  -3.80  -0.01   4.04  16.72 

Coefficients:
             Estimate Std. Error t value      Pr(>|t|)    
(Intercept) -252.1952    69.5868   -3.62       0.00065 ***
logAREA        3.7318     0.4900    7.62 0.00000000046 ***
YRISOL         0.1352     0.0358    3.77       0.00041 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.53 on 53 degrees of freedom
Multiple R-squared:  0.643, Adjusted R-squared:  0.63 
F-statistic: 47.8 on 2 and 53 DF,  p-value: 1.35e-12

Задание

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

Решение: 1) График расстояния Кука

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

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

• Выбросов нет
• Гетерогенность дисперсии?
• Два наблюдения с очень большими предсказанными значениями и большими остатками. Хорошо бы проверить их.

Решение: 3) Графики остатков от предикторов

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

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

res_1 <- gg_resid + aes(x = logAREA)
res_2 <- gg_resid + aes(x = YRISOL)
res_3 <- gg_resid + aes(x = logDIST)
res_4 <- gg_resid + aes(x = logLDIST)
res_5 <- gg_resid + aes(x = ALT)
res_6 <- gg_resid + aes(x = GRAZE)

library(gridExtra)
grid.arrange(res_1, res_2, res_3, res_4,
             res_5, res_6, nrow = 2)

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

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

[1] 11 47

Описываем финальную модель

summary(mod5)

Call:
lm(formula = ABUND ~ logAREA + YRISOL, data = bird)

Residuals:
   Min     1Q Median     3Q    Max 
-14.96  -3.80  -0.01   4.04  16.72 

Coefficients:
             Estimate Std. Error t value      Pr(>|t|)    
(Intercept) -252.1952    69.5868   -3.62       0.00065 ***
logAREA        3.7318     0.4900    7.62 0.00000000046 ***
YRISOL         0.1352     0.0358    3.77       0.00041 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.53 on 53 degrees of freedom
Multiple R-squared:  0.643, Adjusted R-squared:  0.63 
F-statistic: 47.8 on 2 and 53 DF,  p-value: 1.35e-12

График предсказаний модели: один предиктор

\[\widehat{ABUND}_i = -252.20 + 3.73 \cdot logAREA_i + 0.14 \cdot YRISOL_i\]

Задание: постройте график предсказаний модели

Отобразите, как меняются значения зависимой переменной в зависимости от значений одного из предикторов, при условии, что другие предикторы принимают свои средние значения. Дополните этот код:

Код для графика предсказаний модели: один предиктор

\[\widehat{ABUND}_i = -252.20 + 3.73 \cdot logAREA_i + 0.14 \cdot YRISOL_i\]

# Искусственный датафрейм для предсказаний
MyData1 <- data.frame(
  logAREA = seq(min(bird$logAREA), max(bird$logAREA), length.out = 100),
  YRISOL = mean(bird$YRISOL))
# Предсказанные значения
Predictions1 <- predict(mod5, newdata = MyData1, interval = 'confidence')
MyData1 <- data.frame(MyData1, Predictions1)
# Обратная трансформация предиктора, который будем изображать
MyData1$AREA <- exp(MyData1$logAREA)
# График предсказаний модели
Pl_predict1 <- ggplot(MyData1, aes(x = AREA, y = fit)) +
  geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr)) +
  geom_line()
Pl_predict1

График предсказаний модели: два предиктора

\[\widehat{ABUND}_i = -252.20 + 3.73 \cdot logAREA_i + 0.14 \cdot YRISOL_i\]

Код для графика предсказаний модели: два предиктора

\[\widehat{ABUND}_i = -252.20 + 3.73 \cdot logAREA_i + 0.14 \cdot YRISOL_i\]

# Искусственный датафрейм для предсказаний
MyData2 <- expand.grid(
  logAREA = seq(min(bird$logAREA), max(bird$logAREA), length.out = 100),
  YRISOL = round(quantile(bird$YRISOL)))
# Предсказанные значения
Predictions2 <- predict(mod5, newdata = MyData2, interval = 'confidence')
MyData2 <- data.frame(MyData2, Predictions2)
# Обратная трансформация предиктора, который будем изображать
MyData2$AREA <- exp(MyData2$logAREA)
# Делаем год фактором
MyData2$YRISOL <- factor(MyData2$YRISOL)
# График предсказаний модели
Pl_predict2 <- ggplot(MyData2, aes(x = AREA, y = fit, group = YRISOL)) +
  geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr)) +
  geom_line(aes(colour = YRISOL)) +
  scale_colour_brewer(palette = 'Spectral')
Pl_predict2

Недостатки этой модели

\[\widehat{ABUND}_i = -252.20 + 3.73 \cdot logAREA_i + 0.14 \cdot YRISOL_i\]

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

Недостатки и ограничения методологии отбора моделей

• Разные способы отбора оптимальной модели могут приводить к разным результатам. Не важно, какой из способов выбрать, но важно сделать это заранее, до анализа, чтобы не поддаваться соблазну подгонять результаты.
• При отборе моделей приходится делать множество тестов, поэтому увеличивается риск ошибок I рода (подробнее на занятии про дисперсионный анализ). Чтобы хоть как-то себя обезопасить, лучше не включать в модель всё подряд — только самое необходимое.
• Оптимальную модель лучше применять для тестирования гипотез на независимых данных, иначе получаются заниженные стандартные ошибки и завышенные эффекты. Обычно эту проблему игнорируют, но если у вас есть возможность, лучше сделать хотя бы кросс-валидацию.
• Какой-то одной «лучшей» модели, как правило, не существует.

Take-home messages

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

Что почитать

• 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