y_hats <- sapply(fits, function(fit) {
factor(predict(fit, mnist_27$test))
})
accs <- colMeans(y_hats == mnist_27$test$y)
t(t(accs))
| glm | 0.750 |
|---|---|
| lda | 0.750 |
| naive_bayes | 0.795 |
| svmLinear | 0.755 |
| knn | 0.840 |
| gamLoess | 0.845 |
| multinom | 0.750 |
| qda | 0.820 |
| rf | 0.780 |
| adaboost | 0.775 |
pred <- rowMeans(y_hats == 7)
y_hat_ens <- ifelse(pred > 0.5, 7, 2)
y_hat_ens <- factor(y_hat_ens)
confusionMatrix(data = y_hat_ens, reference = mnist_27$test$y)
Confusion Matrix and Statistics
Reference
Prediction 2 7
2 88 20
7 18 74
Accuracy : 0.81
95% CI : (0.7487, 0.8619)
No Information Rate : 0.53
P-Value [Acc > NIR] : <2e-16
Kappa : 0.6182
Mcnemar's Test P-Value : 0.8711
Sensitivity : 0.8302
Specificity : 0.7872
Pos Pred Value : 0.8148
Neg Pred Value : 0.8043
Prevalence : 0.5300
Detection Rate : 0.4400
Detection Prevalence : 0.5400
Balanced Accuracy : 0.8087
'Positive' Class : 2
confusionMatrix(data = y_hat_ens, reference = mnist_27$test$y)$overall["Accuracy"]
We can notice that there are some techniques that achieved individual accuracy better than the ensemble accuracy, with all techniques. These were the knn, gamLoess and qda techniques.
In order to find a better ensemble, we cannot use the accuracy over the test set. We will find the accuracy estimates, obtained from cross validation on the training set.
acc_est <- sapply(fits, function(fit) {
min(fit$results$Accuracy)
})
t(t(acc_est))
| glm | 0.7933559 |
|---|---|
| lda | 0.7887980 |
| naive_bayes | 0.8140575 |
| svmLinear | 0.7951638 |
| knn | 0.8055268 |
| gamLoess | 0.8386924 |
| multinom | 0.7950913 |
| qda | 0.8391139 |
| rf | 0.8050582 |
| adaboost | 0.7967890 |
Now we build a new ensemble only with techniques with estimated accuracy greater or equal to 0.8.
ind <- acc_est >= 0.8
pred <- rowMeans(y_hats[, ind] == 7)
y_hat_ens <- ifelse(pred > 0.5, 7, 2)
y_hat_ens <- factor(y_hat_ens)
confusionMatrix(data = y_hat_ens, reference = mnist_27$test$y)
Confusion Matrix and Statistics
Reference
Prediction 2 7
2 90 18
7 16 76
Accuracy : 0.83
95% CI : (0.7706, 0.8793)
No Information Rate : 0.53
P-Value [Acc > NIR] : <2e-16
Kappa : 0.6584
Mcnemar's Test P-Value : 0.8638
Sensitivity : 0.8491
Specificity : 0.8085
Pos Pred Value : 0.8333
Neg Pred Value : 0.8261
Prevalence : 0.5300
Detection Rate : 0.4500
Detection Prevalence : 0.5400
Balanced Accuracy : 0.8288
'Positive' Class : 2
confusionMatrix(data = y_hat_ens, reference = mnist_27$test$y)$overall["Accuracy"]
In the case of outcomes 2 and 7, the Local Weighted Regression had the highest accuracy. This shows that the models obtained depend on random samples and depending on the case, a particular solution may have better accuracy.
The train set of the full data compromises of 60,000 images, with 784 pixels (predictors) and the test set compromises of 10,000 images with 784 pixels.
Some analysis require more processing time so for these we work with a sample of 10,000 rows from the training set and 1,000 rows from the test set.
We first make an analysis and preprocessing of the data.
ind <- sample(nrow(mnist$train$images), 10000)
x_train <- mnist$train$images[ind,]
y_train <- factor(mnist$train$labels[ind])
ind <- sample(nrow(mnist$test$images), 1000)
x_test <- mnist$test$images[ind,]
y_test <- factor(mnist$test$labels[ind])
colnames(x_train) <- 1:ncol(mnist$train$images)
colnames(x_test) <- colnames(x_train)
We start by looking into the distribution of the pixels variability. Using the standard deviation, we remove predictors (columns) associated with pixels that don't change much and thus can't provide much information for classification.
sds <- colSds(x_train)
qplot(sds, bins = "30", color = I("black"))
The image below shows the variance of the pixels, plotted by location. It is reasonable that higher variability happens in the central area, since rarely the digits are written on the boarders.
options(repr.plot.width=7, repr.plot.height=7)
image(1:28, 1:28, matrix(sds, 28, 28)[, 28:1])
Filtering the pixels with bigger variance, we end up with 322 predictors of the original 784, those on the center.
col_index <- colSds(x_train) > 60
new_x_train <- x_train[ ,colSds(x_train) > 60]
dim(new_x_train)
We can make a similar histogram to see the distribution of all pixel values, separating the pixels with and without ink. We can then binarize the data.
We use a cutoff to zero out low values that likely smudges, since they represent unwritten space and convert all the remaining entries into 1.
qplot(as.vector(x_train), bins = 30, color = I("black"))
bin_x_train <- (new_x_train > 255/2)*1
Below is an example between a random image before and after it was binarized.
ind <- sample(nrow(x_train), 1)
bin_x <- (x_train[ind, ] > 255/2)*1
image(1:28, 1:28, matrix(x_train[ind, ], 28, 28)[, 28:1])
image(1:28, 1:28, matrix(bin_x, 28, 28)[, 28:1])
We will apply the knn and random forests, two very common applied techniques, and study their outcomes.
We use 10-fold cross validation to tune the best k parameter tuning. Since the knn method is very suscetible to neighborhood pixels, we use the training dataset before binarizing.
control <- trainControl(method = "cv", number = 10, p = .9)
train_knn <- train(x = new_x_train, y = y_train,
method = "knn",
tuneGrid = data.frame(k = c(3,5,7,9)),
trControl = control)
train_knn$bestTune
train_knn$finalModel
| k | |
|---|---|
| <dbl> | |
| 1 | 3 |
3-nearest neighbor model Training set outcome distribution: 0 1 2 3 4 5 6 7 8 9 968 1119 987 1011 988 902 957 1086 982 1000
With the optimized k parameter, we can fit the entire dataset for the model.
fit_knn <- knn3(x = new_x_train, y = y_train, k = train_knn$bestTune$k)
y_hat_knn <- predict(fit_knn, x_test[, col_index], type="class")
confusionMatrix(data = y_hat_knn, reference = y_test)
Confusion Matrix and Statistics
Reference
Prediction 0 1 2 3 4 5 6 7 8 9
0 96 0 2 0 1 0 0 0 1 0
1 0 98 1 1 1 0 0 2 0 0
2 0 0 91 0 0 0 0 1 0 0
3 0 0 0 104 0 2 0 0 2 1
4 0 0 1 0 87 0 0 0 2 0
5 0 0 2 2 0 88 0 0 0 0
6 2 1 0 0 0 0 99 0 1 0
7 0 0 2 3 0 0 1 105 2 0
8 0 0 0 1 0 0 0 0 88 0
9 0 1 0 0 4 2 0 0 2 100
Overall Statistics
Accuracy : 0.956
95% CI : (0.9414, 0.9679)
No Information Rate : 0.111
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.9511
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
Sensitivity 0.9796 0.9800 0.9192 0.9369 0.9355 0.9565
Specificity 0.9956 0.9944 0.9989 0.9944 0.9967 0.9956
Pos Pred Value 0.9600 0.9515 0.9891 0.9541 0.9667 0.9565
Neg Pred Value 0.9978 0.9978 0.9912 0.9921 0.9934 0.9956
Prevalence 0.0980 0.1000 0.0990 0.1110 0.0930 0.0920
Detection Rate 0.0960 0.0980 0.0910 0.1040 0.0870 0.0880
Detection Prevalence 0.1000 0.1030 0.0920 0.1090 0.0900 0.0920
Balanced Accuracy 0.9876 0.9872 0.9590 0.9657 0.9661 0.9761
Class: 6 Class: 7 Class: 8 Class: 9
Sensitivity 0.9900 0.9722 0.8980 0.9901
Specificity 0.9956 0.9910 0.9989 0.9900
Pos Pred Value 0.9612 0.9292 0.9888 0.9174
Neg Pred Value 0.9989 0.9966 0.9890 0.9989
Prevalence 0.1000 0.1080 0.0980 0.1010
Detection Rate 0.0990 0.1050 0.0880 0.1000
Detection Prevalence 0.1030 0.1130 0.0890 0.1090
Balanced Accuracy 0.9928 0.9816 0.9484 0.9900
confusionMatrix(data = y_hat_knn, reference = y_test)$overall["Accuracy"]
confusionMatrix(data = y_hat_knn, reference = y_test)$byClass[,1:2]
| Sensitivity | Specificity | |
|---|---|---|
| Class: 0 | 0.9795918 | 0.9955654 |
| Class: 1 | 0.9800000 | 0.9944444 |
| Class: 2 | 0.9191919 | 0.9988901 |
| Class: 3 | 0.9369369 | 0.9943757 |
| Class: 4 | 0.9354839 | 0.9966924 |
| Class: 5 | 0.9565217 | 0.9955947 |
| Class: 6 | 0.9900000 | 0.9955556 |
| Class: 7 | 0.9722222 | 0.9910314 |
| Class: 8 | 0.8979592 | 0.9988914 |
| Class: 9 | 0.9900990 | 0.9899889 |
From the sensitivity and specificity, we can see that with this model, the 8 is the hardest to detect and the most commonly incorrect predicted digit is 9.
ggplot(train_knn, highlight = TRUE)
We show here 12 examples of correct predictions of the test set.
rafalib::mypar(3,4)
for(i in 1:12){
image(matrix(x_test[i,], 28, 28)[, 28:1],
main = paste("Our prediction:", y_hat_knn[i]),
xaxt="n", yaxt="n")
}
And below we have 12 examples of incorrect predictions, with their true classification and prediction.
ind <- which(y_hat_knn != y_test)
ind <- ind[order(y_hat_knn[ind], decreasing = TRUE)]
rafalib::mypar(3,4)
for(i in ind[1:12]){
image(matrix(x_test[i,], 28, 28)[, 28:1],
main = paste0("Pr(",y_test[i],")=",y_hat_knn[i],
" but is a ",y_test[i]),
xaxt="n", yaxt="n")
}
Since the fitting step is the slowest part, rather than predicting, as with knn, we will apply 5-fold cross validation.
We use only a small sample (nSamp) of the data first to build every tree and tune the best parameter.
control <- trainControl(method="cv", number = 5, p = 0.8)
grid <- expand.grid(minNode = c(1,5) , predFixed = c(10, 15, 25, 35, 50))
train_rf <- train(x = new_x_train, y = y_train,
method = "Rborist",
nTree = 50,
trControl = control,
tuneGrid = grid,
nSamp = 5000)
train_rf$bestTune
| predFixed | minNode | |
|---|---|---|
| <dbl> | <dbl> | |
| 2 | 15 | 1 |
Now we can use the best tuned parameter and the entire dataset to fit the model.
fit_rf <- Rborist(x = new_x_train, y = y_train,
nTree = 1000,
minNode = train_rf$bestTune$minNode,
predFixed = train_rf$bestTune$predFixed)
y_hat_rf <- factor(levels(y_train)[predict(fit_rf, x_test[, col_index], type="class")$yPred])
confusionMatrix(data = y_hat_rf, reference = y_test)
Confusion Matrix and Statistics
Reference
Prediction 0 1 2 3 4 5 6 7 8 9
0 96 0 0 0 0 0 1 0 0 0
1 0 98 0 0 0 0 0 1 0 0
2 0 0 95 2 1 0 0 4 0 0
3 0 0 1 104 0 3 0 0 0 1
4 0 0 1 0 86 1 1 0 0 1
5 0 0 0 1 0 88 1 0 1 1
6 2 1 0 0 1 0 97 0 1 0
7 0 1 2 4 0 0 0 100 1 0
8 0 0 0 0 0 0 0 0 93 0
9 0 0 0 0 5 0 0 3 2 98
Overall Statistics
Accuracy : 0.955
95% CI : (0.9402, 0.967)
No Information Rate : 0.111
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.95
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
Sensitivity 0.9796 0.9800 0.9596 0.9369 0.9247 0.9565
Specificity 0.9989 0.9989 0.9922 0.9944 0.9956 0.9956
Pos Pred Value 0.9897 0.9899 0.9314 0.9541 0.9556 0.9565
Neg Pred Value 0.9978 0.9978 0.9955 0.9921 0.9923 0.9956
Prevalence 0.0980 0.1000 0.0990 0.1110 0.0930 0.0920
Detection Rate 0.0960 0.0980 0.0950 0.1040 0.0860 0.0880
Detection Prevalence 0.0970 0.0990 0.1020 0.1090 0.0900 0.0920
Balanced Accuracy 0.9892 0.9894 0.9759 0.9657 0.9602 0.9761
Class: 6 Class: 7 Class: 8 Class: 9
Sensitivity 0.9700 0.9259 0.9490 0.9703
Specificity 0.9944 0.9910 1.0000 0.9889
Pos Pred Value 0.9510 0.9259 1.0000 0.9074
Neg Pred Value 0.9967 0.9910 0.9945 0.9966
Prevalence 0.1000 0.1080 0.0980 0.1010
Detection Rate 0.0970 0.1000 0.0930 0.0980
Detection Prevalence 0.1020 0.1080 0.0930 0.1080
Balanced Accuracy 0.9822 0.9585 0.9745 0.9796
confusionMatrix(data = y_hat_rf, reference = y_test)$overall["Accuracy"]
confusionMatrix(data = y_hat_rf, reference = y_test)$byClass[,1:2]
| Sensitivity | Specificity | |
|---|---|---|
| Class: 0 | 0.9795918 | 0.9988914 |
| Class: 1 | 0.9800000 | 0.9988889 |
| Class: 2 | 0.9595960 | 0.9922309 |
| Class: 3 | 0.9369369 | 0.9943757 |
| Class: 4 | 0.9247312 | 0.9955899 |
| Class: 5 | 0.9565217 | 0.9955947 |
| Class: 6 | 0.9700000 | 0.9944444 |
| Class: 7 | 0.9259259 | 0.9910314 |
| Class: 8 | 0.9489796 | 1.0000000 |
| Class: 9 | 0.9702970 | 0.9888765 |
From the sensitivity and specificity, we can see that with this model, the 4 is the hardest to detect and the most commonly incorrect predicted digit is 9.
ggplot(train_rf, highlight = TRUE)
Below we have 12 examples of incorrect predictions, with their true classification and prediction.
ind <- which(y_hat_rf != y_test)
ind <- ind[order(y_hat_rf[ind], decreasing = TRUE)]
rafalib::mypar(3,4)
for(i in ind[1:12]){
image(matrix(x_test[i,], 28, 28)[, 28:1],
main = paste0("Pr(",y_test[i],")=",y_hat_rf[i],
" but is a ",y_test[i]),
xaxt="n", yaxt="n")
}
In random forests we can investigate the importance of the predictors. By showing the importance of the pixels on the graphic, we can see that the important ones are located in the center, where the digits change more frequently.
rf <- randomForest(x_train, y_train, ntree = 50)
imp <- importance(rf)
image(matrix(imp, 28, 28))
Now we will combine both knn and random forest models in the ensemble, based on the highest average probability.
p_rf <- predict(fit_rf, x_test[, col_index])$census
p_rf <- p_rf / rowSums(p_rf)
p_knn <- predict(fit_knn, x_test[, col_index])
p <- (p_rf + p_knn)/2
y_pred <- factor(apply(p, 1, which.max)-1)
confusionMatrix(data = y_pred, reference = y_test)
Confusion Matrix and Statistics
Reference
Prediction 0 1 2 3 4 5 6 7 8 9
0 96 0 2 0 1 0 0 0 0 0
1 0 98 0 0 1 0 0 2 0 0
2 0 0 92 0 0 0 0 1 0 0
3 0 0 0 106 0 2 0 0 2 1
4 0 0 1 0 87 0 0 0 2 1
5 0 0 2 2 0 89 0 0 0 0
6 2 1 0 0 0 0 99 0 1 0
7 0 0 2 3 0 0 1 105 2 0
8 0 0 0 0 0 0 0 0 90 0
9 0 1 0 0 4 1 0 0 1 99
Overall Statistics
Accuracy : 0.961
95% CI : (0.9471, 0.9721)
No Information Rate : 0.111
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.9566
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
Sensitivity 0.9796 0.9800 0.9293 0.9550 0.9355 0.9674
Specificity 0.9967 0.9967 0.9989 0.9944 0.9956 0.9956
Pos Pred Value 0.9697 0.9703 0.9892 0.9550 0.9560 0.9570
Neg Pred Value 0.9978 0.9978 0.9923 0.9944 0.9934 0.9967
Prevalence 0.0980 0.1000 0.0990 0.1110 0.0930 0.0920
Detection Rate 0.0960 0.0980 0.0920 0.1060 0.0870 0.0890
Detection Prevalence 0.0990 0.1010 0.0930 0.1110 0.0910 0.0930
Balanced Accuracy 0.9881 0.9883 0.9641 0.9747 0.9655 0.9815
Class: 6 Class: 7 Class: 8 Class: 9
Sensitivity 0.9900 0.9722 0.9184 0.9802
Specificity 0.9956 0.9910 1.0000 0.9922
Pos Pred Value 0.9612 0.9292 1.0000 0.9340
Neg Pred Value 0.9989 0.9966 0.9912 0.9978
Prevalence 0.1000 0.1080 0.0980 0.1010
Detection Rate 0.0990 0.1050 0.0900 0.0990
Detection Prevalence 0.1030 0.1130 0.0900 0.1060
Balanced Accuracy 0.9928 0.9816 0.9592 0.9862
confusionMatrix(data = y_pred, reference = y_test)$overall["Accuracy"]
confusionMatrix(data = y_pred, reference = y_test)$byClass[,1:2]
| Sensitivity | Specificity | |
|---|---|---|
| Class: 0 | 0.9795918 | 0.9966741 |
| Class: 1 | 0.9800000 | 0.9966667 |
| Class: 2 | 0.9292929 | 0.9988901 |
| Class: 3 | 0.9549550 | 0.9943757 |
| Class: 4 | 0.9354839 | 0.9955899 |
| Class: 5 | 0.9673913 | 0.9955947 |
| Class: 6 | 0.9900000 | 0.9955556 |
| Class: 7 | 0.9722222 | 0.9910314 |
| Class: 8 | 0.9183673 | 1.0000000 |
| Class: 9 | 0.9801980 | 0.9922136 |
We can see that the ensemble accuracy improved from the two applied individually.
Below are 12 examples of incorrect predictions, with their true classification and prediction.
ind <- which(y_pred != y_test)
ind <- ind[order(y_pred[ind], decreasing = TRUE)]
rafalib::mypar(3,4)
for(i in ind[1:12]){
image(matrix(x_test[i,], 28, 28)[, 28:1],
main = paste0("Pr(",y_test[i],")=",y_pred[i],
" but is a ",y_test[i]),
xaxt="n", yaxt="n")
}