Commit 51708be4 authored by Konrad Zych's avatar Konrad Zych

Merge branch 'dev' of git.embl.de:grp-zeller/SIAMCAT into dev

parents db727a45 e66ec0c9
Pipeline #6328 passed with stages
in 5 minutes and 30 seconds
......@@ -29,7 +29,7 @@ preprocess_data:
- Rscript 01_validate_data.r --metadata_in num_metadata_cancer-vs-healthy_study-pop-I_N141_tax_profile_mocat_bn_genus.tsv --metadata_out valMetaData.tsv --label_in label_cancer-vs-healthy_study-pop-I_N141_tax_profile_mocat_bn_genus.tsv --label_out valLabel.tsv --feat_in feat_cancer-vs-healthy_study-pop-I_N141_tax_profile_mocat_bn_genus.tsv --feat_out valFeat.tsv
- Rscript 02_select_samples.r --metadata_in valMetaData.tsv --metadata_out valMetaData_sel.tsv --label_in valLabel.tsv --label_out valLabel_sel.tsv --feat_in valFeat.tsv --feat_out valFeat_sel.tsv --filter_var="age" --allowed_range="[0,120]"
- Rscript 03_check_for_confounders.r --metadata_in valMetaData_sel.tsv --feat_in valFeat.tsv --plot metaCheck.pdf --label_in valLabel_sel.tsv
- Rscript 04_filter_features.r --feat_in valFeat_sel.tsv --feat_out valFeat_sel_filtered.tsv --method="abundance" --cutoff="0.001" --recomp_prop="FALSE" --rm_unmapped="TRUE"
- Rscript 04_filter_features.r --feat_in valFeat_sel.tsv --feat_out valFeat_sel_filtered.tsv --method="abundance" --cutoff="0.001" --rm_unmapped="TRUE"
- Rscript 05_check_associations.r --feat_in valFeat_sel_filtered.tsv --label_in valLabel_sel.tsv --plot assoc.pdf --col_scheme="RdYlBu" --alpha="0.7" --min_fc="0" --mult_test="fdr" --max_show="50" --detect_limit="1e-06" --plot_type="bean"
- Rscript 06_normalize_features.r --feat_in valFeat_sel_filtered.tsv --feat_out valFeat_sel_norm.tsv --param_out param_out.tsv --method="log.unit" --log_n0="1e-08" --sd_min_quantile="0.2" --vector_norm="2" --norm_margin="1"
- Rscript 07_add_metadata_as_predictor.r --feat_in valFeat_sel_norm.tsv --feat_out valFeat_sel_norm.tsv --metadata_in valMetaData_sel.tsv --pred_names="fobt" --std_meta="TRUE"
......@@ -41,7 +41,7 @@ preprocess_data:
only:
- master
- development
- dev
train_model:
stage: train_model
......@@ -58,7 +58,7 @@ train_model:
only:
- master
- development
- dev
make_predictions:
stage: make_predictions
......@@ -75,7 +75,7 @@ make_predictions:
only:
- master
- development
- dev
evaluate_interpret:
stage: evaluate_interpret
......@@ -83,7 +83,7 @@ evaluate_interpret:
- cd Rscript_flavor
- R CMD INSTALL SIAMCAT_1.0.1.tar.gz
script:
- Rscript 11_evaluate_predictions.r --feat_in valFeat_sel_norm.tsv --label_in valLabel_sel.tsv --plot evalPlot.pdf --pred pred.tsv
- Rscript 11_evaluate_predictions.r --feat_in valFeat_sel_norm.tsv --label_in valLabel_sel.tsv --data_split dataSplit.Rdata --plot evalPlot.pdf --pred pred.tsv
- Rscript 12_interpret_model.r --label_in valLabel_sel.tsv --feat_in valFeat_sel_norm.tsv --origin_feat valFeat_sel.tsv --metadata_in valMetaData_sel.tsv --mlr_models_list models.RData --plot interpretation.pdf --pred pred.tsv --col_scheme="BrBG" --heatmap_type="zscore" --consens_thres="0.5"
artifacts:
......@@ -93,4 +93,4 @@ evaluate_interpret:
only:
- master
- development
- dev
......@@ -49,10 +49,11 @@ Imports:
RColorBrewer,
scales,
stats,
stringr,
utils
License: GPL-3
LazyData: true
RoxygenNote: 6.0.1.9000
RoxygenNote: 6.1.0
biocViews: Metagenomics, Classification, Microbiome, Sequencing, Preprocessing,
Clustering, FeatureExtraction, GeneticVariability, MultipleComparison,
Regression
......
# Generated by roxygen2: do not edit by hand
export("associations<-")
export("data_split<-")
export("eval_data<-")
export("features<-")
export("filt_feat<-")
export("label<-")
export("meta<-")
export("model_list<-")
export("norm_param<-")
export("norm_feat<-")
export("orig_feat<-")
export("physeq<-")
export("pred_matrix<-")
export(accessSlot)
export(add.meta.pred)
export(assoc_param)
export(associations)
export(check.associations)
export(check.confounders)
export(create.data.split)
export(create.label.from.metadata)
export(create.label)
export(data_split)
export(eval_data)
export(evaluate.predictions)
export(features)
export(feature_type)
export(feature_weights)
export(filt_feat)
export(filt_params)
export(filter.features)
export(filter.label)
export(get.features.matrix)
export(get.filt_feat.matrix)
export(get.norm_feat.matrix)
export(get.orig_feat.matrix)
export(label)
export(make.predictions)
......@@ -32,16 +39,14 @@ export(model.interpretation.plot)
export(model_list)
export(model_type)
export(models)
export(norm_param)
export(norm_feat)
export(norm_params)
export(normalize.features)
export(orig_feat)
export(physeq)
export(pred_matrix)
export(read.features)
export(read.labels)
export(read.label)
export(read.lefse)
export(read.meta)
export(reset.features)
export(select.samples)
export(siamcat)
export(siamcat.to.lefse)
......@@ -49,12 +54,14 @@ export(siamcat.to.maaslin)
export(summarize.features)
export(train.model)
export(validate.data)
export(weight_matrix)
exportClasses(associations)
exportClasses(data_split)
exportClasses(eval_data)
exportClasses(filt_feat)
exportClasses(label)
exportClasses(model_list)
exportClasses(norm_param)
exportClasses(orig_feat)
exportClasses(norm_feat)
exportClasses(pred_matrix)
exportClasses(siamcat)
import(LiblineaR)
......@@ -127,6 +134,8 @@ importFrom(stats,rnorm)
importFrom(stats,runif)
importFrom(stats,sd)
importFrom(stats,wilcox.test)
importFrom(stringr,str_extract_all)
importFrom(utils,askYesNo)
importFrom(utils,head)
importFrom(utils,read.csv)
importFrom(utils,read.table)
......
......@@ -7,13 +7,18 @@
#' @title Add metadata as predictors
#' @description This function adds metadata to the feature matrix to be later
#' used as predictors
#' @usage add.meta.pred(siamcat, pred.names = NULL, std.meta =
#' TRUE, verbose = 1)
#' @usage add.meta.pred(siamcat, pred.names,
#' std.meta = TRUE,
#' feature.type='normalized',
#' verbose = 1)
#' @param siamcat object of class \link{siamcat-class}
#' @param pred.names vector of names of the variables within the metadata to be
#' added to the feature matrix as predictors
#' @param std.meta boolean, should added metadata features be standardized?,
#' defaults to \code{TRUE}
#' @param feature.type On which type of features should the function work? Can
#' be either "original", "filtered", or "normalized". Please only change this
#' paramter if you know what you are doing!
#' @param verbose control output: \code{0} for no output at all, \code{1} for
#' only information about progress and success, \code{2} for normal level of
#' information and \code{3} for full debug information, defaults to \code{1}
......@@ -24,17 +29,46 @@
#' @examples
#' data(siamcat_example)
#' # Add the Age of the patients as potential predictor
#' siamcat_age_added <- add.meta.pred(siamcat_example, pred.names=c('age'))
#' siamcat_age_added <- add.meta.pred(siamcat_example, pred.names=c('Age'))
#'
#' # Add Age, BMI, and Gender as potential predictors
#' # Additionally, prevent standardization of the added features
#' siamcat_meta_added <- add.meta.pred(siamcat_example, pred.names=c('age',
#' 'bmi', 'gender'), std.meta=FALSE)
add.meta.pred <- function(siamcat, pred.names = NULL, std.meta = TRUE,
#' siamcat_meta_added <- add.meta.pred(siamcat_example, pred.names=c('Age',
#' 'BMI', 'Gender'), std.meta=FALSE)
add.meta.pred <- function(siamcat, pred.names,
std.meta = TRUE,
feature.type = 'normalized',
verbose = 1) {
if (verbose > 1)
message("+ starting add.meta.pred")
s.time <- proc.time()[3]
if (is.null(meta(siamcat))) {
stop('SIAMCAT object has no metadata. Exiting...')
}
if (!feature.type %in% c('original', 'filtered', 'normalized')){
stop("Unrecognised feature type, exiting...\n")
}
if (feature.type != 'normalized'){
warning('It is recommended to add a meta-predictor only',
' after feature normalization.')
}
# get the right features
if (feature.type == 'original'){
feat <- get.orig_feat.matrix(siamcat)
} else if (feature.type == 'filtered'){
if (is.null(filt_feat(siamcat, verbose=0))){
stop('Features have not yet been filtered, exiting...\n')
}
feat <- get.filt_feat.matrix(siamcat)
} else if (feature.type == 'normalized'){
if (is.null(norm_feat(siamcat, verbose=0))){
stop('Features have not yet been normalized, exiting...\n')
}
feat <- get.norm_feat.matrix(siamcat)
}
### add metadata as predictors to the feature matrix
cnt <- 0
......@@ -47,6 +81,9 @@ add.meta.pred <- function(siamcat, pred.names = NULL, std.meta = TRUE,
if (!p %in% colnames(meta(siamcat))) {
stop("There is no metadata variable called ", p)
}
if (paste0('META_', toupper(p)) %in% rownames(feat)){
stop("This meta-variable has already been added. Exiting...")
}
idx <- which(colnames(meta(siamcat)) == p)
if (length(idx) != 1)
stop(p, "matches multiple columns in the metada")
......@@ -55,7 +92,7 @@ add.meta.pred <- function(siamcat, pred.names = NULL, std.meta = TRUE,
# check if the meta-variable is a factor or numeric
if (class(m) == 'factor'){
message(paste0('WARNING: the meta-variable is a factor and',
warning(paste0('WARNING: meta-variable ', p,' is a factor and',
' not numeric...\n The values will be converted to numerical',
' values with "as.numeric()"\n'))
m <- as.numeric(m)
......@@ -79,12 +116,21 @@ add.meta.pred <- function(siamcat, pred.names = NULL, std.meta = TRUE,
stopifnot(!m.sd == 0)
m <- (m - m.mean)/m.sd
}
features.with.meta <- otu_table(rbind(features(siamcat),m),
features.with.meta <- otu_table(rbind(feat,m),
taxa_are_rows = TRUE)
rownames(features.with.meta)[nrow(features.with.meta)] <- paste(
"META_", toupper(p), sep = "")
features(siamcat) <- features.with.meta
if (feature.type == 'original'){
orig_feat(siamcat) <- features.with.meta
} else if (feature.type == 'filtered'){
filt_feat(siamcat) <- new('filt_feat',
filt.feat=features.with.meta,
filt.param=filt_params(siamcat))
} else if (feature.type == 'normalized'){
norm_feat(siamcat) <- new("norm_feat",
norm.feat=features.with.meta,
norm.param=norm_params(siamcat))
}
cnt <- cnt + 1
}
if (verbose > 1)
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -12,7 +12,7 @@
#' \code{num.resample} times.
#'
#' @usage create.data.split(siamcat, num.folds = 2, num.resample = 1,
#' stratify = TRUE,inseparable = NULL, verbose = 1)
#' stratify = TRUE, inseparable = NULL, verbose = 1)
#'
#' @param siamcat object of class \link{siamcat-class}
#'
......@@ -64,13 +64,9 @@
#' ## # example with a variable which is to be inseparable
#' ## siamcat_split <- create.data.split(siamcat_example, num.folds=10,
#' ## num.resample=5, stratify=FALSE, inseparable='Gender')
create.data.split <-
function(siamcat,
num.folds = 2,
num.resample = 1,
stratify = TRUE,
inseparable = NULL,
verbose = 1) {
create.data.split <- function(siamcat, num.folds = 2, num.resample = 1,
stratify = TRUE, inseparable = NULL, verbose = 1) {
if (verbose > 1)
message("+ starting create.data.split")
s.time <- proc.time()[3]
......
......@@ -3,24 +3,28 @@
### Microbial Communities And host phenoTypes R flavor EMBL
### Heidelberg 2012-2018 GNU GPL 3.0
#' @title create a label object from metadata
#' @title create a label object from metadata or an atomic vector
#'
#' @description This function creates a label object from metadata
#' or an atomic vector
#'
#' @usage create.label.from.metadata(meta, column, case,
#' control=NULL, p.lab = NULL, n.lab = NULL, verbose=1)
#' @usage create.label(label, case,
#' meta=NULL, control=NULL,
#' p.lab = NULL, n.lab = NULL,
#' remove.meta.column=FALSE,
#' verbose=1)
#'
#' @param meta metadata as read by \link{read.meta}
#' of \link[phyloseq]{sample_data-class}
#' @param label named vector to create the label or the name of the metadata
#' column that will be used to create the label
#'
#' @param column name of column that will be used
#' to create the label
#'
#' @param case name of a label that will be used as a positive label. If the
#' @param case name of the group that will be used as a positive label. If the
#' variable is binary, the other label will be used as a negative one. If the
#' variable has multiple values, all the other values will be used a negative
#' label (testing one vs rest).
#'
#' @param meta metadata dataframe object or an object of class
#' \link[phyloseq]{sample_data-class}
#'
#' @param control name of a label or vector with names that will be used as a
#' negative label. All values that are nor equal to case and control will be
#' dropped. Default to NULL in which case: If the variable is binary, the value
......@@ -28,106 +32,150 @@
#' values, all the values not equal to cases will be used a negative label
#' (testing one vs rest).
#'
#' @param p.lab name of the positive label (useful mostly for visualizations).
#' Default to NULL in which case the value of the positive label will be used.
#' @param p.lab name of the positive group (useful mostly for visualizations).
#' Default to NULL in which case the value of the positive group will be used.
#'
#' @param n.lab name of the negative label (useful mostly for visualizations).
#' Default to NULL in which case the value of the negative label will be used
#' @param n.lab name of the negative group (useful mostly for visualizations).
#' Default to NULL in which case the value of the negative group will be used
#' for binary variables and "rest" will be used for variables with multiple
#' values.
#'
#' @param remove.meta.column boolean indicating if the label column in the
#' metadata should be retained. Please note that if this is set to
#' \code{TRUE}, the function will return a list as result. Defaults to
#' \code{FALSE}
#'
#' @param verbose control output: \code{0} for no output at all, \code{1}
#' for only information about progress and success, \code{2} for normal
#' level of information and \code{3} for full debug information,
#' defaults to \code{1}
#'
#' @keywords create.label.from.metadata
#' @keywords create.label
#'
#' @return an object of class \link{label-class}
#' @return an object of class \link{label-class} OR a list with entries
#' \code{meta} and \code{label}, if \code{remove.meta.column} is set to
#' \code{TRUE}
#'
#' @examples
#' data(siamcat_example)
#' label <- create.label.from.metadata(meta(siamcat_example),"fobt",
#' case = 1, control = 0)
#' data('meta_crc_zeller')
#' label <- create.label(label='Group', case='CRC', meta=meta.crc.zeller)
#'
#' @export
create.label.from.metadata <- function(meta, column, case, control = NULL,
verbose=1) {
create.label <- function(label, case, meta=NULL, control = NULL,
p.lab=NULL, n.lab=NULL, remove.meta.column=FALSE, verbose=1) {
if (verbose > 1)
message("+ starting create.label.from.metadata")
message("+ starting create.label")
s.time <- proc.time()[3]
if (!column %in% colnames(meta))
stop("ERROR: Column ", column, " not found in the metadata\n")
#if metadata has been supplied and the label is of length 1
if (!is.null(meta) & length(label) == 1){
if (!label %in% colnames(meta))
stop("ERROR: Column ", label, " not found in the metadata\n")
if (class(meta) == 'sample_data'){
label.vec <- vapply(meta[, label], as.character,
FUN.VALUE = character(nrow(meta)))
} else if (class(meta) == 'data.frame'){
label.vec <- vapply(meta[, label], as.character,
FUN.VALUE = character(1))
} else {
stop(paste0('Please provide the metadata either as a data.frame',
' or a sample_data object!'))
}
names(label.vec) <- rownames(meta)
#if the label is a vector
} else if (is.atomic(label) & length(label) > 1){
label.vec <- label
if (is.null(names(label.vec))){
stop('The label vectors needs to be named!')
}
if (is.factor(label.vec)){
names.old <- names(label.vec)
label.vec <- as.character(label.vec)
names(label.vec) <- names.old
}
} else {
stop(paste0('Could not interpret your input.\n Please provide ',
'either a column name and metadata or a label vector.\nExiting...!'))
}
metaColumn <- vapply(meta[, column], as.character,
FUN.VALUE = character(nrow(meta)))
names(metaColumn) <- rownames(meta)
# remove NAs in the label
if (any(is.na(label.vec))){
if (verbose > 1) message(paste0('+ removing ', sum(is.na(label.vec)),
' instances of NA in the label'))
label.vec <- label.vec[!is.na(label.vec)]
}
labels <- unique(metaColumn)
# get different groups
groups <- unique(label.vec)
### checking case
if(!all(case %in% labels)){
stop("Column ", column, " does not contain values:",
paste(case,collapse=","),"\n")
if(!all(case %in% groups)){
stop("The chosen label does not contain values: ",
paste(case,collapse=","),"\nInstead, contains: ",
paste(groups, collapse=','))
}
### checking control
if (is.null(control)) {
if((length(labels)-length(case))>1){
if((length(groups)-length(case))>1){
control <- "rest"
}else{
control <- setdiff(labels, case)
control <- setdiff(groups, case)
}
}else{
if(!control%in%labels){
stop("Column ", column, " does not contain value:",control,"\n")
if(!control%in%groups){
stop("The chose label does not contain value:",control,
"\nInstead, contains: ", paste(groups, collapse=','))
}
### dropping unused values
if(any(!labels%in%c(case, control))){
metaColumn <- metaColumn[which(metaColumn%in%c(case, control))]
if(any(!groups%in%c(case, control))){
label.vec <- label.vec[which(label.vec%in%c(case, control))]
warning("Dropping values: ",
labels[which(!labels%in%c(case, control))])
groups[which(!groups%in%c(case, control))])
}
}
# message status
if (verbose > 0)
message("Label used as case:\n ",paste(case,collapse=","),
"\nLabel used as control:\n ",paste(control,collapse=","))
label <- list(label = rep(-1, length(metaColumn)))
n.lab <- gsub("[_.-]", " ", control)
if (length(case) > 1) {
p.lab <- "Case"
} else {
p.lab <- gsub("[_.-]", " ", case)
}
# create new label.object
label.new <- list(label = rep(-1, length(label.vec)))
n.lab <- ifelse(is.null(n.lab), gsub("[_.-]", " ", control), n.lab)
p.lab <- ifelse(is.null(p.lab),
ifelse(length(case) > 1, 'Case', gsub("[_.-]", " ", case)), p.lab)
info <- c(-1, 1)
names(info) <- c(n.lab, p.lab)
names(label$label) <- names(metaColumn)
names(label.new$label) <- names(label.vec)
if (length(case) > 1){
label$label[which(metaColumn %in% case)] <- 1
label.new$label[which(label.vec %in% case)] <- 1
} else {
label$label[which(metaColumn == case)] <- 1
label.new$label[which(label.vec == case)] <- 1
}
label$type <- "BINARY"
label$info <- info
label.new$info <- info
label.new$type <- "BINARY"
labelRes <- label(label)
label.new <- label(label.new)
e.time <- proc.time()[3]
e.time <- proc.time()[3]
if (verbose > 0)
message(paste(
"+ finished create.label.from.metadata in",
formatC(e.time - s.time, digits = 3),
"s"
))
return(labelRes)
if (!remove.meta.column){
return(label.new)
} else {
meta.new <- meta[,-which(colnames(meta) == label)]
return(list(label=label.new,
meta=meta.new))
}
}
......@@ -2,15 +2,39 @@
#'
#' Reduced version of the CRC dataset in inst/extdata, containing 100 features
#' (15 associated features at 5\% FDR in the original dataset and 85 random
#' other features) and 141 samples, saved after the complete SIAMCAT
#' other features) and 141 samples, saved after the complete SIAMCAT
#' pipelinehas been run.
#' Therefore, contains entries in every siamcat-object slot,
#' e.g, \code{eval_data} or \code{data_split}. Mainly used for running
#' the examples in the function documentation
#'
#'
#' @name siamcat_example
#'
#'
#' @docType data
#'
#' @keywords data
NULL
#' Documentation for the example feature object in the data folder
#'
#' Feature matrix (as data.frame) for the CRC dataset, containing 141 samples
#' and 1754 bacterial species (features).
#'
#' @name feat.crc.zeller
#'
#' @docType data
#'
#'
#' @keywords data
NULL
#' Documentation for the example metadata object in the data folder
#'
#' Metadata (as data.frame) for the CRC dataset, containing 6 variables (e.g.
#' Age or BMI) for 141 samples.
#'
#' @name meta.crc.zeller
#'
#' @docType data
#'
#' @keywords data
NULL
......@@ -32,23 +32,26 @@
#' the results will be saved in the \code{eval_data}-slot of the
#' supplied \link{siamcat-class}-object. The \code{eval_data}-slot
#' contains a list with several entries: \itemize{
#' \item \code{$roc.average} average ROC-curve across repeats or a
#' \item \code{$roc} average ROC-curve across repeats or a
#' single ROC-curve on complete dataset;
#' \item \code{$auc.average} AUC value for the average ROC-curve;
#' \item \code{$ev.list} list of \code{length(num.folds)}, containing
#' for different decision thresholds the number of false
#' positives, false negatives, true negatives, and
#' true positives;
#' \item \code{$pr.list} list of \code{length(num.folds)}, containing
#' the positive predictive value (precision) and true positive
#' rate (recall) values used to plot the PR curves.}
#' \item \code{$auroc} AUC value for the average ROC-curve;
#' \item \code{$prc} list containing the positive predictive value
#' (precision) and true positive rate (recall) values used
#' to plot the mean PR curve;
#' \item \code{$auprc} AUC value for the mean PR curve;
#' \item \code{$ev} list containing for different decision thresholds
#' the number of false positives, false negatives, true
#' negatives, and true positives.}
#' For the case of repeated cross-validation, the function will
#' additonally return \itemize{
#' \item \code{$roc.all} list of roc objects (see \link[pROC]{roc})
#' for every repeat;
#' \item \code{$aucspr} vector of AUC values for the PR curves
#' \item \code{$auroc.all} vector of AUC values for the ROC curves
#' for every repeat;
#' \item \code{$auc.all} vector of AUC values for the ROC curves
#' \item \code{$prc.all} list of PR curves for every repeat;
#' \item \code{$auprc.all} vector of AUC values for the PR curves
#' for every repeat;
#' \item \code{$ev.all} list of \code{ev} lists (see above)
#' for every repeat.}
#'
#' @export
......@@ -65,8 +68,14 @@
evaluate.predictions <- function(siamcat, verbose = 1) {
if (verbose > 1)
message("+ starting evaluate.predictions")
label <- label(siamcat)
s.time <- proc.time()[3]
label <- label(siamcat)
if (label$type == 'TEST'){
stop('SIAMCAT can not evaluate the predictions on a TEST',
' label. Exiting...')
}
summ.stat = "mean" # TODO make this a possible parameter?
# TODO compare header to label make sure that label and prediction are in
# the same order
m <- match(names(label$label), rownames(pred_matrix(siamcat)))
......@@ -74,80 +83,87 @@ evaluate.predictions <- function(siamcat, verbose = 1) {
pred <- pred_matrix(siamcat)[m, , drop = FALSE]
stopifnot(all(names(label$label) == rownames(pred)))
# ##########################################################################
# ROC curve
if (verbose > 2)
message("+ calculating ROC")
auroc = 0
if (ncol(pred) > 1) {
rocc = list(NULL)
aucs = vector("numeric", ncol(pred))
roc.all = list()
auroc.all = vector("numeric", ncol(pred))
for (c in seq_len(ncol(pred))) {
rocc[c] = list(roc(
roc.all[[c]] = roc(
response = label$label,
predictor = pred[, c],
ci = FALSE
))
aucs[c] = rocc[[c]]$auc
)
auroc.all[c] = roc.all[[c]]$auc
}
l.vec = rep(label$label, ncol(pred))
} else {
l.vec = label$label
}
# average data for plotting one mean prediction curve
if (verbose > 2)
message("+ calculating mean ROC")
summ.stat = "mean"
rocsumm = list(roc(
roc.mean = roc(
response = label$label,
predictor = apply(pred, 1, summ.stat),
ci = TRUE,
of = "se",
sp = seq(0, 1, 0.05)
))
auroc = list(rocsumm[[1]]$auc)
# precision recall curve
pr = list(NULL)
ev = list(NULL)
)
auroc = roc.mean$auc
# ##########################################################################
# PR curve
prc = list()
ev = list()
auprc <- 0
if (ncol(pred) > 1) {