Линейные регрессионные модели

На наборе данных из своего варианта построить модели линейной регрессии с указанными Y и X. Рассмотреть модели с категориальными предикторами, включая их взаимодействие с непрерывными объясняющими переменными. Сгенерировать отчёт по структуре отчёта из практики. Включить в отчёт выводы по каждому из разделов (описание данных, модели, сравнение с kNN).

Описание переменных

Набор данных “Carseats” содержит переменные:
- зависимая переменная: Sales;
- непрерывные объясняющие переменные: Price, Population;
- дискретные объясняющие переменные: Urban

Подготовка данных

Подключаем необходимые библиотеки:

library('GGally')       # графики совместного разброса переменных
library('lmtest')       # тесты остатков регрессионных моделей
library('FNN')          # алгоритм kNN
library('ISLR')         # набор данных
library('caret')

Константы

my.seed <- 12
train.percent <- 0.85

Обучающая и тестовая выборки

set.seed(my.seed)
data("Carseats")
df <- Carseats[, c('Sales', 'Price', 'Population', 'Urban')]
inTrain <- createDataPartition(df$Sales, p = train.percent, list = FALSE)
df.train<- df[inTrain, ]
df.test <- df[-inTrain, -1]

Размерность обучающей выборки:
n=662 строк, p=4 объясняющих переменных. Зависимая переменная – Sales.

ggp <- ggpairs(df.train, upper = list(combo = 'box'))
print(ggp, progress = F)
## Warning in ggmatrix_gtable(x, ...): Please use the 'progress' parameter in
## your ggmatrix-like function call. See ?ggmatrix_progress for a few examples.
## ggmatrix_gtable 'progress' and 'progress_format' will soon be deprecated.TRUE
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Цвета по фактору “Urban”

ggp <- ggpairs(df.train[, c('Population', 'Urban', 'Sales')], aes(color = Urban), upper = list(combo = 'box'))
print(ggp, progress = F)
## Warning in ggmatrix_gtable(x, ...): Please use the 'progress' parameter in
## your ggmatrix-like function call. See ?ggmatrix_progress for a few examples.
## ggmatrix_gtable 'progress' and 'progress_format' will soon be deprecated.TRUE
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Исходя из представленных графиков, можно сделать вывод, что персональные расходы на обучение выше в государственных учебных заведениях, чем в частных.

Модели

model.1 <- lm(Sales ~ . + Urban:Price + Urban:Population, data = df.train)
summary(model.1)
## 
## Call:
## lm(formula = Sales ~ . + Urban:Price + Urban:Population, data = df.train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4754 -1.8469 -0.1737  1.6859  7.5663 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         15.0086905  1.3267085  11.313  < 2e-16 ***
## Price               -0.0662104  0.0106287  -6.229  1.4e-09 ***
## Population           0.0004792  0.0016183   0.296    0.767    
## UrbanYes            -1.7070461  1.6023692  -1.065    0.287    
## Price:UrbanYes       0.0153275  0.0128167   1.196    0.233    
## Population:UrbanYes -0.0002180  0.0019615  -0.111    0.912    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.532 on 335 degrees of freedom
## Multiple R-squared:  0.2111, Adjusted R-squared:  0.1994 
## F-statistic: 17.93 on 5 and 335 DF,  p-value: 9.268e-16
model.2 <- lm(Sales ~ . + Urban:Population,  data = df.train)
summary(model.2)
## 
## Call:
## lm(formula = Sales ~ . + Urban:Population, data = df.train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5132 -1.7384 -0.1203  1.6894  7.4516 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         13.7989199  0.8589271  16.065   <2e-16 ***
## Price               -0.0556696  0.0059434  -9.367   <2e-16 ***
## Population           0.0004933  0.0016193   0.305    0.761    
## UrbanYes             0.0616813  0.6169519   0.100    0.920    
## Population:UrbanYes -0.0002443  0.0019627  -0.124    0.901    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.533 on 336 degrees of freedom
## Multiple R-squared:  0.2078, Adjusted R-squared:  0.1983 
## F-statistic: 22.03 on 4 and 336 DF,  p-value: 3.656e-16
model.3 <- lm(Sales ~ . + Urban:Price,  data = df.train)
summary(model.3)
## 
## Call:
## lm(formula = Sales ~ . + Urban:Price, data = df.train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4812 -1.8485 -0.1808  1.6873  7.5621 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    15.0516089  1.2673712  11.876  < 2e-16 ***
## Price          -0.0662189  0.0106128  -6.240 1.32e-09 ***
## Population      0.0003308  0.0009131   0.362    0.717    
## UrbanYes       -1.7687086  1.5009941  -1.178    0.239    
## Price:UrbanYes  0.0153434  0.0127971   1.199    0.231    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.528 on 336 degrees of freedom
## Multiple R-squared:  0.2111, Adjusted R-squared:  0.2017 
## F-statistic: 22.48 on 4 and 336 DF,  p-value: < 2.2e-16

Остановимся на модели 3. Проверим её остатки.

Проверка остатков

Тест Бройша-Пагана

bptest(model.3)
## 
##  studentized Breusch-Pagan test
## 
## data:  model.3
## BP = 2.0751, df = 4, p-value = 0.7219

Статистика Дарбина-Уотсона

dwtest(model.3)
## 
##  Durbin-Watson test
## 
## data:  model.3
## DW = 1.9421, p-value = 0.2975
## alternative hypothesis: true autocorrelation is greater than 0

Графики остатков

par(mar = c(4.5, 4.5, 2, 1))
par(mfrow = c(1, 3))
plot(model.3, 1)
plot(model.3, 4)
plot(model.3, 5)

Сравнение с kNN

Фактические значения y на тестовой выборке

y.fact <- df[-inTrain,]$Sales
y.model.lm <- predict(model.3, df.test)
MSE.lm <- sum((y.model.lm - y.fact)^2) / length(y.model.lm)

kNN требует на вход только числовые переменные

df.train$Urban = as.numeric(as.factor(df.train$Urban))
df.test$Urban = as.numeric(as.factor(df.test$Urban))
df.train.num <- as.data.frame(df.train)
df.test.num <- as.data.frame(df.test)
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(4.5, 4.5, 1, 1))
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))