practice4
A.Lukyanova
16 марта 2018 г
Построим регрессию по данным из пакета Boston. Для этого загрузим пакет MASS, составим обучающую и тестовую выборки, посчитаем описательные статистики и построим модель.
# константы
my.seed <- 12345
train.percent <- 0.85
# преобразуем категориальные переменные в факторы
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)## zn indus chas nox rm
## Min. : 0.00 Min. : 0.46 0:400 Min. :0.3850 Min. :3.561
## 1st Qu.: 0.00 1st Qu.: 5.13 1: 30 1st Qu.:0.4480 1st Qu.:5.889
## Median : 0.00 Median : 8.56 Median :0.5320 Median :6.214
## Mean :11.52 Mean :10.89 Mean :0.5525 Mean :6.280
## 3rd Qu.:16.25 3rd Qu.:18.10 3rd Qu.:0.6240 3rd Qu.:6.628
## Max. :95.00 Max. :27.74 Max. :0.8710 Max. :8.780
## age dis rad tax
## Min. : 2.90 Min. : 1.130 Min. : 1.000 Min. :187
## 1st Qu.: 42.95 1st Qu.: 2.108 1st Qu.: 4.000 1st Qu.:279
## Median : 76.95 Median : 3.325 Median : 5.000 Median :330
## Mean : 68.06 Mean : 3.854 Mean : 9.367 Mean :403
## 3rd Qu.: 94.08 3rd Qu.: 5.287 3rd Qu.:20.000 3rd Qu.:666
## Max. :100.00 Max. :12.127 Max. :24.000 Max. :711
## ptratio black lstat medv
## Min. :12.60 Min. : 0.32 Min. : 1.730 Min. : 5.00
## 1st Qu.:17.07 1st Qu.:376.06 1st Qu.: 6.723 1st Qu.:17.12
## Median :18.85 Median :391.38 Median :11.035 Median :21.40
## Mean :18.42 Mean :358.18 Mean :12.512 Mean :22.70
## 3rd Qu.:20.20 3rd Qu.:396.23 3rd Qu.:16.635 3rd Qu.:25.27
## Max. :22.00 Max. :396.90 Max. :37.970 Max. :50.00
## crim
## Min. : 0.00632
## 1st Qu.: 0.08204
## Median : 0.25278
## Mean : 3.21937
## 3rd Qu.: 3.43597
## Max. :51.13580# совместный график разброса переменных
ggp <- ggpairs(df.train)
print(ggp, progress = F)
#модели
model.1 <- lm(crim ~indus + age + indus:chas + age:chas,
data = df.train)
summary(model.1)##
## Call:
## lm(formula = crim ~ indus + age + indus:chas + age:chas, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.083 -2.375 -0.446 0.860 43.720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.33522 0.71805 -4.645 4.55e-06 ***
## indus 0.39586 0.05540 7.146 3.91e-12 ***
## age 0.03586 0.01299 2.761 0.00601 **
## indus:chas1 -0.14029 0.26705 -0.525 0.59962
## age:chas1 -0.01350 0.04669 -0.289 0.77260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.756 on 425 degrees of freedom
## Multiple R-squared: 0.2613, Adjusted R-squared: 0.2544
## F-statistic: 37.59 on 4 and 425 DF, p-value: < 2.2e-16model.2 <- lm(crim ~indus + age + indus:chas,
data = df.train)
summary(model.2)##
## Call:
## lm(formula = crim ~ indus + age + indus:chas, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.111 -2.317 -0.444 0.863 43.722
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.33973 0.71711 -4.657 4.29e-06 ***
## indus 0.39888 0.05434 7.340 1.09e-12 ***
## age 0.03534 0.01285 2.751 0.0062 **
## indus:chas1 -0.21415 0.07785 -2.751 0.0062 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.75 on 426 degrees of freedom
## Multiple R-squared: 0.2612, Adjusted R-squared: 0.256
## F-statistic: 50.2 on 3 and 426 DF, p-value: < 2.2e-16Проверим остатки модели 2, так как именно в этой модели все коэффициенты оказались значимы. Остатки на гетероскедастичность проверим с помощью теста Бреша-Пагана, а автокорреляцию - теста Дарбина-Уотсона. Нарисуем также график осатков.
# тест Бройша-Пагана
bptest(model.2)##
## studentized Breusch-Pagan test
##
## data: model.2
## BP = 28.92, df = 3, p-value = 2.327e-06# статистика Дарбина-Уотсона
dwtest(model.2)##
## Durbin-Watson test
##
## data: model.2
## DW = 1.9534, p-value = 0.3139
## alternative hypothesis: true autocorrelation is greater than 0# графики остатков
par(mar = c(4.5, 4.5, 2, 1))
par(mfrow = c(1, 3))
plot(model.2, 1)
plot(model.2, 4)
plot(model.2, 5)
par(mfrow = c(1, 1))В остатках присутствует гетероскедастичность и отсутствует автокорреляция.
Следующим шагом будет построение модели knn-регрессии.
# фактические значения y на тестовой выборке
y.fact <- Boston[-inTrain, 1]
y.model.lm <- predict(model.2, 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))
}
}
min(MSE.knn)## [1] 166.5934# график
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))
Ошибка регрессии на тестовой выборке превысила ошибки метода k ближайших соседей. Минимальна она при k=12 и равна 166,59.
Построенная модель обладает низкими показателями качества, поэтому для дальнейшего прогнозирования она не подходит.