Упражнение №4
Косован Михаил
07 05 2020
На наборе данных из своего варианта построить модели линейной регрессии с указанными Y и X. Рассмотреть модели с категориальными предикторами, включая их взаимодействие с непрерывными объясняющими переменными. Сгенерировать отчёт по структуре отчёта из практики. Включить в отчёт выводы по каждому из разделов (описание данных, модели, сравнение с kNN). Ответить на вопрос, пригодна ли построенная модель регрессии для прогнозирования и почему.
Вариант 10
Данные: Carseats{ISLR};
Зависимая переменная: Sales;
Объясняющие переменные:
Непрерывные: Price, Advertising;
Дискретные (факторы): US.
# загрузка пакетов
library('GGally') # графики совместного разброса переменных
library('lmtest') # тесты остатков регрессионных моделей
library('FNN') # алгоритм kNN
library('ISLR') # данные Carseats
# Загрузка данных Carseats
data('Carseats')
# Отбор необходимых данных для построения моделей
Carseats <- Carseats[,c('Sales', 'Price', 'Advertising', 'US'),drop=FALSE]
# Константы
my.seed <- 12345
train.percent <- 0.85
# Преобразуем категориальную переменную ShelveLoc в фактор
Carseats$US <- as.factor(Carseats$US)
# Обучающая выборка
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– взимаемая плата за автокресла на каждом участке;
Advertising– локальный рекламный бюджет для компании в каждом месте (в тысячах долларов);US- коэффициент с уровнями Нет и Да, чтобы указать, находится ли магазин в США или нет.
Размерность обучающей выборки: \(n = 340\) строк, \(p = 3\) объясняющих переменных. Зависимая переменная – Sales.
## Price Advertising US Sales
## Min. : 49.00 Min. : 0.000 No :130 Min. : 0.000
## 1st Qu.: 99.75 1st Qu.: 0.000 Yes:210 1st Qu.: 5.390
## Median :117.50 Median : 5.000 Median : 7.420
## Mean :115.96 Mean : 6.238 Mean : 7.453
## 3rd Qu.:131.00 3rd Qu.:11.000 3rd Qu.: 9.160
## Max. :191.00 Max. :29.000 Max. :16.270
Судя по коробчатой диаграмме на пересечении Sales и US, средний объем продаж отличается в зависимости от места нахождения магазина: если магазин расположен в US, то объем продаж выше. Нижний правый график показывает, что доли наблюдений с различными значениями признака USв наборе данных имеют следующий вид: наибольшую часть наблюдений отражает коэффициент со значением Yes.
Модели
##
## Call:
## lm(formula = Sales ~ . + Price:US + Advertising:US, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.8853 -1.5866 -0.0768 1.4929 6.3204
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.081302 1.024836 12.764 < 2e-16 ***
## Price -0.054899 0.008805 -6.235 1.37e-09 ***
## Advertising 0.029777 0.123035 0.242 0.809
## USYes -0.119018 1.392011 -0.086 0.932
## Price:USYes 0.001229 0.011505 0.107 0.915
## Advertising:USYes 0.089074 0.126163 0.706 0.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.387 on 334 degrees of freedom
## Multiple R-squared: 0.2785, Adjusted R-squared: 0.2677
## F-statistic: 25.78 on 5 and 334 DF, p-value: < 2.2e-16Совместное влияние Price:US исключаем, т.к. параметры незначимы (наименее занчимы по сравнению с другими незначимыми коэффициентами).
##
## Call:
## lm(formula = Sales ~ . + Advertising:US, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.8887 -1.5806 -0.0698 1.4970 6.3432
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.999431 0.679720 19.125 <2e-16 ***
## Price -0.054179 0.005659 -9.575 <2e-16 ***
## Advertising 0.029107 0.122694 0.237 0.813
## USYes 0.023846 0.387495 0.062 0.951
## Advertising:USYes 0.089603 0.125880 0.712 0.477
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.384 on 335 degrees of freedom
## Multiple R-squared: 0.2784, Adjusted R-squared: 0.2698
## F-statistic: 32.32 on 4 and 335 DF, p-value: < 2.2e-16Взаимодействие Advertising:US также исключаем.
##
## Call:
## lm(formula = Sales ~ Price + Advertising + US, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.8542 -1.6171 -0.0501 1.5037 6.3834
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.977236 0.678506 19.126 < 2e-16 ***
## Price -0.054398 0.005646 -9.635 < 2e-16 ***
## Advertising 0.114278 0.027122 4.214 3.23e-05 ***
## USYes 0.114888 0.365509 0.314 0.753
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.382 on 336 degrees of freedom
## Multiple R-squared: 0.2773, Adjusted R-squared: 0.2709
## F-statistic: 42.98 on 3 and 336 DF, p-value: < 2.2e-16Параметр US является незначимым, поэтому его тоже исключим из уравнения регресии.
##
## Call:
## lm(formula = Sales ~ Price + Advertising, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.8779 -1.5825 -0.0578 1.5017 6.3534
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.996205 0.674913 19.256 < 2e-16 ***
## Price -0.054263 0.005622 -9.651 < 2e-16 ***
## Advertising 0.120117 0.019734 6.087 3.14e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.379 on 337 degrees of freedom
## Multiple R-squared: 0.2771, Adjusted R-squared: 0.2728
## F-statistic: 64.6 on 2 and 337 DF, p-value: < 2.2e-16В данной модели все коэффициенты оказались значимыми, но по характеристикам качества модель недостаточно хороша (\(R^2=0.2728\)). Пробуем добавить взаимодействие Price:US.
##
## Call:
## lm(formula = Sales ~ Price + Advertising + Price:US, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.8505 -1.6175 -0.0394 1.4932 6.3654
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.044788 0.690849 18.882 < 2e-16 ***
## Price -0.055002 0.006037 -9.111 < 2e-16 ***
## Advertising 0.114186 0.026400 4.325 2.01e-05 ***
## Price:USYes 0.001024 0.003023 0.339 0.735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.382 on 336 degrees of freedom
## Multiple R-squared: 0.2774, Adjusted R-squared: 0.2709
## F-statistic: 42.99 on 3 and 336 DF, p-value: < 2.2e-16Очевидно, стоит остановиться на модели без взаимодействий. Проверим её остатки.
Проверка остатков
# тест Бройша-Пагана
bptest(model.4)##
## studentized Breusch-Pagan test
##
## data: model.4
## BP = 2.0049, df = 2, p-value = 0.367# статистика Дарбина-Уотсона
dwtest(model.4)##
## Durbin-Watson test
##
## data: model.4
## DW = 2.1557, p-value = 0.9242
## 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))Судя по графику слева, остатки случайны, и их дисперсия постоянна. В модели есть три влиятельных наблюдения: 51, 368, 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$US <- as.numeric(df.train1$US)
df.test1$US <- as.numeric(df.test1$US)
# 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-nearest neighbors algorithm', 'regression (all factors)'),
lwd = rep(2, 2))
Как можно видеть по графику, ошибка метода k ближайших соседей с k от 2 до 50 на тестовой выборке меньше, чем ошибка регрессии. Ошибка регрессионной модели на тестовой выборке составляет \(\frac {\sqrt{MSE_{TEST}}}{\bar{y}_{TEST}} = 32.5\%\) от среднего значения зависимой переменной. У kNN точность лучше: она ошибается на 31.5% от среднего значения объясняющей переменной. Однако построенная модель регрессии (\(model.4\)) недостаточно пригодна для прогнозирования, так как её параметры недостаточно хорошо объясняют изменение объёма продаж ($R^2=0.2728).