Plot Per-test Type I Error Rate from Null Simulations

Visualizes the per-test type I error rate of differential splicing methods on simulated null data (DID = 0 for all orthologous splice junctions). For each simulation, the per-test type I error rate is defined as the proportion of raw (unadjusted) p-values strictly below alpha. Because no orthologous splice junction has a true effect, every rejection is a false positive and this proportion estimates the per-test type I error under the global null.

Usage

plotSimu(
  params,
  alpha = 0.1,
  nbThreads = 1,
  methodOrder = c("DESeq2", "limma", "limmaCorr", "phylolmBM", "phylolmOU", "EVE",
    "phyloDEBM", "phyloDEOU"),
  ncol = 3,
  textSize = 22,
  breaks = NULL,
  returnData = FALSE
)

Arguments

params

data.frame. Simulation results metadata returned by runSimu, containing columns:

modelProcess: Evolution model ("BM" or "OU")
lambda: Within-species correlation (0 to 1)
thalf: OU half-life, relative to tree height (between-species correlation)
method: Statistical method name
rep: Replicate ID
outputResFile: Path to TSV file with pval and DID columns

alpha

Numeric. Significance threshold applied to raw (unadjusted) p-values. Orthologous splice junctions with pval < alpha are counted as rejections. Under the global null, the expected per-test type I error rate is alpha. Default: 0.1.

nbThreads

Integer. Number of parallel threads for reading result files. Default: 1.

methodOrder

Character vector. Order in which methods are displayed in facets. Default order groups non-phylogenetic methods first, then phylogenetic methods.

ncol

Integer. Number of columns in facet layout. Default: 3.

textSize

Numeric. Base text size for plot theme. Default: 22.

breaks

Numeric vector. Explicit break points for the colour scale. If NULL (default), breaks are chosen automatically by ggplot2::waiver().

returnData

Logical. If TRUE, return the per-replicate data.frame instead of a plot. Default: FALSE.

Value

If returnData = FALSE (default), a ggplot2 object showing a heatmap of the mean per-test type I error rate across lambda (y-axis) and thalf (x-axis), faceted by method. The nominal level alpha is used as the lower bound of the colour scale, so cells below alpha appear white.

If returnData = TRUE, a data.frame with columns:

rep: Replicate ID
modelProcess: Evolution model
lambda: Within-species correlation (factor)
thalf: OU half-life, relative to tree height (factor)
method: Statistical method (factor)
FP: Number of raw p-values < alpha (false positive count)
typeIerrRate: Per-test type I error rate: FP / n_junctions, i.e. the proportion of raw p-values below alpha under the global null

Details

To quantify false positive inflation induced by phylogenetic dependence, the per-test type I error rate is evaluated using raw p-values under the global null. The question addressed here is fundamentally a model misspecification problem: non-phylogenetic methods ignore the covariance structure induced by shared evolutionary history, so their p-values are not uniform under the null when within- or between-species correlations are strong.

The function performs the following steps:

Reads each simulation result file and counts orthologous splice junctions with raw (unadjusted) p-values strictly below alpha as false positives (FP). All data are generated under the global null (DID = 0), so every rejection is a false positive.
Computes the per-test type I error rate as typeIerrRate = FP / n_junctions for each simulation replicate. Under a well-calibrated test this should equal alpha.
Aggregates results by taking the mean per-test type I error rate across replicates for each (lambda, thalf, method) combination.
Generates a heatmap with a log-scaled colour gradient anchored at the nominal level alpha.

Interpreting the Plot

White/light cells: Per-test type I error rate near or below alpha — well-calibrated test.
Red cells: Per-test type I error rate well above alpha — inflated type I error.

Examples

if (FALSE) { # \dontrun{
# Run simulations first
results <- runSimu(data, se, simDir = "sim/", methods = c("limma", "phyloDEOU"))

# Plot per-test type I error rate heatmap
plotSimu(results)

# Get raw data for custom analysis
fp_data <- plotSimu(results, returnData = TRUE)
} # }