resid <- y - with(s,sweep(u[, 1:4], 2, d[1:4], FUN="*") %*% t(v[, 1:4]))
my_image(cor(resid[index,]), zlim = c(-1,1))
title("Residual Y - (PC1+PC2+PC3+PC4)")

Warning message in cor(resid[index, ]):
“the standard deviation is zero”

By looking the first four matrix factorizations, we verify that the first 2 are responsible for most of the pixels (blue and red positive and negative values respectively and white being zero).

Knn method over PCA

Since the first 40 components are responsible for over 80% of the data, we will use these 40 components to fit some classification models, presented in previous study.
We will first fit the model.
We use 10-fold cross validation to tune the best k parameter.

p <- 40
x_train <- pca$x[,1:p]
y_train <- factor(mnist$train$labels)

control <- trainControl(method = "cv", number = 10, p = .9)
train_knn <- train(x = x_train, y = y_train,
                   method = "knn", 
                   tuneGrid = data.frame(k = c(3,5,7,9)),
                   trControl = control)
train_knn$bestTune
train_knn$finalModel

A data.frame: 1 × 1

	k
	<dbl>
1	3

3-nearest neighbor model
Training set outcome distribution:

   0    1    2    3    4    5    6    7    8    9 
5923 6742 5958 6131 5842 5421 5918 6265 5851 5949 

With the best k parameter tuned, we fit the model and transform the test set to look the result.

fit_knn <- knn3(x = x_train, y = y_train, k = train_knn$bestTune$k)

col_means <- colMeans(mnist$test$images)
x_test <- sweep(mnist$test$images, 2, col_means) %*% pca$rotation
x_test <- x_test[,1:p]
y_test <- factor(mnist$test$labels)

y_hat_knn <- predict(fit_knn, x_test, 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  974    0    5    0    1    3    4    0    4    4
         1    1 1130    0    1    0    0    4   16    0    3
         2    1    4 1008    3    0    0    0    6    3    1
         3    0    0    2  973    0    9    0    1   14    6
         4    0    0    0    1  954    1    3    2    2    7
         5    1    0    0   13    1  863    0    0    9    3
         6    2    0    2    0    4    7  945    0    4    1
         7    1    0   13    7    3    1    0  994    4    2
         8    0    0    2    8    1    3    2    0  931    6
         9    0    1    0    4   18    5    0    9    3  976

Overall Statistics
                                          
               Accuracy : 0.9748          
                 95% CI : (0.9715, 0.9778)
    No Information Rate : 0.1135          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.972           
                                          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
Sensitivity            0.9939   0.9956   0.9767   0.9634   0.9715   0.9675
Specificity            0.9977   0.9972   0.9980   0.9964   0.9982   0.9970
Pos Pred Value         0.9789   0.9784   0.9825   0.9682   0.9835   0.9697
Neg Pred Value         0.9993   0.9994   0.9973   0.9959   0.9969   0.9968
Prevalence             0.0980   0.1135   0.1032   0.1010   0.0982   0.0892
Detection Rate         0.0974   0.1130   0.1008   0.0973   0.0954   0.0863
Detection Prevalence   0.0995   0.1155   0.1026   0.1005   0.0970   0.0890
Balanced Accuracy      0.9958   0.9964   0.9874   0.9799   0.9849   0.9823
                     Class: 6 Class: 7 Class: 8 Class: 9
Sensitivity            0.9864   0.9669   0.9559   0.9673
Specificity            0.9978   0.9965   0.9976   0.9956
Pos Pred Value         0.9793   0.9698   0.9769   0.9606
Neg Pred Value         0.9986   0.9962   0.9952   0.9963
Prevalence             0.0958   0.1028   0.0974   0.1009
Detection Rate         0.0945   0.0994   0.0931   0.0976
Detection Prevalence   0.0965   0.1025   0.0953   0.1016
Balanced Accuracy      0.9921   0.9817   0.9767   0.9814

confusionMatrix(data = y_hat_knn, reference = y_test)$byClass[,1:2]

acc_results <- tibble(method = "Knn entire set", acc = 0.956)

knn_pca_acc <- confusionMatrix(data = y_hat_knn, reference = y_test)$overall["Accuracy"]
acc_results <- bind_rows(acc_results,
                         tibble(method="Knn with 40 PCs", acc = knn_pca_acc))
acc_results %>% knitr::kable()

A matrix: 10 × 2 of type dbl

	Sensitivity	Specificity
Class: 0	0.9938776	0.9976718
Class: 1	0.9955947	0.9971799
Class: 2	0.9767442	0.9979929
Class: 3	0.9633663	0.9964405
Class: 4	0.9714868	0.9982258
Class: 5	0.9674888	0.9970356
Class: 6	0.9864301	0.9977881
Class: 7	0.9669261	0.9965448
Class: 8	0.9558522	0.9975626
Class: 9	0.9672944	0.9955511

|method          |    acc|
|:---------------|------:|
|Knn entire set  | 0.9560|
|Knn with 40 PCs | 0.9748|

Below we show 12 examples of correct predictions of the test set.

options(repr.plot.width=7, repr.plot.height=7)
rafalib::mypar(3,4)
for(i in 1:12){
  image(matrix(mnist$test$images[i,], 28, 28)[, 28:1], 
        main = paste("Our prediction:", y_hat_knn[i]),
        xaxt="n", yaxt="n")
}

Following are 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(mnist$test$images[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")
}

The accuracy achieved with knn over 40 principal components was very high, above 95% and in our case, was better than the accuracy achieved in previous knn model, with the entire dataset.

Random Forest Rborist over PCA

We use 5-fold cross validation and a small sample to tune the best parameter.

p <- 40
x_train <- pca$x[,1:p]
y_train <- factor(mnist$train$labels)

control <- trainControl(method="cv", number = 5, p = 0.8)
train_rf <-  train(x = x_train, y = y_train,
                   method = "Rborist",
                   nTree = 50,
                   trControl = control,
                   nSamp = 5000)
train_rf$bestTune

A data.frame: 1 × 2

	predFixed	minNode
	<dbl>	<dbl>
1	2	2

We apply the best parameter tuned and fit the model with the entire training set to look the result.

fit_rf <- Rborist(x = x_train, y = y_train,
                 nTree = 1000,
                 minNode = train_rf$bestTune$minNode,
                 predFixed = train_rf$bestTune$predFixed)

col_means <- colMeans(mnist$test$images)
x_test <- sweep(mnist$test$images, 2, col_means) %*% pca$rotation
x_test <- x_test[,1:p]
y_test <- factor(mnist$test$labels)

y_hat_rf <- factor(levels(y_train)[predict(fit_rf, x_test, 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  967    0    7    2    1    4    7    1    5    7
         1    0 1120    0    1    2    1    3    9    0    7
         2    2    4  979    4    6    1    1   15    6    1
         3    0    3   10  968    0   17    0    0   17   11
         4    0    0    1    1  944    2    3    6   10   22
         5    1    1    0   10    0  854    4    0   14    4
         6    6    4    2    1    6    6  939    0    4    1
         7    1    1   10    8    1    2    0  975    4    8
         8    3    2   22   13    2    3    1    2  908   11
         9    0    0    1    2   20    2    0   20    6  937

Overall Statistics
                                         
               Accuracy : 0.9591         
                 95% CI : (0.955, 0.9629)
    No Information Rate : 0.1135         
    P-Value [Acc > NIR] : < 2.2e-16      
                                         
                  Kappa : 0.9545         
                                         
 Mcnemar's Test P-Value : NA             

Statistics by Class:

                     Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
Sensitivity            0.9867   0.9868   0.9486   0.9584   0.9613   0.9574
Specificity            0.9962   0.9974   0.9955   0.9935   0.9950   0.9963
Pos Pred Value         0.9660   0.9799   0.9607   0.9435   0.9545   0.9617
Neg Pred Value         0.9986   0.9983   0.9941   0.9953   0.9958   0.9958
Prevalence             0.0980   0.1135   0.1032   0.1010   0.0982   0.0892
Detection Rate         0.0967   0.1120   0.0979   0.0968   0.0944   0.0854
Detection Prevalence   0.1001   0.1143   0.1019   0.1026   0.0989   0.0888
Balanced Accuracy      0.9915   0.9921   0.9721   0.9760   0.9782   0.9768
                     Class: 6 Class: 7 Class: 8 Class: 9
Sensitivity            0.9802   0.9484   0.9322   0.9286
Specificity            0.9967   0.9961   0.9935   0.9943
Pos Pred Value         0.9690   0.9653   0.9390   0.9484
Neg Pred Value         0.9979   0.9941   0.9927   0.9920
Prevalence             0.0958   0.1028   0.0974   0.1009
Detection Rate         0.0939   0.0975   0.0908   0.0937
Detection Prevalence   0.0969   0.1010   0.0967   0.0988
Balanced Accuracy      0.9884   0.9723   0.9629   0.9615

confusionMatrix(data = y_hat_rf, reference = y_test)$byClass[,1:2]

acc_results <- bind_rows(acc_results,
                         tibble(method="Rborist entire set", acc = 0.955))
rf_pca_acc <- confusionMatrix(data = y_hat_rf, reference = y_test)$overall["Accuracy"]
acc_results <- bind_rows(acc_results,
                         tibble(method="Rborist with 40 PCs", acc = rf_pca_acc))
acc_results %>% knitr::kable()

A matrix: 10 × 2 of type dbl

	Sensitivity	Specificity
Class: 0	0.9867347	0.9962306
Class: 1	0.9867841	0.9974055
Class: 2	0.9486434	0.9955397
Class: 3	0.9584158	0.9935484
Class: 4	0.9613035	0.9950100
Class: 5	0.9573991	0.9962670
Class: 6	0.9801670	0.9966821
Class: 7	0.9484436	0.9960990
Class: 8	0.9322382	0.9934633
Class: 9	0.9286422	0.9943277

|method              |    acc|
|:-------------------|------:|
|Knn entire set      | 0.9560|
|Knn with 40 PCs     | 0.9748|
|Rborist entire set  | 0.9550|
|Rborist with 40 PCs | 0.9591|

Following we show 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(mnist$test$images[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")
}

Random forest model with only 40 principal components also achieved a high accuracy, similar to the model with entire set. Both got an accuracy above 95%.