inst/doc/reconsiVignette.R

In reconsi: Resampling Collapsed Null Distributions for Simultaneous Inference

## ----installBioConductor, eval = FALSE----------------------------------------
#  library(BiocManager)
#  BiocManager::install("reconsi")

## ----installAndLoadGitHub, eval = FALSE---------------------------------------
#  library(devtools)
#  install_github("CenterForStatistics-UGent/reconsi")

## ----loadReconsi--------------------------------------------------------------
suppressPackageStartupMessages(library(reconsi))
cat("reconsi package version", as.character(packageVersion("reconsi")), "\n")

## ----syntData-----------------------------------------------------------------
#Create some synthetic data with 90% true null hypothesis
 p = 200; n = 50
x = rep(c(0,1), each = n/2)
 mat = cbind(
 matrix(rnorm(n*p/10, mean = 5+x),n,p/10), #DA
 matrix(rnorm(n*p*9/10, mean = 5),n,p*9/10) #Non DA
 )
 #Provide just the matrix and grouping factor, and test using the collapsed null
 fdrRes = reconsi(mat, x)
  #The estimated tail-area false discovery rates.
  estFdr = fdrRes$Fdr

## ----p0-----------------------------------------------------------------------
fdrRes$p0

## ----plotNull-----------------------------------------------------------------
plotNull(fdrRes)

## ----plotApproxCovar----------------------------------------------------------
plotApproxCovar(fdrRes)

## ----customFunction-----------------------------------------------------------
 #With a custom function, here linera regression
fdrResLm = reconsi(mat, x, B = 5e1,
                      test = function(x, y){
fit = lm(y~x)
c(summary(fit)$coef["x","t value"], fit$df.residual)},
distFun = function(q){pt(q = q[1], df = q[2])})

## ----customFunction2----------------------------------------------------------
 #3 groups
 p = 100; n = 60
x = rep(c(0,1,2), each = n/3)
mu0 = 5
 mat = cbind(
 matrix(rnorm(n*p/10, mean = mu0+x),n,p/10), #DA
 matrix(rnorm(n*p*9/10, mean = mu0),n,p*9/10) #Non DA
 )
 #Provide an additional covariate through the 'argList' argument
 z = rpois(n , lambda = 2)
 fdrResLmZ = reconsi(mat, x, B = 5e1,
 test = function(x, y, z){
 fit = lm(y~x+z)
 c(summary(fit)$coef["x","t value"], fit$df.residual)},
distFun = function(q){pt(q = q[1], df = q[2])},
 argList = list(z = z))

## ----kruskal------------------------------------------------------------------
fdrResKruskal = reconsi(mat, x, B = 5e1,
test = function(x, y){kruskal.test(y~x)$statistic}, zValues = FALSE)

## ----resamZvals---------------------------------------------------------------
fdrResKruskalPerm = reconsi(mat, x, B = 5e1,
test = function(x, y){
 kruskal.test(y~x)$statistic}, resamZvals = TRUE)

## ----bootstrap----------------------------------------------------------------
fdrResBootstrap = reconsi(Y = mat, B = 5e1, test = function(y, x, mu){
                                      testRes = t.test(y, mu = mu)
                                      c(testRes$statistic, testRes$parameter)}, argList = list(mu = mu0),
                                  distFun = function(q){pt(q = q[1],
                                                           df = q[2])})

## ----Vandeputte---------------------------------------------------------------
#The grouping and flow cytometry variables are present in the phyloseq object, they only need to be called by their name.
testVanDePutte = testDAA(Vandeputte, groupName = "Health.status", FCname = "absCountFrozen", B = 1e2L)

## ----vandeputteFdr------------------------------------------------------------
FdrVDP = testVanDePutte$Fdr
quantile(FdrVDP)

## ----sessionInfo--------------------------------------------------------------
sessionInfo()