+
+
diff --git a/docs/articles/modelling.html b/docs/articles/modelling.html
index c1754606..917f5f6f 100644
--- a/docs/articles/modelling.html
+++ b/docs/articles/modelling.html
@@ -19,6 +19,8 @@
+
+
@@ -31,18 +33,12 @@
-
-
+
@@ -85,35 +81,35 @@
-
+
Modelling and feature selection
- metabolyseR v0.14.6
- Jasen Finch
+ metabolyseR v0.14.7
+ Jasen Finch
- 17 November, 2021
+ 17 December, 2021
- Source: vignettes/modelling.Rmd
+ Source: vignettes/modelling.Rmd
modelling.Rmd
-
-
-Introduction
+
+Introduction
+
Modelling provides the essential data mining step for extracting biological information and explanatory metabolome features from a data set relating to the experimental conditions. metabolyseR provides a number of both univariate and multivariate methods for data mining.
For an introduction to the usage of metabolyseR for both exploratory and routine analyses, see the introduction vignette using:
-vignette('introduction','metabolyseR')
+vignette('introduction','metabolyseR')
To further supplement this document, a quick start example analysis is also available as a vignette:
-vignette('quick_start','metabolyseR')
+vignette('quick_start','metabolyseR')
To begin, the package can be loaded using:
-library(metabolyseR)
+library(metabolyseR)
#>
#> Attaching package: 'metabolyseR'
#> The following object is masked from 'package:stats':
@@ -122,17 +118,17 @@
#> The following objects are masked from 'package:base':
#>
#> raw, split
-
-
-Example data
-The examples used here will use the abr1 data set from the metaboData package. This is nominal mass flow-injection mass spectrometry (FI-MS) fingerprinting data from a plant-pathogen infection time course experiment. The pipe %>% from the magrittr package will also be used. The example data can be loaded using:
+
+Example data
+
+The examples used here will use the abr1 data set from the metaboData package. This is nominal mass flow-injection mass spectrometry (FI-MS) fingerprinting data from a plant-pathogen infection time course experiment. The pipe %>% from the magrittr package will also be used. The example data can be loaded using:
+library(metaboData)
Only the negative acquisition mode data (abr1$neg) will be used along with the sample meta-information (abr1$fact). Create an AnalysisData class object, assigned to the variable d, using the following:
d <- analysisData(abr1$neg[,1:500],abr1$fact)
As can be seen above the data set contains a total of 120 samples and 500 features.
-
-
-Parallel processing
-The package supports parallel processing using the future package.
-By default, processing by metabolyseR will be done seqentially. However, parallel processing can be activated, prior to analysis, by specifying a parallel implementation using plan(). The following example specifies using the multisession implementation (muliple background R sessions) with two worker processes.
+
+Parallel processing
+
+The package supports parallel processing using the future package.
+By default, processing by metabolyseR will be done seqentially. However, parallel processing can be activated, prior to analysis, by specifying a parallel implementation using plan(). The following example specifies using the multisession implementation (muliple background R sessions) with two worker processes.
-plan(future::multisession,workers = 2)
-See the future package documentation for more information on the types of parallel implementations that are available.
+plan(future::multisession,workers = 2)
+See the future package documentation for more information on the types of parallel implementations that are available.
-
-
-Random Forest
+
+Random Forest
+
Random forest is a versatile ensemble machine learning approach based on forests of decision trees for multivariate data mining. This can include unsupervised analysis, classification of discrete response variables and regression of continuous responses.
-Random forest can be performed in metabolyseR using the randomForest() method. For further details on the arguments for using this function, see ?randomForest. This implementation of random forest in metabolyseR utilises the randomForest package. See ?randomForest::randomForest for more information about that implementation.
-
-
-Unsupervised
+Random forest can be performed in metabolyseR using the randomForest() method. For further details on the arguments for using this function, see ?randomForest. This implementation of random forest in metabolyseR utilises the randomForest package. See ?randomForest::randomForest for more information about that implementation.
+
+Unsupervised
+
The unsupervised random forest approach can be useful starting point for analysis in any experimental context. It can be used to give a general overview of the structure of the data and to identify any possible problems. These could include situations such as the presence of outliers samples or splits in the data caused by the impact of analytical or sample preparation factors. Unsupervised random forest can have advantages in these assessments over other approaches such as Principle Component Analysis (PCA). It is less sensitive to the effect of a single feature that in fact could have little overall impact relative to the other hundreds that could be present in a data set.
The examples below will show the use of unsupervised random forest for assessing the general structure of the example data set and the presence of outlier samples.
Unsupervised random forest can be performed by setting the cls argument of randomForest() to NULL:
-unsupervised_rf <- d %>%
+unsupervised_rf <- d %>%
randomForest(cls = NULL)
The type of random forest that has been performed can be checked using the type method.
Or by printing the results object.
@@ -196,75 +192,75 @@
This shows that it is indeed the experimental factor of interest that is having the greatest impact on the structure of the data. The progression of the experimental time points are obvious across Dimension 1.
The available feature importance metrics for a random forest analysis can be retrieved by:
-importanceMetrics(unsupervised_rf)
+importanceMetrics(unsupervised_rf)
#> [1] "1" "2" "FalsePositiveRate"
#> [4] "MeanDecreaseAccuracy" "MeanDecreaseGini" "SelectionFrequency"
And the importance values of these metrics for each feature can returned using:
-importance(unsupervised_rf)
-#> # A tibble: 3,000 × 3
+importance(unsupervised_rf)
+#> # A tibble: 3,000 × 3
#> Feature Metric Value
-#> <chr> <chr> <dbl>
-#> 1 N1 1 0
-#> 2 N1 2 0
-#> 3 N1 FalsePositiveRate 0.0238
-#> 4 N1 MeanDecreaseAccuracy 0
-#> 5 N1 MeanDecreaseGini 0
-#> 6 N1 SelectionFrequency 0
-#> 7 N10 1 0
-#> 8 N10 2 0
-#> 9 N10 FalsePositiveRate 0.0238
-#> 10 N10 MeanDecreaseAccuracy 0
-#> # … with 2,990 more rows
+#> <chr> <chr> <dbl>
+#> 1 N1 1 0
+#> 2 N1 2 0
+#> 3 N1 FalsePositiveRate 0.0238
+#> 4 N1 MeanDecreaseAccuracy 0
+#> 5 N1 MeanDecreaseGini 0
+#> 6 N1 SelectionFrequency 0
+#> 7 N10 1 0
+#> 8 N10 2 0
+#> 9 N10 FalsePositiveRate 0.0238
+#> 10 N10 MeanDecreaseAccuracy 0
+#> # … with 2,990 more rows
The explanatory features for a given threshold can be extracted for any of the importance metrics. The following will extract the explanatory features below a threshold of 0.05 based on the false positive rate metric.
-unsupervised_rf %>%
- explanatoryFeatures(metric = "FalsePositiveRate",
+unsupervised_rf %>%
+ explanatoryFeatures(metric = "FalsePositiveRate",
threshold = 0.05)
-#> # A tibble: 359 × 3
+#> # A tibble: 359 × 3
#> Feature Metric Value
-#> <chr> <chr> <dbl>
-#> 1 N342 FalsePositiveRate 1.31e-19
-#> 2 N161 FalsePositiveRate 2.34e-16
-#> 3 N341 FalsePositiveRate 6.50e-16
-#> 4 N315 FalsePositiveRate 1.79e-15
-#> 5 N367 FalsePositiveRate 3.47e-14
-#> 6 N173 FalsePositiveRate 9.09e-14
-#> 7 N385 FalsePositiveRate 9.09e-14
-#> 8 N133 FalsePositiveRate 1.52e-12
-#> 9 N439 FalsePositiveRate 1.52e-12
-#> 10 N379 FalsePositiveRate 3.78e-12
-#> # … with 349 more rows
+#> <chr> <chr> <dbl>
+#> 1 N342 FalsePositiveRate 1.31e-19
+#> 2 N161 FalsePositiveRate 2.34e-16
+#> 3 N341 FalsePositiveRate 6.50e-16
+#> 4 N315 FalsePositiveRate 1.79e-15
+#> 5 N367 FalsePositiveRate 3.47e-14
+#> 6 N173 FalsePositiveRate 9.09e-14
+#> 7 N385 FalsePositiveRate 9.09e-14
+#> 8 N133 FalsePositiveRate 1.52e-12
+#> 9 N439 FalsePositiveRate 1.52e-12
+#> 10 N379 FalsePositiveRate 3.78e-12
+#> # … with 349 more rows
In this example there are 359 explanatory features.
The trend of the most highly ranked explanatory feature against the day factor can be plotted using the plotFeature() method.
-unsupervised_rf %>%
+unsupervised_rf %>%
plotFeature(feature = 'N425',
cls = 'day')
-
-
-Classification
+
+Classification
+
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
+
+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()
+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 %>%
+multinomial_rf <- d %>%
randomForest(cls = 'day')
-print(multinomial_rf)
+print(multinomial_rf)
#>
#> Random forest classification
#>
@@ -274,117 +270,117 @@
#> # 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
+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
+#> <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()
+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.
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)
+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',
+multinomial_rf %>%
+ explanatoryFeatures(metric = 'FalsePositiveRate',
threshold = 0.05)
-#> # A tibble: 121 × 5
+#> # 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
+#> <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 %>%
+d %>%
randomForest(cls = 'day',
- comparisons = list(day = c('H~1~2',
- 'H~1~5'))) %>%
- metrics()
-#> # A tibble: 8 × 5
+ 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
+#> <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.
-d %>%
+d %>%
randomForest(cls = 'day',
- comparisons = list(day = c('H~1~2',
- 'H~1~5'))) %>%
+ comparisons = list(day = c('H~1~2',
+ 'H~1~5'))) %>%
{plotMDS(.,cls = 'day') +
plotROC(.) +
- patchwork::plot_layout(ncol = 1)}
+ patchwork::plot_layout(ncol = 1)}
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()
+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
+#> # 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
+#> <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')) %>%
+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",
@@ -393,26 +389,26 @@
#> Unbalanced classes detected. Stratifying sample size to the smallest class size.
-
-
-Binary comparisons
+
+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.
+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')
+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.
-binary_comparisons <- binaryComparisons(d,cls = 'day') %>%
- .[stringr::str_detect(.,'H')]
+binary_comparisons <- binaryComparisons(d,cls = 'day') %>%
+ .[stringr::str_detect(.,'H')]
The binary comparisons can then be performed using the following.
-binary_rf <- d %>%
+binary_rf <- d %>%
randomForest(cls = 'day',
- comparisons = list(day = binary_comparisons))
+ comparisons = list(day = binary_comparisons))
-print(binary_rf)
+print(binary_rf)
#>
#> Random forest classification
#>
@@ -423,59 +419,59 @@
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.
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.
-binary_rf %>%
+binary_rf %>%
plotMetrics()
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',
+binary_rf %>%
+ explanatoryFeatures(metric = 'FalsePositiveRate',
threshold = 0.05)
-#> # A tibble: 251 × 5
+#> # 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
+#> <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'))
+refactor_cls <- clsExtract(binary_rf,
+ cls = 'day') %>%
+ factor(.,levels = c('H','1','2','3','4','5'))
-binary_rf <- clsReplace(binary_rf,
+binary_rf <- clsReplace(binary_rf,
value = refactor_cls,
cls = 'day')
-binary_rf %>%
+binary_rf %>%
plotExplanatoryHeatmap(metric = 'FalsePositiveRate',
threshold = 0.05,
featureNames = TRUE)
-
-
-Regression
+
+Regression
+
Random forest regression can be used to assess the extent of association of the metabolomic data with continuous response variables.
In this example, the extent of association of injection order with the example data will be assessed.
-regression_rf <- d %>%
+regression_rf <- d %>%
randomForest(cls = 'injorder')
-print(regression_rf)
+print(regression_rf)
#>
#> Random forest regression
#>
@@ -484,21 +480,21 @@
#> Response: injorder
The regression model performance metrics, based on the OOB prediction error, can be extracted using the following:
-regression_rf %>%
- metrics()
-#> # A tibble: 5 × 4
+regression_rf %>%
+ metrics()
+#> # A tibble: 5 × 4
#> Response .metric .estimator .estimate
-#> <chr> <chr> <chr> <dbl>
-#> 1 injorder rsq standard 0.476
-#> 2 injorder mae standard 23.5
-#> 3 injorder mape standard 154.
-#> 4 injorder rmse standard 26.5
-#> 5 injorder ccc standard 0.508
+#> <chr> <chr> <chr> <dbl>
+#> 1 injorder rsq standard 0.476
+#> 2 injorder mae standard 23.5
+#> 3 injorder mape standard 154.
+#> 4 injorder rmse standard 26.5
+#> 5 injorder ccc standard 0.508
These regression metrics include R2 (rsq), mean absolute error (mae), mean absolute percentage error (mape), root mean squared error (rmse) and the concordance correlation coefficient (ccc).
The R2 and concordance correlation coefficient metrics suggest that there is some association of features with the injection order, although this is weak. This is in agreement with mean absolute error metric that shows that on average, the injection order could only be predicted to an accuracy of 23 injection order positions.
The MDS plot belows the relative proximities of the samples based on this injection order regression model. This shows that for the most part, there is little correspondence of the sample positions with their injection order. However, there is a small grouping of samples towards the end of the run around sample ~99 to 120. It suggests that there could have been some analytical issues, for certain features, towards the end of the mass spectral analytical run.
-regression_rf %>%
+regression_rf %>%
plotMDS(cls = NULL,
ellipses = FALSE,
label = 'injorder',
@@ -508,51 +504,51 @@
The available feature importance metrics for this regression model can be listed.
-regression_rf %>%
- importanceMetrics()
+regression_rf %>%
+ importanceMetrics()
#> [1] "%IncMSE" "IncNodePurity"
The feature importance metrics can be plotted to give an overview of their distribution. The following will plot the percentage increase in the mean squared error (%IncMSE) importance metric.
-regression_rf %>%
+regression_rf %>%
plotImportance(metric = "%IncMSE",
rank = FALSE)
This shows that there are only a few features that are contributing to the association with injection order. These explanatory features can be extracted with the following, using a threshold of above 5.
-regression_rf %>%
- explanatoryFeatures(metric = '%IncMSE',
+regression_rf %>%
+ explanatoryFeatures(metric = '%IncMSE',
threshold = 5)
-#> # A tibble: 7 × 4
+#> # A tibble: 7 × 4
#> Response Feature Metric Value
-#> <chr> <chr> <chr> <dbl>
-#> 1 injorder N283 %IncMSE 19.9
-#> 2 injorder N135 %IncMSE 8.71
-#> 3 injorder N451 %IncMSE 5.58
-#> 4 injorder N161 %IncMSE 5.51
-#> 5 injorder N306 %IncMSE 5.49
-#> 6 injorder N118 %IncMSE 5.22
-#> 7 injorder N297 %IncMSE 5.07
+#> <chr> <chr> <chr> <dbl>
+#> 1 injorder N283 %IncMSE 19.9
+#> 2 injorder N135 %IncMSE 8.71
+#> 3 injorder N451 %IncMSE 5.58
+#> 4 injorder N161 %IncMSE 5.51
+#> 5 injorder N306 %IncMSE 5.49
+#> 6 injorder N118 %IncMSE 5.22
+#> 7 injorder N297 %IncMSE 5.07
This returned a total of 7 explanatory features above this threshold. The top ranked feature N283 can be plotted to investigate it’s trend in relation to injection order.
-regression_rf %>%
+regression_rf %>%
plotFeature(feature = 'N283',
cls = 'injorder')
This shows an increase in the intensity of that feature for samples above 100 in the injection order which corresponds with the cluster that was seen in the MDS plot above.
-
-
-Univariate analyses
+
+Univariate analyses
+
Univariate methods select features, explanatory for response variables, with features tested on an individual basis. These methods offer simplicity and easy interpretation in their use, however they provide no information as to how features may interact.
The univariate methods currently available in metabolyseR include Welch’s t-test, analysis of variance (ANOVA) and linear regression. The following sections will provide brief examples of the use of each of these methods.
-
-
-Welch’s t-test
+
+Welch’s t-test
+
Welch’s t-test can be used to select explanatory metabolome features for binary comparisons of discrete variables. By default, all the possible binary comparisons for the categories of a response variable will be tested.
Below shows the possible binary comparisons for the day response variable for the example data set.
-binaryComparisons(d,
+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"
@@ -560,11 +556,11 @@
The explanatory features that show a significant difference between the response categories can be extracted as shown below.
-explanatoryFeatures(ttest_analysis,
+explanatoryFeatures(ttest_analysis,
threshold = 0.05)
-#> # A tibble: 73 × 14
+#> # A tibble: 73 × 14
#> Response Comparison Feature estimate estimate1 estimate2 statistic p.value
-#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
-#> 1 day H~5 N163 -735. 19.5 755. -13.8 1.43e-11
-#> 2 day H~5 N341 2445. 2537. 92.6 13.6 2.88e-11
-#> 3 day H~5 N133 1055. 1077. 21.9 13.0 5.44e-11
-#> 4 day H~2 N341 200. 293. 92.6 10.6 1.38e-10
-#> 5 day H~5 N171 62.6 64.7 2.15 11.9 2.62e-10
-#> 6 day H~5 N119 17.2 17.9 0.763 11.0 8.54e-10
-#> 7 day H~5 N342 243. 247. 4.13 10.8 1.42e- 9
-#> 8 day H~5 N343 27.4 28.3 0.961 9.83 5.99e- 9
-#> 9 day H~5 N377 152. 157. 5.05 9.81 6.75e- 9
-#> 10 day H~5 N477 103. 129. 26.1 9.30 1.05e- 8
-#> # … with 63 more rows, and 6 more variables: parameter <dbl>, conf.low <dbl>,
-#> # conf.high <dbl>, method <chr>, alternative <chr>, adjusted.p.value <dbl>
-This will threshold the features based on their adjusted p-value, found in the adjusted.p.value column of the table. The results of all of the features can be returned using the importance() method.
+#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
+#> 1 day H~5 N163 -735. 19.5 755. -13.8 1.43e-11
+#> 2 day H~5 N341 2445. 2537. 92.6 13.6 2.88e-11
+#> 3 day H~5 N133 1055. 1077. 21.9 13.0 5.44e-11
+#> 4 day H~2 N341 200. 293. 92.6 10.6 1.38e-10
+#> 5 day H~5 N171 62.6 64.7 2.15 11.9 2.62e-10
+#> 6 day H~5 N119 17.2 17.9 0.763 11.0 8.54e-10
+#> 7 day H~5 N342 243. 247. 4.13 10.8 1.42e- 9
+#> 8 day H~5 N343 27.4 28.3 0.961 9.83 5.99e- 9
+#> 9 day H~5 N377 152. 157. 5.05 9.81 6.75e- 9
+#> 10 day H~5 N477 103. 129. 26.1 9.30 1.05e- 8
+#> # … with 63 more rows, and 6 more variables: parameter <dbl>, conf.low <dbl>,
+#> # conf.high <dbl>, method <chr>, alternative <chr>, adjusted.p.value <dbl>
+This will threshold the features based on their adjusted p-value, found in the adjusted.p.value column of the table. The results of all of the features can be returned using the importance() method.
A heat map of the explanatory features can be plotted to inspect the relative trends of the explanatory features in relation to the response variable.
plotExplanatoryHeatmap(ttest_analysis)
-
-
-ANOVA
+
+ANOVA
+
ANOVA can be used to select explanatory features for discrete response variables with 3 or more categories. The following example will compare all the categories in the day response variable. However, the comparisons argument can be used to select particular comparisons of interest.
anova_analysis <- anova(d,
cls = 'day')
-print(anova_analysis)
+print(anova_analysis)
#>
#> Univariate ANOVA analysis
#>
@@ -615,22 +611,22 @@
#> # comparisons: 1
The explanatory features that are significantly different between the categories can then be extracted.
-explanatoryFeatures(anova_analysis,
+explanatoryFeatures(anova_analysis,
threshold = 0.05)
-#> # A tibble: 110 × 10
+#> # A tibble: 110 × 10
#> Response Comparison Feature term df sumsq meansq statistic p.value
-#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
-#> 1 day 1~2~3~4~5~H N341 response 5 1.09e8 2.17e7 124. 1.90e-44
-#> 2 day 1~2~3~4~5~H N163 response 5 1.25e7 2.51e6 113. 1.71e-42
-#> 3 day 1~2~3~4~5~H N133 response 5 1.96e7 3.92e6 108. 1.71e-41
-#> 4 day 1~2~3~4~5~H N171 response 5 6.29e4 1.26e4 88.8 1.16e-37
-#> 5 day 1~2~3~4~5~H N342 response 5 1.04e6 2.07e5 85.1 7.61e-37
-#> 6 day 1~2~3~4~5~H N343 response 5 1.19e4 2.38e3 66.1 4.43e-32
-#> 7 day 1~2~3~4~5~H N119 response 5 4.92e3 9.83e2 53.8 2.07e-28
-#> 8 day 1~2~3~4~5~H N497 response 5 1.10e5 2.20e4 49.6 4.83e-27
-#> 9 day 1~2~3~4~5~H N137 response 5 6.32e3 1.26e3 39.9 1.59e-23
-#> 10 day 1~2~3~4~5~H N277 response 5 6.31e4 1.26e4 39.1 3.14e-23
-#> # … with 100 more rows, and 1 more variable: adjusted.p.value <dbl>
+#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
+#> 1 day 1~2~3~4~5~H N341 response 5 1.09e8 2.17e7 124. 1.90e-44
+#> 2 day 1~2~3~4~5~H N163 response 5 1.25e7 2.51e6 113. 1.71e-42
+#> 3 day 1~2~3~4~5~H N133 response 5 1.96e7 3.92e6 108. 1.71e-41
+#> 4 day 1~2~3~4~5~H N171 response 5 6.29e4 1.26e4 88.8 1.16e-37
+#> 5 day 1~2~3~4~5~H N342 response 5 1.04e6 2.07e5 85.1 7.61e-37
+#> 6 day 1~2~3~4~5~H N343 response 5 1.19e4 2.38e3 66.1 4.43e-32
+#> 7 day 1~2~3~4~5~H N119 response 5 4.92e3 9.83e2 53.8 2.07e-28
+#> 8 day 1~2~3~4~5~H N497 response 5 1.10e5 2.20e4 49.6 4.83e-27
+#> 9 day 1~2~3~4~5~H N137 response 5 6.32e3 1.26e3 39.9 1.59e-23
+#> 10 day 1~2~3~4~5~H N277 response 5 6.31e4 1.26e4 39.1 3.14e-23
+#> # … with 100 more rows, and 1 more variable: adjusted.p.value <dbl>
The top ranked explanatory feature N341 can be plotted to inspect it’s trend relative to the day response variable.
plotFeature(anova_analysis,
@@ -638,15 +634,15 @@
cls = 'day')
-
-
-Linear regression
+
+Linear regression
+
Univariate linear regression can be used to associate a continuous response variable with metabolome features. In the example below, the example data will be regressed against injection order to identify any linearly associated metabolome features.
lr_analysis <- linearRegression(d,
cls = 'injorder')
-print(lr_analysis)
+print(lr_analysis)
#>
#> Univariate linear regression analysis
#>
@@ -655,20 +651,20 @@
#> Responses: injorder
The explanatory features can then be extracted.
-explanatoryFeatures(lr_analysis)
-#> # A tibble: 8 × 15
+explanatoryFeatures(lr_analysis)
+#> # A tibble: 8 × 15
#> Response Feature r.squared adj.r.squared sigma statistic p.value df logLik
-#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
-#> 1 injorder N283 0.310 0.304 4.27 53.0 4.10e-11 1 -343.
-#> 2 injorder N135 0.165 0.157 78.7 23.2 4.31e- 6 1 -693.
-#> 3 injorder N221 0.140 0.133 5.87 19.3 2.50e- 5 1 -382.
-#> 4 injorder N473 0.135 0.127 7.24 18.3 3.78e- 5 1 -407.
-#> 5 injorder N335 0.132 0.124 20.1 17.9 4.59e- 5 1 -529.
-#> 6 injorder N452 0.120 0.112 4.00 16.0 1.10e- 4 1 -335.
-#> 7 injorder N255 0.119 0.111 11.1 15.9 1.17e- 4 1 -458.
-#> 8 injorder N267 0.118 0.111 26.4 15.8 1.22e- 4 1 -562.
-#> # … with 6 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>,
-#> # df.residual <int>, nobs <int>, adjusted.p.value <dbl>
+#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
+#> 1 injorder N283 0.310 0.304 4.27 53.0 4.10e-11 1 -343.
+#> 2 injorder N135 0.165 0.157 78.7 23.2 4.31e- 6 1 -693.
+#> 3 injorder N221 0.140 0.133 5.87 19.3 2.50e- 5 1 -382.
+#> 4 injorder N473 0.135 0.127 7.24 18.3 3.78e- 5 1 -407.
+#> 5 injorder N335 0.132 0.124 20.1 17.9 4.59e- 5 1 -529.
+#> 6 injorder N452 0.120 0.112 4.00 16.0 1.10e- 4 1 -335.
+#> 7 injorder N255 0.119 0.111 11.1 15.9 1.17e- 4 1 -458.
+#> 8 injorder N267 0.118 0.111 26.4 15.8 1.22e- 4 1 -562.
+#> # … with 6 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>,
+#> # df.residual <int>, nobs <int>, adjusted.p.value <dbl>
The top ranked explanatory feature N283 can be plotted to inspect inspects it’s association with injection order.
plotFeature(lr_analysis,
@@ -677,36 +673,36 @@
-
-
-Routine analyses
+
+Routine analyses
+
For routine analyses, the initial analysis parameters for pre-treatment of the data and then the modelling can be selected.
-p <- analysisParameters(c('pre-treatment','modelling'))
+p <- analysisParameters(c('pre-treatment','modelling'))
More specific parameters for pre-treatment of the example data can be declared using the following.
-parameters(p,'pre-treatment') <- preTreatmentParameters(
- list(
+parameters(p,'pre-treatment') <- preTreatmentParameters(
+ list(
keep = 'classes',
occupancyFilter = 'maximum',
transform = 'TICnorm'
)
)
-The modellingMethods() function can be used to list the modelling methods that are currently available in metabolyseR.
+The modellingMethods() function can be used to list the modelling methods that are currently available in metabolyseR.
-modellingMethods()
+modellingMethods()
#> [1] "anova" "ttest" "linearRegression" "randomForest"
-The modellingParameters() function can be used to retrieve the default parameters for specific modelling methods. Below, the default modelling parameters for the randomForest and ttest methods are specified.
+The modellingParameters() function can be used to retrieve the default parameters for specific modelling methods. Below, the default modelling parameters for the randomForest and ttest methods are specified.
-parameters(p,'modelling') <- modellingParameters(c('randomForest','ttest'))
+parameters(p,'modelling') <- modellingParameters(c('randomForest','ttest'))
The class parameters can the be universily specified for both the pre-treatment and modelling elements. For this example, the day response variable will be used with just the H and 2 classes.
changeParameter(p,'cls') <- 'day'
-changeParameter(p,'classes') <- c('H','2')
+changeParameter(p,'classes') <- c('H','2')
This gives the following parameters for the analysis.
p
-#> Parameters:
+#> Parameters:
#> pre-treatment
#> keep
#> classes
@@ -736,10 +732,10 @@
#> returnModels = FALSE
The analysis can then be executed.
analysis <- metabolyse(abr1$neg,abr1$fact,p)
-#>
-#> metabolyseR v0.14.6 Wed Nov 17 10:26:39 2021
+#>
[34m
+#> metabolyseR
[39m
[31mv0.14.7
[39m Fri Dec 17 17:57:41 2021
#> ________________________________________________________________________________
-#> Parameters:
+#>
[33m
[33mParameters:
[33m
[39m
#> pre-treatment
#> keep
#> classes
@@ -768,18 +764,18 @@
#> comparisons = list()
#> returnModels = FALSE
#> ________________________________________________________________________________
-#> Pre-treatment …
+#>
[34mPre-treatment
[39m…
-Pre-treatment ✓ [6.9S]
-#> Modelling …
-
-Modelling ✓ [4.7S]
+
[34mPre-treatment
[39m
[32m✓
[39m [6.4S]
+#>
[34mModelling
[39m…
+
[34m
+Modelling
[39m
[32m✓
[39m [4.5S]
#> ________________________________________________________________________________
#>
-#> Complete! [11.7S]
+#>
[32mComplete!
[39m[10.9S]
The results for the modelling can be specifically extracted using the following.
-analysisResults(analysis,'modelling')
+analysisResults(analysis,'modelling')
#> $randomForest
#>
#> Random forest classification
@@ -801,30 +797,30 @@
This returns the results as a list containing the modelling results objects for each specified method.
Alternatively, the modelling results can be assess directly from the Analysis object. Below shows the extraction of the explanatory features, using default parameters for each method, with the results returned in a single table.
-explanatory_features <- analysis %>%
- explanatoryFeatures()
+explanatory_features <- analysis %>%
+ explanatoryFeatures()
-print(explanatory_features)
-#> # A tibble: 100 × 17
+print(explanatory_features)
+#> # A tibble: 100 × 17
#> Method Response Comparison Feature Metric Value estimate estimate1
-#> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
-#> 1 randomForest day 2~H N341 FalsePo… 8.06e-28 NA NA
-#> 2 randomForest day 2~H N377 FalsePo… 5.70e-18 NA NA
-#> 3 randomForest day 2~H N447 FalsePo… 5.70e-18 NA NA
-#> 4 randomForest day 2~H N579 FalsePo… 5.70e-18 NA NA
-#> 5 randomForest day 2~H N1084 FalsePo… 1.19e-16 NA NA
-#> 6 randomForest day 2~H N327 FalsePo… 2.33e-15 NA NA
-#> 7 randomForest day 2~H N580 FalsePo… 4.32e-14 NA NA
-#> 8 randomForest day 2~H N1083 FalsePo… 7.49e-13 NA NA
-#> 9 randomForest day 2~H N1085 FalsePo… 7.49e-13 NA NA
-#> 10 randomForest day 2~H N503 FalsePo… 7.49e-13 NA NA
-#> # … with 90 more rows, and 9 more variables: estimate2 <dbl>, statistic <dbl>,
-#> # p.value <dbl>, parameter <dbl>, conf.low <dbl>, conf.high <dbl>,
-#> # method <chr>, alternative <chr>, adjusted.p.value <dbl>
+#> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
+#> 1 randomForest day 2~H N341 FalsePo… 8.06e-28 NA NA
+#> 2 randomForest day 2~H N377 FalsePo… 5.70e-18 NA NA
+#> 3 randomForest day 2~H N447 FalsePo… 5.70e-18 NA NA
+#> 4 randomForest day 2~H N579 FalsePo… 5.70e-18 NA NA
+#> 5 randomForest day 2~H N1084 FalsePo… 1.19e-16 NA NA
+#> 6 randomForest day 2~H N327 FalsePo… 2.33e-15 NA NA
+#> 7 randomForest day 2~H N580 FalsePo… 4.32e-14 NA NA
+#> 8 randomForest day 2~H N1083 FalsePo… 7.49e-13 NA NA
+#> 9 randomForest day 2~H N1085 FalsePo… 7.49e-13 NA NA
+#> 10 randomForest day 2~H N503 FalsePo… 7.49e-13 NA NA
+#> # … with 90 more rows, and 9 more variables: estimate2 <dbl>, statistic <dbl>,
+#> # p.value <dbl>, parameter <dbl>, conf.low <dbl>, conf.high <dbl>,
+#> # method <chr>, alternative <chr>, adjusted.p.value <dbl>
Heat maps of the explanatory features can also be plotted for both the modelling methods.
-plotExplanatoryHeatmap(analysis) %>%
- patchwork::wrap_plots()
+plotExplanatoryHeatmap(analysis) %>%
+ patchwork::wrap_plots()
@@ -840,11 +836,13 @@
@@ -853,5 +851,7 @@
+
+
-
-
+
@@ -85,35 +81,35 @@
- +----
Modelling and feature selection
-metabolyseR v0.14.6
-Jasen Finch
+metabolyseR v0.14.7
+Jasen Finch
-17 November, 2021
+17 December, 2021
- Source:vignettes/modelling.Rmd+ Source:vignettes/modelling.Rmdmodelling.Rmd--Introduction
++Introduction +
Modelling provides the essential data mining step for extracting biological information and explanatory metabolome features from a data set relating to the experimental conditions.
metabolyseRprovides a number of both univariate and multivariate methods for data mining.For an introduction to the usage of metabolyseR for both exploratory and routine analyses, see the introduction vignette using:
+-
vignette('introduction','metabolyseR')vignette('introduction','metabolyseR')To further supplement this document, a quick start example analysis is also available as a vignette:
+-
vignette('quick_start','metabolyseR')vignette('quick_start','metabolyseR')To begin, the package can be loaded using:
--
library(metabolyseR) +library(metabolyseR) #> #> Attaching package: 'metabolyseR' #> The following object is masked from 'package:stats': @@ -122,17 +118,17 @@#> The following objects are masked from 'package:base': #> #> raw, split
---Example data
-The examples used here will use the
+abr1data set from the metaboData package. This is nominal mass flow-injection mass spectrometry (FI-MS) fingerprinting data from a plant-pathogen infection time course experiment. The pipe%>%from the magrittr package will also be used. The example data can be loaded using:+Example data +
+The examples used here will use the
+abr1data set from the metaboData package. This is nominal mass flow-injection mass spectrometry (FI-MS) fingerprinting data from a plant-pathogen infection time course experiment. The pipe%>%from the magrittr package will also be used. The example data can be loaded using:library(metaboData)Only the negative acquisition mode data (
abr1$neg) will be used along with the sample meta-information (abr1$fact). Create anAnalysisDataclass object, assigned to the variabled, using the following:d <- analysisData(abr1$neg[,1:500],abr1$fact)As can be seen above the data set contains a total of 120 samples and 500 features.
--Parallel processing
-The package supports parallel processing using the future package.
-By default, processing by
+metabolyseRwill be done seqentially. However, parallel processing can be activated, prior to analysis, by specifying a parallel implementation usingplan(). The following example specifies using themultisessionimplementation (muliple background R sessions) with two worker processes.++Parallel processing +
+The package supports parallel processing using the future package.
+By default, processing by
metabolyseRwill be done seqentially. However, parallel processing can be activated, prior to analysis, by specifying a parallel implementation usingplan(). The following example specifies using themultisessionimplementation (muliple background R sessions) with two worker processes.--
plan(future::multisession,workers = 2)See the future package documentation for more information on the types of parallel implementations that are available.
+plan(future::multisession,workers = 2)See the future package documentation for more information on the types of parallel implementations that are available.
---Random Forest
++Random Forest +
Random forest is a versatile ensemble machine learning approach based on forests of decision trees for multivariate data mining. This can include unsupervised analysis, classification of discrete response variables and regression of continuous responses.
-Random forest can be performed in
-metabolyseRusing therandomForest()method. For further details on the arguments for using this function, see?randomForest. This implementation of random forest inmetabolyseRutilises therandomForestpackage. See?randomForest::randomForestfor more information about that implementation.---Unsupervised
+Random forest can be performed in
+metabolyseRusing therandomForest()method. For further details on the arguments for using this function, see?randomForest. This implementation of random forest inmetabolyseRutilises therandomForestpackage. See?randomForest::randomForestfor more information about that implementation.+Unsupervised +
The unsupervised random forest approach can be useful starting point for analysis in any experimental context. It can be used to give a general overview of the structure of the data and to identify any possible problems. These could include situations such as the presence of outliers samples or splits in the data caused by the impact of analytical or sample preparation factors. Unsupervised random forest can have advantages in these assessments over other approaches such as Principle Component Analysis (PCA). It is less sensitive to the effect of a single feature that in fact could have little overall impact relative to the other hundreds that could be present in a data set.
The examples below will show the use of unsupervised random forest for assessing the general structure of the example data set and the presence of outlier samples.
Unsupervised random forest can be performed by setting the
clsargument ofrandomForest()toNULL:-
unsupervised_rf <- d %>% +unsupervised_rf <- d %>% randomForest(cls = NULL)The type of random forest that has been performed can be checked using the
typemethod.Or by printing the results object.
@@ -196,75 +192,75 @@
This shows that it is indeed the experimental factor of interest that is having the greatest impact on the structure of the data. The progression of the experimental time points are obvious across Dimension 1.
The available feature importance metrics for a random forest analysis can be retrieved by:
-
importanceMetrics(unsupervised_rf) +importanceMetrics(unsupervised_rf) #> [1] "1" "2" "FalsePositiveRate" #> [4] "MeanDecreaseAccuracy" "MeanDecreaseGini" "SelectionFrequency"And the importance values of these metrics for each feature can returned using:
+#> <chr> <chr> <dbl> +#> 1 N1 1 0 +#> 2 N1 2 0 +#> 3 N1 FalsePositiveRate 0.0238 +#> 4 N1 MeanDecreaseAccuracy 0 +#> 5 N1 MeanDecreaseGini 0 +#> 6 N1 SelectionFrequency 0 +#> 7 N10 1 0 +#> 8 N10 2 0 +#> 9 N10 FalsePositiveRate 0.0238 +#> 10 N10 MeanDecreaseAccuracy 0 +#> # … with 2,990 more rows-
importance(unsupervised_rf) -#> # A tibble: 3,000 × 3 +importance(unsupervised_rf) +#> # A tibble: 3,000 × 3 #> Feature Metric Value -#> <chr> <chr> <dbl> -#> 1 N1 1 0 -#> 2 N1 2 0 -#> 3 N1 FalsePositiveRate 0.0238 -#> 4 N1 MeanDecreaseAccuracy 0 -#> 5 N1 MeanDecreaseGini 0 -#> 6 N1 SelectionFrequency 0 -#> 7 N10 1 0 -#> 8 N10 2 0 -#> 9 N10 FalsePositiveRate 0.0238 -#> 10 N10 MeanDecreaseAccuracy 0 -#> # … with 2,990 more rowsThe explanatory features for a given threshold can be extracted for any of the importance metrics. The following will extract the explanatory features below a threshold of 0.05 based on the false positive rate metric.
+#> <chr> <chr> <dbl> +#> 1 N342 FalsePositiveRate 1.31e-19 +#> 2 N161 FalsePositiveRate 2.34e-16 +#> 3 N341 FalsePositiveRate 6.50e-16 +#> 4 N315 FalsePositiveRate 1.79e-15 +#> 5 N367 FalsePositiveRate 3.47e-14 +#> 6 N173 FalsePositiveRate 9.09e-14 +#> 7 N385 FalsePositiveRate 9.09e-14 +#> 8 N133 FalsePositiveRate 1.52e-12 +#> 9 N439 FalsePositiveRate 1.52e-12 +#> 10 N379 FalsePositiveRate 3.78e-12 +#> # … with 349 more rows-
unsupervised_rf %>% - explanatoryFeatures(metric = "FalsePositiveRate", +unsupervised_rf %>% + explanatoryFeatures(metric = "FalsePositiveRate", threshold = 0.05) -#> # A tibble: 359 × 3 +#> # A tibble: 359 × 3 #> Feature Metric Value -#> <chr> <chr> <dbl> -#> 1 N342 FalsePositiveRate 1.31e-19 -#> 2 N161 FalsePositiveRate 2.34e-16 -#> 3 N341 FalsePositiveRate 6.50e-16 -#> 4 N315 FalsePositiveRate 1.79e-15 -#> 5 N367 FalsePositiveRate 3.47e-14 -#> 6 N173 FalsePositiveRate 9.09e-14 -#> 7 N385 FalsePositiveRate 9.09e-14 -#> 8 N133 FalsePositiveRate 1.52e-12 -#> 9 N439 FalsePositiveRate 1.52e-12 -#> 10 N379 FalsePositiveRate 3.78e-12 -#> # … with 349 more rowsIn this example there are 359 explanatory features.
The trend of the most highly ranked explanatory feature against the
dayfactor can be plotted using theplotFeature()method.-
unsupervised_rf %>% +unsupervised_rf %>% plotFeature(feature = 'N425', cls = 'day')--Classification
++Classification +
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
++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
dayresponse variable.-
d %>% - clsExtract(cls = 'day') %>% - unique() +d %>% + clsExtract(cls = 'day') %>% + unique() #> [1] 2 3 4 1 H 5 #> Levels: 1 2 3 4 5 HMultinomial classification can be performed by:
-
multinomial_rf <- d %>% +multinomial_rf <- d %>% randomForest(cls = 'day') -print(multinomial_rf) +print(multinomial_rf) #> #> Random forest classification #> @@ -274,117 +270,117 @@#> # 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:
+#> <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-
multinomial_rf %>% - metrics() -#> # A tibble: 4 × 5 +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.146These 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() +d %>% + clsExtract(cls = 'day') %>% + table() #> . #> 1 2 3 4 5 H #> 20 20 20 20 20 20In 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.
This shows that the model is able to discriminate highly between classes such as
5andH. It is less able to discriminate more similar classes such asHand1or4and5whose 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
Hand1classes 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) +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
dayresponse variable.+#> <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-
multinomial_rf %>% - explanatoryFeatures(metric = 'FalsePositiveRate', +multinomial_rf %>% + explanatoryFeatures(metric = 'FalsePositiveRate', threshold = 0.05) -#> # A tibble: 121 × 5 +#> # 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 rowsAs 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
comparisonsargument. This should be specified as a named list, the corresponding to theclsargument. 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 thedayresponse variable and displays the performance metrics.+#> <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-
d %>% +d %>% randomForest(cls = 'day', - comparisons = list(day = c('H~1~2', - 'H~1~5'))) %>% - metrics() -#> # A tibble: 8 × 5 + 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.320The MDS and ROC curve plots can also be plotted simultaneously for the two comparisons.
+ patchwork::plot_layout(ncol = 1)}-
d %>% +d %>% randomForest(cls = 'day', - comparisons = list(day = c('H~1~2', - 'H~1~5'))) %>% + comparisons = list(day = c('H~1~2', + 'H~1~5'))) %>% {plotMDS(.,cls = 'day') + plotROC(.) + - patchwork::plot_layout(ncol = 1)}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
clsargument. In the following, both thenameanddayresponse factors will be analysed and the performance metrics returned in a single table.+#> <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-
d %>% - randomForest(cls = c('name','day')) %>% - metrics() +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 +#> # 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.146The MDS plots can also be returned for both models simultaneously.
-
d %>% - randomForest(cls = c('name','day')) %>% +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", @@ -393,26 +389,26 @@#> Unbalanced classes detected. Stratifying sample size to the smallest class size.
--Binary comparisons
++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 thedayresponse variable.The possible binary comparisons for a given response variable can be displayed using the
binaryComparisons()method. Below shows the 15 comparisons for thedayresponse variable.-
binaryComparisons(d,cls = 'day') +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
Hclass.+-
binary_comparisons <- binaryComparisons(d,cls = 'day') %>% - .[stringr::str_detect(.,'H')]binary_comparisons <- binaryComparisons(d,cls = 'day') %>% + .[stringr::str_detect(.,'H')]The binary comparisons can then be performed using the following.
-
binary_rf <- d %>% +binary_rf <- d %>% randomForest(cls = 'day', - comparisons = list(day = binary_comparisons)) + comparisons = list(day = binary_comparisons)) -print(binary_rf) +print(binary_rf) #> #> Random forest classification #> @@ -423,59 +419,59 @@To run all possible binary comparisons, the
binary = TRUEargument could instead be used.The MDS plots for each comparison can be visualised to inspect the comparisons.
These plots show good separation in all the comparisons except
H~1which 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 theH~1comparison.-
binary_rf %>% +binary_rf %>% plotMetrics()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.
+#> <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-
binary_rf %>% - explanatoryFeatures(metric = 'FalsePositiveRate', +binary_rf %>% + explanatoryFeatures(metric = 'FalsePositiveRate', threshold = 0.05) -#> # A tibble: 251 × 5 +#> # 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 rowsA 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')) +refactor_cls <- clsExtract(binary_rf, + cls = 'day') %>% + factor(.,levels = c('H','1','2','3','4','5')) -binary_rf <- clsReplace(binary_rf, +binary_rf <- clsReplace(binary_rf, value = refactor_cls, cls = 'day') -binary_rf %>% +binary_rf %>% plotExplanatoryHeatmap(metric = 'FalsePositiveRate', threshold = 0.05, featureNames = TRUE)--Regression
++Regression +
Random forest regression can be used to assess the extent of association of the metabolomic data with continuous response variables.
In this example, the extent of association of injection order with the example data will be assessed.
-
regression_rf <- d %>% +regression_rf <- d %>% randomForest(cls = 'injorder') -print(regression_rf) +print(regression_rf) #> #> Random forest regression #> @@ -484,21 +480,21 @@#> Response: injorder
The regression model performance metrics, based on the OOB prediction error, can be extracted using the following:
+#> <chr> <chr> <chr> <dbl> +#> 1 injorder rsq standard 0.476 +#> 2 injorder mae standard 23.5 +#> 3 injorder mape standard 154. +#> 4 injorder rmse standard 26.5 +#> 5 injorder ccc standard 0.508-
regression_rf %>% - metrics() -#> # A tibble: 5 × 4 +regression_rf %>% + metrics() +#> # A tibble: 5 × 4 #> Response .metric .estimator .estimate -#> <chr> <chr> <chr> <dbl> -#> 1 injorder rsq standard 0.476 -#> 2 injorder mae standard 23.5 -#> 3 injorder mape standard 154. -#> 4 injorder rmse standard 26.5 -#> 5 injorder ccc standard 0.508These regression metrics include R2 (
rsq), mean absolute error (mae), mean absolute percentage error (mape), root mean squared error (rmse) and the concordance correlation coefficient (ccc).The R2 and concordance correlation coefficient metrics suggest that there is some association of features with the injection order, although this is weak. This is in agreement with mean absolute error metric that shows that on average, the injection order could only be predicted to an accuracy of 23 injection order positions.
The MDS plot belows the relative proximities of the samples based on this injection order regression model. This shows that for the most part, there is little correspondence of the sample positions with their injection order. However, there is a small grouping of samples towards the end of the run around sample ~99 to 120. It suggests that there could have been some analytical issues, for certain features, towards the end of the mass spectral analytical run.
-
regression_rf %>% +regression_rf %>% plotMDS(cls = NULL, ellipses = FALSE, label = 'injorder', @@ -508,51 +504,51 @@The available feature importance metrics for this regression model can be listed.
-
regression_rf %>% - importanceMetrics() +regression_rf %>% + importanceMetrics() #> [1] "%IncMSE" "IncNodePurity"The feature importance metrics can be plotted to give an overview of their distribution. The following will plot the percentage increase in the mean squared error (
%IncMSE) importance metric.-
regression_rf %>% +regression_rf %>% plotImportance(metric = "%IncMSE", rank = FALSE)This shows that there are only a few features that are contributing to the association with injection order. These explanatory features can be extracted with the following, using a threshold of above 5.
+#> <chr> <chr> <chr> <dbl> +#> 1 injorder N283 %IncMSE 19.9 +#> 2 injorder N135 %IncMSE 8.71 +#> 3 injorder N451 %IncMSE 5.58 +#> 4 injorder N161 %IncMSE 5.51 +#> 5 injorder N306 %IncMSE 5.49 +#> 6 injorder N118 %IncMSE 5.22 +#> 7 injorder N297 %IncMSE 5.07-
regression_rf %>% - explanatoryFeatures(metric = '%IncMSE', +regression_rf %>% + explanatoryFeatures(metric = '%IncMSE', threshold = 5) -#> # A tibble: 7 × 4 +#> # A tibble: 7 × 4 #> Response Feature Metric Value -#> <chr> <chr> <chr> <dbl> -#> 1 injorder N283 %IncMSE 19.9 -#> 2 injorder N135 %IncMSE 8.71 -#> 3 injorder N451 %IncMSE 5.58 -#> 4 injorder N161 %IncMSE 5.51 -#> 5 injorder N306 %IncMSE 5.49 -#> 6 injorder N118 %IncMSE 5.22 -#> 7 injorder N297 %IncMSE 5.07This returned a total of 7 explanatory features above this threshold. The top ranked feature
N283can be plotted to investigate it’s trend in relation to injection order.-
regression_rf %>% +regression_rf %>% plotFeature(feature = 'N283', cls = 'injorder')This shows an increase in the intensity of that feature for samples above 100 in the injection order which corresponds with the cluster that was seen in the MDS plot above.
--Univariate analyses
++Univariate analyses +
Univariate methods select features, explanatory for response variables, with features tested on an individual basis. These methods offer simplicity and easy interpretation in their use, however they provide no information as to how features may interact.
The univariate methods currently available in
-metabolyseRinclude Welch’s t-test, analysis of variance (ANOVA) and linear regression. The following sections will provide brief examples of the use of each of these methods.---Welch’s t-test
+++Welch’s t-test +
Welch’s t-test can be used to select explanatory metabolome features for binary comparisons of discrete variables. By default, all the possible binary comparisons for the categories of a response variable will be tested.
Below shows the possible binary comparisons for the
dayresponse variable for the example data set.@@ -560,11 +556,11 @@-
binaryComparisons(d, +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"The explanatory features that show a significant difference between the response categories can be extracted as shown below.
--
explanatoryFeatures(ttest_analysis, +explanatoryFeatures(ttest_analysis, threshold = 0.05) -#> # A tibble: 73 × 14 +#> # A tibble: 73 × 14 #> Response Comparison Feature estimate estimate1 estimate2 statistic p.value -#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> -#> 1 day H~5 N163 -735. 19.5 755. -13.8 1.43e-11 -#> 2 day H~5 N341 2445. 2537. 92.6 13.6 2.88e-11 -#> 3 day H~5 N133 1055. 1077. 21.9 13.0 5.44e-11 -#> 4 day H~2 N341 200. 293. 92.6 10.6 1.38e-10 -#> 5 day H~5 N171 62.6 64.7 2.15 11.9 2.62e-10 -#> 6 day H~5 N119 17.2 17.9 0.763 11.0 8.54e-10 -#> 7 day H~5 N342 243. 247. 4.13 10.8 1.42e- 9 -#> 8 day H~5 N343 27.4 28.3 0.961 9.83 5.99e- 9 -#> 9 day H~5 N377 152. 157. 5.05 9.81 6.75e- 9 -#> 10 day H~5 N477 103. 129. 26.1 9.30 1.05e- 8 -#> # … with 63 more rows, and 6 more variables: parameter <dbl>, conf.low <dbl>, -#> # conf.high <dbl>, method <chr>, alternative <chr>, adjusted.p.value <dbl>This will threshold the features based on their adjusted p-value, found in the
+#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> +#> 1 day H~5 N163 -735. 19.5 755. -13.8 1.43e-11 +#> 2 day H~5 N341 2445. 2537. 92.6 13.6 2.88e-11 +#> 3 day H~5 N133 1055. 1077. 21.9 13.0 5.44e-11 +#> 4 day H~2 N341 200. 293. 92.6 10.6 1.38e-10 +#> 5 day H~5 N171 62.6 64.7 2.15 11.9 2.62e-10 +#> 6 day H~5 N119 17.2 17.9 0.763 11.0 8.54e-10 +#> 7 day H~5 N342 243. 247. 4.13 10.8 1.42e- 9 +#> 8 day H~5 N343 27.4 28.3 0.961 9.83 5.99e- 9 +#> 9 day H~5 N377 152. 157. 5.05 9.81 6.75e- 9 +#> 10 day H~5 N477 103. 129. 26.1 9.30 1.05e- 8 +#> # … with 63 more rows, and 6 more variables: parameter <dbl>, conf.low <dbl>, +#> # conf.high <dbl>, method <chr>, alternative <chr>, adjusted.p.value <dbl>adjusted.p.valuecolumn of the table. The results of all of the features can be returned using theimportance()method.This will threshold the features based on their adjusted p-value, found in the
adjusted.p.valuecolumn of the table. The results of all of the features can be returned using theimportance()method.A heat map of the explanatory features can be plotted to inspect the relative trends of the explanatory features in relation to the response
variable.plotExplanatoryHeatmap(ttest_analysis)---ANOVA
++ANOVA +
ANOVA can be used to select explanatory features for discrete response variables with 3 or more categories. The following example will compare all the categories in the
dayresponse variable. However, thecomparisonsargument can be used to select particular comparisons of interest.anova_analysis <- anova(d, cls = 'day') -print(anova_analysis) +print(anova_analysis) #> #> Univariate ANOVA analysis #> @@ -615,22 +611,22 @@#> # comparisons: 1
The explanatory features that are significantly different between the categories can then be extracted.
+#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> +#> 1 day 1~2~3~4~5~H N341 response 5 1.09e8 2.17e7 124. 1.90e-44 +#> 2 day 1~2~3~4~5~H N163 response 5 1.25e7 2.51e6 113. 1.71e-42 +#> 3 day 1~2~3~4~5~H N133 response 5 1.96e7 3.92e6 108. 1.71e-41 +#> 4 day 1~2~3~4~5~H N171 response 5 6.29e4 1.26e4 88.8 1.16e-37 +#> 5 day 1~2~3~4~5~H N342 response 5 1.04e6 2.07e5 85.1 7.61e-37 +#> 6 day 1~2~3~4~5~H N343 response 5 1.19e4 2.38e3 66.1 4.43e-32 +#> 7 day 1~2~3~4~5~H N119 response 5 4.92e3 9.83e2 53.8 2.07e-28 +#> 8 day 1~2~3~4~5~H N497 response 5 1.10e5 2.20e4 49.6 4.83e-27 +#> 9 day 1~2~3~4~5~H N137 response 5 6.32e3 1.26e3 39.9 1.59e-23 +#> 10 day 1~2~3~4~5~H N277 response 5 6.31e4 1.26e4 39.1 3.14e-23 +#> # … with 100 more rows, and 1 more variable: adjusted.p.value <dbl>-
explanatoryFeatures(anova_analysis, +explanatoryFeatures(anova_analysis, threshold = 0.05) -#> # A tibble: 110 × 10 +#> # A tibble: 110 × 10 #> Response Comparison Feature term df sumsq meansq statistic p.value -#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> -#> 1 day 1~2~3~4~5~H N341 response 5 1.09e8 2.17e7 124. 1.90e-44 -#> 2 day 1~2~3~4~5~H N163 response 5 1.25e7 2.51e6 113. 1.71e-42 -#> 3 day 1~2~3~4~5~H N133 response 5 1.96e7 3.92e6 108. 1.71e-41 -#> 4 day 1~2~3~4~5~H N171 response 5 6.29e4 1.26e4 88.8 1.16e-37 -#> 5 day 1~2~3~4~5~H N342 response 5 1.04e6 2.07e5 85.1 7.61e-37 -#> 6 day 1~2~3~4~5~H N343 response 5 1.19e4 2.38e3 66.1 4.43e-32 -#> 7 day 1~2~3~4~5~H N119 response 5 4.92e3 9.83e2 53.8 2.07e-28 -#> 8 day 1~2~3~4~5~H N497 response 5 1.10e5 2.20e4 49.6 4.83e-27 -#> 9 day 1~2~3~4~5~H N137 response 5 6.32e3 1.26e3 39.9 1.59e-23 -#> 10 day 1~2~3~4~5~H N277 response 5 6.31e4 1.26e4 39.1 3.14e-23 -#> # … with 100 more rows, and 1 more variable: adjusted.p.value <dbl>The top ranked explanatory feature
N341can be plotted to inspect it’s trend relative to thedayresponse variable.plotFeature(anova_analysis, @@ -638,15 +634,15 @@cls = 'day')
---Linear regression
++Linear regression +
Univariate linear regression can be used to associate a continuous response variable with metabolome features. In the example below, the example data will be regressed against injection order to identify any linearly associated metabolome features.
lr_analysis <- linearRegression(d, cls = 'injorder') -print(lr_analysis) +print(lr_analysis) #> #> Univariate linear regression analysis #> @@ -655,20 +651,20 @@#> Responses: injorder
The explanatory features can then be extracted.
+#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> +#> 1 injorder N283 0.310 0.304 4.27 53.0 4.10e-11 1 -343. +#> 2 injorder N135 0.165 0.157 78.7 23.2 4.31e- 6 1 -693. +#> 3 injorder N221 0.140 0.133 5.87 19.3 2.50e- 5 1 -382. +#> 4 injorder N473 0.135 0.127 7.24 18.3 3.78e- 5 1 -407. +#> 5 injorder N335 0.132 0.124 20.1 17.9 4.59e- 5 1 -529. +#> 6 injorder N452 0.120 0.112 4.00 16.0 1.10e- 4 1 -335. +#> 7 injorder N255 0.119 0.111 11.1 15.9 1.17e- 4 1 -458. +#> 8 injorder N267 0.118 0.111 26.4 15.8 1.22e- 4 1 -562. +#> # … with 6 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>, +#> # df.residual <int>, nobs <int>, adjusted.p.value <dbl>-
explanatoryFeatures(lr_analysis) -#> # A tibble: 8 × 15 +explanatoryFeatures(lr_analysis) +#> # A tibble: 8 × 15 #> Response Feature r.squared adj.r.squared sigma statistic p.value df logLik -#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> -#> 1 injorder N283 0.310 0.304 4.27 53.0 4.10e-11 1 -343. -#> 2 injorder N135 0.165 0.157 78.7 23.2 4.31e- 6 1 -693. -#> 3 injorder N221 0.140 0.133 5.87 19.3 2.50e- 5 1 -382. -#> 4 injorder N473 0.135 0.127 7.24 18.3 3.78e- 5 1 -407. -#> 5 injorder N335 0.132 0.124 20.1 17.9 4.59e- 5 1 -529. -#> 6 injorder N452 0.120 0.112 4.00 16.0 1.10e- 4 1 -335. -#> 7 injorder N255 0.119 0.111 11.1 15.9 1.17e- 4 1 -458. -#> 8 injorder N267 0.118 0.111 26.4 15.8 1.22e- 4 1 -562. -#> # … with 6 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>, -#> # df.residual <int>, nobs <int>, adjusted.p.value <dbl>The top ranked explanatory feature
N283can be plotted to inspect inspects it’s association with injection order.plotFeature(lr_analysis, @@ -677,36 +673,36 @@--Routine analyses
++Routine analyses +
For routine analyses, the initial analysis parameters for pre-treatment of the data and then the modelling can be selected.
+-
p <- analysisParameters(c('pre-treatment','modelling'))p <- analysisParameters(c('pre-treatment','modelling'))More specific parameters for pre-treatment of the example data can be declared using the following.
--
parameters(p,'pre-treatment') <- preTreatmentParameters( - list( +parameters(p,'pre-treatment') <- preTreatmentParameters( + list( keep = 'classes', occupancyFilter = 'maximum', transform = 'TICnorm' ) )The
+modellingMethods()function can be used to list the modelling methods that are currently available inmetabolyseR.The
modellingMethods()function can be used to list the modelling methods that are currently available inmetabolyseR.--
modellingMethods() +modellingMethods() #> [1] "anova" "ttest" "linearRegression" "randomForest"The
+modellingParameters()function can be used to retrieve the default parameters for specific modelling methods. Below, the default modelling parameters for therandomForestandttestmethods are specified.The
modellingParameters()function can be used to retrieve the default parameters for specific modelling methods. Below, the default modelling parameters for therandomForestandttestmethods are specified.+-
parameters(p,'modelling') <- modellingParameters(c('randomForest','ttest'))parameters(p,'modelling') <- modellingParameters(c('randomForest','ttest'))The class parameters can the be universily specified for both the pre-treatment and modelling elements. For this example, the
dayresponse variable will be used with just theHand2classes.+changeParameter(p,'classes') <- c('H','2')changeParameter(p,'cls') <- 'day' -changeParameter(p,'classes') <- c('H','2')This gives the following parameters for the analysis.
p -#> Parameters: +#> Parameters: #> pre-treatment #> keep #> classes @@ -736,10 +732,10 @@#> returnModels = FALSE
The analysis can then be executed.
+#> [32mComplete! [39m[10.9S]analysis <- metabolyse(abr1$neg,abr1$fact,p) -#> -#> metabolyseR v0.14.6 Wed Nov 17 10:26:39 2021 +#> [34m +#> metabolyseR [39m [31mv0.14.7 [39m Fri Dec 17 17:57:41 2021 #> ________________________________________________________________________________ -#> Parameters: +#> [33m [33mParameters: [33m [39m #> pre-treatment #> keep #> classes @@ -768,18 +764,18 @@#> comparisons = list() #> returnModels = FALSE #> ________________________________________________________________________________ -#> Pre-treatment … +#> [34mPre-treatment [39m… -Pre-treatment ✓ [6.9S] -#> Modelling … - -Modelling ✓ [4.7S] + [34mPre-treatment [39m [32m✓ [39m [6.4S] +#> [34mModelling [39m… + [34m +Modelling [39m [32m✓ [39m [4.5S] #> ________________________________________________________________________________ #> -#> Complete! [11.7S]
The results for the modelling can be specifically extracted using the following.
-
analysisResults(analysis,'modelling') +analysisResults(analysis,'modelling') #> $randomForest #> #> Random forest classification @@ -801,30 +797,30 @@This returns the results as a list containing the modelling results objects for each specified method.
Alternatively, the modelling results can be assess directly from the
Analysisobject. Below shows the extraction of the explanatory features, using default parameters for each method, with the results returned in a single table.+#> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> +#> 1 randomForest day 2~H N341 FalsePo… 8.06e-28 NA NA +#> 2 randomForest day 2~H N377 FalsePo… 5.70e-18 NA NA +#> 3 randomForest day 2~H N447 FalsePo… 5.70e-18 NA NA +#> 4 randomForest day 2~H N579 FalsePo… 5.70e-18 NA NA +#> 5 randomForest day 2~H N1084 FalsePo… 1.19e-16 NA NA +#> 6 randomForest day 2~H N327 FalsePo… 2.33e-15 NA NA +#> 7 randomForest day 2~H N580 FalsePo… 4.32e-14 NA NA +#> 8 randomForest day 2~H N1083 FalsePo… 7.49e-13 NA NA +#> 9 randomForest day 2~H N1085 FalsePo… 7.49e-13 NA NA +#> 10 randomForest day 2~H N503 FalsePo… 7.49e-13 NA NA +#> # … with 90 more rows, and 9 more variables: estimate2 <dbl>, statistic <dbl>, +#> # p.value <dbl>, parameter <dbl>, conf.low <dbl>, conf.high <dbl>, +#> # method <chr>, alternative <chr>, adjusted.p.value <dbl>-
explanatory_features <- analysis %>% - explanatoryFeatures() +explanatory_features <- analysis %>% + explanatoryFeatures() -print(explanatory_features) -#> # A tibble: 100 × 17 +print(explanatory_features) +#> # A tibble: 100 × 17 #> Method Response Comparison Feature Metric Value estimate estimate1 -#> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> -#> 1 randomForest day 2~H N341 FalsePo… 8.06e-28 NA NA -#> 2 randomForest day 2~H N377 FalsePo… 5.70e-18 NA NA -#> 3 randomForest day 2~H N447 FalsePo… 5.70e-18 NA NA -#> 4 randomForest day 2~H N579 FalsePo… 5.70e-18 NA NA -#> 5 randomForest day 2~H N1084 FalsePo… 1.19e-16 NA NA -#> 6 randomForest day 2~H N327 FalsePo… 2.33e-15 NA NA -#> 7 randomForest day 2~H N580 FalsePo… 4.32e-14 NA NA -#> 8 randomForest day 2~H N1083 FalsePo… 7.49e-13 NA NA -#> 9 randomForest day 2~H N1085 FalsePo… 7.49e-13 NA NA -#> 10 randomForest day 2~H N503 FalsePo… 7.49e-13 NA NA -#> # … with 90 more rows, and 9 more variables: estimate2 <dbl>, statistic <dbl>, -#> # p.value <dbl>, parameter <dbl>, conf.low <dbl>, conf.high <dbl>, -#> # method <chr>, alternative <chr>, adjusted.p.value <dbl>Heat maps of the explanatory features can also be plotted for both the modelling methods.
+-
plotExplanatoryHeatmap(analysis) %>% - patchwork::wrap_plots()plotExplanatoryHeatmap(analysis) %>% + patchwork::wrap_plots()