1 preparation

we need to install the following packages ¹:

# install needed R packages
remotes::update_packages(c('magrittr', 'dplyr', 'randomForest', 'ROCR', 'GGally'), upgrade = TRUE)

To avoid conflict of function name, in the following code, I will try me best to use pkg::fun() instead of library(pkg). But magrittr is special, we have to library() it.

library(magrittr)

Before we start, let set the random seed to make our results reproducible:

set.seed(0) 

2 generate data set

We use one of R’s built-in data set, iris, Edgar Anderson’s Iris Data set.

The original data set contains observations for four features (sepal length and width, and petal length and width — all in cm) of 150 flowers of three species (each 50).

To make things simple, here we only choose two species, versicolor and virginica. ²

df <- iris %>% tibble::as_tibble() %>% 
    dplyr::filter(Species != 'setosa') %>%
    dplyr::mutate(Species = factor(Species))
df

Remember to drops unused level using factor(), otherwise randomForest::randomForest() would complain. ³

Optional: before we build the model, we may explore the correlation between features ⁴.

3 Split data set

Before we build the model, we need to split the data set into training set and testing set. So we can train our model using data in training set, and evalute the model using data in testing set.

Here we randomly assigns 80 percent samples to the traing set, and the left 20 percent to the testing set ⁵.

nrow_training <- floor(nrow(df) * 0.8)  # Calculate the size of training sets
indexes <- sample(1:nrow(df), nrow_training)  # these rows will be select for training

training <- df[indexes, ] 
testing <- df[-indexes, ]

You may want to have a look at the result ⁶.

4 Build the model

Then we can build a random forest model ⁷.

rf_classifier = randomForest::randomForest(Species ~ ., training)

You can feed randomForest::randomForest() two arguments:

• a formula: here Species is the response variable, . tells that all other variables are features
• a data frame to train the model

5 Evaluate the model - basics

After we build the model, the fundamental use is to make prediction. That’s exactly where the testing set comes in handy:

predicted_class <- predict(rf_classifier, testing[, -ncol(testing)])

predict() needs two arguments: the model and a data frame of features ⁸.

We can compare the predicted class with real class ⁹:

real_class = testing[[ncol(testing)]]

tibble::tibble(predicted_class, real_class) %>% 
    dplyr::mutate(correct = predicted_class == real_class)

To quantify the performance of our model, we can construct a confusion matrix

# we define 'versicolor' as positive class
TP <- sum((predicted_class == 'versicolor') & (real_class == 'versicolor')) #true positive
FP <- sum((predicted_class == 'versicolor') & (real_class != 'versicolor')) #false positive
FN <- sum((predicted_class != 'versicolor') & (real_class == 'versicolor')) #false negative 
TN <- sum((predicted_class != 'versicolor') & (real_class != 'versicolor')) #true negative 

tibble::tribble(
    ~'', ~'Predicted versicolor', ~'Predicted virginica',
    'True versicolor', TP, FP,
    'True virginica',  FN, TN
)

and compute common evaluation indicators

tibble::tribble(
    ~indicator, ~value,
    'sensitivity', TP/(TP + FN),
    'specificity', TN/(TN + FP),
    'accuracy', (TP + TN)/(TP + TN + FP + FN)
)

6 Evaluate the model - advanced

Except for definite class, the model can also tell us the probability that a given sample belongs to each class:

probability <- predict(rf_classifier, testing[, -ncol(testing)], type = "prob")
probability %>% tibble::as_tibble()

For example, 0.998 in the second row means the probability of the second sample in the testing set belonging to versicolor.

Actually, the predicted class we see above is determined by judging whether the probability is greater than a given cutoff (the default is 0.5):

all((probability[ , 1, drop = T] > 0.5) == (predicted_class == 'versicolor'))
## [1] TRUE

Using different curoff, we can draw a ROC curve ¹⁰.

label <- testing[[ncol(testing)]] %>% {. ==  'versicolor'} %>% as.integer

prediction <- ROCR::prediction(probability[, 1], label)
## Error in loadNamespace(name): there is no package called 'ROCR'
prediction %>% ROCR::performance("tpr", "fpr") %>% ROCR::plot(main = "ROC Curve") 
## Error in eval(lhs, parent, parent): object 'prediction' not found

and cauculate the AUC (area under the curve)

ROCR::performance(prediction, 'auc')@y.values[[1]]
## Error in loadNamespace(name): there is no package called 'ROCR'

7 Tips and more

If you want to see Importance of each feature

randomForest::randomForest(Species ~ ., training, importance = TRUE) %>% randomForest::varImpPlot()