library(MASS)
library(psych)
attach(Boston)
library(corrplot)
crim01 <- ifelse(Boston$crim > median(Boston$crim), 1, 0)
mydf <- data.frame(Boston, crim01)
pairs.panels(mydf)
corrplot(cor(mydf), order = "hclust")
На графике представлены корреляции между переменными. В качестве первого набора данных выберем все переменные, кроме crim. Для построения моделей в качестве второго набора выберем те переменные, у которых корреляция с crim01 наибольшая. Это indus, nox, age, dis, rad и tax.
set.seed(1234)
trainNom <- sample(1:nrow(mydf), nrow(mydf)*0.7 , replace=F)
train <- mydf[trainNom,]
test <- mydf[-trainNom,]
glm.fit <- glm(crim01~. - crim, data = train, family = binomial)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(glm.fit)
##
## Call:
## glm(formula = crim01 ~ . - crim, family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7090 -0.1391 -0.0007 0.0018 3.6533
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -30.045086 9.511237 -3.159 0.001584 **
## zn -0.087410 0.043650 -2.003 0.045229 *
## indus -0.084801 0.053913 -1.573 0.115740
## chas 1.689382 0.922114 1.832 0.066940 .
## nox 49.176903 9.179087 5.357 8.44e-08 ***
## rm -0.325560 0.854404 -0.381 0.703176
## age 0.041550 0.016204 2.564 0.010342 *
## dis 0.797573 0.300782 2.652 0.008010 **
## rad 0.644386 0.176346 3.654 0.000258 ***
## tax -0.007404 0.003499 -2.116 0.034351 *
## ptratio 0.589382 0.171444 3.438 0.000587 ***
## black -0.038589 0.017479 -2.208 0.027264 *
## lstat -0.005279 0.061415 -0.086 0.931496
## medv 0.158285 0.075884 2.086 0.036990 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 490.19 on 353 degrees of freedom
## Residual deviance: 142.09 on 340 degrees of freedom
## AIC: 170.09
##
## Number of Fisher Scoring iterations: 9
glm.prob <- predict(glm.fit, test, type = "response")
glm.pred <- ifelse(glm.prob > 0.5, 1,0)
table(glm.pred, test$crim01)
##
## glm.pred 0 1
## 0 64 10
## 1 5 73
mean(glm.pred != test$crim01)
## [1] 0.09868421
Ошибка предсказания для модели логистической регрессии для первого набора данных равна 9,87%.
glm.fit <- glm(crim01~ indus + nox + age + dis + rad + tax, data = train, family = binomial)
summary(glm.fit)
##
## Call:
## glm(formula = crim01 ~ indus + nox + age + dis + rad + tax, family = binomial,
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.90600 -0.27052 -0.04253 0.00717 2.73118
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.24721 4.47290 -4.974 6.57e-07 ***
## indus -0.05260 0.04922 -1.069 0.28525
## nox 37.45549 8.33326 4.495 6.97e-06 ***
## age 0.02092 0.01084 1.930 0.05366 .
## dis 0.14106 0.18028 0.782 0.43393
## rad 0.57337 0.13321 4.304 1.68e-05 ***
## tax -0.00761 0.00291 -2.615 0.00893 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 490.19 on 353 degrees of freedom
## Residual deviance: 175.83 on 347 degrees of freedom
## AIC: 189.83
##
## Number of Fisher Scoring iterations: 8
glm.prob <- predict(glm.fit, test, type = "response")
glm.pred <- ifelse(glm.prob > 0.5, 1,0)
table(glm.pred, test$crim01)
##
## glm.pred 0 1
## 0 66 14
## 1 3 69
mean(glm.pred != test$crim01)
## [1] 0.1118421
Ошибка предсказания для второго набора данных составляет 11,18%.
lda.fit <- lda(crim01~ . -crim, data = train)
lda.pred <- predict(lda.fit, test)
table(lda.pred$class, test$crim01)
##
## 0 1
## 0 66 22
## 1 3 61
mean(lda.pred$class != test$crim01)
## [1] 0.1644737
Ошибка предсказания для модели LDA для первого набора данных равна 16,45%.
lda.fit <- lda(crim01~indus + nox + age + dis + rad + tax, data = train)
lda.pred <- predict(lda.fit, test)
table(lda.pred$class, test$crim01)
##
## 0 1
## 0 67 21
## 1 2 62
mean(lda.pred$class != test$crim01)
## [1] 0.1513158
Ошибка предсказания для второго набора данных составляет 15,13%.
Для полного набора данных
library(class)
train.stan <- as.data.frame(scale(train[,-15]))
test.stan <- as.data.frame(scale(test[,-15]))
train.X <- train.stan[,-1]
test.X <- test.stan[,-1]
set.seed(132)
knn.pred <- knn(train.X, test.X, train$crim01, k=1)
mean(knn.pred != test$crim01)
## [1] 0.07894737
knn.pred <- knn(train.X, test.X, train$crim01, k=10)
mean(knn.pred != test$crim01)
## [1] 0.1513158
knn.pred <- knn(train.X, test.X, train$crim01, k=50)
mean(knn.pred != test$crim01)
## [1] 0.2105263
knn.pred <- knn(train.X, test.X, train$crim01, k=100)
mean(knn.pred != test$crim01)
## [1] 0.2171053
Наименьшая ошибка предсказания составлет 7,8% и достигается при k = 1.
Рассмотрим теперь модель для второго набора параметров.
train.X <- cbind(train.stan$age, train.stan$dis, train.stan$indus, train.stan$nox, train.stan$rad, train.stan$tax)
test.X <- cbind(test.stan$age, test.stan$dis, test.stan$indus, test.stan$nox, test.stan$rad, test.stan$tax)
knn.pred <- knn(train.X, test.X, train$crim01, k=1)
mean(knn.pred != test$crim01)
## [1] 0.07894737
knn.pred <- knn(train.X, test.X, train$crim01, k=10)
mean(knn.pred != test$crim01)
## [1] 0.1052632
knn.pred <- knn(train.X, test.X, train$crim01, k=50)
mean(knn.pred != test$crim01)
## [1] 0.1578947
knn.pred <- knn(train.X, test.X, train$crim01, k=100)
mean(knn.pred != test$crim01)
## [1] 0.1842105
Наименьшая ошибка предсказания составлет 7,8% и достигается при k = 1.
Наилучшей моделью в смысле наименьшей ошибки предсказания для наших данный оказалась KNN model при k=1. Ошибка совпадает для полного и неполного набора предикторов и составляет 7.8%