Тяжелая жизнь статистического аналитика

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

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

Откуда ждать удара?

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

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

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

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

Как справится со сложностями?

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

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

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

1. Дизайн сбора материала

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

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

1. Дизайн сбора материала

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

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

Слева — стандартный график остатков. Выглядит как случайный шум, ничего не заподозришь. Справа — те же остатки, но раскрашенные по ID растения. Видно, что остатки одного растения «кучкуются»: одни растения дают в основном положительные остатки, другие — отрицательные. Это и есть сигнал скрытой группировки.

2. Влиятельные наблюдения

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

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

Слева — данные без выброса, линия регрессии устойчива. Справа — та же выборка плюс одна точка (синий кружок) с экстремальным значением X и большим остатком. Линия регрессии сильно «просела»: угловой коэффициент изменился, хотя все остальные точки остались на месте.

3. Распределение зависимой переменной

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

Слева — линейная регрессия на счётных данных. Справа — остатки этой модели: виден нелинейный паттерн и «раздувание» дисперсии. Линейная модель не подходит: она предполагает непрерывный отклик с постоянной дисперсией, а здесь отклик счётный, дисперсия растёт с x, и предсказания могут быть отрицательными.

4. Характер связи между зависимой переменной и предикторами

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

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

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

Слева — нелинейная связь (парабола) и линейная регрессия. В середине она занижает, по краям — завышает. Справа — остатки этой модели: виден характерный нелинейный паттерн («горб»). Это сигнал, что линейная модель не подходит.

5. Независимость предикторов

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

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

6. Моделирование

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

7. Проверка валидности модели

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

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

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

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

Итог предварительного знакомства с данными

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

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

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

Ищем отскоки

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

Ищем отскоки: точечные диаграммы Кливленда

ggplot(bird, aes(y = 1:nrow(bird), x = ABUND)) +
  geom_point() +
  labs(y = 'Порядковый номер \nв датасете',
       x = 'Значения переменной')

Ищем отскоки: точечные диаграммы Кливленда

ggplot(bird, aes(y = 1:nrow(bird), x = AREA)) +
  geom_point() +
  labs(y = 'Порядковый номер \nв датасете',
       x = 'Значения переменной')

ggplot(bird, aes(y = 1:nrow(bird), x = logAREA)) +
  geom_point() +
  labs(y = 'Порядковый номер \nв датасете',
       x = 'Значения переменной')

Ищем отскоки: диаграммы Кливленда для всех переменных

Ищем отскоки: код для графика

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)

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 + \varepsilon_i\\ \end{aligned}\]

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

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

Вспомним условия применимости линейной регрессии

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

Мультиколлинеарность

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

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

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

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

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

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

Как рассчитывается VIF

Пусть наша модель \(y_i = b_0 + b_1 x_{1i} + b_2 x_{2i} + \ldots + b_{p - 1} x_{p - 1\;i} + \varepsilon_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 + \varepsilon_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} \]

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

Для модели ABUND ~ logAREA + YRISOL с \(n\) наблюдениями:

\[ \begin{pmatrix} ABUND_1 \\ ABUND_2 \\ \vdots \\ ABUND_n \end{pmatrix}_{n \times 1} = \begin{pmatrix} 1 & logAREA_1 & YRISOL_1 \\ 1 & logAREA_2 & YRISOL_2 \\ \vdots & \vdots & \vdots \\ 1 & logAREA_n & YRISOL_n \end{pmatrix}_{n \times 3} \begin{pmatrix} b_0 \\ b_1 \\ b_2 \end{pmatrix}_{3 \times 1} + \begin{pmatrix} e_1 \\ e_2 \\ \vdots \\ e_n \end{pmatrix}_{n \times 1} \]

• Строки — наблюдения (лесные участки)
• Столбцы — параметры: единицы для \(b_0\), затем предикторы
• Первый столбец единиц — чтобы \(b_0\) вошёл в умножение

Модельная матрица в множественной регрессии

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

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

# Обратная трансформация предиктора, который будем изображать
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.}\)).
• При обсуждении значимости отдельных предикторов можно привести таблицу с оценками коэффициентов и тестами их значимости.
• При сравнении влияния отдельных предикторов приводим стандартизованные коэффициенты
• Приводим график предсказаний модели

Take-home messages

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

Что почитать

• Гланц, С., 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