Random forest classification can be used to assess the extent of discrimination (difference) between classes of a discrete response variable. This includes both multinomial (number of classes > 2) and binary (number of classes = 2) comparisons.
In multinomial situations, the suitability of a multinomial comparison versus multiple binary comparisons can depend on the experimental context. For instance, in a treatment/control experiment that includes multiple time points, a multinomial comparison using all available classes could be useful to visualise the general structure of the data. However, it could make any extracted explanatory features difficult to reason about as to how they relate to the individual experimental time point or treatment conditions. An investigator could instead identify the binary comparisons relevant to the biological question and focus the further classification comparisons to better select for explanatory features.
-
-Multinomial comparisons
-
In experiments with more than two classes, multinomial random forest classification can be used to assess the discrimination between the classes and give an overview of the relative structure between classes.
-
The example data set consists of a total of 6 classes for the day response variable.
-
-d %>%
- clsExtract(cls = 'day') %>%
- unique()
-#> [1] 2 3 4 1 H 5
-#> Levels: 1 2 3 4 5 H
-
Multinomial classification can be performed by:
-
-multinomial_rf <- d %>%
- randomForest(cls = 'day')
-
-print(multinomial_rf)
-#>
-#> Random forest classification
-#>
-#> Samples: 120
-#> Features: 500
-#> Response: day
-#> # comparisons: 1
-
The performance of this model can be assessed using metrics based on the success of the out of bag (OOB) predictions. The performance metrics can be extracted using:
-
-multinomial_rf %>%
- metrics()
-#> # A tibble: 4 × 5
-#> Response Comparison .metric .estimator .estimate
-#> <chr> <chr> <chr> <chr> <dbl>
-#> 1 day 1~2~3~4~5~H accuracy multiclass 0.8
-#> 2 day 1~2~3~4~5~H kap multiclass 0.76
-#> 3 day 1~2~3~4~5~H roc_auc hand_till 0.964
-#> 4 day 1~2~3~4~5~H margin <NA> 0.146
-
These metrics include accuracy, Cohen’s kappa (kap), area under the receiver operator characteristic curve (roc_auc, ROC-AUC) and margin. Each metric has both strengths and weaknesses that depend on the context of the classification such as the balance of observations between the classes. As shown below, the class frequencies for this example are balanced with 20 observations per class.
-
-d %>%
- clsExtract(cls = 'day') %>%
- table()
-#> .
-#> 1 2 3 4 5 H
-#> 20 20 20 20 20 20
-
In this context, each of these metrics could be used to assess the predictive performance of the model. The margin metric is the difference between the proportion of votes for the correct class and the maximum proportion of votes for the other classes for a given observation which is then averaged across all the observations. A positive margin value indicates correct classification and values greater than 0.2 can be considered as the models having strong predictive power. The margin also allows the extent of discrimination to be discerned even in very distinct cases above where both the accuracy and ROC-AUC would be registering values of 1.
-
In this example, the values of all the metrics suggest that the model is showing good predictive performance. This can be investigated further by plotting the MDS of observation proximity values.
-
-multinomial_rf %>%
- plotMDS(cls = 'day')
-
-
This shows that the model is able to discriminate highly between classes such as 5 and H. It is less able to discriminate more similar classes such as H and 1 or 4 and 5 whose confidence ellipses show a high degree of overlap. This makes sense in the context of this experiment as these are adjacent time points that are more likely to be similar than time points at each end of the experiment.
-
The ROC curves can also be plotted as shown below.
-
-
-
Classes with their line further from the central dashed line are those that were predicted with the greatest reliability by the model. This plot shows that both the H and 1 classes were least reliably predicted which is a result of their close proximity shown in the MDS plot previously.
-
Importance metrics can be used to identify the metabolome features that contribute most to the class discrimination in the model. The available importance metrics for this model are shown below.
-
-importanceMetrics(multinomial_rf)
-#> [1] "1" "2" "3"
-#> [4] "4" "5" "FalsePositiveRate"
-#> [7] "H" "MeanDecreaseAccuracy" "MeanDecreaseGini"
-#> [10] "SelectionFrequency"
-
Here, we will use the false positive rate metric with a threshold of below 0.05 to identify explanatory features for the day response variable.
-
-multinomial_rf %>%
- explanatoryFeatures(metric = 'FalsePositiveRate',
- threshold = 0.05)
-#> # A tibble: 121 × 5
-#> Response Comparison Feature Metric Value
-#> <chr> <chr> <chr> <chr> <dbl>
-#> 1 day 1~2~3~4~5~H N341 FalsePositiveRate 1.02e-93
-#> 2 day 1~2~3~4~5~H N133 FalsePositiveRate 7.38e-68
-#> 3 day 1~2~3~4~5~H N163 FalsePositiveRate 3.59e-61
-#> 4 day 1~2~3~4~5~H N439 FalsePositiveRate 1.07e-54
-#> 5 day 1~2~3~4~5~H N342 FalsePositiveRate 3.19e-49
-#> 6 day 1~2~3~4~5~H N377 FalsePositiveRate 3.19e-49
-#> 7 day 1~2~3~4~5~H N171 FalsePositiveRate 6.26e-44
-#> 8 day 1~2~3~4~5~H N497 FalsePositiveRate 6.11e-30
-#> 9 day 1~2~3~4~5~H N146 FalsePositiveRate 2.74e-29
-#> 10 day 1~2~3~4~5~H N195 FalsePositiveRate 7.16e-25
-#> # … with 111 more rows
-
As shown above there were a total of 121 explanatory features identified.
-
Within a multinomial experiment, it is also possible to specify the exact class comparisons to include, where it might not be suitable to compare all the classes at once using the comparisons argument. This should be specified as a named list, the corresponding to the cls argument. Each named element should then consist of a vector of comparisons, the classes to compare separated using the ~.
-
The following specifies two comparisons (H~1~2,H~1~5) for the day response variable and displays the performance metrics.
-
-d %>%
- randomForest(cls = 'day',
- comparisons = list(day = c('H~1~2',
- 'H~1~5'))) %>%
- metrics()
-#> # A tibble: 8 × 5
-#> Response Comparison .metric .estimator .estimate
-#> <chr> <chr> <chr> <chr> <dbl>
-#> 1 day H~1~2 accuracy multiclass 0.833
-#> 2 day H~1~2 kap multiclass 0.75
-#> 3 day H~1~5 accuracy multiclass 0.75
-#> 4 day H~1~5 kap multiclass 0.625
-#> 5 day H~1~2 roc_auc hand_till 0.906
-#> 6 day H~1~5 roc_auc hand_till 0.909
-#> 7 day H~1~2 margin <NA> 0.172
-#> 8 day H~1~5 margin <NA> 0.320
-
The MDS and ROC curve plots can also be plotted simultaneously for the two comparisons.
-
-
-
Similarly, it is also possible to model multiple response factors with a single random forest call by specifying a vector of response class information column names to the cls argument. In the following, both the name and day response factors will be analysed and the performance metrics returned in a single table.
-
-d %>%
- randomForest(cls = c('name','day')) %>%
- metrics()
-#> Warning: Classes with < 5 replicates removed: "11_3", "11_4", "11_5", "11_6",
-#> "11_H", "12_1", "12_3", "12_6", "12_H", "13_1", "13_2", "13_3", "13_5", "13_6",
-#> "13_H", "14_2", "14_3", "14_5", "14_6", "14_H", "15_1", "15_2", "15_4", "15_5",
-#> "15_6", "15_H"
-#> Unbalanced classes detected. Stratifying sample size to the smallest class size.
-#> # A tibble: 8 × 5
-#> Response Comparison .metric .estimator .estimate
-#> <chr> <chr> <chr> <chr> <dbl>
-#> 1 name 11_2~12_2~12_4~13_4~14_4~15_3 accuracy multiclass 0.35
-#> 2 name 11_2~12_2~12_4~13_4~14_4~15_3 kap multiclass 0.212
-#> 3 name 11_2~12_2~12_4~13_4~14_4~15_3 roc_auc hand_till 0.753
-#> 4 name 11_2~12_2~12_4~13_4~14_4~15_3 margin <NA> -0.0485
-#> 5 day 1~2~3~4~5~H accuracy multiclass 0.8
-#> 6 day 1~2~3~4~5~H kap multiclass 0.76
-#> 7 day 1~2~3~4~5~H roc_auc hand_till 0.964
-#> 8 day 1~2~3~4~5~H margin <NA> 0.146
-
The MDS plots can also be returned for both models simultaneously.
-
-d %>%
- randomForest(cls = c('name','day')) %>%
- plotMDS()
-#> Warning: Classes with < 5 replicates removed: "11_3", "11_4", "11_5", "11_6",
-#> "11_H", "12_1", "12_3", "12_6", "12_H", "13_1", "13_2", "13_3", "13_5", "13_6",
-#> "13_H", "14_2", "14_3", "14_5", "14_6", "14_H", "15_1", "15_2", "15_4", "15_5",
-#> "15_6", "15_H"
-#> Unbalanced classes detected. Stratifying sample size to the smallest class size.
-
-
-
-Binary comparisons
-
It may in some cases be preferable to analyse class comparisons as multiple binary comparisons.
-
The possible binary comparisons for a given response variable can be displayed using the binaryComparisons() method. Below shows the 15 comparisons for the day response variable.
-
-binaryComparisons(d,cls = 'day')
-#> [1] "1~2" "1~3" "1~4" "1~5" "1~H" "2~3" "2~4" "2~5" "2~H" "3~4" "3~5" "3~H"
-#> [13] "4~5" "4~H" "5~H"
-
For this example we will only use the binary comparisons containing the H class.
-
-
The binary comparisons can then be performed using the following.
-
-binary_rf <- d %>%
- randomForest(cls = 'day',
- comparisons = list(day = binary_comparisons))
-
-print(binary_rf)
-#>
-#> Random forest classification
-#>
-#> Samples: 120
-#> Features: 500
-#> Response: day
-#> # comparisons: 5
-
To run all possible binary comparisons, the binary = TRUE argument could instead be used.
-
The MDS plots for each comparison can be visualised to inspect the comparisons.
-
-binary_rf %>%
- plotMDS(cls = 'day')
-
-
These plots show good separation in all the comparisons except H~1 which is also shown by the plot of the performance metrics below. Each of the comparisons are showing perfect performance for the accuracy, Cohen’s kappa and ROC-AUC metrics as well as very high margin values except for the H~1 comparison.
-
-
-
The explanatory features for these comparisons can be extracted as below using the false positive rate metric and a cut-off threshold of 0.05. This gives a total of 251 explanatory features.
-
-binary_rf %>%
- explanatoryFeatures(metric = 'FalsePositiveRate',
- threshold = 0.05)
-#> # A tibble: 251 × 5
-#> Response Comparison Feature Metric Value
-#> <chr> <chr> <chr> <chr> <dbl>
-#> 1 day 2~H N341 FalsePositiveRate 7.34e-52
-#> 2 day 2~H N439 FalsePositiveRate 1.80e-45
-#> 3 day 3~H N342 FalsePositiveRate 2.71e-39
-#> 4 day 2~H N327 FalsePositiveRate 1.06e-35
-#> 5 day 3~H N439 FalsePositiveRate 1.06e-35
-#> 6 day 2~H N477 FalsePositiveRate 1.60e-34
-#> 7 day 3~H N377 FalsePositiveRate 1.60e-34
-#> 8 day 4~H N477 FalsePositiveRate 7.40e-34
-#> 9 day 2~H N447 FalsePositiveRate 6.48e-30
-#> 10 day 3~H N163 FalsePositiveRate 6.48e-30
-#> # … with 241 more rows
-
A heatmap of these explanatory features can be plotted to show their mean relative intensities across the experiment time points. Here, the classes are also refactored to customise the order of the classes on the x-axis.
-
-refactor_cls <- clsExtract(binary_rf,
- cls = 'day') %>%
- factor(.,levels = c('H','1','2','3','4','5'))
-
-binary_rf <- clsReplace(binary_rf,
- value = refactor_cls,
- cls = 'day')
-binary_rf %>%
- plotExplanatoryHeatmap(metric = 'FalsePositiveRate',
- threshold = 0.05,
- featureNames = TRUE)
-
-
