@@ -3,8 +3,8 @@ Type: Package

 Title: A framework for cross-validated classification problems, with

        applications to differential variability and differential

        distribution testing

-Version: 3.3.13

-Date: 2023-02-13

+Version: 3.3.14

+Date: 2023-02-16

 Authors@R:

     c(

     person(given = "Dario", family = "Strbenac", email = "dario.strbenac@sydney.edu.au", role = c("aut", "cre")),

@@ -72,6 +72,7 @@ Collate:

     'interfaceXGB.R'

     'performancePlot.R'

     'plotFeatureClasses.R'

+    'precisionPathways.R'

     'prepareData.R'

     'previousSelection.R'

     'previousTrained.R'

@@ -76,6 +76,8 @@ exportMethods(models)

 exportMethods(performance)

 exportMethods(performancePlot)

 exportMethods(plotFeatureClasses)

+exportMethods(precisionPathwaysPredict)

+exportMethods(precisionPathwaysTrain)

 exportMethods(predictions)

 exportMethods(prepareData)

 exportMethods(rankingPlot)

@@ -35,6 +35,9 @@ setClassUnion("DataFrameOrNULL", c("DataFrame", "NULL"))

 setClassUnion("tabular", c("data.frame", "DataFrame", "matrix"))

 setClassUnion("tabularOrList", c("tabular", "list"))

+# List-like assay data

+setClassUnion("MultiAssayExperimentOrList", c("MultiAssayExperiment", "list"))

+

 ################################################################################

 #

 # Params

@@ -22,14 +22,18 @@

 #' @param extraParams A list of parameters that will be used to overwrite default settings of transformation, selection, or model-building functions or

 #' parameters which will be passed into the data cleaning function. The names of the list must be one of \code{"prepare"},

 #' \code{"select"}, \code{"train"}, \code{"predict"}. To remove one of the defaults (see the article titled Parameter Tuning Presets for crossValidate and Their Customisation on

-#' the website), specify the list element to be \code{NULL}.

+#' the website), specify the list element to be \code{NULL}. For the valid element names in the \code{"prepare"} list, see \code{?prepareData}.

+#' @param clinicalPredictors Default: \code{NULL}. If \code{measurements} is a \code{MultiAssayExperiment},

+#' a character vector of features to use in modelling. This allows avoidance of things like sample IDs,

+#' sample acquisition dates, etc. which are not relevant for outcome prediction.

 #' @param nFeatures The number of features to be used for classification. If this is a single number, the same number of features will be used for all comparisons

 #' or assays. If a numeric vector these will be optimised over using \code{selectionOptimisation}. If a named vector with the same names of multiple assays, 

 #' a different number of features will be used for each assay. If a named list of vectors, the respective number of features will be optimised over. 

 #' Set to NULL or "all" if all features should be used.

 #' @param selectionMethod Default: \code{"auto"}. A character vector of feature selection methods to compare. If a named character vector with names corresponding to different assays, 

 #' and performing multiview classification, the respective selection methods will be used on each assay. If \code{"auto"}, t-test (two categories) / F-test (three or more categories) ranking

-#' and top \code{nFeatures} optimisation is done. Otherwise, the ranking method is per-feature Cox proportional hazards p-value.

+#' and top \code{nFeatures} optimisation is done. Otherwise, the ranking method is per-feature Cox proportional hazards p-value. \code{"none"} is also a valid value, meaning that no

+#' indepedent feature selection will be performed (but implicit selection might still happen with the classifier).

 #' @param selectionOptimisation A character of "Resubstitution", "Nested CV" or "none" specifying the approach used to optimise \code{nFeatures}.

 #' @param performanceType Default: \code{"auto"}. If \code{"auto"}, then balanced accuracy for classification or C-index for survival. Otherwise, any one of the

 #' options described in \code{\link{calcPerformance}} may otherwise be specified.

@@ -88,12 +92,13 @@

 #' # performancePlot(c(result, resultMerge))

 #' 

 #' @importFrom survival Surv

+#' @usage NULL

 setGeneric("crossValidate", function(measurements, outcome, ...)

     standardGeneric("crossValidate"))

 #' @rdname crossValidate

 #' @export

-setMethod("crossValidate", "DataFrame", 

+setMethod("crossValidate", "DataFrame",

           function(measurements,

                    outcome,

                    nFeatures = 20,

@@ -110,12 +115,15 @@ setMethod("crossValidate", "DataFrame",

           {

               # Check that data is in the right format, if not already done for MultiAssayExperiment input.

-              prepParams <- list(measurements, outcome)

-              if("prepare" %in% names(extraParams))

-                prepParams <- c(prepParams, extraParams[["prepare"]])

-              measurementsAndOutcome <- do.call(prepareData, prepParams)

-              measurements <- measurementsAndOutcome[["measurements"]]

-              outcome <- measurementsAndOutcome[["outcome"]]

+              if(!"assay" %in% colnames(S4Vectors::mcols(measurements))) # Assay is put there by prepareData for MultiAssayExperiment, skip if present. 

+              {

+                prepParams <- list(measurements, outcome)

+                if("prepare" %in% names(extraParams))

+                  prepParams <- c(prepParams, extraParams[["prepare"]])

+                measurementsAndOutcome <- do.call(prepareData, prepParams)

+                measurements <- measurementsAndOutcome[["measurements"]]

+                outcome <- measurementsAndOutcome[["outcome"]]

+              }

               # Ensure performance type is one of the ones that can be calculated by the package.

               if(!performanceType %in% c("auto", .ClassifyRenvir[["performanceTypes"]]))

@@ -305,9 +313,10 @@ setMethod("crossValidate", "DataFrame",

 #' @rdname crossValidate

 #' @export

 # One or more omics data sets, possibly with clinical data.

-setMethod("crossValidate", "MultiAssayExperiment",

+setMethod("crossValidate", "MultiAssayExperimentOrList",

           function(measurements,

-                   outcome, 

+                   outcome,

+                   clinicalPredictors = NULL,

                    nFeatures = 20,

                    selectionMethod = "auto",

                    selectionOptimisation = "Resubstitution",

@@ -321,7 +330,7 @@ setMethod("crossValidate", "MultiAssayExperiment",

                    characteristicsLabel = NULL, extraParams = NULL)

           {

               # Check that data is in the right format, if not already done for MultiAssayExperiment input.

-              prepParams <- list(measurements, outcome)

+              prepParams <- list(measurements, outcome, clinicalPredictors)

               if("prepare" %in% names(extraParams))

                 prepParams <- c(prepParams, extraParams[["prepare"]])

               measurementsAndOutcome <- do.call(prepareData, prepParams)

@@ -408,76 +417,6 @@ setMethod("crossValidate", "matrix", # Matrix of numeric measurements.

                             characteristicsLabel = characteristicsLabel, extraParams = extraParams)

           })

-# This expects that each table is about the same set of samples and thus

-# has the same number of rows as every other table.

-#' @rdname crossValidate                                                            

-#' @export

-setMethod("crossValidate", "list",

-          function(measurements,

-                   outcome, 

-                   nFeatures = 20,

-                   selectionMethod = "auto",

-                   selectionOptimisation = "Resubstitution",

-                   performanceType = "auto",

-                   classifier = "auto",

-                   multiViewMethod = "none",

-                   assayCombinations = "all",

-                   nFolds = 5,

-                   nRepeats = 20,

-                   nCores = 1,

-                   characteristicsLabel = NULL, extraParams = NULL)

-          {

-              # Check data type is valid

-              if (!(all(sapply(measurements, class) %in% c("data.frame", "DataFrame", "matrix")))) {

-                  stop("assays must be of type data.frame, DataFrame or matrix")

-              }

-              

-              # Check the list is named

-              if (is.null(names(measurements))) {

-                  stop("Measurements must be a named list.")

-              }

-              

-              # Check same number of samples for all datasets

-              if (!length(unique(sapply(measurements, nrow))) == 1) {

-                  stop("All datasets must have the same samples.")

-              }

-              

-              # Check the number of outcome is the same

-              if (!all(sapply(measurements, nrow) == length(outcome)) && !is.character(outcome)) {

-                  stop("outcome must have same number of samples as measurements.")

-              }

-              

-              df_list <- sapply(measurements, S4Vectors::DataFrame, check.names = FALSE)

-              

-              df_list <- mapply(function(meas, nam){

-                  S4Vectors::mcols(meas)$assay <- nam

-                  S4Vectors::mcols(meas)$feature <- colnames(meas)

-                  meas

-              }, df_list, names(df_list))

-              

-              

-              combined_df <- do.call("cbind", df_list) 

-              colnames(combined_df) <- S4Vectors::mcols(combined_df)$feature

-

-

-              

-              crossValidate(measurements = combined_df,

-                            outcome = outcome, 

-                            nFeatures = nFeatures,

-                            selectionMethod = selectionMethod,

-                            selectionOptimisation = selectionOptimisation,

-                            performanceType = performanceType,

-                            classifier = classifier,

-                            multiViewMethod = multiViewMethod,

-                            assayCombinations = assayCombinations,

-                            nFolds = nFolds,

-                            nRepeats = nRepeats,

-                            nCores = nCores,

-                            characteristicsLabel = characteristicsLabel, extraParams = extraParams)

-          })

-

-

-

 ######################################

 ######################################

 cleanNFeatures <- function(nFeatures, measurements){

@@ -656,7 +595,6 @@ generateModellingParams <- function(assayIDs,

     }    

     selectionMethod <- unlist(selectionMethod)

-    if(is.null(selectionMethod)) selectionMethod <- "none"

     if(selectionMethod != "none")

     {

@@ -1126,9 +1064,9 @@ train.list <- function(x, outcomeTrain, ...)

 #' @rdname crossValidate

 #' @method train MultiAssayExperiment

 #' @export

-train.MultiAssayExperiment <- function(x, outcome, ...)

+train.MultiAssayExperiment <- function(x, outcome, clinicalPredictors = NULL, ...)

           {

-              prepArgs <- list(x, outcome)

+              prepArgs <- list(x, outcome, clinicalPredictors)

               extraInputs <- list(...)

               prepExtras <- trainExtras <- numeric()

               if(length(extraInputs) > 0)

@@ -1167,7 +1105,7 @@ predict.trainedByClassifyR <- function(object, newData, ...)

     newData <- do.call(cbind, newData)

     } else if(is(newData, "MultiAssayExperiment"))

             {

-              newData <- prepareData(newData, useFeatures = allFeatureNames(object))

+              newData <- prepareData(newData, clinicalPredictors = subset(allFeatureNames(object), assay == "clinical")[, "feature"])

               # Some classifiers dangerously use positional matching rather than column name matching.

               # newData columns are sorted so that the right column ordering is guaranteed.

     }

@@ -6,20 +6,26 @@ elasticNetGLMtrainInterface <- function(measurementsTrain, classesTrain, lambda

     stop("The package 'glmnet' could not be found. Please install it.")

   if(verbose == 3)

     message("Fitting elastic net regularised GLM classifier to data.")

-    

-  measurementsMatrix <- glmnet::makeX(as(measurementsTrain, "data.frame"))

-

-  fitted <- glmnet::glmnet(measurementsMatrix, classesTrain, family = "multinomial", ...)

+  # One-hot encoding needed.    

+  measurementsTrain <- MatrixModels::model.Matrix(~ 0 + ., data = measurementsTrain)

+  fitted <- glmnet::glmnet(measurementsTrain, classesTrain, family = "multinomial", weights = as.numeric(1 / (table(classesTrain)[classesTrain] / length(classesTrain))), ...)

+  # Inverse class size weighting needed to give decent predictions when class imbalance.

+  

   if(is.null(lambda)) # fitted has numerous models for automatically chosen lambda values.

-  { # Pick one lambda based on resubstitution performance.

-    bestLambda <- fitted[["lambda"]][which.min(sapply(fitted[["lambda"]], function(lambda) # Largest Lambda with minimum balanced error rate.

+  { # Pick one lambda based on resubstitution performance. But not the one that makes all variables excluded from model.

+    lambdaConsider <- colSums(as.matrix(fitted[["beta"]][[1]])) != 0

+    bestLambda <- fitted[["lambda"]][lambdaConsider][which.min(sapply(fitted[["lambda"]][lambdaConsider], function(lambda) # Largest Lambda with minimum balanced error rate.

     {

-      classPredictions <- factor(as.character(predict(fitted, measurementsMatrix, s = lambda, type = "class")), levels = fitted[["classnames"]])

+      classPredictions <- factor(as.character(predict(fitted, measurementsTrain, s = lambda, type = "class")), levels = fitted[["classnames"]])

       calcExternalPerformance(classesTrain, classPredictions, "Balanced Error")

     }))[1]]

     attr(fitted, "tune") <- list(lambda = bestLambda)

   }

+  

+  attr(fitted, "featureNames") <- colnames(measurementsTrain)

+  attr(fitted, "featureGroups") <- measurementsTrain@assign

+  

   fitted

 }

 attr(elasticNetGLMtrainInterface, "name") <- "elasticNetGLMtrainInterface"

@@ -28,6 +34,8 @@ attr(elasticNetGLMtrainInterface, "name") <- "elasticNetGLMtrainInterface"

 elasticNetGLMpredictInterface <- function(model, measurementsTest, lambda, ..., returnType = c("both", "class", "score"), verbose = 3)

 { # ... just consumes emitted tuning variables from .doTrain which are unused.

   returnType <- match.arg(returnType)

+  # One-hot encoding needed.    

+  measurementsTest <- MatrixModels::model.Matrix(~ 0 + ., data = measurementsTest)

   # Ensure that testing data has same columns names in same order as training data.

   # Remove those annoying backquotes which glmnet adds if variables have spaces in names.

@@ -41,11 +49,11 @@ elasticNetGLMpredictInterface <- function(model, measurementsTest, lambda, ...,

   if(missing(lambda)) # Tuning parameters are not passed to prediction functions.

     lambda <- attr(model, "tune")[["lambda"]] # Sneak it in as an attribute on the model.

-  testMatrix <- glmnet::makeX(as(measurementsTest, "data.frame"))

-  testMatrix <- testMatrix[, rownames(model[["beta"]][[1]])]

-  classPredictions <- factor(as.character(predict(model, testMatrix, s = lambda, type = "class")), levels = model[["classnames"]])

-  classScores <- predict(model, testMatrix, s = lambda, type = "response")[, , 1]

+  measurementsTest <- measurementsTest[, rownames(model[["beta"]][[1]])]

+  

+  classPredictions <- factor(as.character(predict(model, measurementsTest, s = lambda, type = "class")), levels = model[["classnames"]])

+  classScores <- predict(model, measurementsTest, s = lambda, type = "response")[, , 1]

   if(is.matrix(classScores))

     classScores <- classScores[, model[["classnames"]]]

@@ -61,16 +69,20 @@ elasticNetGLMpredictInterface <- function(model, measurementsTest, lambda, ...,

 #

 # Get selected features (i.e. non-zero model coefficients)

 #

+# Note: Need to convert back to actual features when factors were expanded into

+# multiple columns using indicator variables.

+#

 ################################################################################

 elasticNetFeatures <- function(model)

                       {

-                        inputFeatures <- rownames(model[["beta"]][[1]])            

                         # Floating point numbers test for equality.

                         whichCoefficientColumn <- which(abs(model[["lambda"]] - attr(model, "tune")[["lambda"]]) < 0.00001)[1]

                         coefficientsUsed <- sapply(model[["beta"]], function(classCoefficients) classCoefficients[, whichCoefficientColumn])

                         featureScores <- rowSums(abs(coefficientsUsed))

-                        rankedFeaturesIndices <- order(featureScores, decreasing = TRUE)

-                        selectedFeaturesIndices <- which(featureScores != 0)

+                        featureGroups <- attr(model, "featureGroups")[match(names(featureScores), attr(model, "featureNames"))]

+                        groupScores <- unname(by(featureScores, featureGroups, max))

+                        rankedFeaturesIndices <- order(groupScores, decreasing = TRUE)

+                        selectedFeaturesIndices <- which(groupScores != 0)

                         list(rankedFeaturesIndices, selectedFeaturesIndices)

                       }

\ No newline at end of file

@@ -9,7 +9,7 @@ GLMtrainInterface <- function(measurementsTrain, classesTrain, ..., verbose = 3)

         fitData <- cbind(measurementsTrain, class = classesTrain)

         classesTrain <- "class" # Column name for glm fit.

     } else {fitData <- measurementsTrain}

-  glm(class ~ . + 0, family = binomial, data = fitData, ...)

+  glm(class ~ . + 0, family = binomial, data = fitData, weights = as.numeric(1 / (table(classesTrain)[classesTrain] / length(classesTrain))), ...)

 }

 attr(GLMtrainInterface, "name") <- "GLMtrainInterface"

new file mode 100644

@@ -0,0 +1,180 @@

+#' Precision Pathways for Sample Prediction Based on Prediction Confidence.

+#' 

+#' Precision pathways allows the evaluation of various permutations of multiomics or multiview data.

+#' Samples are predicted by a particular assay if they were consistently predicted as a particular class

+#' during cross-validation. Otherwise, they are passed onto subsequent assays/tiers for prediction. Balanced accuracy

+#' is used to evaluate overall prediction performance and sample-specific accuracy for individual-level evaluation.

+#'

+#' @param measurements Either a \code{\link{MultiAssayExperiment}} or a list of the basic tabular objects containing the data.

+#' @param class Same as \code{measurements} but only training samples. IF \code{measurements} is a \code{list}, may also be

+#' a vector of classes.

+#' @param clinicalPredictors Default: \code{NULL}. Must be a character vector of clinical features to use in modelling. This allows avoidance of things like sample IDs,

+#' sample acquisition dates, etc. which are not relevant for outcome prediction.

+#' @param maxMissingProp Default: 0.0. A proportion less than 1 which is the maximum

+#' tolerated proportion of missingness for a feature to be retained for modelling.

+#' @param topNvariance Default: NULL. An integer number of most variable features per assay to subset to.

+#' Assays with less features won't be reduced in size.

+#' @param fixedAssays A character vector of assay names specifying any assays which must be at the

+#' beginning of the pathway.

+#' @param confidenceCutoff The minimum confidence of predictions for a sample to be predicted by a particular issue

+#' . If a sample was predicted to belong to a particular class a proportion \eqn{p} times, then the confidence is \eqn{2 \times |p - 0.5|}.

+#' @param minAssaySamples An integer specifying the minimum number of samples a tier may have. If a subsequent tier

+#' would have less than this number of samples, the samples are incorporated into the current tier.

+#' @param nFeatures Default: 20. The number of features to consider during feature selection, if feature selection is done.

+#' @param selectionMethod A named character vector of feature selection methods to use for the assays, one for each. The names must correspond to names of \code{measurements}.

+#' @param classifier A named character vector of modelling methods to use for the assays, one for each. The names must correspond to names of \code{measurements}.

+#' @param nFolds A numeric specifying the number of folds to use for cross-validation.

+#' @param nRepeats A numeric specifying the the number of repeats or permutations to use for cross-validation.

+#' @param nCores A numeric specifying the number of cores used if the user wants to use parallelisation.

+#' @param pathways A set of pathways created by \code{precisionPathwaysTrain} to be used for predicting on a new data set.

+#' @rdname precisionPathways

+#' @return An object of class \code{PrecisionPathways} which is basically a named list that other plotting and

+#' tabulating functions can use.

+#' @examples

+#' # To be determined.

+

+#' @usage NULL

+setGeneric("precisionPathwaysTrain", function(measurements, class, ...)

+    standardGeneric("precisionPathwaysTrain"))

+

+#' @rdname precisionPathways

+#' @export

+setMethod("precisionPathwaysTrain", "MultiAssayExperimentOrList", 

+          function(measurements, class, clinicalPredictors = NULL, maxMissingProp = 0.0, topNvariance = NULL,

+                   fixedAssays = "clinical", confidenceCutoff = 0.8, minAssaySamples = 10,

+                   nFeatures = 20, selectionMethod = setNames(c("none", rep("t-test", length(measurements))), c("clinical", names(measurements))),

+                   classifier = setNames(c("elasticNetGLM", rep("randomForest", length(measurements))), c("clinical", names(measurements))),

+                   nFolds = 5, nRepeats = 20, nCores = 1)

+          {

+            if(is.list(measurements)) # Ensure plain list has clinical data.

+            {

+              # One of the tables must be named "clinical".

+              if (!any(names(measurements) == "clinical"))

+                stop("One of the tables must be named \"clinical\".")

+            }

+            prepArgs <- list(measurements, outcomeColumns = class, clinicalPredictors = clinicalPredictors,

+                             maxMissingProp = maxMissingProp, topNvariance = topNvariance)

+            measurementsAndClass <- do.call(prepareData, prepArgs)

+              

+            .precisionPathwaysTrain(measurementsAndClass[["measurements"]], measurementsAndClass[["outcome"]],

+                                   fixedAssays = fixedAssays, confidenceCutoff = confidenceCutoff,

+                                   minAssaySamples = minAssaySamples, nFeatures = nFeatures,

+                                   selectionMethod = selectionMethod, classifier = classifier,

+                                   nFolds = nFolds, nRepeats = nRepeats, nCores = nCores)

+          })

+

+# Internal method which carries out all of the processing, obtaining reformatted data from the

+# MultiAssayExperiment and list (of basic rectangular tables) S4 methods.

+.precisionPathwaysTrain <- function(measurements, class, fixedAssays = "clinical",

+                   confidenceCutoff = 0.8, minAssaySamples = 10,

+                   nFeatures = 20, selectionMethod = setNames(c(NULL, rep("t-test", length(measurements))), c("clinical", names(measurements))),

+                   classifier = setNames(c("elasticNetGLM", rep("randomForest", length(measurements))), c("clinical", names(measurements))),

+                   nFolds = 5, nRepeats = 20, nCores = 1)

+          {

+            # Step 1: Determine all valid permutations of assays, taking into account the

+            # assays to be used and which assays, if any, must be included.

+            assayIDs <- unique(S4Vectors::mcols(measurements)[["assay"]])

+            assaysPermutations <- .permutations(assayIDs, fixed = data.frame(seq_along(fixedAssays), fixedAssays))

+            permutationIDs <- apply(assaysPermutations, 2, function(permutation) paste(permutation, collapse = '-'))

+            

+            # Step 2: Build a classifier for each assay using all of the samples.

+            modelsList <- crossValidate(measurements, class, nFeatures, selectionMethod,

+                                        classifier = classifier, nFolds = nFolds,

+                                        nRepeats = nRepeats, nCores = nCores)

+            modelsList <- lapply(modelsList, calcCVperformance, "Sample Accuracy") # Add sample accuracy, which can be subset later.

+

+            # Step 3: Loop over each pathway and each assay in order to determine which samples are used at that level

+            # and which are passed onwards.

+            precisionPathways <- lapply(as.data.frame(assaysPermutations), function(permutation)

+            {

+              assaysProcessed <- character()

+              samplesUsed <- character()

+              individualsTableAll <- S4Vectors::DataFrame()

+              tierTableAll <- S4Vectors::DataFrame()

+              breakEarly = FALSE

+              for(assay in permutation)

+              {

+                # Step 3a: Identify all samples which are consistently predicted.

+                modelIndex <- match(assay, assayIDs)

+                allPredictions <- predictions(modelsList[[modelIndex]])

+                allSampleIDs <- sampleNames(modelsList[[modelIndex]])

+                predictionsSamplesCounts <- table(allPredictions[, "sample"], allPredictions[, "class"])

+                confidences <- 2 * abs(predictionsSamplesCounts[, 1] / rowSums(predictionsSamplesCounts) - 0.5)

+                sampleIDsUse <- names(confidences)[confidences > confidenceCutoff]

+                

+                # Check if too few samples left for next round. Include them in this round, if so.

+                remainingIDs <- setdiff(allSampleIDs, c(samplesUsed, sampleIDsUse))

+                if(length(remainingIDs) < minAssaySamples)

+                {

+                  sampleIDsUse <- c(sampleIDsUse, remainingIDs)

+                  breakEarly = TRUE

+                }

+                

+                predictionsSamplesCounts <- predictionsSamplesCounts[sampleIDsUse, ]

+                

+                # Step 3b: Individuals predictions and sample-wise accuracy, tier-wise error.

+                maxVotes <- apply(predictionsSamplesCounts, 1, function(sample) which.max(sample))

+                predictedClasses <- factor(colnames(predictionsSamplesCounts)[maxVotes],

+                                           levels = colnames(predictionsSamplesCounts))    

+                individualsTable <- S4Vectors::DataFrame(Tier = assay,

+                                                         `Sample ID` = sampleIDsUse,

+                                                         `Predicted` = predictedClasses,

+                                                         `Accuracy` = performance(modelsList[[modelIndex]])[["Sample Accuracy"]][sampleIDsUse],

+                                                         check.names = FALSE)

+                knownClasses <- actualOutcome(modelsList[[modelIndex]])[match(sampleIDsUse, allSampleIDs)]

+                balancedAccuracy <- calcExternalPerformance(knownClasses, predictedClasses)

+                tierTable <- S4Vectors::DataFrame(Tier = assay,

+                                                  `Balanced Accuracy` = balancedAccuracy, check.names = FALSE)

+                

+                assaysProcessed <- c(assaysProcessed, assay)

+                individualsTableAll <- rbind(individualsTableAll, individualsTable)

+                tierTableAll <- rbind(tierTableAll, tierTable)

+                samplesUsed <- c(samplesUsed, sampleIDsUse)

+                

+                if(breakEarly == TRUE) break

+              }

+              pathwayString <- paste(assaysProcessed, collapse = '-')

+              parameters = list(confidenceCutoff = confidenceCutoff, minAssaySamples = minAssaySamples)

+              list(models = modelsList, parameters = parameters, pathway = pathwayString,

+                  individuals = individualsTableAll, tiers = tierTableAll)

+            })

+            

+            class(precisionPathways) <- "PrecisionPathways"

+            names(precisionPathways) <- sapply(precisionPathways, "[[", "pathway")

+            precisionPathways

+}

+

+# A nice print method to avoid flooding the screen with lots of tables

+# when result is shown in console.

+print.PrecisionPathways <- function(x)

+{

+  cat("An object of class 'PrecisionPathways'.\n")

+  cat("Pathways:\n")

+  cat(paste(names(x), collapse = '\n'))

+}

+

+#' @usage NULL

+setGeneric("precisionPathwaysPredict", function(pathways, measurements, class, ...)

+    standardGeneric("precisionPathwaysPredict"))

+

+#' @rdname precisionPathways

+#' @export

+setMethod("precisionPathwaysPredict", "MultiAssayExperimentOrList", 

+          function(pathways, measurements, class)

+          {

+            if(is.list(measurements)) # Ensure plain list has clinical data.

+            {

+              # One of the tables must be named "clinical".

+              if (!any(names(measurements) == "clinical"))

+                stop("One of the tables must be named \"clinical\".")

+            }

+            prepArgs <- list(measurements, outcomeColumns = class)

+            measurementsAndClass <- do.call(prepareData, prepArgs)

+              

+            .precisionPathwaysPredict(pathways, measurementsAndClass[["measurements"]], measurementsAndClass[["outcome"]])

+          })

+

+.precisionPathwaysPredict <- function(pathways, measurements, class)

+{

+  # To do.

+}

\ No newline at end of file

@@ -18,15 +18,13 @@

 #' @param outcomeColumns If \code{measurements} is a \code{MultiAssayExperiment}, the

 #' names of the column (class) or columns (survival) in the table extracted by \code{colData(data)}

 #' that contain(s) the each individual's outcome to use for prediction.

-#' @param useFeatures If \code{measurements} is a \code{MultiAssayExperiment},

-#' a two-column table of features to use. The first column must have assay names

-#' and the second column must have feature names found for that assay. \code{"clinical"} is

-#' also a valid assay name and refers to the clinical data table. \code{"all"} is a special

-#' keyword that means all features (passing any other filters) of that assay will be used 

-#' for modelling. Otherwise, a character vector of feature names to use suffices.

+#' @param clinicalPredictors If \code{measurements} is a \code{MultiAssayExperiment},

+#' a character vector of features to use in modelling. This allows avoidance of things like sample IDs,

+#' sample acquisition dates, etc. which are not relevant for outcome prediction.

 #' @param maxMissingProp Default: 0.0. A proportion less than 1 which is the maximum

 #' tolerated proportion of missingness for a feature to be retained for modelling.

-#' @param topNvariance Default: NULL. An integer number of most variable features to subset to.

+#' @param topNvariance Default: NULL. An integer number of most variable features per assay to subset to.

+#' Assays with less features won't be reduced in size.

 #' @param ... Variables not used by the \code{matrix} nor the

 #' \code{MultiAssayExperiment} method which are passed into and used by the

 #' \code{DataFrame} method.

@@ -58,18 +56,16 @@ setMethod("prepareData", "data.frame",

 #' @rdname prepareData

 #' @export

 setMethod("prepareData", "DataFrame",

-  function(measurements, outcome, useFeatures = "all", maxMissingProp = 0.0, topNvariance = NULL)

+  function(measurements, outcome, clinicalPredictors = NULL, maxMissingProp = 0.0, topNvariance = NULL)

 {

   if(is.null(rownames(measurements)))

   {

     warning("'measurements' DataFrame must have sample identifiers as its row names. Generating generic ones.")

     rownames(measurements) <- paste("Sample", seq_len(nrow(measurements)))

-  }      

-            

-  if(useFeatures != "all") # Subset to only the desired ones.

-    measurements <- measurements[, useFeatures]

-

+  }

+      

   # Won't ever be true if input data was MultiAssayExperiment because wideFormat already produces valid names.  

+  # Need to check if input data was DataFrame because names might not be valid from user.

   if(!all(colnames(measurements) == make.names(colnames(measurements))))

   {

     warning("Unsafe feature names in input data. Converted into safe names.")

@@ -129,6 +125,19 @@ setMethod("prepareData", "DataFrame",

     else # Three columns. Therefore, counting process data.

       outcome <- survival::Surv(outcome[, 1], outcome[, 2], outcome[, 3])

   }

+  

+  if(!is.null(clinicalPredictors))

+  {

+    if(!is.null(mcols(measurements)$assay))

+    {

+      clinicalIndices <- which(mcols(measurements)$assay == "clinical")

+      usePredictors <- intersect(clinicalIndices, which(mcols(measurements)$feature %in% clinicalPredictors))

+      dropIndices <- setdiff(clinicalIndices, usePredictors)

+      if(length(dropIndices) > 0) measurements <- measurements[, -dropIndices]

+    } else { # The DataFrame is entirely clinical data.

+      measurements <- measurements[, clinicalPredictors]

+    }

+  }  

   # Remove samples with indeterminate outcome.

   dropSamples <- which(is.na(outcome) | is.null(outcome))

@@ -146,11 +155,18 @@ setMethod("prepareData", "DataFrame",

   if(length(dropFeatures) > 0)

     measurements <- measurements[, -dropFeatures]

-  # Use only the most variable features.

+  # Use only the most N variable features per assay.

   if(!is.null(topNvariance))

   {

-    mostVariance <- order(apply(measurements, 2, var, na.rm = TRUE), decreasing = TRUE)[1:topNvariance]

-    measurements <- measurements[, mostVariance]

+    if(is.null(mcols(measurements)$assay)) assays <- rep(1, ncol(measurements)) else assays <- mcols(measurements)$assay

+    do.call(cbind, lapply(unqiue(assays), function(assay)

+    {

+      assayColumns <- which(assays == assay)    

+      if(length(assayColumns) < topNvariance)

+        measurements[, assayColumns]

+      else

+        measurements[, assayColumns][order(apply(measurements[, assayColumns], 2, var, na.rm = TRUE), decreasing = TRUE)[1:topNvariance]]  

+    }))

   }

   list(measurements = measurements, outcome = outcome)

@@ -159,71 +175,75 @@ setMethod("prepareData", "DataFrame",

 #' @rdname prepareData

 #' @export

 setMethod("prepareData", "MultiAssayExperiment",

-  function(measurements, outcomeColumns = NULL, useFeatures = "all", ...)

+  function(measurements, outcomeColumns = NULL, clinicalPredictors = NULL, ...)

 {

-  if(is.character(useFeatures)) useFeatures <- data.frame(assay = names(measurements), feature = "all")

-  omicsTargets <- setdiff(useFeatures[, "assay"], "clinical")

-  if(length(omicsTargets) > 0)

-  {

-    if(any(anyReplicated(measurements[, , omicsTargets])))

-      stop("Data set contains replicates. Please remove or average replicate observations and try again.")

-  }

-  

-  if(!is.null(outcomeColumns) && !all(outcomeColumns %in% colnames(MultiAssayExperiment::colData(measurements))))

+  if(is.null(clinicalPredictors))

+    stop("'clinicalPredictors' must be a vector of informative clinical features (i.e. not sample IDs, sampling dates, etc.) to consider for classification.")      

+

+  if(any(anyReplicated(measurements)))

+    stop("Data set contains replicates. Please remove or average replicate observations and try again.")

+ 

+  if(is.null(outcomeColumns))

+    stop("'outcomeColumns' is a mandatory parameter but was not specified.")

+      

+  if(!all(outcomeColumns %in% colnames(MultiAssayExperiment::colData(measurements))))

     stop("Not all column names specified by 'outcomeColumns' found in clinical table.")  

-  if(!all(useFeatures[, "assay"] %in% c(names(measurements), "clinical")))

-    stop("Some assay names in first column of 'useFeatures' are not assay names in 'measurements' or \"clinical\".")

-  clinicalColumnsDataset <- colnames(MultiAssayExperiment::colData(measurements))

-  if("clinical" %in% useFeatures[, "assay"])

-  {

-    clinicalRows <- useFeatures[, "assay"] == "clinical"      

-    clinicalColumns <- useFeatures[clinicalRows, "feature"]

-    if(length(clinicalColumns) == 1 && clinicalColumns == "all")

-      clinicalColumns <- setdiff(clinicalColumnsDataset, outcomeColumns)

-    useFeatures <- useFeatures[!clinicalRows, ]

-  } else {

-    clinicalColumns <- NULL

-  }

-  

-  if(nrow(useFeatures) > 0)

-  {

-    measurements <- measurements[, , unique(useFeatures[, "assay"])]

-    # Get all desired measurements tables and clinical columns (other than the columns representing outcome).

-    # These form the independent variables to be used for making predictions with.

-    # Variable names will have names like RNA_BRAF for traceability.

-    dataTable <- MultiAssayExperiment::wideFormat(measurements, colDataCols = union(clinicalColumns, outcomeColumns))

-    rownames(dataTable) <- dataTable[, "primary"]

-    S4Vectors::mcols(dataTable)[, "sourceName"] <- gsub("colDataCols", "clinical", S4Vectors::mcols(dataTable)[, "sourceName"])

-    colnames(S4Vectors::mcols(dataTable))[1] <- "assay"

+  # Get all desired measurements tables and clinical columns.

+  # These form the independent variables to be used for making predictions with.

+  # Variable names will have names like RNA_BRAF for traceability.

+  dataTable <- MultiAssayExperiment::wideFormat(measurements, colDataCols = union(clinicalPredictors, outcomeColumns))

+  rownames(dataTable) <- dataTable[, "primary"]

+  S4Vectors::mcols(dataTable)[, "sourceName"] <- gsub("colDataCols", "clinical", S4Vectors::mcols(dataTable)[, "sourceName"])

+  dataTable <- dataTable[, -match("primary", colnames(dataTable))]

+  colnames(S4Vectors::mcols(dataTable))[1] <- "assay"

-    # Sample information variable names not included in column metadata of wide table but only as row names of it.

-    # Create a combined column named "feature" which has feature names of the assays as well as the clinical.

-    S4Vectors::mcols(dataTable)[, "feature"] <- as.character(S4Vectors::mcols(dataTable)[, "rowname"])

-    missingIndices <- is.na(S4Vectors::mcols(dataTable)[, "feature"])

-    S4Vectors::mcols(dataTable)[missingIndices, "feature"] <- colnames(dataTable)[missingIndices]

+  # Sample information variable names not included in column metadata of wide table but only as row names of it.

+  # Create a combined column named "feature" which has feature names of the assays as well as the clinical.

+  S4Vectors::mcols(dataTable)[, "feature"] <- as.character(S4Vectors::mcols(dataTable)[, "rowname"])

+  missingIndices <- is.na(S4Vectors::mcols(dataTable)[, "feature"])

+  S4Vectors::mcols(dataTable)[missingIndices, "feature"] <- colnames(dataTable)[missingIndices]

-    # Finally, a column annotation recording variable name and which table it originated from for all of the source tables.

-    S4Vectors::mcols(dataTable) <- S4Vectors::mcols(dataTable)[, c("assay", "feature")]

+  # Finally, a column annotation recording variable name and which table it originated from for all of the source tables.

+  S4Vectors::mcols(dataTable) <- S4Vectors::mcols(dataTable)[, c("assay", "feature")]

-    # Subset to only the desired features.

-    useFeaturesSubset <- useFeatures[useFeatures[, "feature"] != "all", ]

-    if(nrow(useFeaturesSubset) > 0)

-    {

-      uniqueAssays <- unique(useFeatures[, "assay"])

-      for(filterAssay in uniqueAssays)

-      {

-        dropFeatures <- S4Vectors::mcols(dataTable)[, "assay"] == filterAssay &

-                        !S4Vectors::mcols(dataTable)[, "feature"] %in% useFeatures[useFeatures[, 1] == filterAssay, 2]

-        dataTable <- dataTable[, !dropFeatures]

-      }

-    }

-    dataTable <- dataTable[, -match("primary", colnames(dataTable))]

-  } else { # Must have only been clinical data.

-    dataTable <- MultiAssayExperiment::colData(measurements)

-    S4Vectors::mcols(dataTable) <- DataFrame(assay = "clinical", feature = colnames(dataTable))

-  }

+  # Do other filtering and preparation in DataFrame function.

+  prepareData(dataTable, outcomeColumns, clinicalPredictors = NULL, ...)

+})

+

+#' @rdname prepareData

+#' @export

+setMethod("prepareData", "list",

+  function(measurements, outcome = NULL, clinicalPredictors = NULL, ...)

+{

+  # Check the list is named.

+  if(is.null(names(measurements)))

+    stop("'measurements' must be a named list.")

+

+  # If clinical table is present, features to use must be user-specified.            

+  if("clinical" %in% names(measurements) && is.null(clinicalPredictors))

+    stop("Because one provided table in the list is named \"clinical\", 'clinicalPredictors' must be a vector of informative clinical features (i.e. not sample IDs, sampling dates, etc.) to consider for classification.")

+

+  # Check data type is valid.

+  if(!(all(sapply(measurements, class) %in% c("data.frame", "DataFrame", "matrix"))))

+    stop("assays in the list must be of type data.frame, DataFrame or matrix")

+              

+  # Check same number of samples for all datasets

+  if (!length(unique(sapply(measurements, nrow))) == 1)

+    stop("All datasets must have the same samples.")

+      

+  if("clinical" %in% names(measurements))

+    measurements[["clinical"]] <- measurements[["clinical"]][, clinicalPredictors]

+             

+  allMetadata <- mapply(function(measurementsOne, assayID) {

+                        data.frame(assay = assayID, feature = colnames(measurementsOne))

+                        }, measurements, names(measurements))

+  allMeasurements <- do.call("cbind", measurements)

+  # Different assays e.g. mRNA, protein could have same feature name e.g. BRAF.

+  colnames(allMeasurements) <- paste(allMetadata[, "assay"], allMetadata[, "feature"], sep = '_')

+  allDataFrame <- DataFrame(allMeasurements)

+  S4Vectors::mcols(allMeasurements) <- allMetadata

   # Do other filtering and preparation in DataFrame function.

-  prepareData(dataTable, outcomeColumns, useFeatures = "all", ...)

+  prepareData(dataTable, outcome, clinicalPredictors = NULL, ...)

 })

\ No newline at end of file

@@ -6,8 +6,7 @@

 \alias{crossValidate,DataFrame-method}

 \alias{crossValidate,MultiAssayExperiment-method,}

 \alias{crossValidate,data.frame-method}

-\alias{crossValidate,MultiAssayExperiment-method}

-\alias{crossValidate,list-method}

+\alias{crossValidate,MultiAssayExperimentOrList-method}

 \alias{train.matrix}

 \alias{train.data.frame}

 \alias{train.DataFrame}

@@ -16,8 +15,6 @@

 \alias{predict.trainedByClassifyR}

 \title{Cross-validation to evaluate classification performance.}

 \usage{

-crossValidate(measurements, outcome, ...)

-

 \S4method{crossValidate}{DataFrame}(

   measurements,

   outcome,

@@ -35,9 +32,10 @@ crossValidate(measurements, outcome, ...)

   extraParams = NULL

 )

-\S4method{crossValidate}{MultiAssayExperiment}(

+\S4method{crossValidate}{MultiAssayExperimentOrList}(

   measurements,

   outcome,

+  clinicalPredictors = NULL,

   nFeatures = 20,

   selectionMethod = "auto",

   selectionOptimisation = "Resubstitution",

@@ -86,23 +84,6 @@ crossValidate(measurements, outcome, ...)

   extraParams = NULL

 )

-\S4method{crossValidate}{list}(

-  measurements,

-  outcome,

-  nFeatures = 20,

-  selectionMethod = "auto",

-  selectionOptimisation = "Resubstitution",

-  performanceType = "auto",

-  classifier = "auto",

-  multiViewMethod = "none",

-  assayCombinations = "all",

-  nFolds = 5,

-  nRepeats = 20,

-  nCores = 1,

-  characteristicsLabel = NULL,

-  extraParams = NULL

-)

-

 \method{train}{matrix}(x, outcomeTrain, ...)

 \method{train}{data.frame}(x, outcomeTrain, ...)

@@ -122,7 +103,7 @@ crossValidate(measurements, outcome, ...)

 \method{train}{list}(x, outcomeTrain, ...)

-\method{train}{MultiAssayExperiment}(x, outcome, ...)

+\method{train}{MultiAssayExperiment}(x, outcome, clinicalPredictors = NULL, ...)

 \method{predict}{trainedByClassifyR}(object, newData, ...)

 }

@@ -146,7 +127,8 @@ Set to NULL or "all" if all features should be used.}

 \item{selectionMethod}{Default: \code{"auto"}. A character vector of feature selection methods to compare. If a named character vector with names corresponding to different assays, 

 and performing multiview classification, the respective selection methods will be used on each assay. If \code{"auto"}, t-test (two categories) / F-test (three or more categories) ranking

-and top \code{nFeatures} optimisation is done. Otherwise, the ranking method is per-feature Cox proportional hazards p-value.}

+and top \code{nFeatures} optimisation is done. Otherwise, the ranking method is per-feature Cox proportional hazards p-value. \code{"none"} is also a valid value, meaning that no

+indepedent feature selection will be performed (but implicit selection might still happen with the classifier).}

 \item{selectionOptimisation}{A character of "Resubstitution", "Nested CV" or "none" specifying the approach used to optimise \code{nFeatures}.}

@@ -172,7 +154,11 @@ with each element being a vector of assays to combine. Special value \code{"all"

 \item{extraParams}{A list of parameters that will be used to overwrite default settings of transformation, selection, or model-building functions or

 parameters which will be passed into the data cleaning function. The names of the list must be one of \code{"prepare"},

 \code{"select"}, \code{"train"}, \code{"predict"}. To remove one of the defaults (see the article titled Parameter Tuning Presets for crossValidate and Their Customisation on

-the website), specify the list element to be \code{NULL}.}

+the website), specify the list element to be \code{NULL}. For the valid element names in the \code{"prepare"} list, see \code{?prepareData}.}

+

+\item{clinicalPredictors}{Default: \code{NULL}. If \code{measurements} is a \code{MultiAssayExperiment},

+a character vector of features to use in modelling. This allows avoidance of things like sample IDs,

+sample acquisition dates, etc. which are not relevant for outcome prediction.}

 \item{x}{Same as \code{measurements} but only training samples.}

new file mode 100644

@@ -0,0 +1,80 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/precisionPathways.R

+\name{precisionPathwaysTrain}

+\alias{precisionPathwaysTrain}

+\alias{precisionPathwaysTrain,MultiAssayExperimentOrList-method}

+\alias{precisionPathwaysPredict,MultiAssayExperimentOrList-method}

+\title{Precision Pathways for Sample Prediction Based on Prediction Confidence.}

+\usage{

+\S4method{precisionPathwaysTrain}{MultiAssayExperimentOrList}(

+  measurements,

+  class,

+  clinicalPredictors = NULL,

+  maxMissingProp = 0,

+  topNvariance = NULL,

+  fixedAssays = "clinical",

+  confidenceCutoff = 0.8,

+  minAssaySamples = 10,

+  nFeatures = 20,

+  selectionMethod = setNames(c("none", rep("t-test", length(measurements))),

+    c("clinical", names(measurements))),

+  classifier = setNames(c("elasticNetGLM", rep("randomForest", length(measurements))),

+    c("clinical", names(measurements))),

+  nFolds = 5,

+  nRepeats = 20,

+  nCores = 1

+)

+

+\S4method{precisionPathwaysPredict}{MultiAssayExperimentOrList}(pathways, measurements, class)

+}

+\arguments{

+\item{measurements}{Either a \code{\link{MultiAssayExperiment}} or a list of the basic tabular objects containing the data.}

+

+\item{class}{Same as \code{measurements} but only training samples. IF \code{measurements} is a \code{list}, may also be

+a vector of classes.}

+

+\item{clinicalPredictors}{Default: \code{NULL}. Must be a character vector of clinical features to use in modelling. This allows avoidance of things like sample IDs,

+sample acquisition dates, etc. which are not relevant for outcome prediction.}

+

+\item{maxMissingProp}{Default: 0.0. A proportion less than 1 which is the maximum

+tolerated proportion of missingness for a feature to be retained for modelling.}

+

+\item{topNvariance}{Default: NULL. An integer number of most variable features per assay to subset to.

+Assays with less features won't be reduced in size.}

+

+\item{fixedAssays}{A character vector of assay names specifying any assays which must be at the

+beginning of the pathway.}

+

+\item{confidenceCutoff}{The minimum confidence of predictions for a sample to be predicted by a particular issue

+. If a sample was predicted to belong to a particular class a proportion \eqn{p} times, then the confidence is \eqn{2 \times |p - 0.5|}.}

+

+\item{minAssaySamples}{An integer specifying the minimum number of samples a tier may have. If a subsequent tier

+would have less than this number of samples, the samples are incorporated into the current tier.}

+

+\item{nFeatures}{Default: 20. The number of features to consider during feature selection, if feature selection is done.}

+

+\item{selectionMethod}{A named character vector of feature selection methods to use for the assays, one for each. The names must correspond to names of \code{measurements}.}

+

+\item{classifier}{A named character vector of modelling methods to use for the assays, one for each. The names must correspond to names of \code{measurements}.}

+

+\item{nFolds}{A numeric specifying the number of folds to use for cross-validation.}

+

+\item{nRepeats}{A numeric specifying the the number of repeats or permutations to use for cross-validation.}

+

+\item{nCores}{A numeric specifying the number of cores used if the user wants to use parallelisation.}

+

+\item{pathways}{A set of pathways created by \code{precisionPathwaysTrain} to be used for predicting on a new data set.}

+}

+\value{

+An object of class \code{PrecisionPathways} which is basically a named list that other plotting and