Data Science: инструментарий и жизненный цикл проекта
Глушков Егор Александрович, гр. 20.М04-мм
data(mtcars)
## Переменная Описание Шкала
## 1 mpg Миль/Галлон Отношений
## 2 cyl Число цилиндров Интервальная
## 3 disp Рабочий объем (куб. дюймов) Отношений
## 4 hp Полная мощность (лошадиных сил) Отношений
## 5 drat Передаточное отношение задней оси Отношений
## 6 wt Вес (1000 фунтов) Отношений
## 7 qsec Время 1/4 мили Отношений
## 8 vs Двигатель (0=V-образный, 1=прямой) Номинальная (бинарная)
## 9 am Трансмиссия (0=автоматическая, 1=ручная) Номинальная (бинарная)
## 10 gear Количество передач переднего хода Интервальная
## 11 carb Количество карбюраторов Интервальная
mtcars2 <- within(mtcars, {
vs <- factor(vs, labels = c("V", "S"))
am <- factor(am, labels = c("automatic", "manual"))
cyl <- ordered(cyl)
gear <- ordered(gear)
carb <- ordered(carb)
})
summary(mtcars2)
## mpg cyl disp hp drat
## Min. :10.40 4:11 Min. : 71.1 Min. : 52.0 Min. :2.760
## 1st Qu.:15.43 6: 7 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080
## Median :19.20 8:14 Median :196.3 Median :123.0 Median :3.695
## Mean :20.09 Mean :230.7 Mean :146.7 Mean :3.597
## 3rd Qu.:22.80 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920
## Max. :33.90 Max. :472.0 Max. :335.0 Max. :4.930
## wt qsec vs am gear carb
## Min. :1.513 Min. :14.50 V:18 automatic:19 3:15 1: 7
## 1st Qu.:2.581 1st Qu.:16.89 S:14 manual :13 4:12 2:10
## Median :3.325 Median :17.71 5: 5 3: 3
## Mean :3.217 Mean :17.85 4:10
## 3rd Qu.:3.610 3rd Qu.:18.90 6: 1
## Max. :5.424 Max. :22.90 8: 1
head(mtcars2)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 V manual 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 V manual 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 S manual 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 S automatic 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 V automatic 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 S automatic 3 1
Встроенный набор данных mtcars содержит \(32\) наблюдения по \(11\) признакам. Некоторые признаки перекодируем в соответствующие шкалы (интервальные, бинарные).
Проверим однородность подвыборок (признак Миль/Галлон mt) по фактору Тип двигателя (vs). Построим гистограмму частот с разбиением по указанному фактору.
library(ggplot2)
ggplot(data=mtcars2, aes(x=mpg, fill=vs)) +
geom_histogram(binwidth = 1.25) +
labs(title="Разбиение данных",
x="Миль/Галлон (mpg)", y="Количество") +
scale_fill_discrete(name = "Двигатель (vs)", labels = c("V-образный", "Прямой (S)"))
Проиллюстрируем распределения в подвыборках с помощью “ящика с усами” (boxplot):
mtc <- data.frame(group=mtcars2$vs, X=mtcars2$mpg)
boxplot(X ~ group, mtc, xlab = "Тип двигателя", ylab = "Миль/Галлон")
Проверим нормальность данных в подвыборках:
p.Sh <- with(mtc, tapply(X, group, function(x)shapiro.test(x)$p.value)); p.Sh
## V S
## 0.4491490 0.1665898
Согласие с нормальным распределением для обеих групп не отвергается в соответствии с критерием Шапиро-Уилка с \(p.value = 0.449\) и \(0.167\) при уровне значимости \(\alpha = 0.05\).
p.F <- var.test(X~group, mtc)$p.value; p.F
## [1] 0.1997209
В соответствии с критерием Фишера гипотеза о равенстве дисперсий двух групп не отвергается (\(p.value=0.1997\)), значит, можно использовать критерий Стьюдента (с параметром о равенстве дисперсий).
p.T <- t.test(X~group, mtc, var.equal=TRUE); p.T
##
## Two Sample t-test
##
## data: X by group
## t = -4.8644, df = 30, p-value = 3.416e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -11.274221 -4.606732
## sample estimates:
## mean in group V mean in group S
## 16.61667 24.55714
По критерию Стьюдента гипотеза о равенстве средних не принимается (\(p.value=3.415937 \times 10^{-5}\)) на уровне значимости \(\alpha = 0.05\), то есть различие средних в двух группах (различающихся по типу двигателя) значимо: автомобили с V-образным двигателем (\(16.62\) миль/галлон) действительно проезжают в среднем меньше, чем автомобили с прямым двигателем (\(24.56\)).
Построим модель множественной линейной регресии, зависящей от всех имеющихся переменных:
init.model <- lm(mpg~., data=mtcars)
summary(init.model)
##
## Call:
## lm(formula = mpg ~ ., data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4506 -1.6044 -0.1196 1.2193 4.6271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.30337 18.71788 0.657 0.5181
## cyl -0.11144 1.04502 -0.107 0.9161
## disp 0.01334 0.01786 0.747 0.4635
## hp -0.02148 0.02177 -0.987 0.3350
## drat 0.78711 1.63537 0.481 0.6353
## wt -3.71530 1.89441 -1.961 0.0633 .
## qsec 0.82104 0.73084 1.123 0.2739
## vs 0.31776 2.10451 0.151 0.8814
## am 2.52023 2.05665 1.225 0.2340
## gear 0.65541 1.49326 0.439 0.6652
## carb -0.19942 0.82875 -0.241 0.8122
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.65 on 21 degrees of freedom
## Multiple R-squared: 0.869, Adjusted R-squared: 0.8066
## F-statistic: 13.93 on 10 and 21 DF, p-value: 3.793e-07
Согласно \(p.value=3.793 \times 10^{-7}\) с уровнем значимости \(\alpha = 0.001\) не принимается гипотеза о том, что данная модель не зависит от переменных (предикторов). При этом лишь переменная wt отлична от нуля с уровнем значимости \(0.1\). Коэффициент детерминации равен \(0.869\), скорректированный – \(0.8066\).
Произведем пошаговую подгонку регрессии, когда на каждом шаге добавляется наиболее значимая переменная (в соответствии с информационным критерием Акаике AIC).
library(MASS)
empty_m <- lm(mpg ~ 1, mtcars)
SA.f <- stepAIC(empty_m, direction = 'forward', scope=list(lower=empty_m, upper=~cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb))
## Start: AIC=115.94
## mpg ~ 1
##
## Df Sum of Sq RSS AIC
## + wt 1 847.73 278.32 73.217
## + cyl 1 817.71 308.33 76.494
## + disp 1 808.89 317.16 77.397
## + hp 1 678.37 447.67 88.427
## + drat 1 522.48 603.57 97.988
## + vs 1 496.53 629.52 99.335
## + am 1 405.15 720.90 103.672
## + carb 1 341.78 784.27 106.369
## + gear 1 259.75 866.30 109.552
## + qsec 1 197.39 928.66 111.776
## <none> 1126.05 115.943
##
## Step: AIC=73.22
## mpg ~ wt
##
## Df Sum of Sq RSS AIC
## + cyl 1 87.150 191.17 63.198
## + hp 1 83.274 195.05 63.840
## + qsec 1 82.858 195.46 63.908
## + vs 1 54.228 224.09 68.283
## + carb 1 44.602 233.72 69.628
## + disp 1 31.639 246.68 71.356
## <none> 278.32 73.217
## + drat 1 9.081 269.24 74.156
## + gear 1 1.137 277.19 75.086
## + am 1 0.002 278.32 75.217
##
## Step: AIC=63.2
## mpg ~ wt + cyl
##
## Df Sum of Sq RSS AIC
## + hp 1 14.5514 176.62 62.665
## + carb 1 13.7724 177.40 62.805
## <none> 191.17 63.198
## + qsec 1 10.5674 180.60 63.378
## + gear 1 3.0281 188.14 64.687
## + disp 1 2.6796 188.49 64.746
## + vs 1 0.7059 190.47 65.080
## + am 1 0.1249 191.05 65.177
## + drat 1 0.0010 191.17 65.198
##
## Step: AIC=62.66
## mpg ~ wt + cyl + hp
##
## Df Sum of Sq RSS AIC
## <none> 176.62 62.665
## + am 1 6.6228 170.00 63.442
## + disp 1 6.1762 170.44 63.526
## + carb 1 2.5187 174.10 64.205
## + drat 1 2.2453 174.38 64.255
## + qsec 1 1.4010 175.22 64.410
## + gear 1 0.8558 175.76 64.509
## + vs 1 0.0599 176.56 64.654
summary(SA.f)
##
## Call:
## lm(formula = mpg ~ wt + cyl + hp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9290 -1.5598 -0.5311 1.1850 5.8986
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.75179 1.78686 21.687 < 2e-16 ***
## wt -3.16697 0.74058 -4.276 0.000199 ***
## cyl -0.94162 0.55092 -1.709 0.098480 .
## hp -0.01804 0.01188 -1.519 0.140015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.512 on 28 degrees of freedom
## Multiple R-squared: 0.8431, Adjusted R-squared: 0.8263
## F-statistic: 50.17 on 3 and 28 DF, p-value: 2.184e-11
Итоговая модель построена на трёх переменных: wt, cyl, hp, первые две из которых значимо отличны от нуля (с уровнями значимости \(0.001\), \(0.1\) соответственно). Сама модель также значимо зависит от переменных с \(p.value=2.18 \times 10^{-11}\).
Коэффициент детерминации равен \(0.8431\), скорректированный – \(0.8263\).
Модель с “обратным ходом”, исключающая переменные:
SA.b <- stepAIC(init.model, direction = 'backward')
## Start: AIC=70.9
## mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb
##
## Df Sum of Sq RSS AIC
## - cyl 1 0.0799 147.57 68.915
## - vs 1 0.1601 147.66 68.932
## - carb 1 0.4067 147.90 68.986
## - gear 1 1.3531 148.85 69.190
## - drat 1 1.6270 149.12 69.249
## - disp 1 3.9167 151.41 69.736
## - hp 1 6.8399 154.33 70.348
## - qsec 1 8.8641 156.36 70.765
## <none> 147.49 70.898
## - am 1 10.5467 158.04 71.108
## - wt 1 27.0144 174.51 74.280
##
## Step: AIC=68.92
## mpg ~ disp + hp + drat + wt + qsec + vs + am + gear + carb
##
## Df Sum of Sq RSS AIC
## - vs 1 0.2685 147.84 66.973
## - carb 1 0.5201 148.09 67.028
## - gear 1 1.8211 149.40 67.308
## - drat 1 1.9826 149.56 67.342
## - disp 1 3.9009 151.47 67.750
## - hp 1 7.3632 154.94 68.473
## <none> 147.57 68.915
## - qsec 1 10.0933 157.67 69.032
## - am 1 11.8359 159.41 69.384
## - wt 1 27.0280 174.60 72.297
##
## Step: AIC=66.97
## mpg ~ disp + hp + drat + wt + qsec + am + gear + carb
##
## Df Sum of Sq RSS AIC
## - carb 1 0.6855 148.53 65.121
## - gear 1 2.1437 149.99 65.434
## - drat 1 2.2139 150.06 65.449
## - disp 1 3.6467 151.49 65.753
## - hp 1 7.1060 154.95 66.475
## <none> 147.84 66.973
## - am 1 11.5694 159.41 67.384
## - qsec 1 15.6830 163.53 68.200
## - wt 1 27.3799 175.22 70.410
##
## Step: AIC=65.12
## mpg ~ disp + hp + drat + wt + qsec + am + gear
##
## Df Sum of Sq RSS AIC
## - gear 1 1.565 150.09 63.457
## - drat 1 1.932 150.46 63.535
## <none> 148.53 65.121
## - disp 1 10.110 158.64 65.229
## - am 1 12.323 160.85 65.672
## - hp 1 14.826 163.35 66.166
## - qsec 1 26.408 174.94 68.358
## - wt 1 69.127 217.66 75.350
##
## Step: AIC=63.46
## mpg ~ disp + hp + drat + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - drat 1 3.345 153.44 62.162
## - disp 1 8.545 158.64 63.229
## <none> 150.09 63.457
## - hp 1 13.285 163.38 64.171
## - am 1 20.036 170.13 65.466
## - qsec 1 25.574 175.67 66.491
## - wt 1 67.572 217.66 73.351
##
## Step: AIC=62.16
## mpg ~ disp + hp + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - disp 1 6.629 160.07 61.515
## <none> 153.44 62.162
## - hp 1 12.572 166.01 62.682
## - qsec 1 26.470 179.91 65.255
## - am 1 32.198 185.63 66.258
## - wt 1 69.043 222.48 72.051
##
## Step: AIC=61.52
## mpg ~ hp + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - hp 1 9.219 169.29 61.307
## <none> 160.07 61.515
## - qsec 1 20.225 180.29 63.323
## - am 1 25.993 186.06 64.331
## - wt 1 78.494 238.56 72.284
##
## Step: AIC=61.31
## mpg ~ wt + qsec + am
##
## Df Sum of Sq RSS AIC
## <none> 169.29 61.307
## - am 1 26.178 195.46 63.908
## - qsec 1 109.034 278.32 75.217
## - wt 1 183.347 352.63 82.790
summary(SA.b)
##
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4811 -1.5555 -0.7257 1.4110 4.6610
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.6178 6.9596 1.382 0.177915
## wt -3.9165 0.7112 -5.507 6.95e-06 ***
## qsec 1.2259 0.2887 4.247 0.000216 ***
## am 2.9358 1.4109 2.081 0.046716 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.459 on 28 degrees of freedom
## Multiple R-squared: 0.8497, Adjusted R-squared: 0.8336
## F-statistic: 52.75 on 3 and 28 DF, p-value: 1.21e-11
Модель с “обратным ходом” значимо (\(p.value=1.21 \times 10^{-11}\)) зависит от трех переменных (wt, qsec, am), каждая из которых является ненулевой при уровне значимости \(\alpha = 0.05\). Коэффициент детерминации равен \(0.8497\), скорректированный – \(0.8336\).
Таким образом, последняя модель множественной линейной регрессии, исключающая переменные, показывает наилучший результат (в соответствии со скорректированным коэффициентом детерминации).