Линейная регрессия

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

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

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

Мы рассмотрим

  • Базовые идеи корреляционного анализа
  • Проблему двух статистических подходов: “Тестирование гипотез vs. построение моделей”
  • Разнообразие статистических моделей
  • Основы регрессионного анализа

Вы сможете

  • Оценить взаимосвязь между измеренными величинами
  • Объяснить, что такое линейная модель
  • Формализовать запись модели в виде уравнения
  • Подобрать модель линейной регрессии
  • Протестировать гипотезы о наличии зависимости при помощи t-критерия или F-критерия
  • Оценить предсказательную силу модели

Знакомимся с даными

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

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

Scan_03_11 by bucaorg(Paul_Burnett) on Flickr

Было исследовано 40 человек: 20 женского и 20 мужского пола.

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

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

Пример взят из работы: Willerman, L., Schultz, R., Rutledge, J. N., and Bigler, E. (1991), “In Vivo Brain Size and Intelligence,” Intelligence, 15, 223-228.
Данные представлены в библиотеке “The Data and Story Library” http://lib.stat.cmu.edu/DASL/

Знакомство с данными

Посмотрим на датасет:

brain <- read.csv("data/IQ_brain.csv", header = TRUE)
head(brain)
  Gender FSIQ VIQ PIQ Weight Height MRINACount
1 Female  133 132 124    118   64.5     816932
2   Male  140 150 124     NA   72.5    1001121
3   Male  139 123 150    143   73.3    1038437
4   Male  133 129 128    172   68.8     965353
5 Female  137 132 134    147   65.0     951545
6 Female   99  90 110    146   69.0     928799

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

sum(!complete.cases(brain))
[1] 2

Где пропущенные значения?

Где именно?

sapply(brain, function(x) sum(is.na(x)))
    Gender       FSIQ        VIQ        PIQ     Weight     Height MRINACount 
         0          0          0          0          2          1          0

Что это за случаи?

brain[!complete.cases(brain), ]
   Gender FSIQ VIQ PIQ Weight Height MRINACount
2    Male  140 150 124     NA   72.5    1001121
21   Male   83  83  86     NA     NA     892420

Каков объём выборки?

nrow(brain) ## Это без учета пропущенных значений
[1] 40

Корреляционный анализ

Цель практически любого исследования — поиск взаимосвязи величин и создание базы для предсказания неизвестного на основе имеющихся данных.

Корреляционный анализ

Наличие связи между явлениями не означает, что между ними существует причинно-следственная связь.

Сила связи между явлениями может быть количественно измерена.

Основные типы линейной связи между величинами

Криволинейные связи между величинами

Коэффициент ковариации

Оценивает сонаправленность отклонений двух величин от своих средних значений

\[cov_{x,y} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{n - 1}\]

Коэффициент ковариации варьирует в интервале \(-\infty < cov_{x,y} < +\infty\)

Коэффициент корреляции

Это стандартизованное значение ковариации:

\[r_{x,y} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})} {\sqrt{\sum(x_i-\bar{x})^2}\sqrt{\sum(y_i-\bar{y})^2}} = \frac{cov_{x,y}} {\sigma_x \sigma_y}\]

Коэффициент корреляции варьирует в интервале: \(-1 \le r_{x,y} \le +1\)

Коэффициент корреляции - это еще и параметр двумерного нормального распределения

Коэффициенты корреляции и условия их применимости

Коэффициент Функция Особенности применения
Коэффициент Пирсона cor(x,y,method="pearson") Оценивает связь двух нормально распределенных величин. Выявляет только линейную составляющую взаимосвязи.
Ранговые коэффициенты (коэффициент Спирмена, Кендалла) cor(x,y,method="spearman")
cor(x,y,method="kendall")		 Не зависят от формы распределения. Могут оценивать связь для любых монотонных зависимостей.

Оценка статистической значимости коэффициентов корреляции

  • Коэффициент корреляции — это статистика, значение которой описывает степень взаимосвязи двух сопряженных переменных. Следовательно, применима логика статистического критерия.
  • Нулевая гипотеза \(H_0: r=0\).
  • Бывают двусторонние \(H_a: r\ne 0\) и односторонние критерии \(H_a: r>0\) или \(H_a: r<0\).
  • Ошибка коэффициента Пирсона: \(SE_r=\sqrt{\frac{1-r^2}{n-2}}\).
  • Стандартизованная величина \(t=\frac{r}{SE_r}\) подчиняется распределению Стьюдента с параметром \(df = n-2\).
  • Для ранговых коэффициентов существует проблема “совпадающих рангов” (tied ranks), что приводит к приблизительной оценке \(r\) и приблизительной оценке уровня значимости.
  • Значимость коэффициента корреляции можно оценить пермутационным методом.

Задание

  • Постройте точечную диаграмму, отражающую взаимосвязь между результатами IQ-теста (PIQ) и размером головного мозга (MRINACount)
  • Определите силу и направление связи между этими величинами
  • Оцените сатистическую значимость коэффициента корреляции Пирсона между этими двумя переменными
  • Придумайте способ, как оценить корреляцию между всеми парами исследованных признаков

Hint 1: Обратите внимание на то, что в датафрейме есть пропущенные значения. Изучите, как работают с NA функции, вычисляющие коэффициенты корреляции.

Решение

pl_brain <- ggplot(brain, 
               aes(x = MRINACount, y = PIQ)) + 
  geom_point() + 
  xlab("Brain size") + 
  ylab("IQ test") +
  theme(text = element_text(size = 20))
pl_brain

Решение

cor.test(brain$PIQ, brain$MRINACount, method = "pearson", alternative = "two.sided")

    Pearson's product-moment correlation

data:  brain$PIQ and brain$MRINACount
t = 2.5858, df = 38, p-value = 0.01367
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.08562561 0.62322928
sample estimates:
      cor 
0.3868173

Решение

cor(brain[, 2:6], use = "pairwise.complete.obs")
              FSIQ         VIQ          PIQ       Weight      Height
FSIQ    1.00000000  0.94663878  0.934125147 -0.051482850 -0.08600241
VIQ     0.94663878  1.00000000  0.778135114 -0.076088042 -0.07106813
PIQ     0.93412515  0.77813511  1.000000000  0.002512154 -0.07672318
Weight -0.05148285 -0.07608804  0.002512154  1.000000000  0.69961400
Height -0.08600241 -0.07106813 -0.076723176  0.699614004  1.00000000

Два подхода к исследованию: Тестирование гипотезы VS Построение модели

  • Проведя корреляционный анализ, мы лишь ответили на вопрос “Существует ли статистически значимая связь между величинами?”

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

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

Простейший пример:

  • Между путем, пройденным автомобилем, и временем, проведенным в движении, несомненно есть связь. Хватает ли нам этого знания?
  • Для расчета величины пути в зависимости от времени необходимо построить модель: \(S=Vt\), где \(S\) — зависимая величина, \(t\) — независимая переменная, \(V\) — параметр модели.
  • Зная параметр модели (скорость) и значение независимой переменной (время), мы можем рассчитать (cмоделировать) величину пройденного пути.

Какие бывают модели?

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

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

\[y = b_0 + b_1x\]

\[y = b_0 + b_1x_1 + b_2x_2\]

Нелинейные модели

\[y = b_0 + b_1^x\]

\[y = b_0^{b_1x_1+b_2x_2}\]

Простые и многокомпонентные (множественные) модели

  • Простая модель \[y = b_0 + b_1x\]

  • Множественная модель \[y = b_0 + b_1x_1 + b_2x_2 + b_3x_3 + ... + b_nx_n\]

Детерминистские и стохастические модели

Модель: \(у_i = 2 + 5x_i\)

Два параметра:

угловой коэффициент (slope) \(b_1=5\);

свободный член (intercept) \(b_0=2\).

Чему равен \(y\) при \(x=10\)?

Модель: \(у_i = 2 + 5x_i + \varepsilon_i\)

Появляется дополнительный член \(\varepsilon_i\). Он вводит в модель влияние неучтенных моделью факторов.

Обычно считают, что \(\epsilon \in N(0, \sigma^2)\).

Случайная и фиксированая часть модели

В стохастические модели выделяется две части:

Фиксированная часть: \(у_i = 2 + 5x_i\)
Случайная часть: \(\varepsilon_i\)

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

Модель для зависимости величины IQ от размера головного мозга

Какая из линий “лучше” описывает облако точек?

Найти оптимальную модель позволяет регрессионный анализ

“Essentially, all models are wrong, but some are useful”

(George E. P. Box)

Происхождение термина “регрессия”

Френсис Галтон (Francis Galton)

“the Stature of the adult offspring … [is] … more mediocre than the stature of their Parents” (цит. по Legendre & Legendre, 1998)

Рост регрессирует (возвращается) к популяционной средней.

Угловой коэффициент в зависимости роста потомков от роста родителей — коэффициент регрессии.

Подбор линии регрессии проводится с помощью двух методов

  • С помощью метода наименьших квадратов (Ordinary Least Squares) — используется для простых линейных моделей.
  • Через подбор функции максимального правдоподобия (Maximum Likelihood) — используется для подгонки сложных линейных и нелинейных моделей.

Метод наименьших квадратов

(из кн. Quinn, Keough, 2002, стр. 85)

Остатки (Residuals):
\[\varepsilon_i = y_i - \hat{y_i}\]

Линия регрессии (подобранная модель) — это та линия, у которой \(\sum{\varepsilon_i}^2\) минимальна.

Подбор модели методом наименьших квадратов с помощью функци lm()

fit <- lm(formula, data)

Модель записывается в виде двусторонней формулы:

Модель Формула
Простая линейная регрессия
\(\hat{y_i}=b_0 + b_1x_i\)		 Y ~ X
Y ~ 1 + X
Y ~ X + 1
Простая линейная регрессия
(без \(b_0\), “no intercept”)
\(\hat{y_i}=b_1x_i\)		 Y ~ -1 + X
Y ~ X - 1
Уменьшенная простая линейная регрессия
\(\hat{y_i}=b_0\)		 Y ~ 1
Y ~ 1 - X
Множественная линейная регрессия
\(\hat{y_i}=b_0 + b_1x_i +b_2x_2\)		 Y ~ X1 + X2

Подбор модели методом наименьших квадратов с помощью функци lm()

fit <- lm(formula, data)

Элементы формул для записи множественных моделей

Элемент формулы Значение
: Взаимодействие предикторов
Y ~ X1 + X2 + X1:X2
* Обозначает полную схему взаимодействий
Y ~ X1 * X2 * X3
аналогично
Y ~ X1 + X2 + X3+ X1:X2 + X1:X3 + X2:X3 + X1:X2:X3
. Y ~ .
В правой части формулы записываются все переменные из датафрейма, кроме Y

Подберем модель, наилучшим образом описывающую зависимость результатов IQ-теста от размера головного мозга

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

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

Coefficients:
(Intercept)   MRINACount  
  1.7437570    0.0001203

Как трактовать значения параметров регрессионной модели?

Как трактовать значения параметров регрессионной модели?

  • Угловой коэффициент (slope) показывает на сколько единиц изменяется предсказанное значение \(\hat{y}\) при изменении на одну единицу значения предиктора \((x)\).
  • Свободный член (intercept) — величина во многих случаях не имеющая “смысла”, просто поправочный коэффициент, без которого нельзя вычислить \(\hat{y}\).
    NB! В некоторых линейных моделях он имеет смысл, например, значения \(\hat{y}\) при \(x = 0\).
  • Остатки (residuals) характеризуют влияние неучтенных моделью факторов.

Вопросы:

  1. Чему равны угловой коэффициент и свободный член полученной модели brain_model?
  2. Какое значение IQ-теста предсказывает модель для человека с объемом мозга равным 900000?
  3. Чему равно значение остатка от модели для человека с порядковым номером 10?

Ответы

  1. Чему равны угловой коэффициент и свободный член полученной модели brain_model?

Угловой коэффициент:

coefficients(brain_model) [1]
(Intercept) 
   1.743757

Свободный член:

coefficients(brain_model) [2]
  MRINACount 
0.0001202538

Ответы

Какое значение IQ-теста предсказывает модель для человека с объемом мозга равным 900000?

as.numeric(coefficients(brain_model) [1] + coefficients(brain_model) [2] * 900000)
[1] 109.9722

Ответы

  1. Чему равно значение остатка от модели для человека с порядковым номером 10?
brain$PIQ[10] - fitted(brain_model)[10]
      10 
30.35782 
residuals(brain_model)[10]
      10 
30.35782

Углубляемся в анализ модели: функция 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

Что означают следующие величины?

Estimate
Std. Error
t value
Pr(>|t|)

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

Параметр Оценка Стандартная ошибка
\(\beta_1\)
Slope		 \(b _1 = \frac {\sum _{i=1}^{n} {[(x _i - \bar {x})(y _i - \bar {y})]}}{\sum _{i=1}^{n} {(x _i - \bar x)^2}}\)
или проще
\(b_0 = r\frac{sd_y}{sd_x}\)		 \(SE _{b _1} = \sqrt{\frac{MS _e}{\sum _{i=1}^{n} {(x _i - \bar {x})^2}}}\)
\(\beta_0\)
Intercept		 \(b_0 = \bar y - b_1 \bar{x}\) \(SE _{b _0} = \sqrt{MS _e [\frac{1}{n} + \frac{\bar x}{\sum _{i=1}^{n} {(x _i - \bar x)^2}}]}\)
\(\epsilon _i\) \(e_i = y_i - \hat {y_i}\) \(\approx \sqrt{MS_e}\)

Где \(MS_e\) — Mean Sum of Squares — \(\frac{SS_e}{df_e}\).

Для чего нужны стандартные ошибки?

  • Они нужны, поскольку мы оцениваем параметры по выборке;
  • Позволяют построить доверительные интервалы для параметров;
  • Их используют в статистических тестах.

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

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

pl_brain + geom_smooth(method="lm")

Доверительная зона регрессии. В ней с 95% вероятностью лежит регрессионная прямая, описывающая связь в генеральной совокупности.


Возникает из-за неопределенности оценок коэффициентов модели, вследствие выборочного характера оценок.

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

Симулированный пример: Генеральная совокупность, в которой связь между Y и X, описывается следующей зависимостью

\[ y_i = 10 + 10x_i + \varepsilon_i \\ \varepsilon \in N(0, 20) \]

pop_x <- rnorm(1000, 10, 3)
pop_y <- 10 + 10 * pop_x + rnorm(1000, 0, 20)
population <- data.frame(x = pop_x, y = pop_y)

pop_plot <- ggplot(population, aes(x = x, y = y)) + 
  geom_point(alpha = 0.3, color = "red") + 
  geom_abline(aes(intercept = 10, slope = 10), 
              color = "blue", size = 2) +
  theme(text = element_text(size = 15))
pop_plot

Зависимости, выявленные в нескольких разных выборках

Линии регрессии, полученные для 100 выборок (по 20 объектов в каждой), взятых из одной и той же генеральной совокупности.

Code 
samp_coef <- data.frame(b0 = rep(NA, 100), b1 =  rep(NA, 100))

for(i in 1:100) {
  samp_num <- sample(1:1000, 20)
  samp <- population[samp_num, ]
  fit <- lm(y ~ x, data = samp)
  samp_coef$b0[i] <- coef(fit)[1]
  samp_coef$b1[i] <- coef(fit)[2]
}

ggplot(population, aes(x = x, y = y)) + geom_point(alpha = 0.3, color = "red")+ geom_abline(aes(intercept = b0, slope = b1), data = samp_coef) + geom_abline(aes(intercept = 10, slope = 10), color = "blue", size = 2) +
  theme(text = element_text(size = 18))

Доверительные интервалы для коэффициентов уравнения регрессии

coef(brain_model)
 (Intercept)   MRINACount 
1.7437570401 0.0001202538 
confint(brain_model)
                    2.5 %       97.5 %
(Intercept) -8.407513e+01 8.756265e+01
MRINACount   2.610804e-05 2.143996e-04

Для разных \(\alpha\) можно построить разные доверительные интервалы

Важно!

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

Какое значение IQ можно ожидать у человека с размером головного мозга 900000?

newdata <- data.frame(MRINACount = 9e+05)

predict(brain_model, newdata, interval = "prediction", level = 0.95, se = TRUE)$fit
       fit      lwr      upr
1 109.9722 66.93883 153.0055
  • При размере мозга 900000 среднее значение IQ будет, с вероятностью 95%, находиться в интервале от 67 до 153.

Отражаем на графике область значений, в которую попадут 95% предсказанных величин IQ

Подготавливаем данные:

brain_predicted <- predict(brain_model, interval="prediction")
brain_predicted <- data.frame(brain, brain_predicted)
head(brain_predicted)
  Gender FSIQ VIQ PIQ Weight Height MRINACount       fit      lwr      upr
1 Female  133 132 124    118   64.5     816932  99.98294 56.09762 143.8683
2   Male  140 150 124     NA   72.5    1001121 122.13236 78.23695 166.0278
3   Male  139 123 150    143   73.3    1038437 126.61975 81.89560 171.3439
4   Male  133 129 128    172   68.8     965353 117.83112 74.47698 161.1853
5 Female  137 132 134    147   65.0     951545 116.17066 72.95702 159.3843
6 Female   99  90 110    146   69.0     928799 113.43537 70.36855 156.5022

Отражаем на графике область значений, в которую попадут 95% предсказанных величин IQ

Код для построения графика

pl_brain + 

# 1) Линия регрессии и ее дов. интервал
# Если мы указываем fill внутри aes() и задаем фиксированное значение - 
#  появится соотв. легенда с названием.
# alpha - задает прозрачность
  geom_smooth(method = "lm", aes(fill = "Conf.interval"), alpha = 0.4) +
# 2) Интервал предсказаний создаем при помощи геома ribbon ("лента")
# Данные берем из другого датафрейма - из brain_predicted
# ymin и ymax - эстетики геома ribbon, которые задают нижний и верхний 
# край ленты в точках с заданным x (x = MRINACount было задано в ggplot() 
# при создании pl_brain, поэтому сейчас его указывать не обязательно)
#

geom_ribbon(data = brain_predicted,  aes(ymin = lwr, ymax = upr, fill = "Conf. area for prediction"), alpha = 0.2) +
# 3) Вручную настраиваем цвета заливки при помощи шкалы fill_manual.
# Ее аргумент name - название соотв. легенды, values - вектор цветов
  scale_fill_manual(name = "Intervals", values = c("green", "gray")) +
# 4) Название графика
  ggtitle("Confidence interval \n and confidence area for prediction")

Важно!

Модель “работает” только в том диапазоне значений независимой переменной \((x)\), для которой она построена (интерполяция).


Экстраполяцию надо применять с большой осторожностью.

Итак, что означают следующие величины?

  • Estimate
  • Оценки праметров регрессионной модели
  • Std. Error
  • Стандартная ошибка для оценок
  • Осталось решить, что такое t value, Pr(>|t|)

Summary

  • Модель простой линейной регрессии \(y _i = \beta _0 + \beta _1 x _i + \epsilon _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