Модель линейной регрессии

Описание

Как и корреляция, с которой мы познакомились ранее, регрессия выявляет линейную взаимосвязь между двумя переменными. Но в отличие от корреляции, которая отражает ненаправленную взаимосвязь, регрессия позводяет оценить, как изменяется одна переменная – которая называется зависимой – при изменении другой переменной – которую мы называем независимой. Важно, что выбор зависимой переменнйо (\(y\)) и независимой (\(x\)) осуществляется исследователем, например, на основании теории или предыдущих исследований. Сама по себе регрессия не позволяет выяснить, что является причиной, а что – следствием.

Модель линейной регрессии предполагает, что линейную зависимость между \(y\) и \(x\) можно наилучшим образом приблизить линейной функцией так, чтобы расстояние от каждого наблюдения до оцениваемой функции будет минимальным. Отклонения наблюдений от прямой принято называть отстатками (\(\epsilon\)).

Пускай истинная взаимосвязь между переменными выглядит следующим образом, где \(y_{i}\) и \(x_{i}\) – измеренные характеристики для наблюдения (респондента) \(i\):

\[y_{i}=\alpha + \beta \times x_i + \epsilon_{i}\]

Уравнение прямой, проходящей через облако рассеяния, выгялядит следующим образом, где \(\hat{y_{i}}\) – предсказанное значение зависимой переменной, если бы не существовало остатков:

\[\hat{y_{i}}=\alpha + \beta \times x_i \] Остатки в таком случае будут разностью между наблюдаемым значением зависимой переменной и предсказанным:

\[\epsilon_{i} = y_{i} - \hat{y_{i}}\] Иллюстрация:

Source: Wackerly (2008) Mathematical Statistics with Applications

Source: Wackerly (2008) Mathematical Statistics with Applications

Построение модели регрессии предполагает расчет характеристик \(\alpha\) (константа) и \(\beta\) (угол наклона), наилучшим образом приближающих зависимость между рассматриваемыми переменными, через минимизацию суммы квадратов остатков.

Интерпретация

\[y_{i}=\alpha + \beta \times x_i + \epsilon_{i}\]

Характеристики \(\alpha\) и \(\beta\) называются коэффициентами регрессии. Коэффициент \(\alpha\) – технический, он опускает или поднимает функцию относительно прямой OX и характеризует ожидаемое значение \(y\) при \(x=0\) (но он может быть никогда не равен 0, и тогда коэффициент не будет иметь содержательной интерпретации).

Нас интересует коэффициент \(\beta\) – коэффициент при независимой переменной. Он показывает, на сколько, в среднем, изменяется \(y\) при изменении \(x\) на единицу. Соответственно, знак при коэффициенте будет показывать, связаны ли \(y\) и \(x\) прямой или обратной связью.

Проверка гипотез

Как и любые характеристики, рассчитанные по выборке, наши оценки коэффициентов будут являться случайными величинами. В таких случаях нас интересует, будет ли истинное значение коэффициентов для генеральной совокупности отличаться от нуля при полученных оценках коэффициентов. Иными словами, нас интересует, является ли коэффициент статистически значимым. Если он значим, то мы можем сделать вывод, что влияние \(x\) на \(y\) существует в дейстительности.

Для того, чтобы сделать вывод о значимости или незначимости коэффициента, нам нужно воспользоваться \(t\)-критерием, проверяющим гипотезу о равенстве коэффициента нулю:

\[H_{0}: \beta = 0\]

Нулевая гипотеза утверждает, что коэффициент незначим, и независимая переменная не влияет на зависимую. \[H_{a}: \beta \ne 0\]

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

Аналогично мы можем проверить значимость коэффициента \(\alpha\).

Статистика критерия имеет распределение Стьюдента с \(df=n-2\) и расчитывается следующим образом, где \(se\) – стандартная ошибка коэффициента:

\[t_{\beta} = \frac{\hat{\beta}}{se(\beta)}\] \[t_{\alpha} = \frac{\hat{\alpha}}{se(\alpha)}\]

Качество модели

После оценки модели появляется закономерный вопрос: насколько хорошо она описывает наши данные? Ответить на него помогают (а) коэффициент детерминации \(R^{2}\) и (б) проверка условий Гаусса-Маркова.

Коэффициент детерминации \(R^{2}\) показывает, какую долю вариации зависимой переменной можно объяснить с помощью независимой переменной. Для парной регрессии \(R^{2}\) – это коэффициент корреляции Пирсона, возведенный в квадрат. Как это работает наглядно, можно увидеть здесь. \(R^{2}\) принимает значения от 0 до 1, и чем он выше, тем, считается, модель лучше описывает данные.

Коэффициент детерминации рассчитывается следующим образом: нужно из 1 вычесть отношение суммы квадратов отстатков (Residuals Sum of Squares) к сумме квадратов отклонений значений зависимой переменной от ее среднего (Total Sum of Squares):

\[R^{2}=1-\frac{RSS}{TSS}=1-\frac{\sum (y_{i}-\hat{y}_{i})^2}{\sum(y_{i}-\bar{y})^2}\]

Условия Гаусса-Маркова. Чтобы оценки коэффициентов, которые мы расчитываем, хорошо отражали реальность, нужно, чтобы выполнялось несколько условий.

  1. Существует линейная взаимосвязь между зависимой и независимой переменной.
  2. Математическое ожидание остатков равно нулю: \(E(\epsilon)=0\).
  3. Дисперсия остатков постоянна: \(D(\epsilon)=\sigma^2\).
  4. Случайность остатков.
  5. Нормальное распределение остатков: \(\epsilon \sim N(0,\sigma^2)\).

Если условия не выполняются, то оценки коэффициентов, которые мы получили, не являются самыми лучшими из возможных. Если же условия выполняются, то оценки коэффициентов называются лучшими линейными несмещенными оценками (Best Linear Unbiased Estimators, BLUE). Иными словами, оценкам силы и направления связи, которые мы получили по выборке, а также их статистической проверке на значимость можно доверять.

Оценка линейной регрессии

База данных и описательные статистики

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

df <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Ginzberg.csv")

Мы рассмотрим, как упрощенное видение мира (склонность интерпретировать события в черно-белых тонах) влияет на депрессивность респондента. Мы будем работать с переменными adjsimp – упрощенное видение мира – и adjdep – депрессивность. Эти переменные “очищены” от влияния других факторов, которые потенциально влияют на депрессивность, что позволит посмотреть на чистую взаимосвзяь между переменными.

Посмотрим на описательные статистики для переменных, с которыми мы будем работать, с помощью функции describe из пакета psych:

library(psych)
describe(df$adjsimp)
##    vars  n mean  sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 82    1 0.5   0.97    0.94 0.39 0.24 2.95  2.71 1.22     1.98 0.06
describe(df$adjdep)
##    vars  n mean  sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 82    1 0.5   0.97    0.96 0.49 0.23 2.79  2.55 0.84     0.74 0.06

Описательные статистики позволяют посмотреть на разброс значений и дисперсии переменных. Например, мы видим, что наша независимая переменная adjsimp не принимает значения 0, а значит, константа в модели не будет иметь содержательной интерпретации.

Диаграмма рассеяния и корреляция

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

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

plot(df$adjsimp, df$adjdep, 
     xlab = "Simplicity", ylab = "Depression", 
     col = "red", pch = 16)
abline(lm(df$adjdep~df$adjsimp))

Рассчитаем коэффициент корреляции Пирсона и проверим его на значимость.

cor.test(df$adjdep, df$adjsimp)
## 
##  Pearson's product-moment correlation
## 
## data:  df$adjdep and df$adjsimp
## t = 6.2485, df = 80, p-value = 1.889e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4061632 0.7024061
## sample estimates:
##      cor 
## 0.572693

Коэффициент корреляции равен \(R_{xy}=0.5727\) и статистически отличен от нуля. Если мы возведем это значение в квадрат, то узнаем, что примерно 33% дисперсии депрессии объясняется упрощенным видением мира (и наоборот), так как \(R^2\sim0.33\).

Регрессия

Для оценки линейной регрессии используется команда lm (от linear model – линейная модель). В качестве аргумента прописывается база данных, с которой вы работаете, и вид взаимосвязи: сначала зависимая переменная, тильда (~), независимая переменная.

model <- lm(data = df, adjdep ~ adjsimp)
summary(model)
## 
## Call:
## lm(formula = adjdep ~ adjsimp, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63156 -0.27468 -0.08266  0.19606  1.53530 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.42731    0.10235   4.175 7.53e-05 ***
## adjsimp      0.57269    0.09165   6.248 1.89e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4124 on 80 degrees of freedom
## Multiple R-squared:  0.328,  Adjusted R-squared:  0.3196 
## F-statistic: 39.04 on 1 and 80 DF,  p-value: 1.889e-08

Выдача результатов содежрит три части. Сначала приводятся описательные статистики для остатков модели (Residuals), затем значения коэффициентов и их проверка на значимость (Coefficients), затем проверка качества модели: \(R^2\) и проверка всей модели на значимость с помощью F-критерия (Multiple R-squared и F-statistics).

Рассмотрим выдачу коэффициентов. Она состоит из четырех колонок: оценка коэффициента (Estimate), стандартная ошибка (Std. Error), статистика критерия (t value) и p-value (Pr(>|t|)). Справа от таблцицы звездочками отмечаются значимые коэффициенты. Мы видим, что оба коэффициенты значимы на высоком уровне значимости (\(p-value<0.001\)). Иными словами, гипотеза о равенстве коэффициента нулю отвергается.

Рассмотрим коэффициент при переменной упрощенного видения мира \(\beta = 0.57269\). Он показывает, что при увеличении склонности упрощенно видеть мир на 1, показатель депрессивности растет на \(0.57269\). Его стандартная ошибка равна \(0.09165\). Если мы разделим оценку коэффициента на его стандартную ошибку, мы получим значение статистики критерия \(t_{\beta} = \frac{0.57269}{0.09165}=6.248\).

Коэффициент детерминации \(R^2\) обозначается как (Multiple R-squared) и равен 0.328. Это значит, что примерно треть дисперсии депрессивности объясняется упрощенным видением мира. Значение не очень высоко, но модель имеет объяснительную силу. Об этом же говорит и маленький p-value для F-критерия: гипотеза о незначимости модели отвергается.

Проверка условий Гаусса-Маркова

Чтобы оценить, насколько модель подходит к нашим данным, нужно проверить условия Гаусса-Маркова.

Давайте сохраним в отдельные векторы значения остатков и значения зависимой переменной, предсказанные моделью.

res <- residuals(model)
fitted <- fitted(model)

Нормальность

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

Визуально с помощью графика нормальной вероятностной бумаги:

qqnorm(res)
qqline(res)

Построим гистограмму остатков:

hist(res, freq = F, xlab = "Residuals", col = "plum")
curve(dnorm(x, 
            mean = mean(res),
            sd = sd(res)), add = T)

Наконец, проверим с помощью статистических критериев:

ks.test(res, "pnorm")
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  res
## D = 0.27512, p-value = 5.615e-06
## alternative hypothesis: two-sided
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.92714, p-value = 0.0001742

Как визуальная, так и статистическая проверка говорит о том, что условие о нормальности остатков не выполняется.

Дисперсия остатков

Остальные 3 условия Гаусса-Маркова можно проверить, построив диаграмму рассеяния между предсказанным значением зависимой переменной и остатками. По оси OX мы отложим предсказанные значения депрессивности, а по оси OY – значения остатков.

Что мы хотим увидеть на графике:

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

Иными словами, мы хотели бы, чтобы наши остатки выглядели случайными.

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

Посмотрим на график:

plot(fitted, res, 
     xlab = "Predicted Depression", ylab = "Residuals", 
     col = "red", pch = 16)
abline(h=0)

Иногда бывает удобно подписать точки, чтобы их можно было однозначно идентифицировать внутри базы данных:

plot(fitted, res, 
     xlab = "Predicted Depression", ylab = "Residuals", 
     col = "red", pch = 16)
abline(h=0)
text(fitted, res, labels=df$X, cex= 0.7, pos=3)

График далек от того, что мы бы хотели увидеть. Кажется, что присутствует некоторая нелилейная обратная зависимость между предсказанным значением зависимой переменной и остатками, а значит, мы делаем вывод о невыполнении условий Гаусса-Маркова.

Пример неслучайных остатков

Давайте посмотрим на ситуацию, когда условия Гаусса-Маркова не будут выполняться.

Симиулируруем данные. Создадим два вектора из 300 нвблюдений, взятых из стандартного нормального распределения.

set.seed(123)
s1 <- rnorm(300, mean = 0, sd = 1)
s2 <- rnorm(300, mean = 0, sd = 1)

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

s3<-s1^2+s2

Посмотрим на диаграмму рассеивания так, чтобы s3 лежал по оси Y (он будет нашей зависимой переменной), а s1 по оси X (он будет независимой). Очевидно, что между ними существует нелинейная взаимосвязь, а регрессионная примая не ловит эту взаимосвязь.

plot(s1, s3, col = "red", pch = 16)
abline(lm(s3~s1))

Попробуем оценить линейную регрессию.

model2<-lm(s3~s1)
summary(model2)
## 
## Call:
## lm(formula = s3 ~ s1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5417 -0.9441 -0.0906  0.8439  8.1610 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.89053    0.08694  10.243  < 2e-16 ***
## s1           0.32772    0.09201   3.562 0.000429 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.505 on 298 degrees of freedom
## Multiple R-squared:  0.04083,    Adjusted R-squared:  0.03761 
## F-statistic: 12.69 on 1 and 298 DF,  p-value: 0.0004288

Оценка коэффициента при независимой переменной оказался высоко значим, но можно ли доверять этой оценке? Давайте посмотрим на диаграмму рассеивания остатков.

plot(fitted(model2), residuals(model2), col = "red", pch = 16)
abline(h=0, col = "blue")

В них остался тот же самый паттерн, который мы видели изначально. Мы можем с уверенностью сказать, что остатки не случайны, и условия Гаусса-Маркова не выполняются.