Practice4
Romanova ER
16 марта 2018 г
Математическое моделирование
Практика 4
Линейные регрессионные модели
- оценка параметров модели линейной регрессии
- построение моделей со взаимосдействием переменных
- визуализирование и тестирование остатков модели линейной регрессии Модели: множественная линейная регрессия, kNN. Данные: данные об автомобилях.
Цель: исследовать набор данных Auto с помощью линейной регрессионной модели. Задействовав все возможные регрессоры, сделать вывод о пригодности модели для прогноза. Сравнить с методом k ближайших соседей по MSE на тестовой выборке.
# загрузка пакетов
library('GGally') # графики совместного разброса переменных## Warning: package 'GGally' was built under R version 3.4.3library('lmtest') # тесты остатков регрессионных моделей## Warning: package 'lmtest' was built under R version 3.4.2## Loading required package: zoo## Warning: package 'zoo' was built under R version 3.4.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary("ISLR") ## Warning: package 'ISLR' was built under R version 3.4.3data(Auto)
library('GGally') # графики совместного разброса переменных
library('lmtest') # тесты остатков регрессионных моделей
library('FNN') ## Warning: package 'FNN' was built under R version 3.4.3# константы
my.seed <- 12345
train.percent <- 0.85
str(Auto)## 'data.frame': 392 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...Auto <- Auto[,-4] #исключаем переменные, которые не участвуют в рассчетах
Auto <- Auto[,-6]
Auto <- Auto[,-6]
Auto <- Auto[,-6]
Auto$cylinders <- as.factor(Auto$cylinders)
# обучающая выборка
set.seed(my.seed)
inTrain <- sample(seq_along(Auto$mpg),
nrow(Auto) * train.percent)
df.train <- Auto[inTrain,]
df.test <- Auto[-inTrain,] Описание переменных
Набор данных Auto содержит переменные:
- mpg – миля на галлон;
- cylinders – количество цилиндров от 4 до 8;
- displacement – Объем двигателя;
- weight - Масса автомобиля (кг.);
- acceleration - Время для разгона от 0 до 60 миль / ч (сек.)
Размерность обучающей выборки: n=392 строк, p=4 объясняющих переменных. Зависимая переменная – mpg.
summary(Auto)## mpg cylinders displacement weight acceleration
## Min. : 9.00 3: 4 Min. : 68.0 Min. :1613 Min. : 8.00
## 1st Qu.:17.00 4:199 1st Qu.:105.0 1st Qu.:2225 1st Qu.:13.78
## Median :22.75 5: 3 Median :151.0 Median :2804 Median :15.50
## Mean :23.45 6: 83 Mean :194.4 Mean :2978 Mean :15.54
## 3rd Qu.:29.00 8:103 3rd Qu.:275.8 3rd Qu.:3615 3rd Qu.:17.02
## Max. :46.60 Max. :455.0 Max. :5140 Max. :24.80# совместный график разброса переменных
ggp <- ggpairs(df.train)
print(ggp, progress = F)## `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`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Модели
model.1 <- lm(mpg ~ . + displacement:cylinders + weight:cylinders + acceleration:cylinders,
data = df.train)
summary(model.1)##
## Call:
## lm(formula = mpg ~ . + displacement:cylinders + weight:cylinders +
## acceleration:cylinders, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.4561 -2.2412 -0.0288 1.8979 14.7482
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -38.29971 43.59766 -0.878 0.38035
## cylinders4 85.49236 43.69160 1.957 0.05126 .
## cylinders5 -15.75152 59.62198 -0.264 0.79181
## cylinders6 76.80410 43.87143 1.751 0.08097 .
## cylinders8 55.25845 43.36275 1.274 0.20348
## displacement 0.01977 0.01407 1.404 0.16120
## weight 0.02155 0.01898 1.135 0.25709
## acceleration 0.67784 0.26250 2.582 0.01027 *
## cylinders4:displacement -0.04929 0.02974 -1.657 0.09843 .
## cylinders5:displacement -0.70720 0.27732 -2.550 0.01124 *
## cylinders6:displacement -0.06504 0.02253 -2.887 0.00416 **
## cylinders8:displacement NA NA NA NA
## cylinders4:weight -0.02877 0.01904 -1.511 0.13188
## cylinders5:weight 0.03274 0.03147 1.040 0.29905
## cylinders6:weight -0.02335 0.01909 -1.223 0.22216
## cylinders8:weight -0.02584 0.01904 -1.358 0.17556
## cylinders4:acceleration -0.54152 0.29379 -1.843 0.06623 .
## cylinders5:acceleration NA NA NA NA
## cylinders6:acceleration -0.86653 0.39141 -2.214 0.02755 *
## cylinders8:acceleration NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.968 on 316 degrees of freedom
## Multiple R-squared: 0.7627, Adjusted R-squared: 0.7507
## F-statistic: 63.49 on 16 and 316 DF, p-value: < 2.2e-16Совместное влияние displacement:cylinders и acceleration:cylinders исключаем, т.к. параметры незначимы и недостаточно наблюдений для оценки одного из них.
model.1 <- lm(mpg ~ . + weight:cylinders,
data = df.train)
summary(model.1)##
## Call:
## lm(formula = mpg ~ . + weight:cylinders, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6493 -2.3061 -0.2492 2.0121 14.9319
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -25.141458 44.169736 -0.569 0.570
## cylinders4 72.120890 44.157824 1.633 0.103
## cylinders5 53.202929 50.223675 1.059 0.290
## cylinders6 61.333365 44.404557 1.381 0.168
## cylinders8 53.608265 44.284268 1.211 0.227
## displacement -0.014215 0.009132 -1.557 0.121
## weight 0.019785 0.019372 1.021 0.308
## acceleration 0.156472 0.106510 1.469 0.143
## cylinders4:weight -0.027790 0.019401 -1.432 0.153
## cylinders5:weight -0.020285 0.020868 -0.972 0.332
## cylinders6:weight -0.024750 0.019458 -1.272 0.204
## cylinders8:weight -0.022378 0.019404 -1.153 0.250
##
## Residual standard error: 4.054 on 321 degrees of freedom
## Multiple R-squared: 0.7484, Adjusted R-squared: 0.7398
## F-statistic: 86.81 on 11 and 321 DF, p-value: < 2.2e-16Совместное влияние weight:cylinders тоже исключаем, так как парамеры незначимы
model.1 <- lm(mpg ~ . ,
data = df.train)
summary(model.1)##
## Call:
## lm(formula = mpg ~ ., data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6498 -2.5516 -0.2786 2.4806 15.4552
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.0711064 3.4305481 9.349 < 2e-16 ***
## cylinders4 8.8997965 3.0113810 2.955 0.00335 **
## cylinders5 11.1895234 3.8830294 2.882 0.00422 **
## cylinders6 5.4875247 3.2051543 1.712 0.08783 .
## cylinders8 7.9332194 3.5608808 2.228 0.02657 *
## displacement -0.0142011 0.0092457 -1.536 0.12552
## weight -0.0053782 0.0007964 -6.753 6.67e-11 ***
## acceleration 0.1532485 0.1053885 1.454 0.14688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.15 on 325 degrees of freedom
## Multiple R-squared: 0.7331, Adjusted R-squared: 0.7274
## F-statistic: 127.5 on 7 and 325 DF, p-value: < 2.2e-16Параметр acceleration не значим, исключаем его
model.1 <- lm(mpg ~ cylinders + weight + displacement,
data = df.train)
summary(model.1)##
## Call:
## lm(formula = mpg ~ cylinders + weight + displacement, data = df.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.3050 -2.7124 -0.2882 2.4329 15.6107
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.5795792 3.2755266 10.252 < 2e-16 ***
## cylinders4 9.6282232 2.9744953 3.237 0.00133 **
## cylinders5 12.1074918 3.8379138 3.155 0.00176 **
## cylinders6 6.3470661 3.1555564 2.011 0.04511 *
## cylinders8 8.5691633 3.5399582 2.421 0.01604 *
## weight -0.0050207 0.0007588 -6.617 1.51e-10 ***
## displacement -0.0188922 0.0086794 -2.177 0.03022 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.157 on 326 degrees of freedom
## Multiple R-squared: 0.7314, Adjusted R-squared: 0.7264
## F-statistic: 147.9 on 6 and 326 DF, p-value: < 2.2e-16Все параметры модели значимы, R-квадрат так же достаточно высок, чтобы говорить о прогностической пригодности модели.
Проверка остатков
bptest(model.1)##
## studentized Breusch-Pagan test
##
## data: model.1
## BP = 35.525, df = 6, p-value = 3.408e-06# статистика Дарбина-Уотсона
dwtest(model.1)##
## Durbin-Watson test
##
## data: model.1
## DW = 2.1989, p-value = 0.9658
## 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, 1)
plot(model.1, 4)
plot(model.1, 5)
par(mfrow = c(1, 1))Сравнение с kNN
# фактические значения y на тестовой выборке
y.fact <- Auto[-inTrain, 1]
y.model.lm <- predict(model.1, 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% 'mpg')],
y = df.train.num[, 'mpg'],
test = df.test.num[, !(colnames(df.train.num) %in% 'mpg')], 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('topright', lty = c(1, 2), pch = c(1, NA),
col = c('darkgreen', grey(0.2)),
legend = c('k ближайших соседа', 'регрессия (все факторы)'),
lwd = rep(2, 2))
Как можно видеть по графику, ошибка регрессии на тестовой выборке меньше, чем ошибка метода k ближайших соседей с k от 2 до 7. Далее с увеличением количества соседей точность kNN скачет, и начиная с 25 MSE на тестовой выборке kNN модели уменьшается. В целом, ошибки на тестовой выборке у моделей двух типов невелики сами по себе. Для малого числа соседей лучше использовать регрессию, а для большего - метод kNN