Парная линейная регрессия: практикум

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

surv <- read.csv("http://math-info.hse.ru/f/2018-19/psych-ms/survey01.csv",
                 dec = ",")
surv <- na.omit(surv)

Примечание: в качестве десятичного разделителя используется запятая, поэтому добавлена опция dec=",".

Оставим в датафрейме только строки, где оцененная студентом длина отрезка больше 0:

surv <- surv[surv$length > 0, ]

Добавим столбец с абсолютным отклонением оцененной длины от правильного ответа:

surv$diff_len <- abs(22 - surv$length)

Построим линейную модель, которая будет описывать, каким образом отклонение от правильного ответа (diff_len) зависит от баллов по математике (math):

mod <- lm(data = surv, diff_len ~ math)
summary(mod)
## 
## Call:
## lm(formula = diff_len ~ math, data = surv)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.463 -3.866 -2.823  1.985 22.412 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  7.53033    5.48942   1.372    0.174
## math        -0.02134    0.07424  -0.287    0.775
## 
## Residual standard error: 6.008 on 81 degrees of freedom
## Multiple R-squared:  0.001019,   Adjusted R-squared:  -0.01131 
## F-statistic: 0.08262 on 1 and 81 DF,  p-value: 0.7745

Что у нас с моделью? Модель получилась очень неудачная! Коэффициент при math не является значимым, \(R^2\) очень маленький. Даже если мы бы просто усреднили значения `diff_len, это было бы и то лучше, чем наша модель. Давайте построим другую модель. Модель линейной взаимосвязи, которая описывает, каким образом балл по биологии зависит от балла по математике:

mod2 <- lm(data = surv, bio ~ math)
summary(mod2)
## 
## Call:
## lm(formula = bio ~ math, data = surv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.280  -5.735  -1.492   8.477  17.659 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  51.1897     9.1217   5.612 2.71e-07 ***
## math          0.3939     0.1234   3.193    0.002 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.984 on 81 degrees of freedom
## Multiple R-squared:  0.1118, Adjusted R-squared:  0.1009 
## F-statistic:  10.2 on 1 and 81 DF,  p-value: 0.002003

Проинтерпретируем полученные результаты. Во-первых, судя по p-value для math, на 5% уровне значимости мы можем отвергнуть нулевую гипотезу о равенстве истинного значения коэффициента при math нулю, следовательно, коэффициент значим. Во-вторых, раз коэффициент значим, имеет смысл записать уравнение и содержательно его проинтерпретировать. Уравнение:

\[ \text{bio} = 51.19 + 0.39 \times \text{math}. \] При увеличении балла по математике на 1, балл по биологии, в среднем, увеличивается на 0.39 балла. Средний ожидаемый балл по биологии равен 51.19.

Коэффициент детерминации \(R^2\) равен 0.11, модель объясняет 11% дисперсии зависимой переменной. Посмотрим на остатки (ошибки) модели.

surv$res <- mod2$residuals

Построим диаграмму остатки vs независимая переменная:

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.2
ggplot(data = surv, aes(x = math, y = res)) +
  geom_point() + geom_hline(yintercept = 0, colour = 'red')

Если мы посмотрим на разброс точек относительно линии \(y=0\), то можем заключить следующее:

  1. Точки разбросаны хаотично. В распределении точек нет явных закономерностей.

  2. Дисперсия остатков примерно постоянна при изменениях значений math (если будем двигаться слева направо по графику).

  3. Линия \(y=0\) лежит примерно посередине облака точек, примерно одинаковое число точек находится по разные стороны от этой прямой и при этом они примерно одинаково удалены от нее. Нет оснований считать, что математическое ожидание остатков сильно отличается от нуля.

Можно посчитать выборочное среднее остатков модели:

mean(surv$res)
## [1] 2.182486e-16

Тоже близко к нулю, но, конечно, это еще не гарантия (нужна проверка гипотезы о среднем, оставим пока).

Проверим остатки на нормальность с помощью критерия Шапиро-Уилка:

shapiro.test(surv$res)
## 
##  Shapiro-Wilk normality test
## 
## data:  surv$res
## W = 0.97544, p-value = 0.1126

Судя по p-value, на 5% уровне значимости у нас нет оснований отвергнуть нулевую гипотезу о нормальности распределения. Несмотря на то, что наша модель имеет не слишком высокую предсказательную силу (\(R^2 = 0.1118\), модель объясняет примерно 11% изменчивости зависимой переменной), она применима, удовлетворяет условиям Гаусса-Маркова.