Receiver operator characteristics of simulation results

This function takes the simulation output from simulateDE and computes receiver operator characteristics (ROC) and area under the curve values (AUC) with the help of functions implemented in calculate_performance.

evaluateROC(simRes,
alpha.type=c("adjusted","raw"),
MTC=c('BY', 'BH', 'Storey', 'IHW',
'holm', 'hochberg', 'hommel', 'bonferroni'),
alpha.nominal = 0.1,
target.by=c("lfc", "effectsize"), delta=0,
raw = FALSE, verbose = TRUE)

Arguments

simRes	The result from simulateDE.
alpha.type	A string to represent the way to call DE genes. Available options are "adjusted" i.e. applying multiple testing correction and "raw" i.e. using p-values. Default is "adjusted".
MTC	Multiple testing correction method to use. Available options are 1) see p.adjust.methods, 2) Storey's qvalue see qvalue and 3) Independent Hypothesis Weighting considering mean expression as covariate (see ihw). Default is BY, i.e. Benjamini-Yekutieli FDR correction method.
alpha.nominal	The nomial value of significance. Default is 0.1.
target.by	A string to specify the method to define "biologically important" DE genes. Available options are (1) "lfc": interesting genes are defined by absolute log2 fold changes. (2) "effectsize": interesting genes are defined by absolute log2 fold changes divided by the square root of 1/(mean+dispersion).
delta	A threshold used for defining "biologically important" genes. Genes with absolute log2 fold changes (when target.by is "lfc") or effect sizes (when target.by is "effectsize") greater than this value are deemed DE in error rates calculations. If the default delta=0, then no threshold is applied.
raw	Logical vector. The default FALSE removes intermediate calculations from the output since they can be quite large.
verbose	Logical vector to indicate whether to show progress report of evaluation. Default is TRUE.

Value

A list with the following entries:

Performances

The output of calculate_performance of aspect "fdrtprcurve" calculating the proportions of TN, FN, TP and FP and related rates which uses prediction and performance internally.

TPRvsFDR

The output of calculate_performance of aspect "fdrtpr" calculating the proportions of TN, FN, TP and FP and associated TPR and observed FDR for each nominal level from 0.01 to 1 in steps of 0.01.

Scores

The mean +/- standard error of performance measures and area under the curve. These are calculated once for conservative (extension "conv") and once for liberal nominal (extension "lib") FDR levels. Measures include: accuracy (ACC), F-measure of precision and recall (F1score), Matthew's correlation coefficient (MCC), partial Area under the Curve of TPR vs FDR (TPRvsFDR_pAUC), area under the ROC-curve (TPRvsFPR_AUC) and area under the PR-curve (TPRvsPPV_AUC).

Raw

If raw is TRUE, then the intermediate calculations of the above three list entries is also included in the output.

Settings

Reiterating chosen evaluation parameters.

See also

estimateParam for parameter estimation needed for simulations, Setup for setting up simulation parameters and simulateDE for simulating differential expression.

Examples

if (FALSE) {
# estimate gene parameters
data("SmartSeq2_Gene_Read_Counts")
estparam_gene <- estimateParam(countData = SmartSeq2_Gene_Read_Counts,
                               readData = NULL,
                              batchData = NULL,
                               spikeData = NULL, spikeInfo = NULL,
                               Lengths = NULL, MeanFragLengths = NULL,
                               RNAseq = 'singlecell', Protocol = 'Read',
                               Distribution = 'ZINB', Normalisation = "scran",
                               GeneFilter = 0.1, SampleFilter = 3,
                               sigma = 1.96, NCores = NULL, verbose = TRUE)
# define log2 fold change
p.lfc <- function(x) sample(c(-1,1), size=x,replace=T)*rgamma(x, shape = 1, rate = 2)
# set up simulations
setupres <- Setup(ngenes = 10000, nsims = 10,
                  p.DE = 0.1, pLFC = p.lfc,
                  n1 = c(20,50,100), n2 = c(30,60,120),
                  Thinning = c(1,0.9,0.8), LibSize = 'given',
                  estParamRes = estparam_gene,
                  estSpikeRes = NULL,
                  DropGenes = TRUE,
                  sim.seed = 83596, verbose = TRUE)
# run simulation
simres <- simulateDE(SetupRes = setupres,
                     Prefilter = "FreqFilter", Imputation = NULL,
                     Normalisation = 'scran', Label = 'none',
                     DEmethod = "limma-trend", DEFilter = FALSE,
                     NCores = NULL, verbose = TRUE)
# evaluation
evalrocres <- evaluateROC(simRes = simres,
                          alpha.type = "adjusted",
                          MTC = 'BH', alpha.nominal = 0.05,
                          raw = FALSE)
# plot evaluation
plotEvalROC(evalRes = evalrocres)
}