Установка пакетов и рабочей директории

Установим необходимые для занятия пакеты (это действие необходимо, только если пакеты ещё не установлены на ваш компьютер) (для установки выполните следующую строку без символа # вначале).

# install.packages(c('foreign', 'car', 'ggplot2', 'gvlma'))

Подгружаем загруженные пакеты для их использования в текущей сессии работы R.

library(foreign)
library(car)
library(ggplot2)
library(gvlma)

Устанавливаем рабочую директорию: через меню или командой: Session -> Set working directory -> Choose directory

#setwd('...')  
#(вместо '...' укажите путь к папке, в которой хранится файл)

Открываем файл с данными child_data.sav. Он в формате sav (формат, используемый в SPSS), поэтому его можно открыть с помощью функции read.spss из пакета foreign.

df <- read.spss("child_data.sav", use.value.labels = T, to.data.frame = T, use.missings = T)

Простая линейная регрессия (метод наименьших квадратов, МНК)

Для правильной интерпретации коэффициентов МНК-модели нужно, чтобы ваши данные удовлетворяли ряду требований: • нормальность – значения зависимой переменной нормально распределены при фиксированных значениях независимых переменных; • независимость – значения зависимой переменной независимы друг от друга; • линейность – зависимая переменная линейно связана с независимыми; • гомоскедастичность – дисперсия зависимой переменной постоянна при разных значениях независимых переменных.

Попробуем разобраться на примере. В файле “child_data.sav” содержаться данные о школьниках (AGE - возраст, MEM_SPAN - баллы по тестированию памяти, IQ - коэффициент интеллекта, READ_AB - оценка навыков чтения).

str(df)
## 'data.frame':    20 obs. of  4 variables:
##  $ AGE     : num  6.7 5.9 5.5 6.2 6.4 7.3 5.7 6.15 7.5 6.9 ...
##  $ MEM_SPAN: num  4.4 4 4.1 4.8 5 5.5 3.6 5 5.4 5 ...
##  $ IQ      : num  95 90 105 98 106 100 88 95 96 104 ...
##  $ READ_AB : num  7.2 6 6 6.6 7 7.2 5.3 6.4 6.6 7.3 ...
##  - attr(*, "variable.labels")= Named chr  "age" "short-term memory span" "IQ" "reading ability"
##   ..- attr(*, "names")= chr  "AGE" "MEM_SPAN" "IQ" "READ_AB"

Посмотрим на гистограмму оценки качества чтения

hist(df$READ_AB, breaks = 10)

Видно, что по данной переменной есть разброс - одни дети читают лучше, другие - хуже. От чего это зависит? Можно предположить, что у детей с хорошей памятью лучше получается читать текст, понимать его суть и пересказывать (из этого складывается оценка по чтению). Проверим это предположение с помощью регрессионного анализа.

Построим диаграмму рассеяния для оценки качества чтения и оценки памяти.

plot <- ggplot(df, aes(READ_AB, MEM_SPAN)) + geom_point(size = 4) + theme_bw() + xlab('Чтение') + ylab('Память')
plot

Визуально кажется, что качество чтения связано с качеством памяти. Проведём через эти точки прямую линию.

plot + geom_smooth(method=lm)

Проверим наличие связи с помощью регрессионного анализа методом наименьших квадратов.

model <- lm(READ_AB ~ MEM_SPAN, data=df)
summary(model)
## 
## Call:
## lm(formula = READ_AB ~ MEM_SPAN, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68955 -0.22791 -0.01045  0.21278  1.02209 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.8803     0.7313   2.571   0.0192 *  
## MEM_SPAN      0.9767     0.1604   6.090 9.37e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4037 on 18 degrees of freedom
## Multiple R-squared:  0.6733, Adjusted R-squared:  0.6551 
## F-statistic: 37.09 on 1 and 18 DF,  p-value: 9.374e-06

p-value для F-статистики для модели равняется 9.374e-06, что меньше 0.05, значит модель не является бессмысленной и объясняет какую-то часть вариации READ_AB. Чтобы понять, какую часть вариации умения читать объясняет эта модель, нужно посмотреть на значение коэффициента детерминации Multiple R-squared, которое равно 0.6733. Т.е. наша модель объясняет 67% дисперсии умения читать.

Для переменной MEM_SPAN p-value меньше 0.05, следовательно мы должны отклонить нулевую гипотезу о равенстве нулю коэффициента при MEM_SPAN. Оценка этого коэффициента равна 0.9767. Запишем регрессионное уравнение полностью (уравнение для линии, которую мы провели через точки):

READ_AB = 1.8803 + 0.9767 * MEM_SPAN

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

Стандартную ошибку остатков (Residual standard error: 0.4037 балла) можно интерпретировать как усредненную ошибку предсказания умения читать по оценки памяти с использованием данной модели.

Посмотреть отдельно на коэффициенты

coef(model)
## (Intercept)    MEM_SPAN 
##   1.8803156   0.9767258

Поскольку коэффициент и константа являются оценкой, для них можно расчитать доверительный интервал.

confint(model)
##                 2.5 %   97.5 %
## (Intercept) 0.3439199 3.416711
## MEM_SPAN    0.6397881 1.313664

Посмотреть на остатки

residuals(model)
##           1           2           3           4           5           6 
##  1.02209073  0.21278107  0.11510848  0.03140039  0.23605523 -0.05230769 
##           7           8           9          10          11          12 
## -0.09652860 -0.36394477 -0.55463511  0.53605523 -0.68954635 -0.18256410 
##          13          14          15          16          17          18 
##  0.32441815  0.03605523 -0.38256410 -0.05230769  0.21278107 -0.18256410 
##          19          20 
## -0.38256410  0.21278107

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

fitted(model)
##        1        2        3        4        5        6        7        8 
## 6.177909 5.787219 5.884892 6.568600 6.763945 7.252308 5.396529 6.763945 
##        9       10       11       12       13       14       15       16 
## 7.154635 6.763945 5.689546 5.982564 6.275582 6.763945 5.982564 7.252308 
##       17       18       19       20 
## 5.787219 5.982564 5.982564 5.787219

Тестирование качества модели

Далее необходимо проверить, выполняются ли условия, лежищие в основе регрессионой МНК-модели. Наиболее распространенный подход – применить функцию plot() к объекту, представляющему собой результат действия функции lm(). Выражение par(mfrow=c(2,2)) используется для размещения четырех диаграмм, создаваемых функцией plot(), на одной большой диаграмме 2×2.

par(mfrow=c(2,2))
plot(model)

Второй простой способ - общая проверка выполнения требований, предъявляемых к линейным моделям.

gvmodel <- gvlma(model)
summary(gvmodel)
## 
## Call:
## lm(formula = READ_AB ~ MEM_SPAN, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68955 -0.22791 -0.01045  0.21278  1.02209 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.8803     0.7313   2.571   0.0192 *  
## MEM_SPAN      0.9767     0.1604   6.090 9.37e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4037 on 18 degrees of freedom
## Multiple R-squared:  0.6733, Adjusted R-squared:  0.6551 
## F-statistic: 37.09 on 1 and 18 DF,  p-value: 9.374e-06
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model) 
## 
##                     Value p-value                Decision
## Global Stat        5.8582  0.2100 Assumptions acceptable.
## Skewness           1.0395  0.3079 Assumptions acceptable.
## Kurtosis           0.3692  0.5434 Assumptions acceptable.
## Link Function      2.0686  0.1504 Assumptions acceptable.
## Heteroscedasticity 2.3808  0.1228 Assumptions acceptable.