Описание линейных моделей, проверка их значимости и валидности

Линейные модели…

Вадим Хайтов, Анастасия Лянгузова, Наталья Арбузова, Валентина Домашкина, Юта Тамберг, Марина Варфоломеева

Кафедра Зоологии беспозвоночных, СПбГУ

Описание и проверка значимости линейных моделей

Вы сможете

  • Подобрать линейную модель зависимости переменной-отклика от одного предиктора
  • Протестировать значимость линейной модели в целом и значимость отдельных ее коэффициентов при помощи t и F критериев
  • Проверить условия применимости линейной регрессии при помощи графиков

Вспомним пример из прошлой лекции

Пример: IQ и размеры мозга

Зависит ли уровень интеллекта от размера головного мозга? (Willerman et al. 1991)

Scan_03_11

Scan_03_11 by bucaorg(Paul_Burnett) on Flickr

Было исследовано 20 девушек и 20 молодых людей

У каждого индивида измеряли:

  • вес
  • рост
  • размер головного мозга (количество пикселей на изображении ЯМР сканера)
  • Уровень интеллекта измеряли с помощью IQ тестов

Пример: Willerman, L., Schultz, R., Rutledge, J. N., and Bigler, E. (1991), “In Vivo Brain Size and Intelligence”, Intelligence, 15, p.223–228.
Данные: “The Data and Story Library”
Фото: Scan_03_11 by bucaorg (Paul Burnett) on Flickr

Вспомним, на чем мы остановились

brain_model <- lm(PIQ ~ MRINACount, data = brain)
summary(brain_model)

Call:
lm(formula = PIQ ~ MRINACount, data = brain)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.582 -17.953  -1.601  17.027  42.262 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.744e+00  4.239e+01   0.041   0.9674  
MRINACount  1.202e-04  4.651e-05   2.586   0.0137 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 20.99 on 38 degrees of freedom
Multiple R-squared:  0.1496,    Adjusted R-squared:  0.1272 
F-statistic: 6.686 on 1 and 38 DF,  p-value: 0.01367

Уравнение и график зависимости

\[PIQ_i = 1.744 + 0.0001202 \cdot MRINACount_i\]

Тестирование гипотез о линейных моделях

Способы проверки значимости модели и ее коэффициентов

Два равноправных способа

Значима ли модель целиком?

  • F критерий

Значима ли связь между предиктором и откликом?

  • t-критерий

  • F-критерий

Тестирование гипотез с помощью t-критерия

Тестирование гипотез с помощью t-критерия

Гипотезы

Зависимость есть, если угловой коэффициент \(\beta_k \ne 0\)

Нулевая гипотеза \(H_0: \beta_k = 0\)

Тестирование гипотез с помощью t-критерия

Гипотезы

Зависимость есть, если угловой коэффициент \(\beta_k \ne 0\)

Нулевая гипотеза \(H_0: \beta_k = 0\)

Тестовая статистика

\[t=\frac{b_k - \beta_k}{SE_{b_k}} = \frac{b_k - 0}{SE_{b_k}} = \frac{b_k}{SE_{b_k}}\]

При соблюдении условий применимости, t-статистика должна подчиняться t-распределению с параметром \(df\)

Число степеней свободы: \(df = n − p\),

где

\(n\) — объем выборки,
\(p\) — число параметров модели,
\(b_k\) — конкретный коэффициент регрессии.

Для простой линейной регрессии с одним предиктором \(df = n - 2\).

Зависит ли IQ от размера головного мозга?

\[PIQ_i = 1.744 + 0.0001202 \cdot MRINACount_i\]


Call:
lm(formula = PIQ ~ MRINACount, data = brain)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.582 -17.953  -1.601  17.027  42.262 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.744e+00  4.239e+01   0.041   0.9674  
MRINACount  1.202e-04  4.651e-05   2.586   0.0137 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 20.99 on 38 degrees of freedom
Multiple R-squared:  0.1496,    Adjusted R-squared:  0.1272 
F-statistic: 6.686 on 1 and 38 DF,  p-value: 0.01367

Результаты теста на IQ статистически значимо связаны
с размерами мозга на МРТ ( \(t_{ 0.05, 38}\) = 2.59, \(p\) = 0.014 ).

Тестирование гипотез при помощи F-критерия

Общая изменчивость

Общая изменчивость \(SS_{t}\) — это сумма квадратов отклонений наблюдаемых значений \(y_i\) от общего среднего \(\bar y\)

Структура общей изменчивости

\[SS_t = \color{blue}{SS_r} + \color{green}{SS_e}\] Общую изменчивость зависимой переменной можно разложить на две составляющие:

\(\color{blue}{SS_r}\) - отклонения предсказанных значений от общего среднего (Объясненная изменчивость)

\(\color{green}{SS_e}\) - отклонения наблюдаемых значений от предсказанных (Остаточная изменчивость)

Общая изменчивость

Объясненная изменчивость

Случайная (остаточная) изменчивость

Все изменчивости

От сумм квадратов к дисперсиям

\(MS_t\), полная дисперсия \(\color{blue}{MS_r}\), дисперсия,
объясненная регрессией		 \(\color{green}{MS_e}\), остаточная дисперсия
\(MS_{t} =\frac{SS_{t}}{df_{t}}\) \(\color{blue}{MS_{r}} =\frac{\color{blue}{SS_{r}}}{df_{r}}\) \(\color{green}{MS_{e}} =\frac{\color{green}{SS_{e}}}{df_{e}}\)
\(SS_{t}=\sum{(y_i - \bar{y})^2}\) \(\color{blue}{SS_{r}}=\sum{(\hat{y}-\bar{y})^2}\) \(\color{green}{SS_{e}}=\sum{(y_i - \hat{y})^2}\)
\(df_{t} = n-1\) \(df_{r} = 1\) \(df_{e} = n-2\)

С помощью \(MS_r\) и \(MS_e\) можно тестировать значимость коэффициентов

Если дисперсии остатков для всех значений x равны, то

\(E(\color{blue}{MS_r}) = \color{green}{\sigma^2} + \color{red}{\beta_1}^2\sum(x_i - \bar x)^2 =\)

\(= \color{green}{\sigma^2} + \sigma_{x}^2\)

\(E(\color{green}{MS_e}) = \color{green}{\sigma^2}\)


Если зависимости нет, то \(\color{red}{\beta_1} = 0\), и тогда

\(\color{blue}{MS_r} \approx \color{green}{MS_e}\)

Мы хотим протестировать гипотезы:

  • \(H_0: \color{red}{\beta_1} = 0\)
  • \(H_A: \color{red}{\beta_1} \ne 0\)

Если справедлива \(H_A: \color{red}{\beta_1} \ne 0\), то \(\color{blue}{MS_r} > \color{green}{MS_e}\)

\[F= \frac{\color{blue}{MS _{r}}}{\color{green}{MS_{e}}}\] . . .

Поздравляю, вы только что изобрели F-критерий!

Тестирование значимости коэффициентов регрессии при помощи F-критерия

  • \(H_0: \color{red}{\beta_1} = 0\)
  • \(H_A: \color{red}{\beta_1} \ne 0\)

\[F_{df_r, df_e}= \frac{\color{blue}{MS _{r}}}{\color{green}{MS_{e}}}\]

Для простой линейной регрессии
\(df_{r} = 1\) и \(df_{e} = n - 2\)

  • У F-распределения два параметра: \(df_{r}\) и \(df_{e}\)

  • Логика F-теста аналогична t-критерию: Определяем вероятность того, что при условии справедливости \(H_0\) при повторных выборках появится значение F больше или равное F, полученному в нашем эксперименте.

  • F-тест будет односторонним, т.к. соотношение дисперсий может быть только положительным.

Таблица результатов дисперсионного анализа

Источник изменчивости df SS MS F
Регрессия \(df _r = 1\) \(\color{blue}{SS _r} = \sum{(\hat y _i - \bar y)^2}\) \(\color{blue}{MS _r} = \frac{\color{blue}{SS _r}}{df _r}\) \(F _{df _r, df _e} = \frac{\color{blue}{MS _r}}{\color{green}{MS _e}}\)
Остаточная \(df _e = n - 2\) \(\color{green}{SS _e} = \sum{(y _i - \hat y _i)^2}\) \(\color{green}{MS _e} = \frac{\color{green}{SS _e}}{df _e}\)
Общая \(df _t = n - 1\) \(SS _t = \sum {(y _i - \bar y)^2}\)



Минимальное упоминание результатов в тексте должно содержать \(F _{df _r, df _e}\) и \(p\).

Посмотрим на summary()

summary(brain_model)

Call:
lm(formula = PIQ ~ MRINACount, data = brain)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.582 -17.953  -1.601  17.027  42.262 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.744e+00  4.239e+01   0.041   0.9674  
MRINACount  1.202e-04  4.651e-05   2.586   0.0137 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 20.99 on 38 degrees of freedom
Multiple R-squared:  0.1496,    Adjusted R-squared:  0.1272 
F-statistic: 6.686 on 1 and 38 DF,  p-value: 0.01367

С помощью F-критерия мы показали, что модель целиком значима.

Дисперсионный анализ модели anova()

Analysis of Variance Table

Response: PIQ
           Df  Sum Sq Mean Sq F value  Pr(>F)  
MRINACount  1  2946.6 2946.61  6.6863 0.01367 *
Residuals  38 16746.4  440.69                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

С помощью F-критерия мы показали, что влияние предиктора значимо!

Оценка качества подгонки модели

В чем различие между этими двумя моделями?

–

У этих моделей разный разброс остатков:

  • Модель слева объясняет практически всю изменчивость
  • Модель справа объясняет не очень много изменчивости

Коэффициент детерминации — мера качества подгонки модели

Коэффициент детерминации описывает какую долю дисперсии зависимой переменной объясняет модель

\[R^2 = \frac{\color{blue}{SS_{r}}}{SS_{t}}\]

  • \(0 < R^2 < 1\)
  • \(R^2 = r^2\) — для простой линейной регрессии коэффициент детерминации равен квадрату коэффициента Пирсоновской корреляции

Если в модели много предикторов, нужно внести поправку

Скорректированный коэффициет детерминации (adjusted R-squared)

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

\[ R^2_{adj} = 1- (1-R^2)\frac{n-1}{n-p}\]

\(p\) - количество параметров в модели

Вводится штраф за каждый новый параметр

Еще раз смотрим на результаты регрессионного анализа зависимости IQ от размеров мозга


Call:
lm(formula = PIQ ~ MRINACount, data = brain)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.582 -17.953  -1.601  17.027  42.262 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.744e+00  4.239e+01   0.041   0.9674  
MRINACount  1.202e-04  4.651e-05   2.586   0.0137 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 20.99 on 38 degrees of freedom
Multiple R-squared:  0.1496,    Adjusted R-squared:  0.1272 
F-statistic: 6.686 on 1 and 38 DF,  p-value: 0.01367

Модель зависимости IQ от размеров мозга на МРТ объясняет всего лишь 15 % общей изменчивости IQ в этих данных.

Ценность такой модели сомнительна.

Как записываются результаты регрессионного анализа в тексте статьи?

Мы показали, что связь между результатами теста на IQ и размером головного мозга на МРТ описывается моделью вида
IQ = 1.74 + 0.00012 MRINACount ( \(F_{1,38}\) = 6.686, p = 0.014, \(R^2\) = 0.15)

Пишем статью?

Задание 1 & 2

Выполните задания 1 и 2 в одном из этих файлов:

  • 07_task_assumptions_catsM.R
  • 07_task_assumptions_GAG.R

Зачем нужна диагностика линейных моделей

Если стохастическая модель адекватна, то что мы можем сказать про остатки?

–

  • Остатки — это просто шум. Обычно считается, что \(\varepsilon_i \in N(0, \sigma)\).

  • В остатках не должно быть никаких закономерностей (паттернов).

Зачем нужна диагностика модели? Разве тестов было недостаточно?

Пусть у нас есть некоторая модель, которая показывает зависимость \(V_1 = f(V_2, V_3, V_4, V_5)\)

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

Call:
lm(formula = V1 ~ V2 + V3 + V4 + V5, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0652 -0.9878  0.1419  0.7958  1.9560 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.001194   0.021170   0.056    0.955    
V2          0.986328   0.128630   7.668 2.56e-14 ***
V3          0.972295   0.127506   7.625 3.54e-14 ***
V4          0.861429   0.120431   7.153 1.14e-12 ***
V5          0.927474   0.084432  10.985  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.001 on 2293 degrees of freedom
Multiple R-squared:   0.05, Adjusted R-squared:  0.04834 
F-statistic: 30.17 on 4 and 2293 DF,  p-value: < 2.2e-16

Все значимо!

Пишем статью?

Остатки модели

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

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

Oh, really?

http://www4.stat.ncsu.edu/~stefanski/NSF_Supported/Hidden_Images/stat_res_plots.html

В остатках скрыта информация!

Анализ остатков

Анализ остатков линейных моделей

Этот анализ дает возможность:

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

  2. Проверить соответствует ли построенная модель условиям применимости линейных моделей.

Условия применимости линейных моделей

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

Какие бывают остатки?

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

\[\color{green}{e_i} = y_i - \color{blue}{\hat{y_i}}\]

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

\[s_i = \frac{\color{green}{e_i}}{\sqrt{MS_e(1-\color{purple}{h_{ii}})}}\]

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

Что такое проекционная матрица?

По определению, остатки \(\color{green}{\mathbf{e}} = \mathbf{y} - \color{blue}{\hat{\mathbf{y}}}\).

Тогда \(\color{green}{\mathbf{e}} = \mathbf{y} - \mathbf{X} \mathbf{\beta} = \mathbf{y} - \mathbf{X}[\mathbf{X}'\mathbf{X}]^{-1}\mathbf{X}'\mathbf{y}\).


Обозначим \(\color{purple}{\mathbf{H}} \equiv \mathbf{X}[\mathbf{X}'\mathbf{X}]^{-1}\mathbf{X}'\).

Матрица \(\color{purple}{\mathbf{H}}\) — называется “хэт”-матрица (hat-matrix) или проекционная матрица, т.к. она позволяет получить предсказанные значения из наблюдаемых.

\(\color{blue}{\hat{\mathbf{y}}} = \color{purple}{\mathbf{H}}\mathbf{y}\)


Тогда остатки можно получить как \(\color{green}{\mathbf{e}} = \mathbf{y} - \color{purple}{\mathbf{H}}\mathbf{y} = (\mathbf{I} - \color{purple}{\mathbf{H}})\mathbf{y}\).

Диагональные элементы проекционной матрицы \(\color{purple}{h_{ii}}\) — это мера воздействия точек на ход линии регрессии.

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

Влиятельные наблюдения (influential points)

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

Из кн. Quinn, Keough, 2002

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

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

Воздействие точек (leverage)

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

Из кн. Quinn, Keough, 2002

Weighing Machine by neys fadzil on Flickr

Точки, располагающиеся дальше от \(\bar{x}\), оказывают более сильное влияние на \(\color{blue}{\hat{y_i}}\)

  • \(\color{purple}{h_{ii}}\) варьирует в промежутке от \(1/n\) до 1
  • Если \(\color{purple}{h_{ii}} > 2(p/n)\), то надо внимательно посмотреть на данное значение (\(p\) — число параметров модели, \(n\) — объем выборки)

Расстояние Кука (Cook’s distance)

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

\[D_i = \frac{\sum{( \color{blue}{\hat{y_{j}}} - \color{red}{\hat{y}_{j(i)}})^2}}{p \; MS_e}\]

  • \(\color{blue}{\hat{y_j}}\) — значение предсказанное полной моделью
  • \(\color{red}{\hat{y}_{j(i)}}\) — значение, предсказанное моделью, построенной без учета \(i\)-го значения предиктора
  • \(p\) — количество параметров в модели
  • \(MS_{e}\) — среднеквадратичная ошибка модели ( \(\hat\sigma^2\) )

Расстояние Кука покажет оба аспекта влиятельности

Формулу расстояния Кука для данного \(i\)-того наблюдения можно расписать:

\[D_i = \frac{\color{green}{e_i}^2}{p \; MS_e} \cdot \frac { \color{purple}{h_{ii}} } {(1 - \color{purple}{h_{ii}}) ^ 2}\]

Видно, что \(D_i\) зависит сразу от

  • \(\color{green}{e_i}\) — это величина остатков
  • \(\color{purple}{h_{ii}}\) — это воздействие наблюдений

При каком значении \(D_i\) будем считать наблюдение выбросом?

Два варианта:

  • \(D_i > 1\) — “мягкий” порог (для простоты в этом курсе такой)
  • \(D_i > \frac{4}{(n - p)}\) — “жесткий” порог ( \(n\) — объем выборки, \(p\) — число параметров)

Что делать с наблюдениями-выбросами?

  • Удалить?

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

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

Некоторые виды трансформаций

Трансформация Формула
степень -2 \(1/x^2\)
степень -1 \(1/x\)
степень -0.5 \(1/\sqrt{x}\)
степень 0.5 \(\sqrt{x}\)
логарифмирование \(log(x)\)

Данные для анализа остатков

library(ggplot2)
brain_diag <- fortify(brain_model)
head(brain_diag, 2)
  PIQ MRINACount       .hat   .sigma      .cooksd   .fitted    .resid
1 124     816932 0.06637876 20.87837 0.0498380059  99.98294 24.017065
2 124    1001121 0.06686959 21.27213 0.0003039221 122.13236  1.867637
   .stdresid
1 1.18403934
2 0.09209859
  • .hat — “сила воздействия” данного наблюдения (leverage)
  • .cooksd — расстояние Кука
  • .fitted — предсказанные значения
  • .resid — остатки
  • .stdresid — стандартизованные остатки

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

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

Порядок значений на графике такой же как в исходных данных.

Code 
cook_cutoff <- 4 / (nrow(brain) - length(coef(brain_model)))
  • Есть одно влиятельное наблюдение, которое нужно проверить, но сила его влияния невелика (расстояние Кука \(< 1\), и только одно наблюдение больше \(4/(n-p) = 0.11\))

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

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

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

1. Линейность связи

Нелинейность связи видно на графиках остатков

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

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

Что делать, если связь нелинейна?

  • Добавить неучтенные переменные или взаимодействия;
  • Применить линеаризующее преобразование (Осторожно!);
  • Применить обобщенную линейную модель с другой функцией связи (GLM);
  • Построить аддитивную модель (GAM), если достаточно наблюдений по \(x\);
  • Построить нелинейную модель, если известна форма зависимости.

Примеры линеаризующих преобразований

Логарифмирование

Возведение в степень, и т.д.

Не стоит трансформировать отклик!

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

1.При логарифмировании отклика вы будете изучать поведение мат.ожидания логарифма
\(E(log(y_i)) = b_0 + b_1x_{1i} + ... + e_i\).

2.Трансформация отклика не только линеаризует зависимость, но и затронет величину остатков \(e\).

Вместо трансформации отклика лучше использовать обобщенную линейную модель с подходящей функцией связи, например:

\(log(E(y_i)) = b_0 + b_1x_{1i}+ ... + e_i\)

2. Независимость

Каждое значение \(y_i\) должно быть независимо от любого другого \(y_j\)

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

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

Диагностика нарушений независимости

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

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

Нарушение независимости: Неучтенная переменная

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

Нарушение независимости: Автокорреляция

В данном случае, наблюдения — это временной ряд
или ряд наблюдений в пространстве.

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

Чем опасна автокорреляция (и любая другая зависимость наблюдений друг от друга)

Чем опасна автокорреляция (и любая другая зависимость наблюдений друг от друга)

Построим 1000 моделей на случайных выборках по 10 последовательных наблюдений

Для простого “шума”, когда нет автокорреляций и нет связи между \(y\) и \(x\) (справедлива \(H_0\)) частота ошибок составляет

[1] 0.048

Для автокоррелированных данных, при отсутствии связи между \(y\) и \(x\) (справедлива \(H_0\)) частота ошибок составляет

[1] 0.233

Повышается вероятность ошибок I рода.

Проверка на автокорреляцию

Проверка на автокорреляцию нужна, если данные — это временной ряд, или если известны координаты или время сбора проб. Способы проверки временной автокорреляции (годятся, если наблюдения в ряду расположены через равные интервалы):

  • График автокорреляционной функции остатков (ACF-plot) покажет корреляции с разными лагами.
  • Критерий Дарбина-Уотсона (значимость автокорреляции 1-го порядка).

Для проверки пространственных автокорреляций: вариограммы (большая и сложная тема!)

Что делать, если у вас нарушено условие независимости значений?

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

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

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

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

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

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

3. Нормальное распределение

3. Нормальное распределение \(y\) (для каждого уровня значений \(x\))

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

Проверка нормальности распределения остатков

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

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

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

\[y_i = b_0 + b_1x_{1i} + \cdots + e_i\]

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

\[e_i \sim N(0, \sigma)\]

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

Проверка нормальности распределения остатков

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

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

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

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

По оси Х — квантили теоретического распределения, по оси Y — квантили остатков модели. Если наблюдаемое распределение соответствует теоретическому, то точки должны лечь вдоль прямой по диагонали графика.


Обычные остатки должны подчиняться нормальному распределению \(e \sim N(0, \sigma)\)

Cтьюдентизированные остатки — t-распределению

Code 
library(car)
qqPlot(brain_model, id = FALSE) # из пакета car

Что делать, если остатки распределены не нормально?

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

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

4. Постоянство дисперсии

4. Постоянство дисперсии (= гомогенность дисперсии, гомоскедастичность)

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

Проверка гомогенности дисперсий

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

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

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

Проверка гомогенности дисперсий

Чем опасна гетероскедастичность

Создадим две генеральные совокупности, в которых нет связи между \(y\) и \(x\) (справедлива \(H_0\))

Чем опасна гетероскедастичность

Возьмем 1000 выборок из этих совокупностей по 10 наблюдений и построим много линейных моделей

  • Нет гетероскедастичности и нет связи между \(y\) и \(x\) (справедлива \(H_0\)) частота ошибок составляет
[1] 0.039
  • Для совокупности с гетероскедастичностью, при отсутствии связи между \(y\) и \(x\) (справедлива \(H_0\)) частота ошибок составляет
[1] 0.073

Повышается вероятность ошибок I рода.

Проверка на гетероскедастичность

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

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

Что делать, если найдена гетероскедастичность?

Трансформация может помочь, но лучше изменить модель.

Возможные причины гетероскедастичности

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

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

Самые частые паттерны в остатках

Некоторые частые паттерны на графиках остатков

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

  • Рис. a — Условия применимости соблюдаются, модель хорошая

  • Рис. b — Клиновидный паттерн. Есть гетероскедастичность. Модель плохая

  • Рис. c — Остатки рассеяны равномерно, но нужны дополнительные предикторы

  • Рис. d — Нелинейный паттерн. Линейная модель использована некорректно

Зоопарк нарушений условий применимости регрессионного анализа

Загон №1 (Curva nonlinearis)

Загон №2 (Heteroscedastica horribilis)

Загон №3 (Mensurae atypica)

Загон №4 (Anima amissa)

Загон №5 (Perdidit commercium)

Загон №6 (Data longitudinalis)

Практикум по анализу остатков

Задание

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

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

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

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

Блок 1

Code 
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

Графики

Code 
mod1 <- lm(y1 ~ x1, data = dat1)
mod1_diag <- fortify(mod1)
library(cowplot)
plot_grid(
  ggplot(mod1_diag, aes(x = 1:nrow(mod1_diag), y = .cooksd)) + 
    geom_bar(stat = "identity"),
  ggplot(data = mod1_diag, aes(x = .fitted, y = .stdresid)) + 
    geom_point() + geom_hline(yintercept = 0),
  ggplot(data = mod1_diag, aes(x = x1, y = .stdresid)) + 
    geom_point() + geom_hline(yintercept = 0),
  ggplot(dat1, aes(sample = y1)) + stat_qq() + stat_qq_line(),
  nrow = 2)
  • Выбросов нет, зависимость нелинейна
  • Небольшие отклонения от нормального распределения

Решение, блок 1

Code 
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, id = FALSE)

Блок 2

Code 
  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")

Решение, блок 2

Code 
mod2 <- lm(y2 ~ x2, data = dat2)
mod2_diag <- fortify(mod2)
plot_grid(
  ggplot(mod2_diag, aes(x = 1:nrow(mod2_diag), y = .cooksd)) + 
    geom_bar(stat = "identity"),
  ggplot(data = mod2_diag, aes(x = .fitted, y = .stdresid)) + 
    geom_point() + geom_hline(yintercept = 0),
  ggplot(data = mod2_diag, aes(x = x2, y = .stdresid)) + 
    geom_point() + geom_hline(yintercept = 0),
  ggplot(dat2, aes(sample = y2)) + stat_qq() + stat_qq_line(),
  nrow = 2)
  • Выбросов нет
  • Гетерогенность дисперсий, остатки не подчиняются нормальному распределению

Решение, блок 2

Code 
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, id = FALSE)

Блок 3

Code 
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")

Решение, блок 3

Code 
mod3 <- lm(y3 ~ x3, data = dat3)
mod3_diag <- fortify(mod3)
plot_grid(
  ggplot(mod3_diag, aes(x = 1:nrow(mod3_diag), y = .cooksd)) + 
    geom_bar(stat = "identity"),
  ggplot(data = mod3_diag, aes(x = .fitted, y = .stdresid)) + 
    geom_point() + geom_hline(yintercept = 0),
  ggplot(data = mod3_diag, aes(x = x3, y = .stdresid)) + 
    geom_point() + geom_hline(yintercept = 0),
  ggplot(dat2, aes(sample = y2)) + stat_qq() + stat_qq_line(),
  nrow = 2)
  • 100-е наблюдение сильно влияет на ход регрессии
  • Зависимость может быть нелинейна

Решение, блок 3

Code 
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, id = FALSE)

Задание для самостоятельной работы

Выполните последнее задание в одном из этих файлов: - 07_task_assumptions_catsM.R - 07_task_assumptions_GAG.R

Take-home messages

  • Гипотезы о наличии зависимости можно тестировать при помощи t- или F-теста.

  • Качество подгонки модели можно оценить при помощи коэффициента детерминации \(R^2\).

  • У линейных моделей есть условия применимости, поэтому не спешите описывать результаты — сначала проверьте.

  • Если условия применимости нарушены, то результатам тестов для этой модели нельзя верить (получаются заниженные доверительные вероятности, возрастает вероятность ошибок I рода).

  • Анализ остатков дает разностороннюю информацию о валидности моделей.

Что почитать

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