Performs statistical testing for differential expression using results from a fitted devil model. Supports both standard and robust (clustered) variance estimation, with multiple testing correction and customizable fold change thresholds.

test_de(
  devil.fit,
  contrast,
  pval_adjust_method = "BH",
  max_lfc = 10,
  clusters = NULL,
  parallel.cores = NULL
)

Arguments

devil.fit

A fitted model object from fit_devil(). Must contain beta coefficients, design matrix, and overdispersion parameters.

contrast

Numeric vector or matrix specifying the comparison of interest. Length must match number of coefficients in the model. For example, c(0, 1, -1) tests difference between second and third coefficient.

pval_adjust_method

Character. Method for p-value adjustment. Passed to stats::p.adjust(). Common choices:

  • "BH": Benjamini-Hochberg (default)

  • "bonferroni": Bonferroni correction

  • "holm": Holm's step-down method

max_lfc

Numeric. Maximum absolute log2 fold change to report. Larger values are capped at ±max_lfc. Default: 10

clusters

Numeric vector or factor. Sample cluster assignments for robust variance estimation. Length must match number of samples. Default: NULL

parallel.cores

Integer or NULL. Number of CPU cores for parallel processing. If NULL, uses all available cores. Default: NULL

Value

A tibble with columns:

name

Character. Gene identifiers from input data

pval

Numeric. Raw p-values from statistical tests

adj_pval

Numeric. Adjusted p-values after multiple testing correction

lfc

Numeric. Log2 fold changes, capped at ±max_lfc

Details

The function implements the following analysis pipeline:

  1. Calculates log fold changes using contrast vectors/matrices

  2. Computes test statistics using either standard or robust variance estimation

  3. Calculates p-values using t-distribution with appropriate degrees of freedom

  4. Adjusts p-values for multiple testing

  5. Applies fold change thresholding

The variance estimation can account for sample clustering (e.g., multiple samples from the same patient) using a sandwich estimator for robust inference.

Examples

if (FALSE) { # \dontrun{
# Basic differential expression test
results <- test_de(fit, contrast = c(0, 1, -1))

# With sample clustering and stricter fold change threshold
results <- test_de(fit, contrast = c(0, 1, -1),
                   clusters = patient_ids,
                   max_lfc = 5)
} # }