Простая и множественная регрессия

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

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

Подгружаем загруженные пакеты для их использования в текущей сессии работы 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

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

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

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.

Смешанные регрессионные модели

Используем набор данных Exam из пакета mlmRev

library(mlmRev)
str(Exam)
## 'data.frame':    4059 obs. of  10 variables:
##  $ school  : Factor w/ 65 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ normexam: num  0.261 0.134 -1.724 0.968 0.544 ...
##  $ schgend : Factor w/ 3 levels "mixed","boys",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ schavg  : num  0.166 0.166 0.166 0.166 0.166 ...
##  $ vr      : Factor w/ 3 levels "bottom 25%","mid 50%",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ intake  : Factor w/ 3 levels "bottom 25%","mid 50%",..: 1 2 3 2 2 1 3 2 2 3 ...
##  $ standLRT: num  0.619 0.206 -1.365 0.206 0.371 ...
##  $ sex     : Factor w/ 2 levels "F","M": 1 1 2 1 1 2 2 2 1 2 ...
##  $ type    : Factor w/ 2 levels "Mxd","Sngl": 1 1 1 1 1 1 1 1 1 1 ...
##  $ student : Factor w/ 650 levels "1","2","3","4",..: 143 145 142 141 138 155 158 115 117 113 ...

Посмотрим на связь переменных standLRT и normexam.

ggplot(Exam, aes(standLRT, normexam)) + theme_bw() + geom_point()

Применим простую линейную регрессию (МНК)

model_OLS <- lm(normexam ~ standLRT, data=Exam)
summary(model_OLS)
## 
## Call:
## lm(formula = normexam ~ standLRT, data = Exam)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.65615 -0.51848  0.01265  0.54399  2.97399 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.001191   0.012642  -0.094    0.925    
## standLRT     0.595057   0.012730  46.744   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8054 on 4057 degrees of freedom
## Multiple R-squared:   0.35,  Adjusted R-squared:  0.3499 
## F-statistic:  2185 on 1 and 4057 DF,  p-value: < 2.2e-16

Изобразим регрессионную прямую на основе полученной модели.

Exam$normexam_OLS <- predict(model_OLS)

ggplot(Exam, aes(standLRT, normexam)) + theme_bw() + geom_point() + 
    geom_line(aes(standLRT, normexam_OLS), col='red', size=1)

Подгружаем загруженные пакет lme4.

library(lme4)

Общий синтаксис функции lmer для спецификации смешанных моделей:

lmer(Dep.Var ~ Fixed.Effect + (1 + Fixed.Effect | Random.Effect), data=myDataSet)

Модель со случайным интерсептом

model_mixed_1 <- lmer(normexam ~ standLRT + (1|school), data=Exam)
summary(model_mixed_1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: normexam ~ standLRT + (1 | school)
##    Data: Exam
## 
## REML criterion at convergence: 9368.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7166 -0.6302  0.0294  0.6849  3.2673 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  school   (Intercept) 0.09384  0.3063  
##  Residual             0.56587  0.7522  
## Number of obs: 4059, groups:  school, 65
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.002323   0.040354    0.06
## standLRT    0.563307   0.012468   45.18
## 
## Correlation of Fixed Effects:
##          (Intr)
## standLRT 0.008

p-value для фиксированных эффектов можно посчитать с помощью пакета lmerTest

library(lmerTest)
model_mixed_1 <- lmer(normexam ~ standLRT + (1|school), data=Exam)
summary(model_mixed_1)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: normexam ~ standLRT + (1 | school)
##    Data: Exam
## 
## REML criterion at convergence: 9368.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7166 -0.6302  0.0294  0.6849  3.2673 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  school   (Intercept) 0.09384  0.3063  
##  Residual             0.56587  0.7522  
## Number of obs: 4059, groups:  school, 65
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 2.323e-03  4.035e-02 6.100e+01   0.058    0.954    
## standLRT    5.633e-01  1.247e-02 4.050e+03  45.180   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## standLRT 0.008

Визуализируем полученную модель.

Exam$normexam_mixed_1 <- predict(model_mixed_1)

ggplot(Exam, aes(standLRT, normexam)) + theme_bw() + geom_point() + 
    geom_line(aes(standLRT, normexam_mixed_1, col=school))

Модель со случайным интерсептом и случайным наклоном

model_mixed_2 <- lmer(normexam ~ standLRT + (1 + standLRT|school), data=Exam)
summary(model_mixed_2)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: normexam ~ standLRT + (1 + standLRT | school)
##    Data: Exam
## 
## REML criterion at convergence: 9327.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8323 -0.6317  0.0339  0.6834  3.4562 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  school   (Intercept) 0.09212  0.3035       
##           standLRT    0.01497  0.1223   0.49
##  Residual             0.55364  0.7441       
## Number of obs: 4059, groups:  school, 65
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) -0.01165    0.04011 60.65000   -0.29    0.772    
## standLRT     0.55653    0.02011 56.34000   27.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## standLRT 0.365

Визуализируем полученную модель.

Exam$normexam_mixed_2 <- predict(model_mixed_2)

ggplot(Exam, aes(standLRT, normexam)) + theme_bw() + geom_point() + 
    geom_line(aes(standLRT, normexam_mixed_2, col=school))

Модель со случайным наклоном (но без случайного интерсепта)

model_mixed_3 <- lmer(normexam ~ standLRT + (0 + standLRT|school), data=Exam)
summary(model_mixed_3)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: normexam ~ standLRT + (0 + standLRT | school)
##    Data: Exam
## 
## REML criterion at convergence: 9700.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2793 -0.6501  0.0241  0.6833  3.9119 
## 
## Random effects:
##  Groups   Name     Variance Std.Dev.
##  school   standLRT 0.0259   0.1609  
##  Residual          0.6251   0.7906  
## Number of obs: 4059, groups:  school, 65
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   -0.01463    0.01284 4019.00000  -1.139    0.255    
## standLRT       0.58782    0.02423   58.00000  24.256   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## standLRT 0.000

Визуализируем полученную модель.

Exam$normexam_mixed_3 <- predict(model_mixed_3)

ggplot(Exam, aes(standLRT, normexam)) + theme_bw() + geom_point() + 
    geom_line(aes(standLRT, normexam_mixed_3, col=school))

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

library(sjPlot)
sjt.lmer(model_mixed_1, model_mixed_2, model_mixed_3)
    normexam   normexam   normexam
    B CI p   B CI p   B CI p
Fixed Parts
(Intercept)   0.00 -0.08 – 0.08 .954   -0.01 -0.09 – 0.07 .772   -0.01 -0.04 – 0.01 .255
standLRT   0.56 0.54 – 0.59 <.001   0.56 0.52 – 0.60 <.001   0.59 0.54 – 0.64 <.001
Random Parts
σ2   0.566   0.554   0.625
τ00, school   0.094   0.092   0.026
ρ01       0.494    
Nschool   65   65   65
ICCschool   0.142   0.143   0.040
Observations   4059   4059   4059
R2 / Ω02   .441 / .441   .458 / .457   .381 / .380