08.2_model_selection_code.R

# ---
# title: "Сравнение линейных моделей"
# author: Марина Варфоломеева, Вадим Хайтов

# # Модель множественной линейной регрессии с прошлой лекции

# ## Пример: птицы в лесах Австралии
# От каких характеристик лесного участка зависит
# обилие птиц в лесах юго-западной Виктории,
# Австралия (Loyn, 1987)
# Переменных много, мы хотим из них выбрать
# __оптимальный небольшой__ набор.
# 56 лесных участков:
# - ABUND - обилие птиц
# - AREA - площадь участка
# - YRISOL - год изоляции участка
# - DIST - расстояние до ближайшего леса
# - LDIST - расстояние до ближайшего большого леса
# - GRAZE - пастбищная нагрузка (1-5)
# - ALT - высота над уровнем моря
# Пример из кн. 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)
# ABUND_i = - 226.00 + 3.69 logAREA_i + 0.12 YRISOL_i -
#                    - 0.10 logDIST_i - 0.33 logLDIST_i +
#                    + 0.03 ALT_i     + e_i


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




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

# ## Задание 1  ------------------------------------------------
#
# Для тренировки запишем вложенные модели для данной полной модели
#
# (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
# (2)
# (3)
# и т.д.


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



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

# F_{ALT | logAREA, YRISOL, logDIST, logLDIST}
mod2_reduced <- update(mod2, . ~ . - ALT)
anova(mod2_reduced, mod2)


# ## Результаты F- и t-тестов эквивалентны
# t^2 = F
F_val <- anova(mod2, mod2_reduced)[2, 'F']
F_val
t_val <- coef(summary(mod2))['ALT', 't value']
t_val^2


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

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


# # I и II типы сумм квадратов
# `Y ~ A + B + C`
# ### I тип сумм квадратов (SS type I)
# - SS_{A}
# - SS_{B|A}
# - SS_{C|A, B}
# ### II тип сумм квадратов (SS type II)
# - SS_{A|B, C}
# - SS_{B|A, C}
# - SS_{C|A, B}

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

# ## Осторожно, при последовательном тестировании результат зависит от порядка предикторов
mod2_reordered <- lm(ABUND ~ ALT + logAREA + YRISOL + logDIST + logLDIST,
                     data = bird)
anova(mod2_reordered)

# ## Поочередное тестирование `Anova(модель)` или `drop1()`
library(car)
Anova(mod2)
# drop1(mod2, test = 'F')

# ## При поочередном тестировании результат НЕ зависит от порядка предикторов
Anova(mod2_reordered)



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


# ## Что происходит на каждом шаге обратного пошагового алгоритма?
# `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)




# ##  Частный F-критерий при помощи `drop1()` (шаг 1)
# Тестируем значимость всех предикторов за один раз. Затем выбираем предиктор, удаление которого меньше всего ухудшает модель.
drop1(mod2, test = "F")

# ## Тестируем предикторы (шаг 2)
# Удаляем выбранный предиктор и повторяем тесты снова. И т.д.
# Убрали logDIST
mod3 <- update(mod2, . ~ . - logDIST)
drop1(mod3, test = "F")
# Нужно убрать logLDIST

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

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

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



# ## Задание 2 -----------------------------------------------
#
# Проверьте финальную модель на выполнение условий
# применимости. Дополните код
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



# ## Описываем финальную модель
summary(mod5)



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

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

# Искуственный датафрейм для предсказаний
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


# ## График предсказаний модели: два предиктора
# Искуственный датафрейм для предсказаний
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


# ## В чем недостатки этой модели?



## Задание 4 ------------------------------------------------
# В датасете cystfibr из пакета ISwR находится
# информация характеристиках дыхательной системы
# больных муковисцидозом (в возрасте от 7 до 23
# лет).
# Первые два действия вы уже сделали на прошлом занятии.
# Продолжите с 3-го шага.
# ---
# 1) Постройте модель описывающую связь между усилием мышц,
# осуществляющих выдох (`pemax`) и следующими переменными:
# - `age`    --- Возраст
# - `height` --- Рост (см)
# - `weight` --- Вес (кг)
# - `bmp`    --- Отклонения в весе от нормы (% от нормы)
# - `fev1`   --- Объем наполненных легких
# - `rv`     --- Остаточный объем легких
# - `frc`    --- Функциональная остаточная емкость легких
# - `tlc`    --- Общая емкость легких
# Не учитывайте влияние пола `sex` (мы пока не
# умеем работать с дискретными переменными, о них
# позже)
# 2) Исключите из модели коллинеарные предикторы.
# ---
# 3) Проведите оставшуюся диагностику модели (это
# нужно сделать уже сейчас, чтобы понять, стоит ли
# дальше возиться с этой моделью)
# 4) Попробуйте упростить модель при помощи частного F-критерия.
# 5) Проведите диагностику упрощенной модели.
# 6) Нарисуйте график предсказаний финальной модели.

library(ISwR)
data(cystfibr)