Getting Started with ClassifyR

Case Study: Diagnosing Asthma

To demonstrate some key features of ClassifyR, a data set consisting of the 2000 most variably expressed genes and 190 people will be used to quickly obtain results. The journal article corresponding to the data set was published in Scientific Reports in 2018 and is titled A Nasal Brush-based Classifier of Asthma Identified by Machine Learning Analysis of Nasal RNA Sequence Data.

Load the package.

library(ClassifyR)

A glimpse at the RNA measurements and sample classes.

data(asthma) # Contains measurements and classes variables.
measurements[1:5, 1:5]

##            HBB BPIFA1  XIST FCGR3B HBA2
## Sample 1  9.72  14.06 12.28  11.42 7.83
## Sample 2 11.98  13.89  6.35  13.25 9.42
## Sample 3 12.15  17.44 10.21   7.87 9.68
## Sample 4 10.60  11.87  6.27  14.75 8.96
## Sample 5  8.18  15.01 11.21   6.77 6.43

head(classes)

## [1] No  No  No  No  Yes No 
## Levels: No Yes

The numeric matrix variable measurements stores the normalised values of the RNA gene abundances for each sample and the factor vector classes identifies which class the samples belong to. The measurements were normalised using DESeq2’s varianceStabilizingTransformation function, which produces log₂-like data.

For more complex data sets with multiple kinds of experiments (e.g. DNA methylation, copy number, gene expression on the same set of samples) a MultiAssayExperiment is recommended for data storage and supported by ClassifyR’s methods.

Data Integration with crossValidate

crossValidate also allows data from multiple sources to be integrated into a single model. The integration method can be specified with multiViewMethod argument. In this example, suppose the first 10 variables in the asthma data set are from a certain source and the remaining 1990 variables are from a second source. To integrate multiple data sets, each variable must be labeled with the data set it came from. This is done in a different manner depending on the data type of measurements.

If using Bioconductor’s DataFrame, this can be specified using mcols. In the column metadata, each feature must have an assay and a feature name.

measurementsDF <- DataFrame(measurements)
mcols(measurementsDF) <- data.frame(
  assay = rep(c("assay_1", "assay_2"), times = c(10, 1990)),
  feature = colnames(measurementsDF)
)

result <- crossValidate(measurementsDF, classes, classifier = "SVM", nFolds = 5,
                        nRepeats = 3, multiViewMethod = "merge")

performancePlot(result, characteristicsList = list(x = "Assay Name"))

## Warning in .local(results, ...): Balanced Accuracy not found in all elements of results. Calculating it now.

If using a list of data.frames, the name of each element in the list will be used as the assay name.

# Assigns first 10 variables to dataset_1, and the rest to dataset_2
measurementsList <- list(
  (measurements |> as.data.frame())[1:10],
  (measurements |> as.data.frame())[11:2000]
)
names(measurementsList) <- c("assay_1", "assay_2")

result <- crossValidate(measurementsList, classes, classifier = "SVM", nFolds = 5,
                        nRepeats = 3, multiViewMethod = "merge")

performancePlot(result, characteristicsList = list(x = "Assay Name"))

## Warning in .local(results, ...): Balanced Accuracy not found in all elements of results. Calculating it now.

Comparison to Existing Classification Frameworks

There are a few other frameworks for classification in R. The table below provides a comparison of which features they offer.

Provided Functionality

Although being a cross-validation framework, a number of popular feature selection and classification functions are provided by the package which meet the requirements of functions to be used by it (see the last section).

Fine-grained Cross-validation and Modelling Using runTests

For more control over the finer aspects of cross-validation of a single data set, runTests may be employed in place of crossValidate. For the variety of cross-validation, the parameters are specified by a CrossValParams object. The default setting is for 100 permutations and five folds and parameter tuning is done by resubstitution. It is also recommended to specify a parallelParams setting. On Linux and MacOS operating systems, it should be MulticoreParam and on Windows computers it should be SnowParam. Note that each of these have an option RNGseed and this needs to be set by the user because some classifiers or feature selection functions will have some element of randomisation. One example that works on all operating systems, but is best-suited to Windows is:

CVparams <- CrossValParams(parallelParams = SnowParam(16, RNGseed = 123))
CVparams

For the actual operations to do to the data to build a model of it, each of the stages should be specified by an object of class ModellingParams. This controls how class imbalance is handled (default is to downsample to the smallest class), any transformation that needs to be done inside of cross-validation (i.e. involving a computed value from the training set), any feature selection and the training and prediction functions to be used. The default is to do an ordinary t-test (two groups) or ANOVA (three or more groups) and classification using diagonal LDA.

ModellingParams()

## An object of class "ModellingParams"
## Slot "balancing":
## [1] "downsample"
## 
## Slot "transformParams":
## NULL
## 
## Slot "selectParams":
## An object of class 'SelectParams'.
## Selection Name: Difference in Means.
## 
## Slot "trainParams":
## An object of class 'TrainParams'.
## Classifier Name: Diagonal LDA.
## 
## Slot "predictParams":
## An object of class 'PredictParams'.
## 
## Slot "doImportance":
## [1] FALSE

runTests Driver Function of Cross-validated Classification

runTests is the main function in ClassifyR which handles the sample splitting and parallelisation, if used, of cross-validation. To begin with, a simple classifier will be demonstrated. It uses a t-test or ANOVA ranking (depending on the number of classes) for feature ranking and DLDA for classification. This classifier relies on differences in means between classes. No parameters need to be specified, because this is the default classification of runTests. By default, the number of features is tuned by resubstitution on the training set.

crossValParams <- CrossValParams(permutations = 5)
DMresults <- runTests(measurements, classes, crossValParams, verbose = 1)

Here, 5 permutations (non-default) and 5 folds cross-validation (default) is specified. For computers with more than 1 CPU, the number of cores to use can be given to runTests by using the argument parallelParams. The parameter seed is important to set for result reproducibility when doing a cross-validation such as this, because it employs randomisation to partition the samples into folds. Also, RNGseed is highly recommended to be set to the back-end specified to BPPARAM if doing parallel processing. The first seed mentioned does not work for parallel processes. For more details about runTests and the parameter classes used by it, consult the help pages of such functions.

Comparison of Different Classifications

The samplesMetricMap function allows the visual comparison of sample-wise error rate or accuracy measures from different ClassifyResult objects. Firstly, a classifier will be run that uses Kullback-Leibler divergence ranking and resubstitution error as a feature selection heuristic and a naive Bayes classifier for classification. This classification will use features that have either a change in location or in scale between classes.

modellingParamsDD <- ModellingParams(selectParams = SelectParams("KL"),
                                     trainParams = TrainParams("naiveBayes"),
                                     predictParams = NULL)
DDresults <- runTests(measurements, classes, crossValParams, modellingParamsDD, verbose = 1)
DDresults

## An object of class 'ClassifyResult'.
## Characteristics:
##    characteristic                       value
##    Selection Name Kullback-Leibler Divergence
##   Classifier Name          Naive Bayes Kernel
##  Cross-validation     5 Permutations, 5 Folds
## Features: List of length 25 of feature identifiers.
## Predictions: A data frame of 950 rows.
## Performance Measures: None calculated yet.

The naive Bayes kernel classifier by default uses the vertical distance between class densities but it can instead use the horizontal distance to the nearest non-zero density cross-over point to confidently classify samples in the tails of the densities.

The per-sample classification accuracy is automatically calculated for both the differential means and differential distribution classifiers and plotted with samplesMetricMap.

resultsList <- list(Abundance = DMresults, Distribution = DDresults)
samplesMetricMap(resultsList, showXtickLabels = FALSE)

## Warning in .local(results, ...): Sample Accuracy not found in all elements of results. Calculating it now.

## Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_tile()`).

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name                grob
## 1 1 (2-2,1-1) arrange      gtable[layout]
## 2 2 (1-1,1-1) arrange text[GRID.text.581]

The benefit of this plot is that it allows the easy identification of samples which are hard to classify and could be explained by considering additional information about them. Differential distribution class prediction appears to be biased to the majority class (No Asthma).

More traditionally, the distribution of performance values of each complete cross-validation can be visualised by performancePlot by providing them as a list to the function. The default is to draw box plots, but violin plots could also be made. The default performance metric to plot is balanced accuracy. If it’s not already calculated for all classifications, as in this case for DD, it will be done automatically.

performancePlot(resultsList)

## Warning in .local(results, ...): Balanced Accuracy not found in all elements of results. Calculating it now.

We can observe that the spread of balanced accuracy rates is small, but slightly wider for the differential distribution classifier.

The features being ranked and selected in the feature selection stage can be compared within and between classifiers by the plotting functions rankingPlot and selectionPlot. Consider the task of visually representing how consistent the feature rankings of the top 100 different features were for the differential distribution classifier for all 5 folds in the 5 cross-validations.

rankingPlot(DDresults, topRanked = 1:100, xLabelPositions = c(1, seq(10, 100, 10)))

The top-ranked features are fairly similar between all pairs of the 20 cross-validations.

For a large cross-validation scheme, such as leave-2-out cross-validation, or when results contains many classifications, there are many feature set comparisons to make. Note that rankingPlot and selectionPlot have a parallelParams options which allows for the calculation of feature set overlaps to be done on multiple processors.

Generating a ROC Plot

Some classifiers can output scores or probabilities representing how likely a sample is to be from one of the classes, instead of, or as well as, class labels. This enables different score thresholds to be tried, to generate pairs of false positive and false negative rates. The naive Bayes classifier used previously by default has its returnType parameter set to “both”, so class predictions and scores are both stored in the classification result. So does diagonal LDA. In this case, a data frame with class predictions and scores for each class is returned by the classifier to the cross-validation framework. Setting returnType to “score” for a classifier which has such an option is also sufficient to generate a ROC plot. Many existing classifiers in other R packages also have an option that allows a score or probability to be calculated.

By default, scores from different iterations of prediction are merged and one line is drawn per classification. Alternatively, setting mode = “average” will consider each iteration of prediction separately, average them and also calculate and draw confidence intervals. The default interval is a 95% interval and is customisable by setting interval.

ROCplot(resultsList, fontSizes = c(24, 12, 12, 12, 12))

This ROC plot shows the classifiability of the asthma data set is high. Some examples of functions which output scores are fisherDiscriminant, DLDApredictInterface, and SVMpredictInterface.

Using an Independent Test Set

Sometimes, cross-validation is unnecessary. This happens when studies have large sample sizes and are designed such that a large number of samples is prespecified to form a test set. The classifier is only trained on the training sample set, and makes predictions only on the test sample set. This can be achieved by using the function runTest directly. See its documentation for required inputs.

Cross-validating Selected Features on a Different Data Set

Once a cross-validated classification is complete, the usefulness of the features selected may be explored in another dataset. previousSelection is a function which takes an existing ClassifyResult object and returns the features selected at the equivalent iteration which is currently being processed. This is necessary, because the models trained on one data set are not directly transferrable to a new dataset; the classifier training (e.g. choosing thresholds, fitting model coefficients) is redone. Of course, the features in the new dataset should have the same naming system as the ones in the old dataset.

Parameter Tuning

Some feature ranking methods or classifiers allow the choosing of tuning parameters, which controls some aspect of their model learning. An example of doing parameter tuning with a linear SVM is presented. This particular SVM has a single tuning parameter, the cost. Higher values of this parameter penalise misclassifications more. Moreover, feature selection happens by using a feature ranking function and then trying a range of top-ranked features to see which gives the best performance, the range being specified by a list element named nFeatures and the performance type (e.g. Balanced Accuracy) specified by a list element named performanceType. Therefore, some kind of parameter tuning always happens, even if the feature ranking or classifier function does not have any explicit tuning parameters.

Tuning is achieved in ClassifyR by providing a variable called tuneParams to the SelectParams or TrainParams constructor. tuneParams is a named list, with the names being the names of the tuning variables, except for one which is named “performanceType” and specifies the performance metric to use for picking the parameter values. Any of the non-sample-specific performance metrics which calcCVperformance calculates can be optimised.

tuneList <- list(cost = c(0.01, 0.1, 1, 10))
SVMparams <- ModellingParams(trainParams = TrainParams("SVM", kernel = "linear", tuneParams = tuneList),
                             predictParams = PredictParams("SVM"))
SVMresults <- runTests(measurements, classes, crossValParams, SVMparams)

The index of chosen of the parameters, as well as all combinations of parameters and their associated performance metric, are stored for every validation, and can be accessed with the tunedParameters function.

length(tunedParameters(SVMresults))

## [1] 25

tunedParameters(SVMresults)[1:5]

## [[1]]
## [[1]]$tuneCombinations
##    topN  cost Balanced Accuracy
## 1    10  0.01         0.8022680
## 2    20  0.01         0.8167524
## 3    30  0.01         0.8268534
## 4    40  0.01         0.8268534
## 5    50  0.01         0.8268534
## 6    60  0.01         0.8362874
## 7    70  0.01         0.8413379
## 8    80  0.01         0.8369545
## 9    90  0.01         0.8463884
## 10  100  0.01         0.8369545
## 11   10  0.10         0.7935010
## 12   20  0.10         0.8514389
## 13   30  0.10         0.8746903
## 14   40  0.10         0.8847913
## 15   50  0.10         0.8905089
## 16   60  0.10         0.9049933
## 17   70  0.10         0.9093768
## 18   80  0.10         0.9231942
## 19   90  0.10         0.9332952
## 20  100  0.10         0.9194778
## 21   10  1.00         0.7796836
## 22   20  1.00         0.8709739
## 23   30  1.00         0.8847913
## 24   40  1.00         0.8861254
## 25   50  1.00         0.8999428
## 26   60  1.00         0.9383457
## 27   70  1.00         0.9760816
## 28   80  1.00         0.9811321
## 29   90  1.00         0.9811321
## 30  100  1.00         0.9811321
## 31   10 10.00         0.7753002
## 32   20 10.00         0.8766914
## 33   30 10.00         0.9049933
## 34   40 10.00         0.9238613
## 35   50 10.00         0.8396226
## 36   60 10.00         1.0000000
## 37   70 10.00         1.0000000
## 38   80 10.00         1.0000000
## 39   90 10.00         1.0000000
## 40  100 10.00         1.0000000
## 
## [[1]]$bestIndex
## [1] 36
## 
## 
## [[2]]
## [[2]]$tuneCombinations
##    topN  cost Balanced Accuracy
## 1    10  0.01         0.8022680
## 2    20  0.01         0.8016009
## 3    30  0.01         0.8160854
## 4    40  0.01         0.8312369
## 5    50  0.01         0.8211359
## 6    60  0.01         0.8312369
## 7    70  0.01         0.8218029
## 8    80  0.01         0.8312369
## 9    90  0.01         0.8645893
## 10  100  0.01         0.8645893
## 11   10  0.10         0.7891176
## 12   20  0.10         0.8369545
## 13   30  0.10         0.8514389
## 14   40  0.10         0.9036592
## 15   50  0.10         0.9036592
## 16   60  0.10         0.9238613
## 17   70  0.10         0.9087097
## 18   80  0.10         0.9188107
## 19   90  0.10         0.9477797
## 20  100  0.10         0.9572136
## 21   10  1.00         0.8413379
## 22   20  1.00         0.8652563
## 23   30  1.00         0.8716409
## 24   40  1.00         0.9427292
## 25   50  1.00         0.9383457
## 26   60  1.00         0.9427292
## 27   70  1.00         0.9383457
## 28   80  1.00         0.9666476
## 29   90  1.00         0.9905660
## 30  100  1.00         1.0000000
## 31   10 10.00         0.8652563
## 32   20 10.00         0.8659234
## 33   30 10.00         0.9376787
## 34   40 10.00         0.9804650
## 35   50 10.00         1.0000000
## 36   60 10.00         1.0000000
## 37   70 10.00         1.0000000
## 38   80 10.00         1.0000000
## 39   90 10.00         1.0000000
## 40  100 10.00         1.0000000
## 
## [[2]]$bestIndex
## [1] 30
## 
## 
## [[3]]
## [[3]]$tuneCombinations
##    topN  cost Balanced Accuracy
## 1    10  0.01         0.8261864
## 2    20  0.01         0.8167524
## 3    30  0.01         0.8073185
## 4    40  0.01         0.8079855
## 5    50  0.01         0.8551553
## 6    60  0.01         0.8645893
## 7    70  0.01         0.8551553
## 8    80  0.01         0.8696398
## 9    90  0.01         0.8797408
## 10  100  0.01         0.8847913
## 11   10  0.10         0.8703068
## 12   20  0.10         0.8558224
## 13   30  0.10         0.8804079
## 14   40  0.10         0.8753573
## 15   50  0.10         0.8847913
## 16   60  0.10         0.8891748
## 17   70  0.10         0.9043263
## 18   80  0.10         0.9231942
## 19   90  0.10         0.9332952
## 20  100  0.10         0.9231942
## 21   10  1.00         0.8615399
## 22   20  1.00         0.8841243
## 23   30  1.00         0.8898418
## 24   40  1.00         0.9188107
## 25   50  1.00         0.9420621
## 26   60  1.00         0.9477797
## 27   70  1.00         0.9666476
## 28   80  1.00         0.9905660
## 29   90  1.00         1.0000000
## 30  100  1.00         1.0000000
## 31   10 10.00         0.8804079
## 32   20 10.00         0.8709739
## 33   30 10.00         0.8955594
## 34   40 10.00         0.9716981
## 35   50 10.00         0.9905660
## 36   60 10.00         1.0000000
## 37   70 10.00         1.0000000
## 38   80 10.00         1.0000000
## 39   90 10.00         1.0000000
## 40  100 10.00         1.0000000
## 
## [[3]]$bestIndex
## [1] 29
## 
## 
## [[4]]
## [[4]]$tuneCombinations
##    topN  cost Balanced Accuracy
## 1    10  0.01         0.8305698
## 2    20  0.01         0.8261864
## 3    30  0.01         0.8740233
## 4    40  0.01         0.8979417
## 5    50  0.01         0.8979417
## 6    60  0.01         0.8885077
## 7    70  0.01         0.8986087
## 8    80  0.01         0.8986087
## 9    90  0.01         0.9225272
## 10  100  0.01         0.9130932
## 11   10  0.10         0.8827902
## 12   20  0.10         0.8878407
## 13   30  0.10         0.9269106
## 14   40  0.10         0.9124261
## 15   50  0.10         0.9124261
## 16   60  0.10         0.9130932
## 17   70  0.10         0.9326282
## 18   80  0.10         0.9332952
## 19   90  0.10         0.9332952
## 20  100  0.10         0.9282447
## 21   10  1.00         0.8935582
## 22   20  1.00         0.8935582
## 23   30  1.00         0.9457785
## 24   40  1.00         0.9370116
## 25   50  1.00         0.9464456
## 26   60  1.00         0.9609301
## 27   70  1.00         0.9659806
## 28   80  1.00         0.9905660
## 29   90  1.00         1.0000000
## 30  100  1.00         1.0000000
## 31   10 10.00         0.8942253
## 32   20 10.00         0.9319611
## 33   30 10.00         0.9797980
## 34   40 10.00         0.9848485
## 35   50 10.00         1.0000000
## 36   60 10.00         1.0000000
## 37   70 10.00         1.0000000
## 38   80 10.00         1.0000000
## 39   90 10.00         1.0000000
## 40  100 10.00         1.0000000
## 
## [[4]]$bestIndex
## [1] 29
## 
## 
## [[5]]
## [[5]]$tuneCombinations
##    topN  cost Balanced Accuracy
## 1    10  0.01         0.8153846
## 2    20  0.01         0.8107692
## 3    30  0.01         0.8542308
## 4    40  0.01         0.8738462
## 5    50  0.01         0.8784615
## 6    60  0.01         0.8784615
## 7    70  0.01         0.8930769
## 8    80  0.01         0.9126923
## 9    90  0.01         0.9223077
## 10  100  0.01         0.9223077
## 11   10  0.10         0.8446154
## 12   20  0.10         0.8880769
## 13   30  0.10         0.9126923
## 14   40  0.10         0.9273077
## 15   50  0.10         0.9369231
## 16   60  0.10         0.9469231
## 17   70  0.10         0.9419231
## 18   80  0.10         0.9419231
## 19   90  0.10         0.9419231
## 20  100  0.10         0.9615385
## 21   10  1.00         0.8734615
## 22   20  1.00         0.8934615
## 23   30  1.00         0.9657692
## 24   40  1.00         0.9807692
## 25   50  1.00         0.9807692
## 26   60  1.00         0.9903846
## 27   70  1.00         0.9903846
## 28   80  1.00         1.0000000
## 29   90  1.00         1.0000000
## 30  100  1.00         1.0000000
## 31   10 10.00         0.8926923
## 32   20 10.00         0.9273077
## 33   30 10.00         0.9803846
## 34   40 10.00         1.0000000
## 35   50 10.00         1.0000000
## 36   60 10.00         1.0000000
## 37   70 10.00         1.0000000
## 38   80 10.00         1.0000000
## 39   90 10.00         1.0000000
## 40  100 10.00         1.0000000
## 
## [[5]]$bestIndex
## [1] 28

The cost value of 1 or 10 appears to often be chosen.