Упражнение_4
Нестерова А.И.
16 03 2020
Математическое моделирование
Линейные регрессионные модели
На наборе данных из своего варианта построить модели линейной регрессии с указанными Y и X. Рассмотреть модели с категориальными предикторами, включая их взаимодействие с непрерывными объясняющими переменными. Сгенерировать отчёт по структуре отчёта из практики. Включить в отчёт выводы по каждому из разделов (описание данных, модели, сравнение с kNN). Ответить на вопрос, пригодна ли построенная модель регрессии для прогнозирования и почему. Код отчёта в файле .Rmd разместить на github.com, отчёт в формате html – на rpubs.com. Ссылки на репозиторий и отчёт на rpubs выслать на почту.
Вариант - 13
# загрузка пакетов
library('GGally') # графики совместного разброса переменных
library('lmtest') # тесты остатков регрессионных моделей
library('FNN') # алгоритм kNN
library('ISLR')
# загрузка данных Carseats
data('Carseats')
# отбор необходимых данных для построения моделей
Carseats <- Carseats[,c('Sales', 'Price', 'Population', 'ShelveLoc'),drop=FALSE]
# константы
my.seed <- 12345
train.percent <- 0.85
# преобразуем категориальную переменную ShelveLoc в фактор
Carseats$ShelveLoc <- as.factor(Carseats$ShelveLoc)
# обучающая выборка
set.seed(my.seed)
inTrain <- sample(seq_along(Carseats$Price),
nrow(Carseats) * train.percent)
df.train <- Carseats[inTrain, c(colnames(Carseats)[-1], colnames(Carseats)[1])]
df.test <- Carseats[-inTrain,-1]Описание переменных
Набор данных Carseats содержит переменные:
Sales- объем продаж в каждом месте (в тысячах);Price– взимаемая плата за автокресла на каждом участке;
Population– численность населения в регионе (в тысячах);ShelveLoc– коэффициент, указывающий на качество места для размещения автомобильных сидений на каждом участке: 1 – плохое (Bad); 2 – хорошее (Good); 3 - среднее (Medium).
Размерность обучающей выборки: \(n = 340\) строк, \(p = 3\) объясняющих переменных. Зависимая переменная – Sales.
## Price Population ShelveLoc Sales
## Min. : 49.00 Min. : 10.0 Bad : 83 Min. : 0.000
## 1st Qu.: 99.75 1st Qu.:131.0 Good : 72 1st Qu.: 5.390
## Median :117.50 Median :266.0 Medium:185 Median : 7.420
## Mean :115.96 Mean :259.4 Mean : 7.453
## 3rd Qu.:131.00 3rd Qu.:384.0 3rd Qu.: 9.160
## Max. :191.00 Max. :509.0 Max. :16.270
Судя по коробчатой диаграмме на пересечении Sales и ShelveLoc, средний объем продаж отличается зависимости от места размещения автомобильного кресла: чем лучше расположено автомобильное кресло, тем больше объём продаж. Нижний правый график показывает, что доли наблюдений с различными значениями признака ShelveLocв наборе данных имеют следующий вид: наибольшую часть наблюдений отражает коэффициент со значением Medium, а c Bad и Good примерно равны (с коэффициентом Bad чуть больше).
Модели
##
## Call:
## lm(formula = Sales ~ . + Price:ShelveLoc + Population:ShelveLoc,
## data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7337 -1.3521 -0.0655 1.2537 4.9755
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.983881 1.156665 8.632 2.58e-16 ***
## Price -0.048566 0.008656 -5.611 4.26e-08 ***
## Population 0.004149 0.001422 2.918 0.00377 **
## ShelveLocGood 7.924355 1.632048 4.855 1.86e-06 ***
## ShelveLocMedium 3.440151 1.394985 2.466 0.01417 *
## Price:ShelveLocGood -0.016531 0.012235 -1.351 0.17757
## Price:ShelveLocMedium -0.005061 0.010739 -0.471 0.63774
## Population:ShelveLocGood -0.004119 0.002207 -1.867 0.06284 .
## Population:ShelveLocMedium -0.004072 0.001674 -2.432 0.01555 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.849 on 331 degrees of freedom
## Multiple R-squared: 0.5709, Adjusted R-squared: 0.5606
## F-statistic: 55.06 on 8 and 331 DF, p-value: < 2.2e-16Совместное влияние Price:ShelveLoc исключаем, т.к. параметры незначимы (наименее занчимы по сравнению с другими незначимыми коэффициентами).
##
## Call:
## lm(formula = Sales ~ . + Population:ShelveLoc, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7586 -1.2998 -0.0519 1.3353 5.0581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.821872 0.690813 15.665 < 2e-16 ***
## Price -0.055295 0.004407 -12.546 < 2e-16 ***
## Population 0.003896 0.001394 2.795 0.0055 **
## ShelveLocGood 5.906933 0.647657 9.120 < 2e-16 ***
## ShelveLocMedium 2.794785 0.497060 5.623 3.98e-08 ***
## Population:ShelveLocGood -0.003751 0.002183 -1.718 0.0867 .
## Population:ShelveLocMedium -0.003814 0.001652 -2.309 0.0215 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.849 on 333 degrees of freedom
## Multiple R-squared: 0.5684, Adjusted R-squared: 0.5606
## F-statistic: 73.09 on 6 and 333 DF, p-value: < 2.2e-16Взаимодействие Population:ShelveLoc также исключаем.
##
## Call:
## lm(formula = Sales ~ Price + Population + ShelveLoc, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5466 -1.3152 -0.0686 1.3418 4.9608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.7117532 0.5836817 20.065 < 2e-16 ***
## Price -0.0563697 0.0044058 -12.794 < 2e-16 ***
## Population 0.0010052 0.0006861 1.465 0.144
## ShelveLocGood 4.9172335 0.2996406 16.410 < 2e-16 ***
## ShelveLocMedium 1.7936158 0.2457789 7.298 2.14e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.859 on 335 degrees of freedom
## Multiple R-squared: 0.5611, Adjusted R-squared: 0.5558
## F-statistic: 107.1 on 4 and 335 DF, p-value: < 2.2e-16Параметр Population является незначимым, поэтому его тоже исключим из уравнения регресии.
##
## Call:
## lm(formula = Sales ~ Price + ShelveLoc, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4731 -1.3552 -0.0733 1.2653 5.1448
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.023705 0.544391 22.087 < 2e-16 ***
## Price -0.056765 0.004405 -12.887 < 2e-16 ***
## ShelveLocGood 4.915635 0.300150 16.377 < 2e-16 ***
## ShelveLocMedium 1.784519 0.246119 7.251 2.88e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.862 on 336 degrees of freedom
## Multiple R-squared: 0.5583, Adjusted R-squared: 0.5543
## F-statistic: 141.5 on 3 and 336 DF, p-value: < 2.2e-16В данной модели все коэффициенты оказались значимыми, но по характеристикам качества модель недостаточно хороша (\(R^2=0.5583\)). Пробуем добавить взаимодействие Price:ShelveLoc.
##
## Call:
## lm(formula = Sales ~ Price + ShelveLoc + Price:ShelveLoc, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4781 -1.3598 -0.0816 1.2427 4.9877
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.7470469 0.9943558 11.814 < 2e-16 ***
## Price -0.0543499 0.0084953 -6.398 5.33e-10 ***
## ShelveLocGood 6.1701700 1.4415001 4.280 2.44e-05 ***
## ShelveLocMedium 1.6956329 1.2493941 1.357 0.176
## Price:ShelveLocGood -0.0107554 0.0121608 -0.884 0.377
## Price:ShelveLocMedium 0.0007331 0.0106409 0.069 0.945
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.864 on 334 degrees of freedom
## Multiple R-squared: 0.5599, Adjusted R-squared: 0.5533
## F-statistic: 84.98 on 5 and 334 DF, p-value: < 2.2e-16Очевидно, стоит остановиться на модели без взаимодействий. Проверим её остатки.
Проверка остатков
# тест Бройша-Пагана
bptest(model.4)##
## studentized Breusch-Pagan test
##
## data: model.4
## BP = 0.35539, df = 3, p-value = 0.9493# статистика Дарбина-Уотсона
dwtest(model.4)##
## Durbin-Watson test
##
## data: model.4
## DW = 2.0934, p-value = 0.805
## alternative hypothesis: true autocorrelation is greater than 0# графики остатков
par(mar = c(4.5, 4.5, 2, 1))
par(mfrow = c(1, 3))
plot(model.4, 1)
plot(model.4, 4)
plot(model.4, 5)
par(mfrow = c(1, 1))Судя по графику слева, остатки случайны, и их дисперсия постоянна. В модели есть три влиятельных наблюдения: 311, 357, 377, – которые, однако, не выходят за пределы доверительных границ на третьем графике. Графики остатков не заставляют усомниться в том, что остатки удовлетворяют условиям Гаусса-Маркова.
Сравнение с kNN
# фактические значения y на тестовой выборке
y.fact <- Carseats[-inTrain, 1]
y.model.lm <- predict(model.4, df.test)
MSE.lm <- sum((y.model.lm - y.fact)^2) / length(y.model.lm)
df.train1 <- df.train
df.test1 <- df.test
df.train1$ShelveLoc <- as.numeric(df.train1$ShelveLoc)
df.test1$ShelveLoc <- as.numeric(df.test1$ShelveLoc)
# kNN требует на вход только числовые переменные
df.train.num <- as.data.frame(apply(df.train1, 2, as.numeric))
df.test.num <- as.data.frame(apply(df.test1, 2, as.numeric))
for (i in 2:50){
model.knn <- knn.reg(train = df.train.num[, !(colnames(df.train.num) %in% 'Sales')],
y = df.train.num[, 'Sales'],
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(3, 3, 1, 1))
# ошибки kNN
# ошибка регрессии
plot(2:50, MSE.knn, ylim = c(4,9), 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 до 50. Ошибка регрессионной модели на тестовой выборке не очень велика и составляет \(\frac {\sqrt{MSE_{TEST}}}{\bar{y}_{TEST}} = 28.5\%\) от среднего значения зависимой переменной. У лучшей модели kNN точность ещё хуже: она ошибается на 34.9% от среднего значения объясняющей переменной.
Однако построенная модель регрессии (\(model.6\)) недостаточно пригодна для прогнозирования, так как её параметры недостатчно хорошо объясняют изменение объёма продаж (\(R^2=0.5583\)).