От каких характеристик лесного участка зависит обилие птиц в лесах юго-западной Виктории, Австралия (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.1 ## 2 50 2127.9 1 78.19 1.8372 0.1814
Если высота над уровнем моря — единственный предиктор в модели, то его влияние статистически значимо (\(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 6337.9 ## 2 54 5394.4 1 943.52 9.445 0.003316 ** ## --- ## 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.0 3471.0 81.5580 4.42e-12 *** ## YRISOL 1 607.4 607.4 14.2709 0.0004221 *** ## logDIST 1 28.5 28.5 0.6687 0.4173708 ## logLDIST 1 25.0 25.0 0.5878 0.4468656 ## ALT 1 78.2 78.2 1.8372 0.1813659 ## Residuals 50 2127.9 42.6 ## --- ## 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 943.52 943.52 22.1701 0.0000201703257 *** ## logAREA 1 2755.14 2755.14 64.7378 0.0000000001412 *** ## YRISOL 1 502.48 502.48 11.8069 0.001196 ** ## logDIST 1 3.68 3.68 0.0866 0.769786 ## logLDIST 1 5.17 5.17 0.1214 0.728963 ## Residuals 50 2127.92 42.56 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(модель)library(car)
## Загрузка требуемого пакета: carData
## ## Присоединяю пакет: 'car'
## Следующий объект скрыт от 'package:dplyr': ## ## recode
Anova(mod2)
## Anova Table (Type II tests) ## ## Response: ABUND ## Sum Sq Df F value Pr(>F) ## logAREA 1613.84 1 37.9206 0.0000001243 *** ## YRISOL 436.58 1 10.2584 0.002368 ** ## logDIST 0.33 1 0.0077 0.930318 ## logLDIST 5.17 1 0.1214 0.728963 ## ALT 78.19 1 1.8372 0.181366 ## Residuals 2127.92 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.19 1 1.8372 0.181366 ## logAREA 1613.84 1 37.9206 0.0000001243 *** ## YRISOL 436.58 1 10.2584 0.002368 ** ## logDIST 0.33 1 0.0077 0.930318 ## logLDIST 5.17 1 0.1214 0.728963 ## Residuals 2127.92 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 1613.84 1 37.9206 0.0000001243 *** ## YRISOL 436.58 1 10.2584 0.002368 ** ## logDIST 0.33 1 0.0077 0.930318 ## logLDIST 5.17 1 0.1214 0.728963 ## ALT 78.19 1 1.8372 0.181366 ## Residuals 2127.92 50 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Если не все предикторы влияют, то может можно удалить лишние, то есть упростить модель?
Постепенно удаляем предикторы. Модели обязательно должны быть вложенными! *
summary(mod2)
## ## Call: ## lm(formula = ABUND ~ logAREA + YRISOL + logDIST + logLDIST + ## ALT, data = bird) ## ## Residuals: ## Min 1Q Median 3Q Max ## -15.4386 -3.6686 0.3315 2.7647 15.7594 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -226.00065 74.25257 -3.044 0.00372 ** ## logAREA 3.68816 0.59892 6.158 0.000000124 *** ## YRISOL 0.12073 0.03770 3.203 0.00237 ** ## logDIST -0.10335 1.17593 -0.088 0.93032 ## logLDIST -0.32814 0.94171 -0.348 0.72896 ## ALT 0.03180 0.02346 1.355 0.18137 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 6.524 on 50 degrees of freedom ## Multiple R-squared: 0.6643, Adjusted R-squared: 0.6307 ## F-statistic: 19.78 on 5 and 50 DF, p-value: 7.968e-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> 2127.9 215.70 ## logAREA 1 1613.84 3741.8 245.31 37.9206 0.0000001243 *** ## YRISOL 1 436.58 2564.5 224.15 10.2584 0.002368 ** ## logDIST 1 0.33 2128.3 213.71 0.0077 0.930318 ## logLDIST 1 5.17 2133.1 213.84 0.1214 0.728963 ## ALT 1 78.19 2206.1 215.72 1.8372 0.181366 ## --- ## 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.3 213.71 ## logAREA 1 1628.32 3756.6 243.53 39.0200 0.0000000841 *** ## YRISOL 1 437.22 2565.5 222.18 10.4772 0.002126 ** ## logLDIST 1 8.52 2136.8 211.94 0.2043 0.653227 ## ALT 1 80.79 2209.0 213.80 1.9360 0.170141 ## --- ## 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> 2136.8 211.94 ## logAREA 1 2111.53 4248.3 248.42 51.3856 0.000000002667 *** ## YRISOL 1 502.48 2639.3 221.76 12.2283 0.0009725 *** ## ALT 1 122.82 2259.6 213.06 2.9888 0.0897756 . ## --- ## 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> 2259.6 213.06 ## logAREA 1 2472.50 4732.1 252.46 57.994 0.000000000462 *** ## YRISOL 1 607.35 2866.9 224.40 14.246 0.0004067 *** ## --- ## 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.9560 -3.7974 -0.0105 4.0383 16.7191 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -252.19516 69.58676 -3.624 0.000650 *** ## logAREA 3.73177 0.49003 7.615 0.000000000462 *** ## YRISOL 0.13525 0.03583 3.774 0.000407 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 6.529 on 53 degrees of freedom ## Multiple R-squared: 0.6435, Adjusted R-squared: 0.63 ## F-statistic: 47.83 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.9560 -3.7974 -0.0105 4.0383 16.7191 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -252.19516 69.58676 -3.624 0.000650 *** ## logAREA 3.73177 0.49003 7.615 0.000000000462 *** ## YRISOL 0.13525 0.03583 3.774 0.000407 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 6.529 on 53 degrees of freedom ## Multiple R-squared: 0.6435, Adjusted R-squared: 0.63 ## F-statistic: 47.83 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