Permutation-Based FDR Point and Confidence Interval Estimation

Description

FDR functions for permutation-based estimators, including pi0 as well as FDR confidence intervals. The confidence intervals account for dependencies between tests by the incorporation of an overdispersion parameter, which is estimated from the permuted data. The package also includes a parametric analog of the same approach.

Details

This method is designed to compute FDR when a permutation-based approach has been utilized. The objective here is to identify a subset of positive tests that have corresponding statistics with a more exteme distribution than the permuted results, which are assumed to represent the null. The significance of the subset is described in terms of the FDR and uncertainty in the FDR estimate by computing a confidence interval. Say a set of p-values(or simply a set of test statistics) were recorded for a set of hypothesis tests, and data were permuted B times with test results generated for each permutation. The function fdr_od() can be used to estimate FDR and and a confidence interval along with pi0, the proportion of true null hypotheses, given a selected significance threshold. The function fdrTbl()uses fdr_od() to create a table of results over a sequence of possible significance thresholds. Finally, the function FDRplot will plot results from fdrTbl(), facilitating the selection of a final significance threshold.

Author(s)

Joshua Millstein

Maintainer: Joshua Millstein <joshua.millstein@usc.edu> Joshua Millstein

References

Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.

Millstein, J., Battaglin, F., Arai, H., Zhang, W., Jayachandran, P., Soni, S., Parikh, A.R., Mancao, C. and Lenz, H.J., 2022. fdrci: FDR confidence interval selection and adjustment for large-scale hypothesis testing. Bioinformatics advances, 2(1), p.vbac047.

Plot results of FDR table generated by fdrTbl()

Description

This function plots FDR point and CI estimates over a sequence of possible significance thresholds. Results from fdrTbl() can be plotted directly as input to FDRplot.

Usage

FDRplot(
  plotdat,
  lowerbound,
  upperbound,
  ymax = 1,
  annot = "",
  xpos = 0.8,
  ypos = 0.8
)

Arguments

Value

ggplot2 object

Author(s)

Joshua Millstein, joshua.millstein@usc.edu

Joshua Millstein

References

Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.

Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.

Examples

ss = 100
nvar = 100
X = as.data.frame(matrix(rnorm(ss*nvar),nrow=ss,ncol=nvar))
e = as.data.frame(matrix(rnorm(ss*nvar),nrow=ss,ncol=nvar))
Y = .1*X + e
nperm = 10

myanalysis = function(X,Y){
	ntests = ncol(X)
	rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2))
	names(rslts) = c("ID","pvalue")
	rslts[,"ID"] = 1:ntests
	for(i in 1:ntests){
		fit = cor.test(X[,i],Y[,i],na.action="na.exclude",
			alternative="two.sided",method="pearson")
		rslts[i,"pvalue"] = fit$p.value
	}
	return(rslts)
} # End myanalysis

# Generate observed results
obs = myanalysis(X,Y)

# Generate permuted results
perml = vector('list',nperm)
for(perm in 1:nperm){
	X1 = X[order(runif(nvar)),]
	perml[[perm]] = myanalysis(X1,Y)
}

# FDR results table
myfdrtbl = fdrTbl(obs$pvalue,perml,"pvalue",nvar,0,3)
# Plot results
FDRplot(myfdrtbl,0,3,annot="A. An Example")

MV q-values and confidence intervals

Description

q-values with confidence intervals are generated, based in the Millstein and Volfson (MV) estimators.

Usage

MV_q(obsp, permp, pnm, ntests, cl = 0.95, c1 = NA)

Arguments

Details

Millstein and Volfson (2013) FDR is based on the idea that FDR is estimated at a level specified by the investigator. Storey and Tibshirani (2003) developed the q-value concept, where FDR is estimated at each observed p-value. However, Millstein and Volfson argued that in order to be informative, uncertainty in the estimate should be quantified, thus the development of confidence intervals for FDR. The MV FDR estimator is less conservative than the BH estimator.

Value

A dataframe which includes:

Author(s)

Joshua Millstein

References

Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.

Storey, John D., and Robert Tibshirani. "Statistical significance for genomewide studies." Proceedings of the National Academy of Sciences 100.16 (2003): 9440-9445.

Examples

ss=100
nvar=100
X = as.data.frame(matrix(rnorm(ss*nvar),nrow=ss,ncol=nvar))
e = as.data.frame(matrix(rnorm(ss*nvar),nrow=ss,ncol=nvar))
Y = .1*X + e
nperm = 10

myanalysis = function(X,Y){
	ntests = ncol(X)
	rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2))
	names(rslts) = c("ID","pvalue")
	rslts[,"ID"] = 1:ntests
	for(i in 1:ntests){
		fit = cor.test(X[,i],Y[,i],na.action="na.exclude",
			alternative="two.sided",method="pearson")
		rslts[i,"pvalue"] = fit$p.value
	}
	return(rslts)
} # End myanalysis

# Generate observed results
obs = myanalysis(X,Y)

# Generate permuted results
perml = vector('list',nperm)
for(p_ in 1:nperm){
	X1 = X[order(runif(nvar)),]
	perml[[p_]] = myanalysis(X1,Y)
}

q.values.MV = MV_q(obs$pvalue,perml,"pvalue",nvar)

FDR Estimate and Confidence Interval Sequence Table

Description

Computes FDR estimates and confidence intervals for a sequence of potential significance thresholds.

Usage

fdrTbl(
  obs.vec,
  perm.list = NULL,
  pname,
  ntests,
  lowerbound,
  upperbound,
  incr = 0.1,
  cl = 0.95,
  c1 = NA,
  correct = "none",
  meff = TRUE,
  seff = TRUE,
  mymat,
  nperms = 5
)

Arguments

Details

fdrTbl calls fdr_od for a series of discovery thresholds. Output from fdrTbl() can be used for FDRplot() input.

If correct = "BH", then confidence intervals will be corrected according to the thresholds specified by lowerbound, upperbound, and incr. Thresholds will be selected if FDR is determined to be significantly different than 1. First a Z-score test is conducted using the Millstein & Volfson standard error estimate. Then BH FDR is computed according to the Benjamini and Yekutieli (2005) approach. CIs for selected thresholds will be adjusted to account for multiple CI estimation. For thresholds that are not selected, NA values are returned.

Value

A dataframe is returned where rows correspond to p-value thresholds in the sequence from lowerbound to upperbound and columns are:

If permutation: c("threshold","fdr","ll","ul","pi0","odp","S","Sp")

If parametric: c("threshold","fdr","ll","ul","M","M.eff","S","S.eff")

Author(s)

Joshua Millstein, Eric S. Kawaguchi

References

Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.

Benjamini, Yoav, and Daniel Yekutieli. "False discovery rate adjusted multiple confidence intervals for selected parameters." Journal of the American Statistical Association 100.469 (2005): 71-81.

Examples

n.row=100
n.col=100
X = as.data.frame(matrix(rnorm(n.row*n.col),nrow=n.row,ncol=n.col))
e = as.data.frame(matrix(rnorm(n.row*n.col),nrow=n.row,ncol=n.col))
Y = .1*X + e
nperm = 10

myanalysis = function(X,Y){
	ntests = ncol(X)
	rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2))
	names(rslts) = c("ID","pvalue")
	rslts[,"ID"] = 1:ntests
	for(i in 1:ntests){
		fit = cor.test(X[,i],Y[,i],na.action="na.exclude",
			alternative="two.sided",method="pearson")
		rslts[i,"pvalue"] = fit$p.value
	}
	return(rslts)
} # End myanalysis

## Generate observed results
obs = myanalysis(X,Y)

## Generate permuted results
perml = vector('list',nperm)
for(perm in 1:nperm){
	X1 = X[order(runif(n.col)),]
	perml[[perm]] = myanalysis(X1,Y)
}

## FDR results table
fdrTbl(obs$pvalue,perml,"pvalue",n.col,1,2)
fdrTbl(obs$pvalue,perml,"pvalue",n.col,1,2,correct="BH")

## FDR results table (parametric)
fdrTbl(obs$pvalue, NULL, "pvalue",n.col,1,2,meff = TRUE, seff = TRUE, mymat = X, nperms = 5)

Tail-area FDR and Confidence Interval

Description

This function can be used to estimate FDR, corresponding confidence interval, and pi0, the proportion of true null hypotheses, given a selected significance threshold, and results from permuted data.

Usage

fdr_od(
  obsp,
  permp = NULL,
  pnm,
  ntests,
  thres,
  cl = 0.95,
  c1 = NA,
  meff = TRUE,
  seff = TRUE,
  mymat,
  nperms = 5
)

Arguments

Details

If a very large number of tests are conducted, it may be useful to filter results, that is, save only results of those tests that meet some relaxed nominal significance threshold. This alleviates the need to record results for tests that are clearly non-significant. Results from fdr_od() are valid as long as thres < the relaxed nomimal significance threshold for both observed and permuted results. It is not necessary for the input to fdr_od() to be p-values, however, fdr_od() is designed for statistics in which smaller values are more extreme than larger values as is the case for p-values. Therefore, if raw statistics are used, then a transformation may be necessary to insure that smaller values are more likely associated with false null hypotheses than larger values. In certain situations, for instance when a large proportion of tests meet the significance threshold, pi0 is estimated to be very small, and thus has a large influence on the FDR estimate. To limit this influence, pi0 is constrained to be .5 or greater, resulting in a more conservative estimate under these conditions.

Value

For the permutation-based estimator: A list which includes:

For the parametric-based estimator: A list which includes:

Author(s)

Joshua Millstein, Eric S. Kawaguchi

References

Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.

Examples

ss=100
nvar=100
X = as.data.frame(matrix(rnorm(ss*nvar),nrow=ss,ncol=nvar))
e = as.data.frame(matrix(rnorm(ss*nvar),nrow=ss,ncol=nvar))
Y = .1*X + e
nperm = 10

myanalysis = function(X,Y){
	ntests = ncol(X)
	rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2))
	names(rslts) = c("ID","pvalue")
	rslts[,"ID"] = 1:ntests
	for(i in 1:ntests){
		fit = cor.test(X[,i],Y[,i],na.action="na.exclude",
			alternative="two.sided",method="pearson")
		rslts[i,"pvalue"] = fit$p.value
	}
	return(rslts)
} # End myanalysis

# Generate observed results
obs = myanalysis(X,Y)

## Generate permuted results
perml = vector('list',nperm)
for(p_ in 1:nperm){
	X1 = X[order(runif(nvar)),]
	perml[[p_]] = myanalysis(X1,Y)
}

## FDR results (permutation)
fdr_od(obs$pvalue,perml,"pvalue",nvar, thres = .05)

## FDR results (parametric)
fdr_od(obs$pvalue, permp = NULL, thres = 0.05, meff = FALSE, seff = FALSE, mymat = X)

Estimate the Effective Number of Tests

Description

Estimate the effective number of tests using a permutation-based approach.

Usage

meff.jm(mydat, B = 1)

Arguments

Details

The effective no of independent vars is the sum of the scaled eigenvalues for truly independent vars (n) minus the variance inflation of the top eigenvalues, that is, (observed - permuted) The sum of the eigenvalues of the correlation matrix = trace = sum of diag = n. The eigenvalues of variables X can be thought of as normalized Dependencies between variables will increase the magnitude of top eigenvalues, therefore, the effective number of independent variables is equal to the proportion of the sum of eigenvalues attributed to independence minus the proportion attributed to dependencies.

Value

Numeric. Returns the estimated effective number of tests averaged over B permutations.

Author(s)

Joshua Millstein, Eric S. Kawaguchi

References

Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.

Examples

# Independent
ss=300
nvar=100
X = as.data.frame(matrix(rnorm(ss * nvar), nrow = ss, ncol = nvar))
meff.jm(X, B = 5)

# High correlation
S = matrix(0.9, nvar, nvar)
diag(S) = 1
X = as.matrix(X) %*% chol(S)
meff.jm(X, B = 5)