rotateStat {randRotation} | R Documentation

## Generate data rotations and calculate statistics on it

### Description

This function generates rotations of data and calculates the provided `statistic` on each rotation and the non-rotated (original) data. This is the central function of the package.

### Usage

    rotateStat(
      initialised.obj,
      R = 10,
      statistic,
      ...,
      parallel = FALSE,
      BPPARAM = BiocParallel::bpparam()
    )

### Arguments

`initialised.obj` | An initialised random rotation object as returned by `initRandrot` and `initBatchRandrot`.
`R` | The number of resamples/rotations. Single `numeric` larger than 1.
`statistic` | A function which takes a data matrix (same dimensions as `Y` - see also `initRandrot`) as first argument and returns a statistic of interest. Any further arguments are passed to it by `...`. We highly recommend using pivotal quantities as `statistic` if possible (see also `Details` in `pFdr`). Note that `pFdr` considers larger values of statistics as more significant, so one-tailed tests may require reversal of the sign and two-tailed tests may require taking absolute values, see `Examples`. The results of `statistic` for each resample are finally combined with `as.matrix` and `cbind`, so ensure that `statistic` returns either a vector or a matrix. Results with multiple columns are possible and handled adequately in subsequent functions (e.g. `pFdr`).
`...` | Further named arguments for `statistic` which are passed unchanged each time it is called. Avoid partial matching to arguments of `rotateStat`. See also the `Examples`.
`parallel` | `logical` if parallel computation should be performed, see details for use of parallel computing.
`BPPARAM` | An optional `BiocParallelParam` instance, see documentation of `BiocParallel` package of Bioconductor.

### Details

The function takes an initialised randrot object (`initRandrot`) and a function that calculates a statistic on the data. The statistic function thereby takes the a matrix `Y` as first argument. Any further arguments are passed to it by `...`.

Together with `pFdr`, this function implements the `workflow` described in (Hettegger et al. 2021).

Be aware that only data is rotated (see also `randrot`), so any additional information including `weights`, `X` etc. need to be provided to `statistic`. See also package vignette and `Examples`.

Parallel processing is implemented with the BiocParallel package of Bioconductor. The default argument `BiocParallel::bpparam()` for `BPPARAM` returns the registered default backend. See package documentation for further information and usage options. If `parallel = TRUE` the function calls in `statistic` need to be called explicitly with package name and "::". So e.g. calling `lmFit` from the `limma` package is done with `limma::lmFit(...)`, see also the examples in the package vignette.

### Value

An object of class `rotateStat`.

### Author(s)

Peter Hettegger

### References

Hettegger P, Vierlinger K, Weinhaeusel A (2021). “Random rotation for identifying differentially expressed genes with linear models following batch effect correction.” Bioinformatics. ISSN 1367-4803, doi: 10.1093/bioinformatics/btab063, https://doi.org/10.1093/bioinformatics/btab063.

### Examples

    #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.2.2 Index]