R Notebook
Математическое моделирование
Практика 4
Линейные регрессионные модели
В практических примерах ниже показано:
- как оценить параметры модели линейной регрессии
- как строить модель со взаимосдействием переменных
- как визуализировать и тестировать остатки модели линейной регрессии
Модели: множественная линейная регрессия, kNN.
Данные: статистика по зарплатам в г.Москва за 2012 г., российский мониторинг экономического положения и здоровья населения НИУ-ВШЭ http://www.hse.ru/rlms.
Цель: исследовать набор данных wages.ru с помощью линейной регрессионной модели. Задействовав все возможные регрессоры, сделать вывод о пригодности модели для прогноза. Сравнить с методом k ближайших соседей по MSE на тестовой выборке.
# загрузка пакетов
library('GGally') # графики совместного разброса переменных
library('lmtest') # тесты остатков регрессионных моделей
library('FNN') # алгоритм kNN
# константы
my.seed <- 12345
train.percent <- 0.85
# загрузка данных
fileURL <- 'https://sites.google.com/a/kiber-guu.ru/msep/mag-econ/salary_data.csv?attredirects=0&d=1'
# преобразуем категориальные переменные в факторы
wages.ru <- read.csv(fileURL, row.names = 1, sep = ';', as.is = T)
wages.ru$male <- as.factor(wages.ru$male)
wages.ru$educ <- as.factor(wages.ru$educ)
wages.ru$forlang <- as.factor(wages.ru$forlang)
# обучающая выборка
set.seed(my.seed)
inTrain <- sample(seq_along(wages.ru$salary),
nrow(wages.ru) * train.percent)
df.train <- wages.ru[inTrain, c(colnames(wages.ru)[-1], colnames(wages.ru)[1])]
df.test <- wages.ru[-inTrain, -1]Описание переменных
Набор данных wages содержит переменные:
salary– среднемесячная зарплата после вычета налогов за последние 12 месяцев (рублей);
male– пол: 1 -- мужчина, 0 -- женщина;educ– уровень образования:
1 – 0-6 классов,
2 – незаконченное среднее (7-8 классов),
3 - незаконченное среднее плюс что-то ещё,
4 – законченное среднее,
5 – законченное среднее специальное, 6 – законченное высшее образование и выше;
forlang- иност. язык: 1 -- владеет, 0 – нет;exper– официальный стаж c 1.01.2002 (лет).
Размерность обучающей выборки: \(n = 127\) строк, \(p = 4\) объясняющих переменных. Зависимая переменная -- salary.
# описательные статистики по переменным
summary(df.train)## male educ forlang exper salary
## 0:63 2: 3 0:77 Min. : 1.000 Min. : 8000
## 1:64 3: 6 1:50 1st Qu.:10.000 1st Qu.: 20000
## 4:33 Median :10.000 Median : 30000
## 5:25 Mean : 8.858 Mean : 39322
## 6:60 3rd Qu.:10.000 3rd Qu.: 45000
## Max. :11.000 Max. :280000# совместный график разброса переменных
ggp <- ggpairs(df.train, upper = list(combo = 'box'))
print(ggp, progress = F)
# цвета по фактору male
ggp <- ggpairs(df.train[, c('exper', 'male', 'salary')],
aes(color = male), upper = list(combo = 'box'))
print(ggp, progress = F)
Судя по коробчатой диаграмме на пересечении sale и male, средний уровень зарплат у мужчин и женщин отличается: мужчины в среднем получают больше. Центральный график показывает, что доли наблюдений с различными значениями признака male в наборе данных примерно равны (мужчин чуть больше).
# цвета по фактору educ
ggp <- ggpairs(df.train[, c('exper', 'educ', 'salary')],
aes(color = educ), upper = list(combo = 'box'))
print(ggp, progress = F)
Коробчатые диаграммы на пересечении salary и educ показывают, что самому низкому уровню образования соответствует самая низкая средняя зарплата, и дисперсия зарплат здесь тоже наименьшая. У респондентов с самым высоким уровнем образования дисперсия зарплаты наибольшая. Судя по графику в центре, наблюдения распределены по значениям переменой educ неравномерно: группа с самым высоким уровнем образования самая многочисленная.
# цвета по фактору forlang
ggp <- ggpairs(df.train[, c('exper', 'forlang', 'salary')],
aes(color = forlang), upper = list(combo = 'box'))
print(ggp, progress = F)
Судя по графику, у сотрудников со знанием иностранного языка (forlang = 1) в среднем зарплата выше. Доли наблюдений с forlang = 1 и forlang = 0 сопоставимы.
Модели
model.1 <- lm(salary ~ . + exper:educ + exper:forlang + exper:male,
data = df.train)
summary(model.1)##
## Call:
## lm(formula = salary ~ . + exper:educ + exper:forlang + exper:male,
## data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -50045 -14603 -4430 9827 209955
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17659.57 31052.60 -0.569 0.5707
## male1 9008.51 20297.11 0.444 0.6580
## educ3 24476.15 53505.44 0.457 0.6482
## educ4 41594.50 29809.06 1.395 0.1656
## educ5 49783.80 32187.93 1.547 0.1247
## educ6 32269.82 18758.06 1.720 0.0881 .
## forlang1 -5659.91 20440.87 -0.277 0.7824
## exper 1875.17 2592.31 0.723 0.4709
## educ3:exper -40.13 5373.94 -0.007 0.9941
## educ4:exper -2393.94 2460.23 -0.973 0.3326
## educ5:exper -3613.94 2835.59 -1.274 0.2051
## educ6:exper NA NA NA NA
## forlang1:exper 2110.16 2221.56 0.950 0.3442
## male1:exper 1223.27 2181.54 0.561 0.5761
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30800 on 114 degrees of freedom
## Multiple R-squared: 0.2251, Adjusted R-squared: 0.1435
## F-statistic: 2.76 on 12 and 114 DF, p-value: 0.002509Совместное влияние exper:educ исключаем, т.к. параметры незначимы и недостаточно наблюдений для оценки одного из них.
model.2 <- lm(salary ~ . + exper:forlang + exper:male,
data = df.train)
summary(model.2)##
## Call:
## lm(formula = salary ~ . + exper:forlang + exper:male, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49068 -14149 -2670 8831 210932
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2385.4 25855.2 0.092 0.927
## male1 6028.6 19295.1 0.312 0.755
## educ3 21998.5 21941.9 1.003 0.318
## educ4 18520.2 18794.5 0.985 0.326
## educ5 16399.1 19216.2 0.853 0.395
## educ6 30615.5 18624.8 1.644 0.103
## forlang1 -8973.3 19520.7 -0.460 0.647
## exper -123.5 1954.4 -0.063 0.950
## forlang1:exper 2509.3 2126.6 1.180 0.240
## male1:exper 1515.5 2091.7 0.725 0.470
##
## Residual standard error: 30670 on 117 degrees of freedom
## Multiple R-squared: 0.2115, Adjusted R-squared: 0.1509
## F-statistic: 3.488 on 9 and 117 DF, p-value: 0.0007519Взаимодействие male1:exper также исключаем.
model.3 <- lm(salary ~ . + exper:forlang,
data = df.train)
summary(model.3)##
## Call:
## lm(formula = salary ~ . + exper:forlang, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47861 -14025 -3421 7880 212139
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4419.7 24040.2 -0.184 0.854449
## male1 19411.5 5565.1 3.488 0.000685 ***
## educ3 20700.0 21824.5 0.948 0.344828
## educ4 18221.9 18752.2 0.972 0.333177
## educ5 15516.9 19138.9 0.811 0.419141
## educ6 29765.4 18550.3 1.605 0.111260
## forlang1 -6770.1 19243.6 -0.352 0.725607
## exper 734.2 1551.9 0.473 0.637045
## forlang1:exper 2253.2 2092.8 1.077 0.283820
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30610 on 118 degrees of freedom
## Multiple R-squared: 0.208, Adjusted R-squared: 0.1543
## F-statistic: 3.874 on 8 and 118 DF, p-value: 0.000438Коэффициент при forlang1 наименее значимый, и его знак не соответствует здравому смыслу.
model.4 <- lm(salary ~ male + educ + exper,
data = df.train)
summary(model.4)##
## Call:
## lm(formula = salary ~ male + educ + exper, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44025 -15112 -4058 10440 220975
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10433 21403 -0.487 0.62682
## male1 18119 5626 3.221 0.00165 **
## educ3 23749 22092 1.075 0.28453
## educ4 20345 19024 1.069 0.28702
## educ5 19843 19351 1.025 0.30723
## educ6 36691 18607 1.972 0.05092 .
## exper 1465 1038 1.411 0.16081
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31090 on 120 degrees of freedom
## Multiple R-squared: 0.1691, Adjusted R-squared: 0.1276
## F-statistic: 4.072 on 6 and 120 DF, p-value: 0.0009425В модели практичски нет значимых объясняющих переменных. Вероятно, это из-за того, что подвыборки по уровням фактора educ очень маленькие. Попробуем сделать educ дискретной количественной переменной.
df.train$educ <- as.numeric(df.train$educ)
df.test$educ <- as.numeric(df.test$educ)
model.6 <- lm(salary ~ .,
data = df.train)
summary(model.6)##
## Call:
## lm(formula = salary ~ ., data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44990 -15724 -3497 7773 215010
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17011 13955 -1.219 0.225219
## male1 19797 5456 3.628 0.000418 ***
## educ 5659 2690 2.103 0.037477 *
## forlang1 13379 5968 2.242 0.026794 *
## exper 2053 1036 1.982 0.049773 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30530 on 122 degrees of freedom
## Multiple R-squared: 0.1852, Adjusted R-squared: 0.1585
## F-statistic: 6.932 on 4 and 122 DF, p-value: 4.613e-05Эта модель лучше, но по характеристикам качества очень слабая. Пробуем добавить взаимодействие exper:male.
df.train$educ <- as.numeric(df.train$educ)
model.7 <- lm(salary ~ . + exper:male,
data = df.train)
summary(model.7)##
## Call:
## lm(formula = salary ~ . + exper:male, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45656 -15591 -3293 8149 214344
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13331 15875 -0.840 0.4027
## male1 10935 18838 0.580 0.5627
## educ 5767 2708 2.130 0.0352 *
## forlang1 13408 5987 2.239 0.0270 *
## exper 1579 1418 1.113 0.2678
## male1:exper 1002 2038 0.492 0.6239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30630 on 121 degrees of freedom
## Multiple R-squared: 0.1868, Adjusted R-squared: 0.1532
## F-statistic: 5.56 on 5 and 121 DF, p-value: 0.0001205Очевидно, стоит остановиться на модели без взаимодействий. Проверим её остатки.
Проверка остатков
# тест Бройша-Пагана
bptest(model.6)##
## studentized Breusch-Pagan test
##
## data: model.6
## BP = 8.6687, df = 4, p-value = 0.06994# статистика Дарбина-Уотсона
dwtest(model.6)##
## Durbin-Watson test
##
## data: model.6
## DW = 1.9505, p-value = 0.39
## alternative hypothesis: true autocorrelation is greater than 0# графики остатков
par(mar = c(4.5, 4.5, 2, 1))
par(mfrow = c(1, 3))
plot(model.7, 1)
plot(model.7, 4)
plot(model.7, 5)
par(mfrow = c(1, 1))Судя по графику слева, остатки не случайны, и их дисперсия непостоянна. В модели есть три влиятельных наблюдения: 104, 105, 114, -- которые, однако, не выходят за пределы доверительных границ на третьем графике. Графики остатков заставляют усомниться в том, что остатки удовлетворяют условиям Гаусса-Маркова.
Сравнение с kNN
# фактические значения y на тестовой выборке
y.fact <- wages.ru[-inTrain, ]$salary
y.model.lm <- predict(model.6, df.test)
MSE.lm <- sum((y.model.lm - y.fact)^2) / length(y.model.lm)
# kNN требует на вход только числовые переменные
df.train.num <- as.data.frame(apply(df.train, 2, as.numeric))
df.test.num <- as.data.frame(apply(df.test, 2, as.numeric))
for (i in 2:50){
model.knn <- knn.reg(train = df.train.num[, !(colnames(df.train.num) %in% 'salary')],
y = df.train.num[, 'salary'],
test = df.test.num, k = i)
y.model.knn <- model.knn$pred
if (i == 2){
MSE.knn <- sum((y.model.knn - y.fact)^2) / length(y.model.knn)
} else {
MSE.knn <- c(MSE.knn,
sum((y.model.knn - y.fact)^2) / length(y.model.knn))
}
}
# график
par(mar = c(4.5, 4.5, 1, 1))
# ошибки kNN
plot(2:50, MSE.knn, type = 'b', col = 'darkgreen',
xlab = 'значение k', ylab = 'MSE на тестовой выборке')
# ошибка регрессии
lines(2:50, rep(MSE.lm, 49), lwd = 2, col = grey(0.2), lty = 2)
legend('bottomright', lty = c(1, 2), pch = c(1, NA),
col = c('darkgreen', grey(0.2)),
legend = c('k ближайших соседа', 'регрессия (все факторы)'),
lwd = rep(2, 2))
Как можно видеть по графику, ошибка регрессии на тестовой выборке больше, чем ошибка метода k ближайших соседей с k от 2 до 30. Далее с увеличением количества соседей точность kNN падает. Ошибка регрессионной модели на тестовой выборке очень велика и составляет \(\frac {\sqrt{MSE_{TEST}}}{\bar{y}_{TEST}} = 68.2\%\) от среднего значения зависимой переменной. Для модели регрессии это может означать отсутствие важного объясняющего фактора. Отсутствие важной объясняющей переменной подтверждается тем, что у лучшей модели kNN также низкая точность: она ошибается на 41.5% от среднего значения объясняющей переменной.
Упражнение 4
На наборе данных из своего варианта построить модели линейной регрессии с указанными Y и X. Рассмотреть модели с категориальными предикторами, включая их взаимодействие с непрерывными объясняющими переменными. Сгенерировать отчёт по структуре отчёта из практики. Включить в отчёт выводы по каждому из разделов (описание данных, модели, сравнение с kNN). Ответить на вопрос, пригодна ли построенная модель регрессии для прогнозирования и почему.
Код отчёта в файле .Rmd разместить на github.com, отчёт в формате html – на rpubs.com. Ссылки на репозиторий и отчёт на rpubs выслать на почту.
Варианты
| Номер варианта | Данные | Зависимая переменная | Объясняющие переменные | |
| непрерывные | дискретные (факторы) | |||
| 1 | Boston {MASS} | crim | zn, rm | chas |
| 2 | Boston {MASS} | crim | tax, black, indus | chas |
| 3 | Boston {MASS} | crim | indus, age | chas |
| 4 | Auto {ISLR} | mpg | weight, acceleration, year | cylinders |
| 5 | Auto {ISLR} | mpg | weight, displacement, horsepower | cylinders |
| 6 | Auto {ISLR} | mpg | weight, displacement, acceleration | cylinders |
| 7 | Auto {ISLR} | mpg | weight, acceleration, year | drv |
| 8 | Auto {ISLR} | mpg | weight, displacement, horsepower | drv |
| 9 | Auto {ISLR} | mpg | weight, displacement, acceleration | drv |
| 10 | Carseats {ISLR} | Sales | Price, Advertising | US |
| 11 | Carseats {ISLR} | Sales | Price, Advertising | ShelveLoc |
| 12 | Carseats {ISLR} | Sales | Price, Population | US |
| 13 | Carseats {ISLR} | Sales | Price, Population | ShelveLoc |
| 14 | Carseats {ISLR} | Sales | Price, Population | Urban |
| 15 | Carseats {ISLR} | Sales | Price, Advertising | Urban |
| 16 | College {ISLR} | Accept | Books, Personal, PhD | Private |
| 17 | College {ISLR} | Accept | BGrad.Rate, PhD, Expend | Private |
| 18 | College {ISLR} | Accept | Top25perc, PhD, Room.Board | Private |
Как открыть данные
На примере первого варианта: Boston {MASS} означает, что это набор данных Boston из пакета MASS.
library(‘MASS’) # загружаем пакет
data(Boston) # открываем данные
?Boston # справка по даннымКак опубликовать отчёт на RPubs
Зарегистрироваться на rpubs.com. Логин и пароль понадобятся при публикации отчёта.
Связать отчёт (кнопка «Knit» над окном скрипта).
Опубликовать отчёт (кнопка «Publish» в правом верхнем углу окна с html-документом).
Рис.1. Генерация отчёта
Рис. 2. Публикация отчёта
Источники
- James G., Witten D., Hastie T. and Tibshirani R. An Introduction to Statistical Learning with Applications in R. URL: http://www-bcf.usc.edu/~gareth/ISL/ISLR%20First%20Printing.pdf