Task 2. Chaper 4, ex 13

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,]

Logistic Regression model

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 model

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%.

KNN model

Для полного набора данных

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%