• Зачем нужен анализ остатков?
• Влиятельные наблюдения
• Типы остатков
• Проверка на влиятельные наблюдения в R
• Проверка условий применимости линейных моделей
  1. Линейность связи
  2. Независимость
  3. Нормальное распределение
  4. Постоянство дисперсии
• Тренинг по анализу остатков

dat <- read.table('data/orly_owl_Lin_4p_5_flat.txt')
fit <- lm(V1 ~ V2 + V3 + V4 + V5 - 1, data = dat)
coef(summary(fit))

##    Estimate Std. Error t value Pr(>|t|)
## V2    0.986     0.1280    7.70 1.99e-14
## V3    0.971     0.1266    7.67 2.50e-14
## V4    0.861     0.1196    7.20 8.30e-13
## V5    0.927     0.0833   11.13 4.78e-28

Все достоверно? Пишем статью?

Постройте график зависимости остатков от предсказанных значений при помощи этого кода

library(car)
residualPlot(fit, pch = ".")

library(car)
residualPlot(fit, pch = ".")

http://www4.stat.ncsu.edu/~stefanski/NSF\_Supported/Hidden\_Images/stat\_res\_plots.html

1) Проверка на наличие влиятельных наблюдений

2) Проверка условий применимости линейных моделей

1. Линейная связь
2. Независимость
3. Нормальное распределение
4. Гомогенность дисерсий
5. Отсутствие коллинеарности предикторов (для можественной регрессии)

Как зависит размер сердца от размера тела у котов? (Fisher 1947; Venables, Ripley 1994)

97 котов (самцы) весом больше 2 кг. Про каждого кота известно:
- Sex — пол - Bwt — вес тела в кг - Hwt — вес сердца в г

## 
## Call:
## lm(formula = Hwt ~ Bwt, data = cats)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.773 -1.048 -0.298  0.984  4.865 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -1.184      0.998   -1.19     0.24    
## Bwt            4.313      0.340   12.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.56 on 95 degrees of freedom
## Multiple R-squared:  0.629,  Adjusted R-squared:  0.625 
## F-statistic:  161 on 1 and 95 DF,  p-value: <2e-16

\(Hwt = - 1.2 + 4.3Bwt\)

“Сырые” остатки

\[\varepsilon_i = y_i - \hat{y_i}\]

Пирсоновские остатки

\[p_i = \frac{\varepsilon_i}{\sqrt{Var(\hat{y_i})}}\] - легко сравнивать (cтандартизованы)

Стандартизованные (внутренне-стьюдентизированные) остатки

\[s_i = \frac{p_i}{\sqrt{1 - h_{ii}}} = \frac{\varepsilon_i}{\sqrt{Var(\hat{y_i})(1-h_{ii})}}\] - легко сравнивать (стандартизованы), учитывают силу влияния наблюдений

• \(\sqrt{Var(\hat{y_i})}\) — cтандартное отклонение предсказанных значений
• \(h_{ii}\) — “сила воздействия” отдельных наблюдений (leverage, рычаг проекционной матрицы)

Влиятельные наблюдения — это наблюдения, которые вносят слишком большой вклад в оценку парметров (коэффициентов) модели.

Из кн. Quinn, Keugh, 2002

Учет каких из этих точек повлияет на ход регрессии и почему?

• Точка 1 почти не повлияет, т.к. у нее маленький остаток, хоть и большой \(X\)
• Точка 2 почти не повлияет, т.к. ее \(X\) близок к среднему, хоть и большой остаток
• Точка 3 повлияет сильно, т.к. у нее не только большой остаток, но и большой \(X\)

показывает силу влияния значений \(x_i\) на ход линии регрессии, то есть на \(\hat{y_i}\)

Из кн. Quinn, Keough, 2002

Weighing Machine by neys fadzil on Flickr

Из кн. Quinn, Keough, 2002

• Точки, располагающиеся дальше от \(\bar{x}\), оказывают более сильное влияние на \(\hat{y_i}\)
  
• Эта величина, в норме, варьирует в промежутке от \(1/n\) до 1
  
• Если \(h_{ii} > 2(p/n)\), то надо внимательно посмотреть на данное значение (p — число параметров, n — объем выборки)

описывает, как повлияет на модель удаление данного наблюдения

\[D_i = \frac{\sum{(\hat{y_j}-\hat{y}_{j(i)})^2}}{p \cdot MSe} \large( \frac {h_{ii}} {1 - h_{ii}} \large)\]

• \(\hat{y_j}\) - значение предсказанное полной моделью
• \(\hat{y}_{j(i)}\) - значение, предказанное моделью, построенной без учета \(i\)-го значения предиктора
• \(p\) - количество параметров в модели
• \(MSe\) - среднеквадратичная ошибка модели (\(\hat\sigma^2\))
• \(h_{ii}\) — “сила воздействия” отдельных наблюдений (leverage)

• Зависит одновременно от величины остатков и “силы воздействия” наблюдений.

• Условное пороговое значение. Наблюдение является выбросом (outlier), если:
• \(D_i > 1\)
• \(D_i > 4/(n − p)\) (\(n\) - объем выборки)

• Удалить?

Осторожно! Только очевидные ошибки в наблюдениях можно удалять. Лучше найти причины.

• Трансформировать? Это не всегда поможет.

library(ggplot2)
cat_diag <- fortify(cat_model)
head(cat_diag, 2)

##   Hwt Bwt   .hat .sigma .cooksd .fitted .resid .stdresid
## 1 6.5   2 0.0489   1.56 0.00988    7.44 -0.941     -0.62
## 2 6.5   2 0.0489   1.56 0.00988    7.44 -0.941     -0.62

• .hat — “сила воздействия” данного наблюдения (leverage)
• .cooksd — расстояние Кука
• .fitted — предсказанные значения
• .resid — остатки
• .stdresid — стандартизованные остатки

Проверяем наличие влиятельных наблюдений в cat_model.

Значения на графике расстояния Кука приведены в том же порядке, что и в исходных данных.

# График расстояния Кука
ggplot(cat_diag, aes(x = 1:nrow(cat_diag), y = .cooksd)) + 
  geom_bar(stat = "identity")

• Есть одно влиятельное наблюдение, которое нужно проверить, но сила его влияния невелика (расстояние Кука < 1, и даже меньше 4/(n-k-1) = 0.04)

Большую часть того, что нужно знать про остатки вы увидите на этом графике. А сейчас давайте научимся читать такой график.

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

Проверка на линейность связи

• График зависимости \(Y\) от \(x\) (и от других переменных, не включенных в модель)
• График остатков от предсказанных значений

То, что мы сможем сделать в ходе этого курса:
- Добавить неучтенные переменные или взаимодействия - Применить линеаризующее преобразование (Осторожно!)

То, что еще может помочь:
- Применить обобщенную линейную модель с другой функцией связи (GLM) - Построить аддитивную модель (GAM), если достаточно наблюдений по \(x\) - Построить нелинейную модель, если известна форма зависимости

Осторожно! Вы рискуете изучить не то, что хотели

Это нужно контролировать на этапе планирования сбора материала

• Наиболее частые источники зависимостей:
  • псевдоповторности (повторно измеренные объекты)
  • неучтенные переменные
  • временные автокорреляции (если данные - временной ряд)
  • пространственные автокорреляции (если пробы взяты в разных местах)
  • и т.п.

Взаимозависимости можно заметить на графиках остатков

• остатки vs. предсказанные значения
• остатки vs. переменные в модели
• остатки vs. переменные не в модели

• Слева: Если в модели не учтена переменная \(X2\), внешне все нормально, но величина остатков зависит от \(X2\)
• Справа: Если \(X2\) учесть, то зависимость остатков от \(X2\) исчезает

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

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

Выбор зависит от обстоятельств. Вот несколько возможных вариантов.

• псевдоповторности
  • избавляемся от псевдоповторностей, вычислив среднее
  • подбираем модель со случайным фактором
• неучтенные переменные
  • включаем в модель (если возможно)
• автокорреляции
  • моделируем автокорреляцию
  • подбираем модель со случайным фактором (= random effects model, mixed model)

Графики зависимости остатков от предикторов в модели

# Полный код
ggplot(data = cat_diag, aes(x =  Bwt, y = .stdresid)) + 
  geom_point() + geom_hline(yintercept = 0)

# То же самое с использованием ранее созданного gg_resid
gg_resid + aes(x =  Bwt)

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

В данном случае их нет

Это условие невозможно проверить “влоб”, т.к. обычно каждому \(X\) сообветствует лишь небольшое число \(Y\)

Если \(Y\) это нормально распределенная случайная величина

\[Y_i \sim N(\mu_{y_i}, \sigma^2)\]

и мы моделируем ее как

\[Y_i = b_0 + b_1x_{1i} + \cdots + \varepsilon_i\]

то остатки от этой модели — тоже нормально распределенная случайная величина

\[\varepsilon_i \sim N(0, \sigma^2)\]

Т.е. выполнение этого условия можно оценить по поведению случайной части модели.

Есть формальные тесты, но:

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

Лучший способ проверки — квантильный график остатков.

Квантиль — значение, которое заданная случайная величина не превышает с фиксированной вероятностью

Если точки — это реализации случайной величины из \(N(0, \sigma^2)\), то они должны лечь вдоль прямой \(Y=X\). Если это стьюдентизированные остатки — то используются квантили t-распределения

library(car)
qqPlot(cat_model) # из пакета car

## [1] 93 97

Зависит от причины

• Нелинейная связь?
  • Построить аддитивную модель (если достаточно наблюдений по \(x\))
  • Построить нелинейную модель (если известна форма зависимости)
• Неучтенные переменные?
  • добавляем в модель
• Зависимая переменная распределена по-другому?
  • трансформируем данные (неудобно)
  • подбираем модель с другим распределением остатков (обобщенную линейную модель)

Это самое важное условие, поскольку многие тесты чувствительны к гетероскедастичности.

Есть формальные тесты (тест Бройша-Пагана, тест Кокрана), но:

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

Лучший способ проверки на гомогенность дисперсий — график остатков.

Этот график у нас уже есть

gg_resid

• Гетерогенность дисперсий не выражена.

Трансформация может помочь…

Даже если трансформация может помочь, лучше поискать причину гетерогенности дисперсий

• Неучтенные переменные
  • добавляем в модель
• Зависимая переменная распределена по-другому
  • трансформируем данные (неудобно)
  • подбираем модель с другим распределением остатков (обобщенную линейную модель)
• Моделируем гетерогенность дисперсии.

Из кн. Logan, 2010, стр. 174

• Рис. a — Условия применимости соблюдаются, модель хорошая
• Рис. b — Клиновидный паттерн. Есть гетероскедастичность. Модель плохая
• Рис. c — Остатки рассеяны равномерно, но нужны дополнительные предикторы
• Рис. d — Нелинейный паттерн. Линейная модель использована некорректно

Выполните три блока кода

Какие нарушения условий применимости линейных моделей здесь наблюдаются?

Вам понадобятся

1. График расстояния Кука
2. График остатков от предсказанных значений
3. Графики остатков от предикторов в модели и не в модели
4. Квантильный график остатков

set.seed(90829)
x1 <- seq(1, 100, 1)
y1 <-  diffinv(rnorm(99))  + rnorm(100, 0.2, 4)
dat1 = data.frame(x1, y1)
ggplot(dat1, aes(x = x1, y =  y1)) + geom_point()+ 
  geom_smooth(method="lm", alpha = 0.7)

Графики

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

## [1] 14 81

mod1 <- lm(y1 ~ x1, data = dat1)

# Данные для графиков остатков
mod1_diag <- fortify(mod1)

# 1) График расстояния Кука
ggplot(mod1_diag, aes(x = 1:nrow(mod1_diag), y = .cooksd)) + 
  geom_bar(stat = "identity")

# 2) График остатков от предсказанных значений
gg_resid <- ggplot(data = mod1_diag, aes(x = .fitted, y = .stdresid)) + 
  geom_point() + geom_hline(yintercept = 0)
gg_resid

# 3) Графики остатков от предикторов в модели и не в модели
gg_resid + aes(x = x1)

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

  set.seed(7657674)
  x2 <- runif(1000, 1, 100)
  b_0 <- 100;  b_1 <- -20
  h <- function(x) x^0.5
  eps <- rnorm(1000, 0, h(x2))
  y2 <- b_0 + b_1 * x2 + eps
  dat2 <- data.frame(x2, y2)
  ggplot(dat2, aes(x = x2, y = y2)) + geom_point() + geom_smooth(method = "lm")

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

## [1] 887 928

mod2 <- lm(y2 ~ x2, data = dat2)

# Данные для графиков остатков
mod2_diag <- fortify(mod2)

# 1) График расстояния Кука
ggplot(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) Графики остатков от предикторов в модели и не в модели
gg_resid + aes(x = x2)

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

set.seed(9283)
x3 <- rnorm(25, 50, 10)
b_0 <- 20; b_1 <- 20; eps <- rnorm(50, 0, 100)
y3 <- b_0 + b_1*x3 + eps
y3[100] <- 1000; x3[100] <- 95; y3[99] <- 1300; x3[99] <- 90; y3[98] <- 1500; x3[98] <- 80
dat3 <- data.frame(x3, y3)
ggplot(dat3, aes(x=x3, y=y3)) + geom_point() + geom_smooth(method="lm")

• 100-е наблюдение сильно влияет на ход регрессии
• Зависимость нелинейна

## [1]   2 100

mod3 <- lm(y3 ~ x3, data = dat3)

# Данные для графиков остатков
mod3_diag <- fortify(mod3)

# 1) График расстояния Кука
ggplot(mod3_diag, aes(x = 1:nrow(mod3_diag), y = .cooksd)) + 
  geom_bar(stat = "identity")

# 2) График остатков от предсказанных значений
gg_resid <- ggplot(data = mod3_diag, aes(x = .fitted, y = .stdresid)) + 
  geom_point() + geom_hline(yintercept = 0)
gg_resid

# 3) Графики остатков от предикторов в модели и не в модели
gg_resid + aes(x = x3)

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

• У линейных моделей есть условия применимости
• Если условия применимости нарушены, то такой моделью нельзя пользоваться, даже если анализ показывает достоверную зависимость
• Анализ остатков дает разностороннюю информацию о валидности моделей

• Кабаков Р.И. 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