Вы сможете

• Подобрать модель множественной линейной регрессии
• Протестировать ее статистическую значимость и валидность

Фрагментация лесных местообитаний - одна из важнейших проблем Австралии. От каких характеристик лесного участка зависит обилие птиц во фрагментированных лесных массивах? (Loyn, 1987)

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

bird <- read.csv("data/loyn.csv")

Все ли правильно открылось?

str(bird)

## 'data.frame':    56 obs. of  7 variables:
##  $ ABUND : num  5.3 2 1.5 17.1 13.8 14.1 3.8 2.2 3.3 3 ...
##  $ AREA  : num  0.1 0.5 0.5 1 1 1 1 1 1 1 ...
##  $ YRISOL: int  1968 1920 1900 1966 1918 1965 1955 1920 1965 1900 ...
##  $ DIST  : int  39 234 104 66 246 234 467 284 156 311 ...
##  $ LDIST : int  39 234 311 66 246 285 467 1829 156 571 ...
##  $ GRAZE : int  2 5 5 3 5 3 5 5 4 5 ...
##  $ ALT   : int  160 60 140 160 140 130 90 60 130 130 ...

Есть ли пропущенные значения?

colSums(is.na(bird))

##  ABUND   AREA YRISOL   DIST  LDIST  GRAZE    ALT 
##      0      0      0      0      0      0      0

cor(bird)

##          ABUND     AREA   YRISOL   DIST   LDIST  GRAZE    ALT
## ABUND   1.0000  0.25597  0.50336  0.236  0.0872 -0.683  0.386
## AREA    0.2560  1.00000 -0.00149  0.108  0.0346 -0.310  0.388
## YRISOL  0.5034 -0.00149  1.00000  0.113 -0.0833 -0.636  0.233
## DIST    0.2361  0.10834  0.11322  1.000  0.3172 -0.256 -0.110
## LDIST   0.0872  0.03458 -0.08332  0.317  1.0000 -0.028 -0.306
## GRAZE  -0.6825 -0.31040 -0.63557 -0.256 -0.0280  1.000 -0.407
## ALT     0.3858  0.38775  0.23272 -0.110 -0.3060 -0.407  1.000

• Нет

• Обычная корреляция не учитывает, что взаимосвязь между переменными может находиться под контролем других переменных и их взаимодействий.
• Множественные тесты. При тестировании значимости множества коэффициентов корреляции нужно вводить поправку для уровня значимости. Лучше было бы учесть все в одном анализе.

\[y_i = b_0 + b_1x_{1i} + b_2x_{2i} + b_3x_{3i} + ... + b_{p - 1}x_{p - 1\;i} + e_i\]

• \(y_i\) — значение зависимой переменной для \(i\)-того наблюдения
• \(b_0\) — свободный член (intercept). Значение \(y\) при \(x_1 = x_2 = x_3= \ldots = x_{p - 1} = 0\)
• \(b_1\) — частный угловой коэффициент для зависимости \(y\) от \(x_1\). Показывает, на сколько единиц изменяется \(y\) при изменении \(x_1\) на одну единицу и при условии, что все остальные предикторы не изменяются \(b_2\), \(b_3\), …., \(b_{p - 1}\) — аналогично
• \(e_i\) — варьирование \(y\), не объясненное данной моделью
• \(p\) — число параметров модели

Для случая с одним предиктором \(y_i = b_0 + b_1x_{1i} + e_i\) — линия регрессии

## `geom_smooth()` using formula 'y ~ x'

Для случая с двумя предикторами \(y_i = b_0 + b_1x_{1i} + b_2x_{2i} + e_i\) — плоскость в трехмерном пространстве

Для случая с большим количеством предикторов

\[y_i = b_0 + b_1x_{1i} + b_2x_{2i} + b_3x_{3i} + ... + b_{p - 1}x_{p - 1\;i} + e_i\]

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

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

Главный принцип — “garbage in, garbage out”.

Сбор данных часто проходит по принципу “ДДПР” (“Давай-давай, потом разберемся”), аналитик может не знать подводных камней сбора материала, но он должен их выявить.

Данные могут не отвечать многочисленным ограничениям регрессионного анализа.

Характер данных и связей между ними может требовать доработки данных (преобразования, стандартизации и т.п.).

Но! Самое главное гипотеза должна быть сформулирована до начала анализа, точнее, до начала сбора материала.

Часть осмысления данных необходимо провести до построения модели (разведочный анализ).

Часть — после того как модель построена (анализ валидности модели).

Давайте подробнее рассмотрим тактику анализа данных.

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

• Нет ли скрытых группирующих (случайных) факторов
• Нет ли временных и пространственных автокорреляций

Нужно проверить, нет ли отскакивающих значений.

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

Простая линейная регрессия часто не подходит для моделирования некоторых величин (например для счетных данных).

Связь между откликом и предикторами может выглядеть нелинейной.

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

Не все зависимости можно моделировать с помощью простых регрессионных моделей.

Нужно проверить, нет ли взаимозависимости (коллинеарности) между предикторами.

О том, почему это важно и как это проверять, поговорим позже.

Строим модель, соответствующую нашей гипотезе.

Для моделей с нормальным распределением отклика:

1. Независимы ли наблюдения друг от друга (проверяем еще раз, т.к. речь идет уже не о дизайне, а о модели).
2. Присутствует ли гетерогенность дисперсии
   
3. Соответствует ли распределение остатков модели нормальному распределению

Если нарушений условий применимости не выявлено, то можно начать осмыслять результаты построения модели.

Если выявлены нарушения условий применимости, то надо задуматься о том, верно ли подобран тип модели и все ли хорошо с данными. Возвращаемся к пункту 1.

• Большая часть значений AREA, DIST, LDISТ сгруппирована в начале области определения. Связь между AREA и откликом выглядит нелинейной. Нужно логарифмировать эти переменные.

• Переменная GRAZE — это уровень выпаса скота в баллах, ее лучше было бы анализировать как дискретную переменную. Технически, ее можно анализировать и как непрерывную, но нужно помнить, это предполагает одинаковые различия между разными соседними уровнями выпаса скота. Это нереалистично.

Трансформируем переменные

bird$logAREA <- log(bird$AREA)
bird$logDIST <- log(bird$DIST)
bird$logLDIST <- log(bird$LDIST)

Отскоки (выбросы, outliers) - наблюдения, которые имеют более высокие (или низкие) значения относительно большинства других наблюдений.

Код для графика

gg_dot <- ggplot(bird, aes(y = 1:nrow(bird))) + geom_point() + ylab('index')
Pl1 <- gg_dot + aes(x = ABUND)
Pl2 <- gg_dot + aes(x = YRISOL)
Pl3 <- gg_dot + aes(x = logAREA)
Pl4 <- gg_dot + aes(x = logDIST)
Pl5 <- gg_dot + aes(x = logLDIST)
Pl6 <- gg_dot + aes(x = ALT)
Pl7 <- gg_dot + aes(x = GRAZE)

library(cowplot) # пакет для группировки графиков
theme_set(theme_bw())
plot_grid(Pl1, Pl2, Pl3, Pl4, Pl5, Pl6, 
          Pl7, ncol = 3, nrow = 3)

\[\begin{aligned}{ABUND}_i &= b_0 + b_1 \cdot logAREA_i + b_2 \cdot YRISOL_i + \\ &+ b_3 \cdot logDIST_i + b_4 \cdot logLDIST_i + b_5 \cdot GRAZE_i + b_6 \cdot ALT_i + e_i\\ \end{aligned}\]

Возможно, что эту модель придется изменить.

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

• Линейная связь между зависимой переменной (\(Y\)) и предикторами (\(X\))
• Независимость значений \(Y\) друг от друга
• Нормальное распределение \(Y\) для каждого уровня значений \(X\)
• Гомогенность дисперсий \(Y\) для каждого уровня значений \(X\)
• Отсутствие коллинеарности предикторов (для можественной регрессии)

Мультиколлинеарность — наличие линейной зависимости между независимыми переменными (предикторами) в регрессионной модели.

При наличии мультиколлинеарности оценки параметров неточны, а значит сложно интерпретировать влияние предикторов на отклик.

Косвенные признаки мультиколлинеарности:

• Большие ошибки оценок параметров
  
• Большинство значений параметров модели незначимо отличается от нуля, но F критерий говорит, что модель в целом значима.

Проверка на мультиколлинеарность

• Коэффициент раздутия дисперсии (Variance inflation factor, VIF).

Пусть наша модель \(y_i = b_0 + b_1 x_{1i} + b_2 x_{2i} + \ldots + b_{p - 1} x_{p - 1\;i} + e_i\).

Нужно оценить какую долю изменчивости конкретного предиктора могут объяснить другие предикторы (т.е. насколько предикторы независимы).

Для каждого предиктора:

1. Строим регрессионную модель данного предиктора от всех остальных: \[x_1 = c_0 + c_1x_2 +c_2x_3 + .... + c_{p - 2}x_{p - 1}\]
2. Находим \(R^2\) этой модели.
3. Вычисляем коэффициент раздутия дисперсии: \[VIF = \frac{1}{1-R^2}\]

В случае наличия мультиколлинеарности:

• Оценки коэффициентов модели нестабильны (даже могут менять знак при небольших изменениях модели или исходных данных).
• Стандартные ошибки оценок параметров увеличатся в \(\sqrt{VIF}\) раз.
• В результате меньше шансов заметить влияние предиктора, т.к. уровень значимости (p-value) в тестах будет выше.

• Можно последовательно удалить из модели избыточные предикторы с VIF > 2
  1. подбираем модель
  2. считаем VIF
  3. удаляем предиктор с самым большим VIF
  4. повторяем 1-3
• Можно заменить исходные предикторы новыми независимыми друг от друга переменными, сконструированными методом главных компонент (Principal component analysis, PCA).

• Постройте множественную линейную регрессию для зависимости обилия птиц (ABUND) от других переменных (logAREA, YRISOL, logDIST, logLDIST, GRAZE, ALT)

\[\begin{aligned}{ABUND}_i &= b_0 + b_1 \cdot logAREA_i + b_2 \cdot YRISOL_i + \\ &+ b_3 \cdot logDIST_i + b_4 \cdot logLDIST_i + b_5 \cdot GRAZE_i + b_6 \cdot ALT_i + e_i\\ \end{aligned}\]

• Используйте функцию vif(), чтобы проверить, коллинеарны ли предикторы.

Дополните код:

mod1 <- lm(formula = , data = )
vif()

# Строим модель
mod1 <- lm(formula = ABUND ~  logAREA + YRISOL + logDIST + logLDIST + GRAZE + ALT, 
           data = bird)
# Проверяем, есть ли коллинеарность?
vif(mod1)

##  logAREA   YRISOL  logDIST logLDIST    GRAZE      ALT 
##     1.91     1.80     1.65     2.01     2.52     1.47

В нашей модели сильной мультиколлинеарности нет.

Однако, возможно GRAZE — избыточный предиктор.

mod2 <- update(mod1, . ~ . -GRAZE)
vif(mod2)

##  logAREA   YRISOL  logDIST logLDIST      ALT 
##     1.62     1.20     1.62     2.01     1.35

Теперь мультиколлинеарности нет. В модели осталось пять предикторов (и шесть параметров).

\[\begin{aligned}{ABUND}_i &= b_0 + b_1 \cdot logAREA_i + b_2 \cdot YRISOL_i + \\ &+ b_3 \cdot logDIST_i + b_4 \cdot logLDIST_i + b_5 \cdot ALT_i + e_i\\ \end{aligned}\]

\[\begin{aligned}{ABUND}_i &= b_0 + b_1 \cdot logAREA_i + b_2 \cdot YRISOL_i + \\ &+ b_3 \cdot logDIST_i + b_4 \cdot logLDIST_i + b_5 \cdot ALT_i + e_i\\ \end{aligned}\]

Мы подобрали коэффициенты и можем записать уравнение модели.

coef(mod2)

## (Intercept)     logAREA      YRISOL     logDIST    logLDIST 
##   -226.0006      3.6882      0.1207     -0.1033     -0.3281 
##         ALT 
##      0.0318

\[\begin{aligned}{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 + e_i\\ \end{aligned}\]

Аналогично простой регрессии

\[ \mathbf{Y} = \mathbf{Xb} + \mathbf{e} \]

Отличие лишь в форме модельной матрицы

X <- model.matrix(mod2)

head(X)

##   (Intercept) logAREA YRISOL logDIST logLDIST ALT
## 1           1  -2.303   1968    3.66     3.66 160
## 2           1  -0.693   1920    5.46     5.46  60
## 3           1  -0.693   1900    4.64     5.74 140
## 4           1   0.000   1966    4.19     4.19 160
## 5           1   0.000   1918    5.51     5.51 140
## 6           1   0.000   1965    5.46     5.65 130

• Оставить все как есть. Если значение коэффициента при предикторе не значимо отличается от нуля, значит, этот предиктор не влияет на обилие птиц.
• Провести пошаговый подбор оптимальной модели (Об этом на следующей лекции).

Сейчас оставим все как есть.

Теперь давайте проверим, выполняются ли оставшиеся условия применимости, а после этого попытаемся выяснить, какие предикторы влияют сильнее всего.

Проверьте, выполняются ли условия применимости для модели mod2. Дополните код:

library()
mod2_diag <- data.frame(fortify(), $GRAZE)
# 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(gridExtra)
grid.arrange(res_1, res_2, nrow = 2)
# 4) Квантильный график остатков
library(car)
qq

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

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

library(ggplot2)
mod2_diag <- data.frame(fortify(mod2), GRAZE = bird$GRAZE)

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

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

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

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

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 = GRAZE)
res_6 <- gg_resid + aes(x = ALT)

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

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

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

library(car)
qqPlot(mod2)

## [1] 18 47

Для ответа на этот вопрос надо “уравнять” шкалы, всех предикторов, то есть стандартизовать их.

Коэффициенты при стандартизованных предикторах покажут, насколько сильно меняется отклик при изменении предиктора на одно стандартное отклонение.

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

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

coef(summary(mod2_scaled))

##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)      19.5143      0.872 22.3848 1.03e-27
## scale(logAREA)    6.8992      1.120  6.1580 1.24e-07
## scale(YRISOL)     3.0885      0.964  3.2029 2.37e-03
## scale(logDIST)   -0.0985      1.120 -0.0879 9.30e-01
## scale(logLDIST)  -0.4344      1.247 -0.3485 7.29e-01
## scale(ALT)        1.3842      1.021  1.3554 1.81e-01

• Сильнее всего на обилие птиц влияют логарифм площади леса и продолжительность изоляции
• При изменении логарифма площади на 1 стандартное отклонение, обилие птиц изменяется на 6.9
• При изменении продолжительности изоляции на 1 стандартное отклонение, обилие птиц изменяется на 3.09

Для простой линейной регрессии легко нарисовать график на плоскости, поскольку есть только две переменные: отклик и предиктор.

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

Самое частое решение — построить график \(y\) от наиболее важного предиктора при средних значениях всех остальных предикторов. Но можно рассмотреть и любые другие интересные вам сценарии.

\[\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}\]

coef(mod2_scaled)

##     (Intercept)  scale(logAREA)   scale(YRISOL)  scale(logDIST) 
##         19.5143          6.8992          3.0885         -0.0985 
## scale(logLDIST)      scale(ALT) 
##         -0.4344          1.3842

Судя по стандартизованным коэффициентам, самая важная здесь переменная — это логарифм площади леса.

Построим график предсказанных значений обилия птиц для лесов разной площади при средних значениях всех остальных предикторов.

\[\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}\]

В нашей модели предикторы трансформированы и ABUND от них зависит линейно.

# Искусственный датафрейм для предсказаний
MyData <- data.frame(
  logAREA = seq(min(bird$logAREA), max(bird$logAREA), length.out = 100),
  YRISOL = mean(bird$YRISOL),
  logDIST = mean(bird$logDIST),
  logLDIST = mean(bird$logLDIST),
  ALT = mean(bird$ALT))
# Предсказанные значения
Predictions <- predict(mod2, newdata = MyData,  interval = 'confidence')
MyData <- data.frame(MyData, Predictions)
# График предсказаний модели
Pl_predict <- ggplot(MyData, aes(x = logAREA, y = fit)) +
  geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr)) +
  geom_line()
Pl_predict

\[\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}\]

Для удобства восприятия лучше сделать обратную трансформацию предикторов.

# Обратная трансформация предиктора, который будем изображать
MyData$AREA <- exp(MyData$logAREA)
# График предсказаний модели
Pl_predict <- ggplot(MyData, aes(x = AREA, y = fit)) +
  geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr)) +
  geom_line() + xlab("Площадь леса") + ylab("Обилие птиц")
Pl_predict

Имеет ли смысл на таком графике изображать исходные наблюдения?

• Обычно это делают, чтобы (1) посмотреть на величину остатков, (2) посмотреть на диапазон исходных наблюдений.

• Нет, нет смысла. Наш график показывает предсказания только для определенных наблюдений — со средними значениями нескольких предикторов. Если мы нарисуем значения наблюдений, то их остатки не будут иметь смысл, т.к. у этих наблюдений в реальности совсем не обязательно средние значения этих предикторов.

• Записываем уравнение модели.
• Общую значимость модели оцениваем при помощи F-критерия.
• Качество подгонки модели описываем при помощи коэффициента детерминации с поправкой (\(R^2_{adj.}\)).
• При обсуждении значимости отдельных предикторов можно привести таблицу с оценками коэффициентов и тестами их значимости.
• При сравнении влияния отдельных предикторов приводим стандартизованные коэффициенты
• Приводим график предсказаний модели

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

• Гланц, С., 1998. Медико-биологическая статистика. М., Практика
• Кабаков Р.И. R в действии. Анализ и визуализация данных на языке R. М.: ДМК Пресс, 2014
• Diez, D.M., Barr, C.D. and Çetinkaya-Rundel, M., 2015. OpenIntro Statistics. OpenIntro.
• 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