pFdr {randRotation} | R Documentation

## Calculate resampling based p-values and FDRs

### Description

This function calculates either (1) resampling based p-values with subsequent p-value adjustment using `stats::p.adjust` or (2) resampling based false-discovery-rates (FDRs) for rotated statistics from a `rotateStat` object.

### Usage

    pFdr(obj, method = "none", pooled = TRUE, na.rm = FALSE, beta = 0.05)

### Arguments

`obj` | A `rotateStat` object as returned by `rotateStat`.
`method` | Can be either `"none"` (default), `"fdr.q"`, `"fdr.qu"` or any term that can be passed as `method` argument to `stats::p.adjust`, see `Details`. If `method = "none"`, resampling based p-values without further adjustment are calculated.
`pooled` | `logical`. `TRUE` (default) if marginal distributions are exchangeable for all features so that rotated stats can be pooled, see `Details`.
`na.rm` | `logical`. `NA` values are ignored if set `TRUE`. `NA` values should be avoided and could e.g. be removed by imputation in original data or by removing features that contain `NA` values. Few `NA` values do not have a large effect, but many `NA` values can lead to wrong estimations of p-values and FDRs. We highly recommend avoiding `NA` values.
`beta` | `numeric` between `0` and `1`. Corresponds to beta in (Yekutieli and Benjamini 1999).

### Details

Larger values of obj$s0 are considered more significant when compared to the empirical distribution. E.g. for calculation of resampling based p-values (with `pooled = FALSE`) we in principle use `p.val <- (rowSums(obj$stats >= obj$s0)+1)/(ncol(obj$stats)+1)` according to (Phipson and Smyth 2010).

`method = "fdr.q"` and `method = "fdr.qu"` are resampling based fdr estimates and can only be used with `pooled = TRUE`. `method = "fdr.q"` is the FDR local estimator and `method = "fdr.qu"` is the FDR upper limit, see (Reiner et al. 2003; Yekutieli and Benjamini 1999). For all other `method` arguments resampling based p-values are calculated and passed to `stats::p.adjust` for p-value adjustment. So these methods provide resampling based p-values with (non-resampling based) p-value adjustment. `method = "fdr.q"` and `method = "fdr.qu"` were adapted from package `fdrame` (Benjamini et al. 2019; Reiner et al. 2003).

When `pooled = TRUE`, marginal distributions of the test statistics are considered exchangeable for all features. The resampling based p-values of each feature are then calculated from all rotated statistics (all features, all rotations). For these cases, if the number of features is reasonably large, usually only few resamples (argument `R` in `rotateStat`) are required. We want to emphasize that in order for the marginal distributions to be exchangeable, the statistics must be a pivotal quantity (i.e. it must be scale independent). Pivotal quantities are e.g. t values. Using e.g. linear models with `coef` as statistics is questionable if the different features are measured on different scales. The resampled coefficients then have different variances and `pooled = TRUE` is not applicable. We thus highly recommend using pivotal quantities as `statistics` in `rotateStat` if possible.

When `pooled = FALSE` the resampling based p-values are calculcated for each feature separately. This is required if one expects the resampling based statistics to be distributed differently for individual features. For most common applications this should not be the case and the marginal distribution are exchangeable for all features, hence `pooled = TRUE` by default.

If `method = "fdr.q"` or `method = "fdr.qu"` and `weights` were specified when initialising the random rotation object (see parameter `initialised.obj` in `rotateStat`), a warning is displayed. The correlation structure (dependence structure) of linear model coefficients between different features is not generally preserved if different weights are used for different features. Methods `fdr.q` and `fdr.qu` rely on preserved correlation structure of dependent statistics and thus should not be used if statistics based on model coefficients (e.g. t statistics of model coefficients) are used in combination with different weights.

P-values and FDRs are calculated for each column of `obj$s0` separately.

### Value

A `numeric` matrix of corrected p-values or FDRs with dimension `dim(obj$s0)`.

### Author(s)

Peter Hettegger

### References

Benjamini Y, Kenigsberg E, Reiner A, Yekutieli D (2019). fdrame: FDR adjustments of Microarray Experiments (FDR-AME). R package version 1.56.0.
Phipson B, Smyth GK (2010). “Permutation P-values should never be zero: Calculating exact P-values when permutations are randomly drawn.” Statistical Applications in Genetics and Molecular Biology, 9(1). ISSN 15446115, doi: 10.2202/1544-6115.1585, 1603.05766, http://www.ncbi.nlm.nih.gov/pubmed/21044043.
Reiner A, Yekutieli D, Benjamini Y (2003). “Identifying differentially expressed genes using false discovery rate controlling procedures.” Bioinformatics, 19(3), 368–375. ISSN 13674803, doi: 10.1093/bioinformatics/btf877, http://www.ncbi.nlm.nih.gov/pubmed/12584122.
Yekutieli D, Benjamini Y (1999). “Resampling-based false discovery rate controlling multiple test procedures for correlated test statistics.” Journal of Statistical Planning and Inference, 82(1-2), 171–196. ISSN 03783758, doi: 10.1016/S0378-3758(99)00041-5, https://doi.org/10.1016/S0378-3758(99)00041-5.

### See Also

`rotateStat`

### Examples

    # See also '?rotateStat':

    #set.seed(0)

    # Dataframe of phenotype data (sample information)
    # We simulate 2 sample classes processed in 3 batches
    pdata <- data.frame(batch = rep(1:3, c(10,10,10)),
                       phenotype = rep(c("Control", "Cancer"), c(5,5)))
    features <- 100

    # Matrix with random gene expression data
    edata <- matrix(rnorm(features * nrow(pdata)), features)
    rownames(edata) <- paste("feature", 1:nrow(edata))

    mod1 <- model.matrix(~phenotype, pdata)

    # Initialisation of the random rotation class
    init1 <- initBatchRandrot(Y = edata, X = mod1, coef.h = 2, batch = pdata$batch)
    init1

    # Definition of the batch effect correction procedure with subsequent calculation
    # of two-sided test statistics
    statistic <- function(., batch, mod, coef){

      # The "capture.output" and "suppressMessages" simply suppress any output
      capture.output(suppressMessages(
        Y.tmp <- sva::ComBat(., batch = batch, mod)
      ))

      fit1 <- lm.fit(mod, t(Y.tmp))
      abs(coef(fit1)[coef,])
    }

    # We calculate test statistics for the second coefficient

    res1 <- rotateStat(initialised.obj = init1,
                        R = 10,
                        statistic = statistic,
                        batch = pdata$batch, mod = mod1, coef = 2)

    hist(pFdr(res1))

[Package randRotation version 1.6.0 Index]