От каких характеристик лесного участка зависит обилие птиц в лесах юго-западной Виктории, Австралия (Loyn, 1987)
Переменных много, мы хотим из них выбрать оптимальный небольшой набор.
Mystic Forest - Warburton, Victoria by ¡kuba! on flickr
56 лесных участков:
Пример из кн. 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}\]
Две модели являются вложенными, если одну из них можно получить из другой путем удаления некоторых предикторов.
Удаление предиктора - коэффициент при данном предикторе равен нулю.
М1: \(y _i = b _0 + b _1 x _{1i} + b _2 x _{2i} + e _i\)
М2: \(y _i = b _0 + b _1 x _{1i} + e _i\)
М3: \(y _i = b _0 + b _2 x _{2i} + e _i\)
M2 вложена в M1
M3 вложена в M1
M2 и M3 не вложены друг в друга
\(y _i = b _0 + e _i\)
Для тренировки запишем вложенные модели для данной полной модели
(1)\(y _i = b _0 + b _1 x _{1i} + b _2 x _{2i} + b _3 x _{3i} + e _i\)
Для тренировки запишем вложенные модели для данной полной модели
(1)\(y _i = b _0 + b _1 x _{1i} + b _2 x _{2i} + b _3 x _{3i} + e _i\)
Модели:
Вложенность:
\(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\)
\(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 = \frac {(SS_{e, reduced} - SS _{e, full})/(df _{e, reduced} - df _{e, full})}{SS _{e, full} / df _{e, full}}\]
Если модель значимо ухудшается от удаления предиктора, то влияние этого предиктора значимо.
Чтобы проверить значимость предиктора ALT, нам нужно сравнить две модели
ABUND ~ logAREA + YRISOL + logDIST + logLDIST + ALTABUND ~ logAREA + YRISOL + logDIST + logLDIST
Это можно обозначить
\(F_{(\text{ALT | logAREA, YRISOL, logDIST, logLDIST})}\)
А разницу остаточных сумм квадратов \(SS_{e, reduced} - SS _{e, full}\), которую мы тестируем, можно обозначить
\(SS_{(\text{ALT | logAREA, YRISOL, logDIST, logLDIST})}\)
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
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
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
Если не все предикторы влияют, то может можно удалить лишние, то есть упростить модель?
“Essentially, all models are wrong,
but some are useful”
Georg E. P. Box
“Entia non sunt multiplicanda praeter necessitatem”
Gulielmus Occamus
Since all models are wrong the scientist cannot obtain a “correct” one by excessive elaboration. On the contrary following William of Occam he should seek an economical description of natural phenomena. […] overelaboration and overparameterization is often the mark of mediocrity (Box, 1976).
Поскольку все модели ошибочны, ученый не может получить «правильную» модель даже если очень постарается. Напротив, вслед за Уильямом Оккамом он должен искать экономичное описание природы. […] чрезмерное усложнение модели часто является признаком посредственности. (Box, 1976).
While a model can never be “truth,” a model might be ranked from very useful, to useful, to somewhat useful to, finally, essentially useless (Burnham & Anderson, 2002).
Хотя модель никогда не может быть «правдой», модели можно ранжировать от очень полезных к полезным, до некоторой степени полезным и , наконец, к абсолютно бесполезным (Burnham & Anderson, 2002).
На этих графиках показаны предсказания loess-регрессии для одних и тех же исходных данных.
Cложность модели — в общем случае, это число параметров. Для loess-регрессии сложность модели отражает степень сглаживания: у более сложных моделей маленькая степень сглаживания.
На новых данных предсказания моделей не будут идеальными.
Ошибка предсказаний на новых данных практически всегда больше, чем на исходных данных. Более сложные модели лучше описывают существующие данные, но на новых данных их предсказания хуже.
Обычно при усложнении модели:
Предсказания, сделанные на новых данных, будут отличаться от истинных значений не только из-за погрешности или неточности. Еще один источник отличий — это так называемая неснижаемая ошибка.
\[Полная~ошибка = Дисперсия + (Погрешность)^2 + Неснижаемая~ошибка\] При увеличении сложности модели снижается погрешность предсказаний, но возрастает их разброс. Поэтому общая ошибка предсказаний велика у недообученных или переобученных моделей, а у моделей средней сложности она будет минимальной.
Объяснение закономерностей, описание функциональной зависимости
Предсказание значений зависимой переменной
Хорошая модель должна соответствовать условиям применимости, иначе вы не сможете доверять результатам тестов.
Другие соображения: разумность, целесообразность модели, простота, ценность выводов, важность предикторов.
Постепенно удаляем предикторы. Модели обязательно должны быть вложенными! *
summary(mod2)
## ## Call: ## lm(formula = ABUND ~ logAREA + YRISOL + logDIST + logLDIST + ## ALT, data = bird) ## ## Residuals: ## Min 1Q Median 3Q Max ## -15.439 -3.669 0.332 2.765 15.759 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -226.0006 74.2526 -3.04 0.0037 ** ## logAREA 3.6882 0.5989 6.16 0.00000012 *** ## YRISOL 0.1207 0.0377 3.20 0.0024 ** ## logDIST -0.1033 1.1759 -0.09 0.9303 ## logLDIST -0.3281 0.9417 -0.35 0.7290 ## ALT 0.0318 0.0235 1.36 0.1814 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 6.52 on 50 degrees of freedom ## Multiple R-squared: 0.664, Adjusted R-squared: 0.631 ## F-statistic: 19.8 on 5 and 50 DF, p-value: 7.97e-11
anova(модель_1, модель_2)
Подбираем полную модель, удаляем предикторы по-одному, тестируем ухудшилась ли модель. Для окончательного удаления на этом шаге выбираем предиктор, удаление которого меньше всего ее ухудшает.
mod3 <- update(mod2, . ~ . - logAREA) anova(mod2, mod3) mod4 <- update(mod2, . ~ . - YRISOL) anova(mod2, mod4) mod5 <- update(mod2, . ~ . - logDIST) anova(mod2, mod5) mod6 <- update(mod2, . ~ . - logLDIST) anova(mod2, mod6) mod7 <- update(mod2, . ~ . - ALT) anova(mod2, mod7)
drop1()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
Удаляем выбранный предиктор и повторяем тесты снова. И т.д.
# Убрали 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
# Убрали 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
# Убрали 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
Проверьте финальную модель на выполнение условий применимости. Дополните код
library() mod_diag <- data.frame(fortify(), bird[, c()]) # 1) График расстояния Кука ggplot(data = , aes(x = 1:, y = .cooksd)) + geom_bar(stat = "") # 2) График остатков от предсказанных значений gg_resid <- ggplot(data = , aes(x = , y = )) + geom_point() + geom_hline() gg_resid # 3) Графики остатков от предикторов в модели и нет res_1 <- gg_resid + aes(x = logAREA) res_1 res_2 <- gg_resid res_3 <- gg_resid res_4 <- gg_resid res_5 <- gg_resid res_6 <- gg_resid # все графики вместе library() grid.arrange(res_1, res_2, nrow = 2) # 4) Квантильный график остатков library() qq
mod5_diag <- data.frame(
fortify(mod5),
bird[, c("logDIST", "logLDIST", "GRAZE", "ALT")]
)
ggplot(mod5_diag, aes(x = 1:nrow(mod5_diag), y = .cooksd)) +
geom_bar(stat = "identity")
gg_resid <- ggplot(data = mod5_diag, aes(x = .fitted, y = .stdresid)) + geom_point() + geom_hline(yintercept = 0) gg_resid
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)
library(car) qqPlot(mod5)
## [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\]
Постройте график предсказаний модели. Отобразите, как меняются значения зависимой переменной в зависимости от значений одного из предикторов, при условии, что другие предикторы принимают свои средние значения. Дополните этот код:
# Искусственный датафрейм для предсказаний MyData1 <- data.( = seq(min( ), max( ), ), YRISOL = ) # Предсказанные значения Predictions1 <- predict( , newdata = , interval = ) MyData1 <- data.frame(MyData1, ) # Обратная трансформация предиктора, который будем изображать MyData1$AREA <- # График предсказаний модели Pl_predict1 <- ggplot( , aes(x = AREA, y = )) + geom_ (alpha = 0.2, aes(ymin = , ymax = )) + geom_ () Pl_predict1
\[\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 рода (подробнее на занятии про дисперсионный анализ). Чтобы хоть как-то себя обезопасить, лучше не включать в модель все подряд — только самое необходимое.
Оптимальную модель лучше применять для тестирования гипотез на независимых данных, иначе получаются заниженные стандартные ошибки и завышенные эффекты. Обычно эту проблему игнорируют, но если у вас есть возможность, лучше сделать хотя бы кросс-валидацию.
Какой-то одной “лучшей” модели, как правило, не существует.
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