Project 1 - Classifying Iris Flowers

a. Discriminant Analysis

It seems fitting (excuse the pun) to start with the method that Fisher himself used, and then I will attempt to improve upon it using other machine learning algorithms. LDAs are useful for classification problems where there are more than 2 values for the predictor variable which is the case here, otherwise we could use logistic regression. Using a LDA is predicated on the normality assumption, and from the bell-curve shape of the density plots, its reasonable to proceed.

model_LDA <- lda(Species ~ ., train)
model_LDA

## Call:
## lda(Species ~ ., data = train)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa           4.9725      3.3925       1.4500      0.2375
## versicolor       5.9375      2.7675       4.2675      1.3375
## virginica        6.6450      3.0100       5.6075      2.0450
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.8776422  0.08459848
## Sepal.Width   1.3427090  2.27738505
## Petal.Length -2.2499175 -0.91449073
## Petal.Width  -2.6969658  2.67386351
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9905 0.0095

From the model outputs, I can use the coefficients of each linear discriminant to show the decision rules (to 2 decimal places): \(LD1 = 0.88 * Sepal.Length + 1.34 * Sepal.Width - 2.25 * Petal.Length - 2.70 * Petal.Width\) \(LD2 = 0.08 * Sepal.Length + 2.28 * Sepal.Width - 0.91 * Petal.Length + 2.67 * Petal.Width\)

The proportion of trace tells us how well each discriminant distinguishes between the species, and given the very high size of the first discriminant (0.9905), we see that LD1 explains 99% of the variance and so LD2 is contributing very little to the distinction of species, as can be seen from the plot below.

dataforplot <- data.frame(species = train$Species, lda = predict(model_LDA)$x)
ggscatter(dataforplot, x="lda.LD1", y="lda.LD2", shape = "species", color = "species", palette = "lancet", size = 2, alpha = 0.6)

And now on to the predictions: is this model good at predicting Iris species based on 4 measurements?

predict_LDA <- predict(model_LDA, test)
confusion_LDA <- confusionMatrix(predict_LDA$class, test$Species)
confusion_LDA

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         10          0         0
##   versicolor      0          9         0
##   virginica       0          1        10
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9667          
##                  95% CI : (0.8278, 0.9992)
##     No Information Rate : 0.3333          
##     P-Value [Acc > NIR] : 2.963e-13       
##                                           
##                   Kappa : 0.95            
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                 1.0000            0.9000           1.0000
## Specificity                 1.0000            1.0000           0.9500
## Pos Pred Value              1.0000            1.0000           0.9091
## Neg Pred Value              1.0000            0.9524           1.0000
## Prevalence                  0.3333            0.3333           0.3333
## Detection Rate              0.3333            0.3000           0.3333
## Detection Prevalence        0.3333            0.3000           0.3667
## Balanced Accuracy           1.0000            0.9500           0.9750

96.7% accuracy has been achieved with only one mistaken prediction, so this is a high benchmark, but lets see if we can better it. I’ll store this result for later.

Models_Accuracies <- tibble(model = "LDA", accuracy = confusion_LDA$overall['Accuracy'])

b. Decision Tree

I have chosen this algorithm next as the outputs are generally easily interpretable (especially when the number of predictor variables is low) and being non-parametric, they are good at handling a mixed bag of variable types.

library(rpart)
library(rpart.plot)
iristree <- rpart(formula = Species ~ ., data = train, control = rpart.control(minsplit = 5, minbucket = 2), method = "class") # method = "class" directs the function to treat Species as a categorical variable.
print(iristree)

## n= 120 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 120 80 setosa (0.33333333 0.33333333 0.33333333)  
##    2) Petal.Length< 2.45 40  0 setosa (1.00000000 0.00000000 0.00000000) *
##    3) Petal.Length>=2.45 80 40 versicolor (0.00000000 0.50000000 0.50000000)  
##      6) Petal.Length< 5.05 44  4 versicolor (0.00000000 0.90909091 0.09090909)  
##       12) Petal.Width< 1.75 39  0 versicolor (0.00000000 1.00000000 0.00000000) *
##       13) Petal.Width>=1.75 5  1 virginica (0.00000000 0.20000000 0.80000000) *
##      7) Petal.Length>=5.05 36  0 virginica (0.00000000 0.00000000 1.00000000) *

rpart.plot(iristree)

The plot itself acts as a useful set of rules to apply to further datasets to predict species, and so the allocation rules would be:

1. If Petal.Length < 2.5cm then species = Setosa

2. If Petal.Length > 5.1cm OR If Petal.Length between 2.5cm and 5.1cm and Petal.Width > 1.8cm then species = Virginica

3. If Petal.Length between 2.5cm and 5.1cm and Petal.Width < 1.8cm then species = Versicolor

Its interesting to note that this model does not accommodate the two sepal variables. I suspect that this will hamper the model’s accuracy.

Let’s look at how well the model performs:

predict_tree <- predict(iristree, test, type="class")
confusion_tree <- confusionMatrix(predict_tree, test$Species)
Models_Accuracies <- add_row(Models_Accuracies, model = "Classification Tree", accuracy = confusion_tree$overall['Accuracy'])
confusion_tree

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         10          0         0
##   versicolor      0          9         2
##   virginica       0          1         8
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9             
##                  95% CI : (0.7347, 0.9789)
##     No Information Rate : 0.3333          
##     P-Value [Acc > NIR] : 1.665e-10       
##                                           
##                   Kappa : 0.85            
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                 1.0000            0.9000           0.8000
## Specificity                 1.0000            0.9000           0.9500
## Pos Pred Value              1.0000            0.8182           0.8889
## Neg Pred Value              1.0000            0.9474           0.9048
## Prevalence                  0.3333            0.3333           0.3333
## Detection Rate              0.3333            0.3000           0.2667
## Detection Prevalence        0.3333            0.3667           0.3000
## Balanced Accuracy           1.0000            0.9000           0.8750

We have 3 flowers misclassified, and so we have an accuracy rate of 3 / 30 = 90%, therefore the model does not perform as well as the LDA. The next logical step then would be to combine multiple trees, so on to random forest.

c. Random Forests

Unlike regression models, we need to select our hyperparameters before the model is fit. But first I will run a simple baseline random forest from which we can then tune the hyperparameters to improve accuracy.

library(randomForest)

rf_baseline <- randomForest(Species ~ ., train, importance=TRUE, proximity=TRUE)
rf_baseline

## 
## Call:
##  randomForest(formula = Species ~ ., data = train, importance = TRUE,      proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 3.33%
## Confusion matrix:
##            setosa versicolor virginica class.error
## setosa         40          0         0        0.00
## versicolor      0         38         2        0.05
## virginica       0          2        38        0.05

The Out-of-bag error for this model is 3.33% which gives us an accuracy of 96.7%, and the error breaks down into 2 misclassifications of versicolor when it should be virginica, and 2 misclassifications of virginica when it should be versicolor.

varImpPlot(rf_baseline)

The importance information and plots show that the most important predictors in the model are petal length and petal width. No surprise given what we’ve seen already in our scatterplots of length v width, and also the decision tree earlier on that completely disregarded sepal length and width.

Now let’s see how the trained model performs on the test set.

predict_rf_baseline <- predict(rf_baseline, test)
confusion_RF <- confusionMatrix(predict_rf_baseline, test$Species)
Models_Accuracies <- add_row(Models_Accuracies, model = "Random Forest baseline", accuracy = confusion_RF$overall['Accuracy'])
confusion_RF

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         10          0         0
##   versicolor      0          9         2
##   virginica       0          1         8
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9             
##                  95% CI : (0.7347, 0.9789)
##     No Information Rate : 0.3333          
##     P-Value [Acc > NIR] : 1.665e-10       
##                                           
##                   Kappa : 0.85            
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                 1.0000            0.9000           0.8000
## Specificity                 1.0000            0.9000           0.9500
## Pos Pred Value              1.0000            0.8182           0.8889
## Neg Pred Value              1.0000            0.9474           0.9048
## Prevalence                  0.3333            0.3333           0.3333
## Detection Rate              0.3333            0.3000           0.2667
## Detection Prevalence        0.3333            0.3667           0.3000
## Balanced Accuracy           1.0000            0.9000           0.8750

The test accuracy here is not as strong as the out of bag accuracy and the class errors seem to be pretty similar (no difficulty in classifying the setosa flowers, the misclassification comes from the other two species).

I will now complexify the random forest a little by increasing the number of trees and variables to try at each split.

rf_2 <- randomForest(formula = Species ~ ., data = train, importance = TRUE, proximity = TRUE, ntree=1000, mtry=4) 
rf_2

## 
## Call:
##  randomForest(formula = Species ~ ., data = train, importance = TRUE,      proximity = TRUE, ntree = 1000, mtry = 4) 
##                Type of random forest: classification
##                      Number of trees: 1000
## No. of variables tried at each split: 4
## 
##         OOB estimate of  error rate: 2.5%
## Confusion matrix:
##            setosa versicolor virginica class.error
## setosa         40          0         0       0.000
## versicolor      0         38         2       0.050
## virginica       0          1        39       0.025

There is a small uptick in accuracy thanks to decreasing the number of misclassifications from 4 to 3.

varImpPlot(rf_2)

Petal length is now even more important in this model and the importance of sepal length and width has fallen even more.

predict_rf_2 <- predict(rf_2, test)
confusion_rf2 <- confusionMatrix(predict_rf_2, test$Species)
Models_Accuracies <- add_row(Models_Accuracies, model = "Random Forest complex", accuracy = confusion_rf2$overall['Accuracy'])
confusion_rf2

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         10          0         0
##   versicolor      0          9         2
##   virginica       0          1         8
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9             
##                  95% CI : (0.7347, 0.9789)
##     No Information Rate : 0.3333          
##     P-Value [Acc > NIR] : 1.665e-10       
##                                           
##                   Kappa : 0.85            
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                 1.0000            0.9000           0.8000
## Specificity                 1.0000            0.9000           0.9500
## Pos Pred Value              1.0000            0.8182           0.8889
## Neg Pred Value              1.0000            0.9474           0.9048
## Prevalence                  0.3333            0.3333           0.3333
## Detection Rate              0.3333            0.3000           0.2667
## Detection Prevalence        0.3333            0.3667           0.3000
## Balanced Accuracy           1.0000            0.9000           0.8750

Unfortunately, although the accuracy improved on the training data, it did not alter the accuracy on the test set: the model still misclassifies 3 out of 30 training set cases.

d. K-nearest Neighbour

Finally, I will use K-nearest neighbours as an example of an instance-based algorithm which is the only one used here that doesn’t learn a model, so it will be interesting to see how it performs. I will create 10 models, where k lies between 1 and 10 inclusive.

library(class)
model_knn <- list()
accuracy_knn <- numeric()
for (i in 1:10) {
 model_knn[[i]] <- knn(train[,-5], test[,-5], train$Species, k=i, prob=TRUE) 
 accuracy_knn[i] <- sum(model_knn[[i]]==test$Species)/length(test$Species)*100
}
accuracy_knn

##  [1] 90.00000 90.00000 90.00000 93.33333 93.33333 93.33333 93.33333
##  [8] 93.33333 93.33333 93.33333

The accuracy of the model improves a little when k>3 but the ceiling of the model seems to be around the 93.3% mark. I think taking k any higher than 10 for such a small dataset would lead to overfitting.

model_knn4 <- knn(train[,-5], test[,-5], train$Species, k=4, prob=TRUE)
confusion_knn <- table(test$Species, model_knn4)
confusion_knn

##             model_knn4
##              setosa versicolor virginica
##   setosa         10          0         0
##   versicolor      0          9         1
##   virginica       0          1         9

Models_Accuracies <- add_row(Models_Accuracies, model = "KNN", accuracy = accuracy_knn[4]/100)