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

В данной задаче выполняются следующие пункты:

Модели: множественная линейная регрессия, kNN.
Данные: статистика стоимости жилья в пригороде Бостона.

Цель: исследовать набор данных Boston с помощью линейной регрессионной модели. Задействовав все возможные регрессоры, сделать вывод о пригодности модели для прогноза. Сравнить с методом k ближайших соседей по MSE на тестовой выборке.

# загрузка пакетов
library('GGally')       # графики совместного разброса переменных
library('lmtest')       # тесты остатков регрессионных моделей
library('FNN')          # алгоритм kNN
library('mlbench')
library('MASS')
# константы
my.seed <- 12345
train.percent <- 0.85
data(Boston)            # открываем данные
?Boston
Boston <- Boston[,-c(2,5,6,8,9,10,11,12,13,14)]
Boston$chas <- as.factor(Boston$chas)
# обучающая выборка
set.seed(my.seed)
inTrain <- sample(seq_along(Boston$crim), 
                  nrow(Boston) * train.percent)
df.train <- Boston[inTrain, c(colnames(Boston)[-1], colnames(Boston)[1])]
df.test <- Boston[-inTrain, -1]

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

# описательные статистики по переменным
summary(df.train)
##      indus       chas         age              crim         
##  Min.   : 0.46   0:400   Min.   :  2.90   Min.   : 0.00632  
##  1st Qu.: 5.13   1: 30   1st Qu.: 42.95   1st Qu.: 0.08204  
##  Median : 8.56           Median : 76.95   Median : 0.25278  
##  Mean   :10.89           Mean   : 68.06   Mean   : 3.21937  
##  3rd Qu.:18.10           3rd Qu.: 94.08   3rd Qu.: 3.43597  
##  Max.   :27.74           Max.   :100.00   Max.   :51.13580
# совместный график разброса переменных
ggp <- ggpairs(df.train)
print(ggp, progress = F)

# цвета по фактору chas
ggp <- ggpairs(df.train, mapping = ggplot2::aes(color = chas))
print(ggp, progress = F)

Коробчатые диаграммы на пересечении chas & indus и chas & age показывают, что если район не граничит с рекой (chas = 0), то доля акров некоммерческого бизнеса в городе и доля занятых собственниками, построенных до 1940 года, будут выше. Однако и уровень преступности в такой местности будет значительно превосходить, чем в местности с установленным трактом.

Модели

model <- lm(crim ~ age + indus + chas,
              data = df.train)
summary(model)
## 
## Call:
## lm(formula = crim ~ age + indus + chas, data = df.train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.917  -2.425  -0.506   0.901  43.808 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.22099    0.71799  -4.486 9.34e-06 ***
## age          0.03530    0.01288   2.741  0.00638 ** 
## indus        0.38776    0.05420   7.154 3.71e-12 ***
## chas1       -2.67537    1.09532  -2.443  0.01499 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.76 on 426 degrees of freedom
## Multiple R-squared:  0.2585, Adjusted R-squared:  0.2532 
## F-statistic: 49.49 on 3 and 426 DF,  p-value: < 2.2e-16

В модели все объясняющие переменные ялвяются значимыми. Проверим её остатки.

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

# тест Бройша-Пагана
bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 27.303, df = 3, p-value = 5.085e-06
# статистика Дарбина-Уотсона
dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.9621, p-value = 0.347
## alternative hypothesis: true autocorrelation is greater than 0
# графики остатков
par(mar = c(4.5, 4.5, 2, 1))
par(mfrow = c(1, 3))
plot(model, 1)
plot(model, 4)
plot(model, 5)

par(mfrow = c(1, 1))

Судя по графику слева, остатки не случайны, и их дисперсия непостоянна. В модели есть три влиятельных наблюдения: 405, 411, 415 – которые, однако, не выходят за пределы доверительных границ на третьем графике. Графики остатков заставляют усомниться в том, что остатки удовлетворяют условиям Гаусса-Маркова.

Сравнение с kNN

y.fact <- Boston[-inTrain, 1]
y.model.lm <- predict(model, 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% 'crim')], 
                       y = df.train.num[, 'crim'], 
                       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, ylim = c(0,1000),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 ближайших соседей. С увеличением количества соседей точность kNN практически не изменяется, что говорит о том, что данная регрессионная модель пригодна для прогнозирования.Ошибка регрессионной модели на тестовой выборке составляет \(\frac {\sqrt{MSE_{TEST}}}{\bar{y}_{TEST}} = 251\%\) от среднего значения зависимой переменной.