| Title: | Resampling-based multiple hypothesis testing |
|---|---|
| Description: | Non-parametric bootstrap and permutation resampling-based multiple testing procedures (including empirical Bayes methods) for controlling the family-wise error rate (FWER), generalized family-wise error rate (gFWER), tail probability of the proportion of false positives (TPPFP), and false discovery rate (FDR). Several choices of bootstrap-based null distribution are implemented (centered, centered and scaled, quantile-transformed). Single-step and step-wise methods are available. Tests based on a variety of t- and F-statistics (including t-statistics based on regression parameters from linear and survival models as well as those based on correlation parameters) are included. When probing hypotheses with t-statistics, users may also select a potentially faster null distribution which is multivariate normal with mean zero and variance covariance matrix derived from the vector influence function. Results are reported in terms of adjusted p-values, confidence regions and test statistic cutoffs. The procedures are directly applicable to identifying differentially expressed genes in DNA microarray experiments. |
| Authors: | Katherine S. Pollard, Houston N. Gilbert, Yongchao Ge, Sandra Taylor, Sandrine Dudoit |
| Maintainer: | Katherine S. Pollard <[email protected]> |
| License: | LGPL |
| Version: | 2.63.0 |
| Built: | 2024-10-30 08:57:42 UTC |
| Source: | https://github.com/bioc/multtest |
A user-level function to perform empirical Bayes multiple testing procedures (EBMTP). A variety of t- and F-tests, including robust versions of most tests, are implemented. A common-cutoff method is used to control the chosen type I error rate (FWER, gFWER, TPPFP, or FDR). Bootstrap-based null distributions are available. Additionally, for t-statistics, one may wish to sample from an appropriate multivariate normal distribution with mean zero and correlation matrix derived from the vector influence function. In EBMTP, realizations of local q-values, obtained via density estimation, are used to partition null and observed test statistics into guessed sets of true and false null hypotheses at each round of (re)sampling. In this manner, parameters of any type I error rate which can be expressed as a function the number of false positives and true positives can be estimated. Arguments are provided for user control of output. Gene selection in microarray experiments is one application.
EBMTP(X, W = NULL, Y = NULL, Z = NULL, Z.incl = NULL, Z.test = NULL, na.rm = TRUE, test = "t.twosamp.unequalvar", robust = FALSE, standardize = TRUE, alternative = "two.sided", typeone = "fwer", method = "common.cutoff", k = 0, q = 0.1, alpha = 0.05, smooth.null = FALSE, nulldist = "boot.cs", B = 1000, psi0 = 0, marg.null = NULL, marg.par = NULL, ncp = NULL, perm.mat = NULL, ic.quant.trans = FALSE, MVN.method = "mvrnorm", penalty = 1e-06, prior = "conservative", bw = "nrd", kernel = "gaussian", seed = NULL, cluster = 1, type = NULL, dispatch = NULL, keep.nulldist = TRUE, keep.rawdist = FALSE, keep.falsepos = FALSE, keep.truepos = FALSE, keep.errormat = FALSE, keep.Hsets=FALSE, keep.margpar = TRUE, keep.index = FALSE, keep.label = FALSE)EBMTP(X, W = NULL, Y = NULL, Z = NULL, Z.incl = NULL, Z.test = NULL, na.rm = TRUE, test = "t.twosamp.unequalvar", robust = FALSE, standardize = TRUE, alternative = "two.sided", typeone = "fwer", method = "common.cutoff", k = 0, q = 0.1, alpha = 0.05, smooth.null = FALSE, nulldist = "boot.cs", B = 1000, psi0 = 0, marg.null = NULL, marg.par = NULL, ncp = NULL, perm.mat = NULL, ic.quant.trans = FALSE, MVN.method = "mvrnorm", penalty = 1e-06, prior = "conservative", bw = "nrd", kernel = "gaussian", seed = NULL, cluster = 1, type = NULL, dispatch = NULL, keep.nulldist = TRUE, keep.rawdist = FALSE, keep.falsepos = FALSE, keep.truepos = FALSE, keep.errormat = FALSE, keep.Hsets=FALSE, keep.margpar = TRUE, keep.index = FALSE, keep.label = FALSE)
For brevity, the presentation of arguments below will highlight those which differ significantly from arguments in the other main-level user function MTP. See MTP for further details.
typeone |
Character string indicating which type I error rate to control, by default family-wise error rate ('fwer'). Other options include generalized family-wise error rate ('gfwer'), with parameter |
method |
Character string indicating the EBMTP method. Currently only 'common.cutoff' is implemented. This method is most similar to 'ss.maxT' in |
nulldist |
Character string indicating which resampling method to use for estimating the joint test statistics null distribution, by default the non-parametric bootstrap with centering and scaling ('boot.cs'). The old default 'boot' will still compile and will correspond to 'boot.cs'. Other null distribution options include 'boot.ctr', 'boot.qt', and 'ic', corresponding to the centered-only bootstrap distribution, quantile-transformed bootstrap distribution, and influence curve multivariate normal joint null distribution, respectively. The permutation distribution is not available. |
prior |
Character string indicating which choice of prior probability to use for estimating local q-values (i.e., the posterior probabilities of a null hypothesis being true given the value of its corresponding test statistic). Default is 'conservative', in which case the prior is set to its most conservative value of 1, meaning that all hypotheses are assumed to belong to the set of true null hypotheses. Other options include 'ABH' for the adaptive Benjamini-Hochberg estimator of the number/proportion of true null hypotheses, and 'EBLQV' for the empirical Bayes local q-value value estimator of the number/proportion of true null hypotheses. If 'EBLQV', the estimator of the prior probability is taken to be the sum of the estimated local q-values divided by the number of tests. Relaxing the prior may result in more rejections, albeit at a cost of type I error control under certain conditions. See details and references. |
bw |
A character string argument to |
kernel |
A character string argument to |
keep.falsepos |
A logical indicating whether or not to store the matrix of guessed false positives at each round of (re)sampling. The matrix has rows equal to the number of cut-offs (observed test statistics) and columns equal to the |
keep.truepos |
A logical indicating whether or not to store the matrix of guessed true positives at each round of (re)sampling. The matrix has rows equal to the number of cut-offs (observed test statistics) and columns equal to the |
keep.errormat |
A logical indicating whether or not to store the matrix of type I error rate values at each round of (re)sampling. The matrix has rows equal to the number of cut-offs (observed test statistics) and columns equal to the |
keep.Hsets |
A logical indicating whether or not to return the matrix of indicators which partition the hypotheses into guessed sets of true and false null hypotheses at each round of (re)sampling. Default is 'FALSE'. |
X, W, Y, Z, Z.incl, Z.test, na.rm, test, robust, standardize, alternative, k, q, alpha, smooth.null, B, psi0, marg.null, marg.par, ncp, perm.mat, ic.quant.trans, MVN.method, penalty, seed, cluster, type, dispatch, keep.nulldist, keep.rawdist, keep.margpar, keep.index, keep.label
|
These arguments are all similarly used by the |
The EBMTP begins with a marginal nonparametric mixture model for estimating local q-values. By definition, q-values are 'the opposite' of traditional p-values. That is, q-values represent the probability of null hypothesis being true given the value of its corresponding test statistic. If the test statistics Tn have marginal distribution f = pi*f_0 + (1-pi)f_1, where pi is the prior probability of a true null hypothesis and f_0 and f_1 represent the marginal null and alternative densities, respectively, then the local q-value function is given by pi*f_0(Tn)/f(Tn).
One can estimate both the null density f_0 and full density f by applying kernel density estimation over the matrix of null test statistics and the vector of observed test statistics, respectively. Practically, this step in EBMTP also ensures that sidedness is correctly accounted for among the test statistics and their estimated null distribution. The prior probability pi can be set to its most conservative value of 1 or estimated by some other means, e.g., using the adaptive Benjamini Hochberg ('ABH') estimator or by summing up the estimated local q-values themselves ('EBLQV')and dividing by the number of tests. Bounding these estimated probabilities by one provides a vector of estimated local q-values with length equal to the number of hypotheses. Bernoulli 0/1 realizations of the posterior probabilities indicate which hypotheses are guessed as belonging to the true set of null hypotheses given the value of their test statistics. Once this partitioning has been achieved, one can count the numbers of guessed false positives and guessed true positives at each round of (re)sampling that are obtained when using the value of an observed test statistic as a cut-off.
EBMTPs use function closures to represent type I error rates in terms of their defining features. Restricting the choice of type I error rate to 'fwer', 'gfwer', 'tppfp', and 'fdr', means that these features include whether to control the number of false positives or the proportion of false positives among the number of rejetions made (i.e., the false discovery proportion), whether we are controlling a tail probability or expected value error rate, and, in the case of tail probability error rates, what bound we are placing on the random variable defining the type I error rate (e.g., k for 'gfwer' or 'q' for 'tppfp'). Averaging the type I error results over B (bootstrap or multivariate normal) samples provides an estimator of the evidence against the null hypothesis (adjusted p-values) with respect to the choice of type I error rate. Finally, a monotonicity constraint is placed on the adjusted p-values before being returned as output.
As detailed in the references, relaxing the prior may result in a more powerful multiple testing procedure, albeit sometimes at the cost of type I error control. Additionally, when the proportion of true null hypotheses is close to one, type I error control may also become an issue, even when using the most conservative prior probability of one. This feature is known to occur with some other procedures which rely on the marginal nonparametric mixture model for estimating (local) q-values. The slot EB.h0M returned by objects of class EBMTP is the sum of the local q-values estimated via kernel density estimation (divided by the total number of tests). If this value is close to one (>0.9-0.95), the user will probably not want relax the prior, as even the conservative EBMTP might be approaching a performance bound with respect to type I error control. The user is advised to begin by using the most 'conservative' prior, assess the estimated proportion of true null hypotheses, and then decide if relaxing the prior might be desired. Gains in power over other multiple testing procedures have been observed even when using the most conservative prior of one.
Situations of moderate-high to high levels of correlation may also affect the results of multiple testing methods which use the same mixture model for generating q-values. Microarray analysis represents a scenario in which dependence structures are typically weak enough to mitigate this concern. On the other hand, the analysis of densely sampled SNPs, for example, may present problems.
An object of class EBMTP. Again, for brevity, the values below represent slots which distinguish objects of class EBMTP from those of class MTP.
falsepos |
A matrix with rows equal to the number of hypotheses and columns the number of samples of null test statistics ( |
truepos |
A matrix with rows equal to the number of hypotheses and columns the number of samples of null test statistics ( |
errormat |
The matrix obtained after applying to type I error rate function closure to the matrices in |
EB.h0M |
The sum of the local q-values obtained after density estimation. This number serves as an estimate of the proportion of true null hypotheses. Values close to one indicate situations in which type I error control may not be guaranteed by the EBMTP. When |
prior |
The numeric value of the prior 'pi' used when evaluating the local q-value function. |
prior.type |
Character string returning the value of |
lqv |
A numeric vector of length the number of hypotheses with the estimated local q-values used for generating guessed sets of true null hypotheses. |
Hsets |
A numeric matrix with the same dimension as |
Houston N. Gilbert, based on the original MTP code written by Katherine S. Pollard
H.N. Gilbert, K.S. Pollard, M.J. van der Laan, and S. Dudoit (2009). Resampling-based multiple
hypothesis testing with applications to genomics: New developments in R/Bioconductor
package multtest. Journal of Statistical Software (submitted). Temporary URL: http://www.stat.berkeley.edu/~houston/JSSNullDistEBMTP.pdf.
Y. Benjamini and Y. Hochberg (2000). On the adaptive control of the false
discovery rate in multiple testing with independent statistics. J. Behav.
Educ. Statist. Vol 25: 60-83.
Y. Benjamini, A. M. Krieger and D. Yekutieli (2006). Adaptive linear step-up
procedures that control the false discovery rate. Biometrika.
Vol. 93: 491-507.
M.J. van der Laan, M.D. Birkner, and A.E. Hubbard (2005). Empirical Bayes and Resampling Based Multiple Testing Procedure Controlling the Tail Probability of the Proportion of False Positives. Statistical Applications in Genetics and Molecular Biology, 4(1).
http://www.bepress.com/sagmb/vol4/iss1/art29/
S. Dudoit and M.J. van der Laan. Multiple Testing Procedures and Applications to Genomics. Springer Series in Statistics. Springer, New York, 2008.
S. Dudoit, H.N. Gilbert, and M J. van der Laan (2008).
Resampling-based empirical Bayes multiple testing procedures for controlling
generalized tail probability and expected value error rates: Focus on the false
discovery rate and simulation study. Biometrical Journal, 50(5):716-44. http://www.stat.berkeley.edu/~houston/BJMCPSupp/BJMCPSupp.html.
H.N. Gilbert, M.J. van der Laan, and S. Dudoit. Joint multiple testing procedures for
graphical model selection with applications to biological networks. Technical report,
U.C. Berkeley Division of Biostatistics Working Paper Series, April 2009. URL http://www.bepress.com/ucbbiostat/paper245.
MTP, EBMTP-class, EBMTP-methods, Hsets
set.seed(99) data<-matrix(rnorm(90),nr=9) group<-c(rep(1,5),rep(0,5)) #EB fwer control with centered and scaled bootstrap null distribution #(B=100 for speed) eb.m1<-EBMTP(X=data,Y=group,alternative="less",B=100,method="common.cutoff") print(eb.m1) summary(eb.m1) par(mfrow=c(2,2)) plot(eb.m1,top=9)set.seed(99) data<-matrix(rnorm(90),nr=9) group<-c(rep(1,5),rep(0,5)) #EB fwer control with centered and scaled bootstrap null distribution #(B=100 for speed) eb.m1<-EBMTP(X=data,Y=group,alternative="less",B=100,method="common.cutoff") print(eb.m1) summary(eb.m1) par(mfrow=c(2,2)) plot(eb.m1,top=9)
An object of class EBMTP is the output of a particular multiple testing procedure, as generated by the function EBMTP. The object has slots for the various data used to make multiple testing decisions, in particular adjusted p-values.
Objects can be created by calls of the form
new('MTP',
statistic = ...., object of class numeric
estimate = ...., object of class numeric
sampsize = ...., object of class numeric
rawp = ...., object of class numeric
adjp = ...., object of class numeric
reject = ...., object of class matrix
rawdist = ...., object of class matrix
nulldist = ...., object of class matrix
nulldist.type = ...., object of class character
marg.null = ...., object of class character
marg.par = ...., object of class matrix
label = ...., object of class numeric
falsepos = ...., object of class matrix
truepos = ...., object of class matrix
errormat = ...., object of class matrix
EB.h0M = ...., object of class numeric
prior = ...., object of class numeric
prior.type= ...., object of class character
lqv = ...., object of class numeric
Hsets = ...., object of class matrix
index = ...., object of class matrix
call = ...., object of class call
seed = ...., object of class integer
)
statisticObject of class numeric, observed test statistics for each hypothesis, specified by the values of the MTP arguments test, robust, standardize, and psi0.
estimateFor the test of single-parameter null hypotheses using t-statistics (i.e., not the F-tests), the numeric vector of estimated parameters corresponding to each hypothesis, e.g. means, differences in means, regression parameters.
sampsizeObject of class numeric, number of columns (i.e. observations) in the input data set.
rawpObject of class numeric, unadjusted, marginal p-values for each hypothesis.
adjpObject of class numeric, adjusted (for multiple testing) p-values for each hypothesis (computed only if the get.adjp argument is TRUE).
rejectObject of class 'matrix', rejection indicators (TRUE for a rejected null hypothesis), for each value of the nominal Type I error rate alpha.
rawdistThe numeric matrix for the estimated nonparametric non-null test statistics distribution (returned only if keep.rawdist=TRUE and if nulldist is one of 'boot.ctr', 'boot.cs', or 'boot.qt'). This slot must not be empty if one wishes to call update to change choice of bootstrap-based null distribution.
nulldistThe numeric matrix for the estimated test statistics null distribution (returned only if keep.nulldist=TRUE). By default (i.e., for nulldist='boot.cs'), the entries of nulldist are the null value shifted and scaled bootstrap test statistics, with one null test statistic value for each hypothesis (rows) and bootstrap iteration (columns).
nulldist.typeCharacter value describing which choice of null distribution was used to generate the MTP results. Takes on one of the values of the original nulldist argument in the call to MTP, i.e., 'boot.cs', 'boot.ctr', 'boot.qt', or 'ic'.
marg.nullIf nulldist='boot.qt', a character value returning which choice of marginal null distribution was used by the MTP. Can be used to check default values or to ensure manual settings were correctly applied.
marg.parIf nulldist='boot.qt', a numeric matrix returning the parameters of the marginal null distribution(s) used by the MTP. Can be used to check default values or to ensure manual settings were correctly applied.
falseposA matrix with rows equal to the number of hypotheses and columns the number of samples of null test statistics (B) indicating the number of guessed false positives when using the corresponding value of the observed test statistic as a cut-off. Not returned unless keep.falsepos=TRUE.
trueposA matrix with rows equal to the number of hypotheses and columns the number of samples of null test statistics (B) indicating the number of guessed true positives when using the corresponding value of the observed test statistic as a cut-off. Not returned unless keep.truepos=TRUE.
errormatThe matrix obtained after applying to type I error rate function closure to the matrices in falsepos, and, if applicable, truepos. Not returned unless keep.errormat=TRUE.
EB.h0MThe sum of the local q-values obtained after density estimation. This number serves as an estimate of the proportion of true null hypotheses. Values close to one indicate situations in which type I error control may not be guaranteed by the EBMTP. When prior='EBLQV', this value is used as the prior 'pi' during evaluation of the local q-value function.
priorThe numeric value of the prior 'pi' used when evaluating the local q-value function.
prior.typeCharacter string returning the value of prior in the original call to EBMTP. One of 'conservative', 'ABH', or 'EBLQV'.
lqvA numeric vector of length the number of hypotheses with the estimated local q-values used for generating guessed sets of true null hypotheses.
HsetsA numeric matrix with the same dimension as nulldist, containing the Bernoulli realizations of the estimated local q-values stored in lqv which were used to partition the hypotheses into guessed sets of true and false null hypotheses at each round of (re)sampling. Not returned unless keep.Hsets=TRUE.
labelIf keep.label=TRUE, a vector storing the values used in the argument Y. Storing this object is particularly important when one wishes to update EBMTP objects with F-statistics using default marg.null and marg.par settings when nulldist='boot.qt'.
indexFor tests of correlation parameters a matrix corresponding to t(combn(p,2)), where p is the number of variables in X. This matrix gives the indices of the variables considered in each pairwise correlation. For all other tests, this slot is empty, as the indices are in the same order as the rows of X.
callObject of class call, the call to the MTP function.
seedAn integer or vector for specifying the state of the random number generator used to create the resampled datasets. The seed can be reused for reproducibility in a repeat call to MTP. This argument is currently used only for the bootstrap null distribution (i.e., for nulldist="boot.xx"). See ?set.seed for details.
signature(x = "EBMTP")
: Subsetting method for EBMTP class, which operates selectively on each slot of an EBMTP instance to retain only the data related to the specified hypotheses.
: Converts an object of class EBMTP to an object of class list, with an entry for each slot.
: plot methods for EBMTP class, produces the following graphical summaries of the results of a EBMTP. The type of display may be specified via the which argument.
1. Scatterplot of number of rejected hypotheses vs. nominal Type I error rate.
2. Plot of ordered adjusted p-values; can be viewed as a plot of Type I error rate vs. number of rejected hypotheses.
3. Scatterplot of adjusted p-values vs. test statistics (also known as "volcano plot").
4. Plot of unordered adjusted p-values.
The plot method for objects of class EBMTP does not return the plots associated with which=5 (using confidence regions) or with which=6 (pertaining to cut-offs) as it does for objects of class MTP. This is because the function EBMTP currently only returns adjusted p-values. The argument logscale (by default equal to FALSE) allows one to use the negative decimal logarithms of the adjusted p-values in the second, third, and fourth graphical displays. The arguments caption and sub.caption allow one to change the titles and subtitles for each of the plots (default subtitle is the MTP function call). Note that some of these plots are implemented in the older function mt.plot.
: print method for EBMTP class, returns a description of an object of class EBMTP, including sample size, number of tested hypotheses, type of test performed (value of argument test), Type I error rate (value of argument typeone), nominal level of the test (value of argument alpha), name of the EBMTP (value of argument method), call to the function EBMTP.
In addition, this method produces a table with the class, mode, length, and dimension of each slot of the EBMTP instance.
: summary method for EBMTP class, provides numerical summaries of the results of an EBMTP and returns a list with the following three components.
1. rejections: A data.frame with the number(s) of rejected hypotheses for the nominal Type I error rate(s) specified by the alpha argument of the function MTP.
2. index: A numeric vector of indices for ordering the hypotheses according to first adjp, then rawp, and finally the absolute value of statistic (not printed in the summary).
3. summaries: When applicable (i.e., when the corresponding quantities are returned by MTP), a table with six number summaries of the distributions of the adjusted p-values, unadjusted p-values, test statistics, and parameter estimates.
: update method for EBMTP class, provides a mechanism to re-run the MTP with different choices of the following arguments - nulldist, alternative, typeone, k, q, alpha, smooth.null, bw, kernel, prior, keep.nulldist, keep.rawdist, keep.falsepos, keep.truepos, keep.errormat, keep.margpar. When evaluate is 'TRUE', a new object of class EBMTP is returned. Else, the updated call is returned. The EBMTP object passed to the update method must have either a non-empty rawdist slot or a non-empty nulldist slot (i.e., must have been called with either 'keep.rawdist=TRUE' or 'keep.nulldist=TRUE').
Additionally, when calling EBupdate for any Type I error rate other than FWER, the typeone argument must be specified (even if the original object did not control FWER). For example,
typeone="fdr", would always have to be specified, even if the original object also controlled the FDR. In other words, for all function arguments, it is safest to always assume that you
are updating from the EBMTP default function settings, regardless of the original call to the EBMTP function. Currently, the main advantage of the EBupdate method is that it prevents the need for repeated estimation of the test statistics null distribution.
To save on memory, if one knows ahead of time that one will want to compare different choices of bootstrap-based null distribution, then it is both necessary and sufficient to specify 'keep.rawdist=TRUE', as there is no other means of moving between null distributions other than through the non-transformed non-parametric bootstrap distribution. In this case, 'keep.nulldist=FALSE' may be used. Specifically, if an object of class EBMTP contains a non-empty rawdist slot and an empty nulldist slot, then a new null distribution will be generated according to the values of the nulldist= argument in the original call to EBMTP or any additional specifications in the call to update. On the other hand, if one knows that one wishes to only update an EBMTP object in ways which do not involve choice of null distribution, then 'keep.nulldist=TRUE' will suffice and 'keep.rawdist' can be set to FALSE (default settings). The original null distribution object will then be used for all subsequent calls to update.
N.B.: Note that keep.rawdist=TRUE is only available for the bootstrap-based resampling methods. The non-null distribution does not exist for the permutation or influence curve multivariate normal null distributions.
: coersion method for converting objects of class EBMTP to objects of class MTP. Slots common to both objects are taken from the object of class EBMTP and used to create a new object of class MTP. Once an object of class MTP is created, one may use the method update to perform resampling-based multiple testing (as would have been done with calls to MTP) without the need for repeated resampling.
Houston N. Gilbert, based on the original MTP class and method definitions written by Katherine S. Pollard
H.N. Gilbert, K.S. Pollard, M.J. van der Laan, and S. Dudoit (2009). Resampling-based multiple
hypothesis testing with applications to genomics: New developments in R/Bioconductor
package multtest. Journal of Statistical Software (submitted). Temporary URL: http://www.stat.berkeley.edu/~houston/JSSNullDistEBMTP.pdf.
Y. Benjamini and Y. Hochberg (2000). On the adaptive control of the false
discovery rate in multiple testing with independent statistics. J. Behav. Educ. Statist. Vol 25: 60-83.
Y. Benjamini, A. M. Krieger and D. Yekutieli (2006). Adaptive linear step-up
procedures that control the false discovery rate. Biometrika.
Vol. 93: 491-507.
M.J. van der Laan, M.D. Birkner, and A.E. Hubbard (2005). Empirical Bayes and Resampling Based Multiple Testing Procedure Controlling the Tail Probability of the Proportion of False Positives. Statistical Applications in Genetics and Molecular Biology, 4(1).
http://www.bepress.com/sagmb/vol4/iss1/art29/
S. Dudoit and M.J. van der Laan. Multiple Testing Procedures and Applications to Genomics. Springer Series in Statistics. Springer, New York, 2008.
S. Dudoit, H. N. Gilbert, and M. J. van der Laan (2008).
Resampling-based empirical Bayes multiple testing procedures for controlling
generalized tail probability and expected value error rates: Focus on the false discovery rate and simulation study. Biometrical Journal, 50(5):716-44. http://www.stat.berkeley.edu/~houston/BJMCPSupp/BJMCPSupp.html.
H.N. Gilbert, M.J. van der Laan, and S. Dudoit. Joint multiple testing procedures for
graphical model selection with applications to biological networks. Technical report,
U.C. Berkeley Division of Biostatistics Working Paper Series, April 2009. URL http://www.bepress.com/ucbbiostat/paper245.
EBMTP, EBMTP-methods,
MTP, MTP-methods,
[-methods, as.list-methods, print-methods, plot-methods, summary-methods, mtp2ebmtp,
ebmtp2mtp
## See EBMTP function: ? EBMTP## See EBMTP function: ? EBMTP
Gene expression data (3051 genes and 38 tumor mRNA samples) from the leukemia microarray study of Golub et al. (1999). Pre-processing was done as described in Dudoit et al. (2002). The R code for pre-processing is available in the file ../doc/golub.R.
data(golub)data(golub)
golub |
matrix of gene expression levels for the 38 tumor mRNA samples, rows correspond to genes (3051 genes) and columns to mRNA samples. |
golub.cl |
numeric vector indicating the tumor class, 27 acute lymphoblastic leukemia (ALL) cases (code 0) and 11 acute myeloid leukemia (AML) cases (code 1). |
golub.gnames |
a matrix containing the names of the 3051 genes for the expression matrix |
Golub et al. (1999). Molecular classification of cancer: class
discovery and class prediction by gene expression
monitoring, Science, Vol. 286:531-537.
http://www-genome.wi.mit.edu/MPR/
.
S. Dudoit, J. Fridlyand, and T. P. Speed (2002). Comparison of discrimination methods for the classification of tumors using gene expression data. Journal of the American Statistical Association, Vol. 97, No. 457, p. 77–87.
These functions are called internally by the main user-level function EBMTP. They are used for estimating local q-values, generating guessed sets of true null hypotheses, and applying these results to function closures defining the choice of type I error rate (FWER, gFWER, TPPFP, and FDR).
Hsets(Tn, nullmat, bw, kernel, prior, B, rawp) ABH.h0(rawp) G.VS(V, S = NULL, tp = TRUE, bound)Hsets(Tn, nullmat, bw, kernel, prior, B, rawp) ABH.h0(rawp) G.VS(V, S = NULL, tp = TRUE, bound)
Tn |
The vector of observed test statistics. |
nullmat |
The matrix of null test statistics obtained either through null transformation of the bootstrap distribution or by sampling from an appropriate multivariate normal distribution (when |
bw |
A character string argument to |
kernel |
A character string argument to |
prior |
Character string indicating which choice of prior probability to use for estimating local q-values (i.e., the posterior probabilities of a null hypothesis being true given the value of its corresponding test statistic). Default is 'conservative', in which case the prior is set to its most conservative value of 1, meaning that all hypotheses are assumed to belong to the set of true null hypotheses. Other options include 'ABH' for the adaptive Benjamini-Hochberg estimator of the number/proportion of true null hypotheses, and 'EBLQV' for the empirical Bayes local q-value value estimator of the number/proportion of true null hypotheses. If 'EBLQV', the estimator of the prior probability is taken to be the sum of the estimated local q-values divided by the number of tests. Relaxing the prior may result in more rejections, albeit at a cost of type I error control under certain conditions. See references. |
B |
The number of bootstrap iterations (i.e. how many resampled data sets) or the number of samples from the multivariate normal distribution (if |
rawp |
A vector of raw (unadjusted) p-values obtained bootstrap-based or influence curve null distribution. |
V |
A matrix of the numbers of guessed false positives for each cut-off, i.e., observed value of a test statistic, within each sample in |
S |
A matrix of the numbers of guessed true positives for each cut-off, i.e., observed value of a test statistic, within each sample in |
tp |
Logical indicator which is TRUE if type I error rate is a tail probability error rate and FALSE is if it is an expected value error rate. |
bound |
If a tail probability error rate, the bound to be placed on function of guessed false positives and guessed true positives. For, 'fwer', equal to 0; 'gfwer', equal to 'k'; and tppfp, equal to 'q'. |
The most important object to be returned from the function Hsets is a matrix of indicators, i.e., Bernoulli realizations of the estimated local q-values, taking the value of 1 if the hypothesis is guessed as belonging to the set of true null hypotheses and 0 otherwise (guessed true alternative). Realizations of these probabilities are generated with a call to rbinom, meaning that this function will set the RNG seed forward another B*(the number of hypotheses) places. This matrix, with rows equal to the number of hypotheses and columns the number of (bootstrap or multivariate normal) samples is used to subset the matrix of null test statistics and the vector of observed test statistics at each round of (re)sampling into samples of statistics guessed as belonging to the sets of true null and true alternative hypotheses, respectively. Using the values of the observed test statistics themselves as cut-offs, the numbers of guessed false positives and (if applicable) guessed true positives can be counted and eventually used to estimate the distribution of a type I error rate characterized by the closure returned from G.VS. Counting of guessed false positives and guessed true positives is performed in C through the function VScount.
For the function Hsets, a list with the following elements:
Hsets.mat |
A matrix of numeric indicators with rows equal to the number of test (hypotheses, typically |
EB.h0M |
The estimated proportion of true null hypotheses as determined by nonparametric density estimation. This value is the sum of the estimated local q-values divided by the total number of tests (hypotheses). |
prior |
The value of the prior applied to the local q-value function. If 'conservative', the prior is set to one. Otherwise, the prior is the value obtained from the estimator of the adaptive Benjamini-Hochberg procedure (if |
pn.out |
The vector of estimated local q-values. This vector is returned in the |
For the function ABH.h0, the estimated number of true null hypotheses using the estimator from the linear step-up adaptive Benjamini-Hochberg procedure.
For the function G.VS, a closure which accepts as arguments the matrices of guessed false positive and true positives (if applicable) and applies the appropriate function defining the desired type I error rate.
Houston N. Gilbert
H.N. Gilbert, K.S. Pollard, M.J. van der Laan, and S. Dudoit (2009). Resampling-based multiple
hypothesis testing with applications to genomics: New developments in R/Bioconductor
package multtest. Journal of Statistical Software (submitted). Temporary URL: http://www.stat.berkeley.edu/~houston/JSSNullDistEBMTP.pdf.
Y. Benjamini and Y. Hochberg (2000). On the adaptive control of the false
discovery rate in multiple testing with independent statistics. J. Behav.
Educ. Statist. Vol 25: 60-83.
Y. Benjamini, A.M. Krieger and D. Yekutieli (2006). Adaptive linear step-up
procedures that control the false discovery rate. Biometrika.
Vol. 93: 491-507.
M.J. van der Laan, M.D. Birkner, and A.E. Hubbard (2005). Empirical Bayes and Resampling Based Multiple Testing Procedure Controlling the Tail Probability of the Proportion of False Positives. Statistical Applications in Genetics and Molecular Biology, 4(1).
http://www.bepress.com/sagmb/vol4/iss1/art29/
S. Dudoit and M.J. van der Laan. Multiple Testing Procedures and Applications to Genomics. Springer Series in Statistics. Springer, New York, 2008.
S. Dudoit, H.N. Gilbert, and M.J. van der Laan (2008).
Resampling-based empirical Bayes multiple testing procedures for controlling
generalized tail probability and expected value error rates: Focus on the false
discovery rate and simulation study. Biometrical Journal, 50(5):716-44. http://www.stat.berkeley.edu/~houston/BJMCPSupp/BJMCPSupp.html.
H.N. Gilbert, M.J. van der Laan, and S. Dudoit. Joint multiple testing procedures for
graphical model selection with applications to biological networks. Technical report,
U.C. Berkeley Division of Biostatistics Working Paper Series, April 2009. URL http://www.bepress.com/ucbbiostat/paper245.
EBMTP, EBMTP-class, EBMTP-methods
set.seed(99) data<-matrix(rnorm(90),nr=9) group<-c(rep(1,5),rep(0,5)) #EB fwer control with centered and scaled bootstrap null distribution #(B=100 for speed) eb.m1<-EBMTP(X=data,Y=group,alternative="less",B=100,method="common.cutoff") print(eb.m1) summary(eb.m1) par(mfrow=c(2,2)) plot(eb.m1,top=9) abh <- ABH.h0(eb.m1@rawp) abh eb.m2 <- EBupdate(eb.m1,prior="ABH") eb.m2@priorset.seed(99) data<-matrix(rnorm(90),nr=9) group<-c(rep(1,5),rep(0,5)) #EB fwer control with centered and scaled bootstrap null distribution #(B=100 for speed) eb.m1<-EBMTP(X=data,Y=group,alternative="less",B=100,method="common.cutoff") print(eb.m1) summary(eb.m1) par(mfrow=c(2,2)) plot(eb.m1,top=9) abh <- ABH.h0(eb.m1@rawp) abh eb.m2 <- EBupdate(eb.m1,prior="ABH") eb.m2@prior
These functions compute permutation adjusted -values for step-down multiple testing procedures described in Westfall & Young (1993).
mt.maxT(X,classlabel,test="t",side="abs",fixed.seed.sampling="y",B=10000,na=.mt.naNUM,nonpara="n") mt.minP(X,classlabel,test="t",side="abs",fixed.seed.sampling="y",B=10000,na=.mt.naNUM,nonpara="n")mt.maxT(X,classlabel,test="t",side="abs",fixed.seed.sampling="y",B=10000,na=.mt.naNUM,nonpara="n") mt.minP(X,classlabel,test="t",side="abs",fixed.seed.sampling="y",B=10000,na=.mt.naNUM,nonpara="n")
X |
A data frame or matrix, with |
classlabel |
A vector of integers corresponding to observation (column)
class labels. For |
test |
A character string specifying the statistic to be
used to test the null hypothesis of no association between the
variables and the class labels. |
side |
A character string specifying the type of rejection region. |
fixed.seed.sampling |
If |
B |
The number of permutations. For a complete
enumeration, |
na |
Code for missing values (the default is |
nonpara |
If |
These functions compute permutation adjusted -values for the step-down maxT and minP multiple testing procedures, which provide strong control of the family-wise Type I error rate (FWER). The adjusted -values for the minP procedure are defined in equation (2.10) p. 66 of Westfall & Young (1993), and the maxT procedure is discussed p. 50 and 114. The permutation algorithms for estimating the adjusted -values are given in Ge et al. (In preparation). The procedures are for the simultaneous test of null hypotheses, namely, the null hypotheses of no association between the variables corresponding to the rows of the data frame X and the class labels classlabel. For gene expression data, the null hypotheses correspond to no differential gene expression across mRNA samples.
A data frame with components
index |
Vector of row indices, between 1 and |
teststat |
Vector of test statistics, ordered according to |
rawp |
Vector of raw (unadjusted) |
adjp |
Vector of adjusted |
plower |
For |
Yongchao Ge, [email protected],
Sandrine Dudoit, http://www.stat.berkeley.edu/~sandrine.
S. Dudoit, J. P. Shaffer, and J. C. Boldrick (Submitted). Multiple hypothesis testing in microarray experiments.
Y. Ge, S. Dudoit, and T. P. Speed. Resampling-based multiple testing for microarray data hypothesis, Technical Report \#633 of UCB Stat. http://www.stat.berkeley.edu/~gyc
P. H. Westfall and S. S. Young (1993). Resampling-based
multiple testing: Examples and methods for -value adjustment. John Wiley \& Sons.
mt.plot, mt.rawp2adjp, mt.reject, mt.sample.teststat, mt.teststat, golub.
# Gene expression data from Golub et al. (1999) # To reduce computation time and for illustrative purposes, we condider only # the first 100 genes and use the default of B=10,000 permutations. # In general, one would need a much larger number of permutations # for microarray data. data(golub) smallgd<-golub[1:100,] classlabel<-golub.cl # Permutation unadjusted p-values and adjusted p-values # for maxT and minP procedures with Welch t-statistics resT<-mt.maxT(smallgd,classlabel) resP<-mt.minP(smallgd,classlabel) rawp<-resT$rawp[order(resT$index)] teststat<-resT$teststat[order(resT$index)] # Plot results and compare to Bonferroni procedure bonf<-mt.rawp2adjp(rawp, proc=c("Bonferroni")) allp<-cbind(rawp, bonf$adjp[order(bonf$index),2], resT$adjp[order(resT$index)],resP$adjp[order(resP$index)]) mt.plot(allp, teststat, plottype="rvsa", proc=c("rawp","Bonferroni","maxT","minP"),leg=c(0.7,50),lty=1,col=1:4,lwd=2) mt.plot(allp, teststat, plottype="pvsr", proc=c("rawp","Bonferroni","maxT","minP"),leg=c(60,0.2),lty=1,col=1:4,lwd=2) mt.plot(allp, teststat, plottype="pvst", proc=c("rawp","Bonferroni","maxT","minP"),leg=c(-6,0.6),pch=16,col=1:4) # Permutation adjusted p-values for minP procedure with F-statistics (like equal variance t-statistics) mt.minP(smallgd,classlabel,test="f",fixed.seed.sampling="n") # Note that the test statistics used in the examples below are not appropriate # for the Golub et al. data. The sole purpose of these examples is to # demonstrate the use of the mt.maxT and mt.minP functions. # Permutation adjusted p-values for maxT procedure with paired t-statistics classlabel<-rep(c(0,1),19) mt.maxT(smallgd,classlabel,test="pairt") # Permutation adjusted p-values for maxT procedure with block F-statistics classlabel<-rep(0:18,2) mt.maxT(smallgd,classlabel,test="blockf",side="upper")# Gene expression data from Golub et al. (1999) # To reduce computation time and for illustrative purposes, we condider only # the first 100 genes and use the default of B=10,000 permutations. # In general, one would need a much larger number of permutations # for microarray data. data(golub) smallgd<-golub[1:100,] classlabel<-golub.cl # Permutation unadjusted p-values and adjusted p-values # for maxT and minP procedures with Welch t-statistics resT<-mt.maxT(smallgd,classlabel) resP<-mt.minP(smallgd,classlabel) rawp<-resT$rawp[order(resT$index)] teststat<-resT$teststat[order(resT$index)] # Plot results and compare to Bonferroni procedure bonf<-mt.rawp2adjp(rawp, proc=c("Bonferroni")) allp<-cbind(rawp, bonf$adjp[order(bonf$index),2], resT$adjp[order(resT$index)],resP$adjp[order(resP$index)]) mt.plot(allp, teststat, plottype="rvsa", proc=c("rawp","Bonferroni","maxT","minP"),leg=c(0.7,50),lty=1,col=1:4,lwd=2) mt.plot(allp, teststat, plottype="pvsr", proc=c("rawp","Bonferroni","maxT","minP"),leg=c(60,0.2),lty=1,col=1:4,lwd=2) mt.plot(allp, teststat, plottype="pvst", proc=c("rawp","Bonferroni","maxT","minP"),leg=c(-6,0.6),pch=16,col=1:4) # Permutation adjusted p-values for minP procedure with F-statistics (like equal variance t-statistics) mt.minP(smallgd,classlabel,test="f",fixed.seed.sampling="n") # Note that the test statistics used in the examples below are not appropriate # for the Golub et al. data. The sole purpose of these examples is to # demonstrate the use of the mt.maxT and mt.minP functions. # Permutation adjusted p-values for maxT procedure with paired t-statistics classlabel<-rep(c(0,1),19) mt.maxT(smallgd,classlabel,test="pairt") # Permutation adjusted p-values for maxT procedure with block F-statistics classlabel<-rep(0:18,2) mt.maxT(smallgd,classlabel,test="blockf",side="upper")
This function produces a number of graphical summaries
for the results of multiple testing procedures and their corresponding
adjusted -values.
mt.plot(adjp, teststat, plottype="rvsa", logscale=FALSE, alpha=seq(0, 1, length = 100), proc, leg=c(0, 0), ...)mt.plot(adjp, teststat, plottype="rvsa", logscale=FALSE, alpha=seq(0, 1, length = 100), proc, leg=c(0, 0), ...)
adjp |
A matrix of adjusted p-values, with rows
corresponding to hypotheses (genes) and columns to multiple testing
procedures. This matrix could be obtained from the functions
|
teststat |
A vector of test statistics for each of the hypotheses. This vector could be obtained from the functions |
plottype |
A character string specifying the type of graphical
summary for the results of the multiple testing procedures. |
logscale |
A logical variable for the |
alpha |
A vector of nominal Type I error rates for the |
proc |
A vector of character strings containing the names of the multiple testing procedures, to be used in the legend. |
... |
Graphical parameters such as |
leg |
A vector of coordinates for the legend. |
Sandrine Dudoit, http://www.stat.berkeley.edu/~sandrine,
Yongchao Ge, [email protected].
S. Dudoit, J. P. Shaffer, and J. C. Boldrick (Submitted). Multiple hypothesis testing in microarray experiments.
Y. Ge, S. Dudoit, and T. P. Speed. Resampling-based multiple testing for microarray data hypothesis, Technical Report \#633 of UCB Stat. http://www.stat.berkeley.edu/~gyc
mt.maxT, mt.minP, mt.rawp2adjp, mt.reject, mt.teststat, golub.
# Gene expression data from Golub et al. (1999) # To reduce computation time and for illustrative purposes, we condider only # the first 100 genes and use the default of B=10,000 permutations. # In general, one would need a much larger number of permutations # for microarray data. data(golub) smallgd<-golub[1:100,] classlabel<-golub.cl # Permutation unadjusted p-values and adjusted p-values for maxT procedure res1<-mt.maxT(smallgd,classlabel) rawp<-res1$rawp[order(res1$index)] teststat<-res1$teststat[order(res1$index)] # Permutation adjusted p-values for simple multiple testing procedures procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY") res2<-mt.rawp2adjp(rawp,procs) # Plot results from all multiple testing procedures allp<-cbind(res2$adjp[order(res2$index),],res1$adjp[order(res1$index)]) dimnames(allp)[[2]][9]<-"maxT" procs<-dimnames(allp)[[2]] procs[7:9]<-c("maxT","BH","BY") allp<-allp[,procs] cols<-c(1:4,"orange","brown","purple",5:6) ltypes<-c(3,rep(1,6),rep(2,2)) # Ordered adjusted p-values mt.plot(allp,teststat,plottype="pvsr",proc=procs,leg=c(80,0.4),lty=ltypes,col=cols,lwd=2) # Adjusted p-values in original data order mt.plot(allp,teststat,plottype="pvsi",proc=procs,leg=c(80,0.4),lty=ltypes,col=cols,lwd=2) # Number of rejected hypotheses vs. level of the test mt.plot(allp,teststat,plottype="rvsa",proc=procs,leg=c(0.05,100),lty=ltypes,col=cols,lwd=2) # Adjusted p-values vs. test statistics mt.plot(allp,teststat,plottype="pvst",logscale=TRUE,proc=procs,leg=c(0,4),pch=ltypes,col=cols)# Gene expression data from Golub et al. (1999) # To reduce computation time and for illustrative purposes, we condider only # the first 100 genes and use the default of B=10,000 permutations. # In general, one would need a much larger number of permutations # for microarray data. data(golub) smallgd<-golub[1:100,] classlabel<-golub.cl # Permutation unadjusted p-values and adjusted p-values for maxT procedure res1<-mt.maxT(smallgd,classlabel) rawp<-res1$rawp[order(res1$index)] teststat<-res1$teststat[order(res1$index)] # Permutation adjusted p-values for simple multiple testing procedures procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY") res2<-mt.rawp2adjp(rawp,procs) # Plot results from all multiple testing procedures allp<-cbind(res2$adjp[order(res2$index),],res1$adjp[order(res1$index)]) dimnames(allp)[[2]][9]<-"maxT" procs<-dimnames(allp)[[2]] procs[7:9]<-c("maxT","BH","BY") allp<-allp[,procs] cols<-c(1:4,"orange","brown","purple",5:6) ltypes<-c(3,rep(1,6),rep(2,2)) # Ordered adjusted p-values mt.plot(allp,teststat,plottype="pvsr",proc=procs,leg=c(80,0.4),lty=ltypes,col=cols,lwd=2) # Adjusted p-values in original data order mt.plot(allp,teststat,plottype="pvsi",proc=procs,leg=c(80,0.4),lty=ltypes,col=cols,lwd=2) # Number of rejected hypotheses vs. level of the test mt.plot(allp,teststat,plottype="rvsa",proc=procs,leg=c(0.05,100),lty=ltypes,col=cols,lwd=2) # Adjusted p-values vs. test statistics mt.plot(allp,teststat,plottype="pvst",logscale=TRUE,proc=procs,leg=c(0,4),pch=ltypes,col=cols)
This function computes adjusted -values for simple
multiple testing procedures from a vector of raw (unadjusted)
-values. The procedures include the Bonferroni, Holm (1979),
Hochberg (1988), and Sidak procedures for strong control of the
family-wise Type I error rate (FWER), and the Benjamini & Hochberg
(1995) and Benjamini & Yekutieli (2001) procedures for (strong)
control of the false discovery rate (FDR). The less conservative
adaptive Benjamini & Hochberg (2000) and two-stage Benjamini & Hochberg
(2006) FDR-controlling procedures are also included.
mt.rawp2adjp(rawp, proc=c("Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY","ABH","TSBH"), alpha = 0.05, na.rm = FALSE)mt.rawp2adjp(rawp, proc=c("Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY","ABH","TSBH"), alpha = 0.05, na.rm = FALSE)
rawp |
A vector of raw (unadjusted) |
proc |
A vector of character strings containing the names of the
multiple testing procedures for which adjusted Adjusted
|
alpha |
A nominal type I error rate, or a vector of error
rates, used for estimating the number of true null hypotheses in the
two-stage Benjamini & Hochberg procedure ( |
na.rm |
An option for handling |
A list with components:
adjp |
A matrix of adjusted |
index |
A vector of row indices, between 1 and
|
h0.ABH |
The estimate of the number of true null hypotheses as proposed
by Benjamini & Hochberg (2000) used when computing adjusted |
h0.TSBH |
The estimate (or vector of estimates) of the number of true
null hypotheses as proposed by Benjamini et al. (2006) when computing adjusted
|
Sandrine Dudoit, http://www.stat.berkeley.edu/~sandrine,
Yongchao Ge, [email protected],
Houston Gilbert, http://www.stat.berkeley.edu/~houston.
Y. Benjamini and Y. Hochberg (1995). Controlling the false discovery
rate: a practical and powerful approach to multiple
testing. J. R. Statist. Soc. B. Vol. 57: 289-300.
Y. Benjamini and Y. Hochberg (2000). On the adaptive control of the false discovery rate in multiple testing with independent statistics. J. Behav. Educ. Statist. Vol 25: 60-83.
Y. Benjamini and D. Yekutieli (2001). The control of the false discovery rate in multiple hypothesis testing under dependency. Annals of Statistics. Vol. 29: 1165-88.
Y. Benjamini, A. M. Krieger and D. Yekutieli (2006). Adaptive linear step-up procedures that control the false discovery rate. Biometrika. Vol. 93: 491-507.
S. Dudoit, J. P. Shaffer, and J. C. Boldrick (2003). Multiple
hypothesis testing in microarray experiments. Statistical Science. Vol. 18: 71-103.
S. Dudoit, H. N. Gilbert, and M. J. van der Laan (2008).
Resampling-based empirical Bayes multiple testing procedures for controlling generalized tail probability and expected value error rates: Focus on the false discovery rate and simulation study. Biometrical Journal, 50(5):716-44. http://www.stat.berkeley.edu/~houston/BJMCPSupp/BJMCPSupp.html.
Y. Ge, S. Dudoit, and T. P. Speed (2003). Resampling-based multiple testing for microarray data analysis. TEST. Vol. 12: 1-44 (plus discussion p. 44-77).
Y. Hochberg (1988). A sharper Bonferroni procedure for multiple tests of significance, Biometrika. Vol. 75: 800-802.
S. Holm (1979). A simple sequentially rejective multiple test procedure. Scand. J. Statist.. Vol. 6: 65-70.
mt.maxT, mt.minP,
mt.plot, mt.reject, golub.
# Gene expression data from Golub et al. (1999) # To reduce computation time and for illustrative purposes, we condider only # the first 100 genes and use the default of B=10,000 permutations. # In general, one would need a much larger number of permutations # for microarray data. data(golub) smallgd<-golub[1:100,] classlabel<-golub.cl # Permutation unadjusted p-values and adjusted p-values for maxT procedure res1<-mt.maxT(smallgd,classlabel) rawp<-res1$rawp[order(res1$index)] # Permutation adjusted p-values for simple multiple testing procedures procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY","ABH","TSBH") res2<-mt.rawp2adjp(rawp,procs)# Gene expression data from Golub et al. (1999) # To reduce computation time and for illustrative purposes, we condider only # the first 100 genes and use the default of B=10,000 permutations. # In general, one would need a much larger number of permutations # for microarray data. data(golub) smallgd<-golub[1:100,] classlabel<-golub.cl # Permutation unadjusted p-values and adjusted p-values for maxT procedure res1<-mt.maxT(smallgd,classlabel) rawp<-res1$rawp[order(res1$index)] # Permutation adjusted p-values for simple multiple testing procedures procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY","ABH","TSBH") res2<-mt.rawp2adjp(rawp,procs)
This function returns the identity and number of rejected hypotheses for several multiple testing procedures and different nominal Type I error rates.
mt.reject(adjp, alpha)mt.reject(adjp, alpha)
adjp |
A matrix of adjusted p-values, with rows
corresponding to hypotheses and columns to multiple testing
procedures. This matrix could be obtained from the function
|
alpha |
A vector of nominal Type I error rates. |
A list with components
r |
A matrix containing the number of rejected hypotheses for several multiple testing procedures and different nominal Type I error rates. Rows correspond to Type I error rates and columns to multiple testing procedures. |
which |
A matrix of indicators for the rejection of individual hypotheses by different multiple testing procedures for a nominal Type I error rate |
Sandrine Dudoit, http://www.stat.berkeley.edu/~sandrine,
Yongchao Ge, [email protected].
mt.maxT, mt.minP, mt.rawp2adjp, golub.
# Gene expression data from Golub et al. (1999) # To reduce computation time and for illustrative purposes, we condider only # the first 100 genes and use the default of B=10,000 permutations. # In general, one would need a much larger number of permutations # for microarray data. data(golub) smallgd<-golub[1:100,] classlabel<-golub.cl # Permutation unadjusted p-values and adjusted p-values for maxT procedure res<-mt.maxT(smallgd,classlabel) mt.reject(cbind(res$rawp,res$adjp),seq(0,1,0.1))$r# Gene expression data from Golub et al. (1999) # To reduce computation time and for illustrative purposes, we condider only # the first 100 genes and use the default of B=10,000 permutations. # In general, one would need a much larger number of permutations # for microarray data. data(golub) smallgd<-golub[1:100,] classlabel<-golub.cl # Permutation unadjusted p-values and adjusted p-values for maxT procedure res<-mt.maxT(smallgd,classlabel) mt.reject(cbind(res$rawp,res$adjp),seq(0,1,0.1))$r
These functions provide tools to investigate the permutation distribution
of test statistics, raw (unadjusted) -values, and class labels.
mt.sample.teststat(V,classlabel,test="t",fixed.seed.sampling="y",B=10000,na=.mt.naNUM,nonpara="n") mt.sample.rawp(V,classlabel,test="t",side="abs",fixed.seed.sampling="y",B=10000,na=.mt.naNUM,nonpara="n") mt.sample.label(classlabel,test="t",fixed.seed.sampling="y",B=10000)mt.sample.teststat(V,classlabel,test="t",fixed.seed.sampling="y",B=10000,na=.mt.naNUM,nonpara="n") mt.sample.rawp(V,classlabel,test="t",side="abs",fixed.seed.sampling="y",B=10000,na=.mt.naNUM,nonpara="n") mt.sample.label(classlabel,test="t",fixed.seed.sampling="y",B=10000)
V |
A numeric vector containing the data for one of the variables (genes). |
classlabel |
A vector of integers corresponding to observation (column)
class labels. For |
test |
A character string specifying the statistic to be
used to test the null hypothesis of no association between the
variables and the class labels. |
side |
A character string specifying the type of rejection region. |
fixed.seed.sampling |
If |
B |
The number of permutations. For a complete
enumeration, |
na |
Code for missing values (the default is |
nonpara |
If |
For mt.sample.teststat, a vector containing B permutation test statistics.
For mt.sample.rawp, a vector containing B permutation unadjusted -values.
For mt.sample.label, a matrix containing B
sets of permuted class labels. Each row corresponds to one permutation.
Yongchao Ge, [email protected],
Sandrine Dudoit, http://www.stat.berkeley.edu/~sandrine.
# Gene expression data from Golub et al. (1999) data(golub) mt.sample.label(golub.cl,B=10) permt<-mt.sample.teststat(golub[1,],golub.cl,B=1000) qqnorm(permt) qqline(permt) permt<-mt.sample.teststat(golub[50,],golub.cl,B=1000) qqnorm(permt) qqline(permt) permp<-mt.sample.rawp(golub[1,],golub.cl,B=1000) hist(permp)# Gene expression data from Golub et al. (1999) data(golub) mt.sample.label(golub.cl,B=10) permt<-mt.sample.teststat(golub[1,],golub.cl,B=1000) qqnorm(permt) qqline(permt) permt<-mt.sample.teststat(golub[50,],golub.cl,B=1000) qqnorm(permt) qqline(permt) permp<-mt.sample.rawp(golub[1,],golub.cl,B=1000) hist(permp)
These functions provide a convenient way to compute test statistics, e.g., two-sample Welch t-statistics, Wilcoxon statistics, F-statistics, paired t-statistics, block F-statistics, for each row of a data frame.
mt.teststat(X,classlabel,test="t",na=.mt.naNUM,nonpara="n") mt.teststat.num.denum(X,classlabel,test="t",na=.mt.naNUM,nonpara="n")mt.teststat(X,classlabel,test="t",na=.mt.naNUM,nonpara="n") mt.teststat.num.denum(X,classlabel,test="t",na=.mt.naNUM,nonpara="n")
X |
A data frame or matrix, with |
classlabel |
A vector of integers corresponding to observation (column)
class labels. For |
test |
A character string specifying the statistic to be
used to test the null hypothesis of no association between the
variables and the class labels. |
na |
Code for missing values (the default is |
nonpara |
If |
For mt.teststat, a vector of test statistics for each row (gene).
For mt.teststat.num.denum, a data frame with
teststat.num |
the numerator of the test statistics for each row, depending on the
specific |
teststat.denum |
the denominator of the test statistics for each row, depending on the
specific |
Yongchao Ge, [email protected],
Sandrine Dudoit, http://www.stat.berkeley.edu/~sandrine.
# Gene expression data from Golub et al. (1999) data(golub) teststat<-mt.teststat(golub,golub.cl) qqnorm(teststat) qqline(teststat) tmp<-mt.teststat.num.denum(golub,golub.cl,test="t") num<-tmp$teststat.num denum<-tmp$teststat.denum plot(sqrt(denum),num) tmp<-mt.teststat.num.denum(golub,golub.cl,test="f")# Gene expression data from Golub et al. (1999) data(golub) teststat<-mt.teststat(golub,golub.cl) qqnorm(teststat) qqline(teststat) tmp<-mt.teststat.num.denum(golub,golub.cl,test="t") num<-tmp$teststat.num denum<-tmp$teststat.denum plot(sqrt(denum),num) tmp<-mt.teststat.num.denum(golub,golub.cl,test="f")
A user-level function to perform multiple testing procedures (MTP). A variety of t- and F-tests, including robust versions of most tests, are implemented. Single-step and step-down minP and maxT methods are used to control the chosen type I error rate (FWER, gFWER, TPPFP, or FDR). Bootstrap and permutation null distributions are available. Additionally, for t-statistics, one may wish to sample from an appropriate multivariate normal distribution with mean zero and correlation matrix derived from the vector influence function. Arguments are provided for user control of output. Gene selection in microarray experiments is one application.
MTP(X, W = NULL, Y = NULL, Z = NULL, Z.incl = NULL, Z.test = NULL, na.rm = TRUE, test = "t.twosamp.unequalvar", robust = FALSE, standardize = TRUE, alternative = "two.sided", psi0 = 0, typeone = "fwer", k = 0, q = 0.1, fdr.method = "conservative", alpha = 0.05, smooth.null = FALSE, nulldist = "boot.cs", B = 1000, ic.quant.trans = FALSE, MVN.method = "mvrnorm", penalty = 1e-06, method = "ss.maxT", get.cr = FALSE, get.cutoff = FALSE, get.adjp = TRUE, keep.nulldist = TRUE, keep.rawdist = FALSE, seed = NULL, cluster = 1, type = NULL, dispatch = NULL, marg.null = NULL, marg.par = NULL, keep.margpar = TRUE, ncp = NULL, perm.mat = NULL, keep.index = FALSE, keep.label = FALSE)MTP(X, W = NULL, Y = NULL, Z = NULL, Z.incl = NULL, Z.test = NULL, na.rm = TRUE, test = "t.twosamp.unequalvar", robust = FALSE, standardize = TRUE, alternative = "two.sided", psi0 = 0, typeone = "fwer", k = 0, q = 0.1, fdr.method = "conservative", alpha = 0.05, smooth.null = FALSE, nulldist = "boot.cs", B = 1000, ic.quant.trans = FALSE, MVN.method = "mvrnorm", penalty = 1e-06, method = "ss.maxT", get.cr = FALSE, get.cutoff = FALSE, get.adjp = TRUE, keep.nulldist = TRUE, keep.rawdist = FALSE, seed = NULL, cluster = 1, type = NULL, dispatch = NULL, marg.null = NULL, marg.par = NULL, keep.margpar = TRUE, ncp = NULL, perm.mat = NULL, keep.index = FALSE, keep.label = FALSE)
X |
A matrix, data.frame or ExpressionSet containing the raw data. In the case of an ExpressionSet, |
W |
A vector or matrix containing non-negative weights to be used in computing the test statistics. If a matrix, |
Y |
A vector, factor, or |
Z |
A vector, factor, or matrix containing covariate data to be used in the regression (linear and Cox) models. Each variable should be in one column, so that |
Z.incl |
The indices of the columns of |
Z.test |
The index or names of the column of |
na.rm |
Logical indicating whether to remove observations with an NA. Default is 'TRUE'. |
test |
Character string specifying the test statistics to use, by default 't.twosamp.unequalvar'. See details (below) for a list of tests. |
robust |
Logical indicating whether to use the robust version of the chosen test, e.g. Wilcoxon singed rank test for robust one-sample t-test or |
standardize |
Logical indicating whether to use the standardized version of the test statistics (usual t-statistics are standardized). Default is 'TRUE'. |
alternative |
Character string indicating the alternative hypotheses, by default 'two.sided'. For one-sided tests, use 'less' or 'greater' for null hypotheses of 'greater than or equal' (i.e. alternative is 'less') and 'less than or equal', respectively. |
psi0 |
The hypothesized null value, typically zero (default). Currently, this should be a single value, which is used for all hypotheses. |
typeone |
Character string indicating which type I error rate to control, by default family-wise error rate ('fwer'). Other options include generalized family-wise error rate ('gfwer'), with parameter |
k |
The allowed number of false positives for gFWER control. Default is 0 (FWER). |
q |
The allowed proportion of false positives for TPPFP control. Default is 0.1. |
fdr.method |
Character string indicating which FDR controlling method should be used when |
alpha |
The target nominal type I error rate, which may be a vector of error rates. Default is 0.05. |
smooth.null |
Indicator of whether to use a kernel density estimate for the tail of the null distributon for computing raw pvalues close to zero. Only used if 'rawp' would be zero without smoothing. Default is 'FALSE'. |
nulldist |
Character string indicating which resampling method to use for estimating the joint test statistics null distribution, by default the non-parametric bootstrap with centering and scaling ('boot.cs'). The old default 'boot' will still compile and will correspond to 'boot.cs'. Other null distribution options include 'perm', 'boot.ctr', 'boot.qt', and 'ic', corresponding to the permutation distribution, centered-only bootstrap distribution, quantile-transformed bootstrap distribution, and influence curve multivariate normal joint null distribution, respectively. More details below. |
B |
The number of bootstrap iterations (i.e. how many resampled data sets), the number of permutations (if |
ic.quant.trans |
If |
MVN.method |
If |
penalty |
If |
method |
The multiple testing procedure to use. Options are single-step maxT ('ss.maxT', default), single-step minP ('ss.minP'), step-down maxT ('sd.maxT'), and step-down minP ('sd.minP'). |
get.cr |
Logical indicating whether to compute confidence intervals for the estimates. Not available for F-tests. Default is 'FALSE'. |
get.cutoff |
Logical indicating whether to compute thresholds for the test statistics. Default is 'FALSE'. |
get.adjp |
Logical indicating whether to compute adjusted p-values. Default is 'TRUE'. |
keep.nulldist |
Logical indicating whether to return the computed bootstrap or influence curve null distribution, by default 'TRUE'. Not available for |
keep.rawdist |
Logical indicating whether to return the computed non-null (raw) bootstrap distribution, by default 'FALSE'. Not available when using |
seed |
Integer or vector of integers to be used as argument to |
cluster |
Integer for number of nodes to create or a cluster object created through the package snow. With |
type |
Interface system to use for computer cluster. See |
dispatch |
The number or percentage of bootstrap iterations to dispatch at a time to each node of the cluster if a computer cluster is used. If dispatch is a percentage, |
marg.null |
If |
marg.par |
If |
keep.margpar |
If |
ncp |
If |
perm.mat |
If |
keep.index |
If |
keep.label |
Default is 'FALSE'. A logical indicating whether or not the label in |
A multiple testing procedure (MTP) is defined by choices of test statistics, type I error rate, null distribution and method for error rate control. Each component is described here. For two-sample t-tests, the group with the smaller-valued label is substracted from the group with the larger-valued label. That is, differences in means are calculated as "mean of group 2 - mean of group 1" or "mean of group B - mean of group A". For paired t-tests, the arrangement of group indices does not matter, as long as the columns are arranged in the same corresponding order between groups. For example, if group 1 is coded as 0, and group 2 is coded as 1, for 3 pairs of data, it does not matter if the label Y is coded as "0,0,0,1,1,1", "1,1,1,0,0,0" "0,1,0,1,0,1" or "1,0,1,0,1,0", the paired differences between groups will be calculated as "group 2 - group 1". See references for more detail.
Test statistics are determined by the values of test:
one-sample t-statistic for tests of means;
equal variance two-sample t-statistic for tests of differences in means (two-sample t-statistic);
unequal variance two-sample t-statistic for tests of differences in means (two-sample Welch t-statistic);
two-sample paired t-statistic for tests of differences in means;
multi-sample F-statistic for tests of equality of population means (assumes constant variance across groups, but not normality);
multi-sample F-statistic for tests of equality of population means in a block design (assumes constant variance across groups, but not normality). This test is not available with the bootstrap null distribution;
multi-sample F-statistic for tests of equality of population means in a block design (assumes constant variance across groups, but not normality). Differs from f.block in requiring multiple observations per group*block combintation. This test uses the means of each group*block combination as response variable and test for group main effects assuming a randomized block design;
t-statistic for tests of regression coefficients for variable Z.test in linear models, each with a row of X as outcome, possibly adjusted by covariates Z.incl from the matrix Z (in the case of no covariates, one recovers the one-sample t-statistic, t.onesamp);
t-statistic for tests of regression coefficients in linear models, with outcome Y and each row of X as covariate of interest, with possibly other covariates Z.incl from the matrix Z;
t-statistic for tests of regression coefficients in Cox proportional hazards survival models, with outcome Y and each row of X as covariate of interest, with possibly other covariates Z.incl from the matrix Z.
t-statistics for tests of pairwise correlation parameters for all variables in X. Note that the number of hypotheses can become quite large very fast. This test is only available with the influence curve null distribution.
Fisher's z-statistics for tests of pairwise correlation parameters for all variables in X. Note that the number of hypotheses can become quite large very fast. This test is only available with the influence curve null distribution.
When robust=TRUE, non-parametric versions of each test are performed. For the linear models, this means rlm is used instead of lm. There is not currently a robust version of test=coxph.YvsXZ. For the t- and F-tests, data values are simply replaced by their ranks. This is equivalent to performing the following familiar named rank-based tests. The conversion after each test is the formula to convert from the MTP test to the statistic reported by the listed R function (where num is the numerator of the MTP test statistics, n is total sample size, nk is group k sample size, K is total number of groups or treatments, and rk are the ranks in group k).
Wilcoxon signed rank, wilcox.test with y=NULL or paired=TRUE,
conversion: num/n
Wilcoxon rank sum or Mann-Whitney, wilcox.test,
conversion: n2*(num+mean(r1)) - n2*(n2+1)/2
Kruskal-Wallis rank sum, kruskal.test,
conversion: num*12/(n*(n-1))
Friedman rank sum, friedman.test,
conversion: num*12/(K*(K+1))
Friedman rank sum, friedman.test,
conversion: num*12/(K*(K+1))
The implemented MTPs are based on control of the family-wise error rate, defined as the probability of any false positives. Let Vn denote the (unobserved) number of false positives. Then, control of FWER at level alpha means that Pr(Vn>0)<=alpha. The set of rejected hypotheses under a FWER controlling procedure can be augmented to increase the number of rejections, while controlling other error rates. The generalized family-wise error rate is defined as Pr(Vn>k)<=alpha, and it is clear that one can simply take an FWER controlling procedure, reject k more hypotheses and have control of gFWER at level alpha. The tail probability of the proportion of false positives depends on both the number of false postives (Vn) and the number of rejections (Rn). Control of TPPFP at level alpha means Pr(Vn/Rn>q)<=alpha, for some proportion q. Control of the false discovery rate refers to the expected proportion of false positives (rather than a tail probability). Control of FDR at level alpha means E(Vn/Rn)<=alpha.
In practice, one must choose a method for estimating the test statistics null distribution. We have implemented several versions of an ordinary non-parametric bootstrap estimator and a permutation estimator (which makes sense in certain settings, see references). The non-parametric bootstrap estimator (default) provides asymptotic control of the type I error rate for any data generating distribution, whereas the permutation estimator requires the subset pivotality assumption. One draw back of both methods is the discreteness of the estimated null distribution when the sample size is small. Furthermore, when the sample size is small enough, it is possible that ties will lead to a very small variance estimate. Using standardize=FALSE allows one to avoid these unusually small test statistic denominators. Parametric bootstrap estimators are another option (not yet implemented). For asymptotically linear estimators, such as those commonly probed using t-statistics, another choice of null distribution is provided when sampling from a multivariate normal distribution with mean zero and correlation matrix derived from the vector influence function. Sampling from a multivariate normal may alleviate the discreteness of the bootstrap and permutation distributions, although accuracy in estimation of the test statistics correlation matrix will be of course also affected by sample size.
For the nonparametric bootstrap distribution with marginal null quantile transformation, the following defaults for marg.null and marg.par are available based on choice of test statistics, sample size 'n', and various other parameters:
t-distribution with df=n-1;
t-distribution with df=n-2;
N(0,1);
t-distribution with df=n-1, where n is the number of unique samples, i.e., the number of observed differences between paired samples;
F-distribution with df1=k-1, df2=n-k, for k groups;
NA. Only available with permutation distribution;
F-distribution with df1=k-1,df2=n-k*l, for k groups and l blocks;
N(0,1);
N(0,1);
N(0,1);
t-distribution with df=n-2;
N(0,1).
The above defaults, however, can be overridden by manually setting values of marg.null and marg.par. In the case of nulldist='ic', and ic.quant.trans=TRUE, the defaults are the same as above except that 'lm.XvsZ' and 'lm.YvsXZ' are replaced with t-distributions with df=n-p.
Given observed test statistics, a type I error rate (with nominal level), and a test statistics null distribution, MTPs provide adjusted p-values, cutoffs for test statistics, and possibly confidence regions for estimates. Four methods are implemented, based on minima of p-values and maxima of test statistics. Only the step down methods are currently available with the permutation null distribution.
Computation times using a bootstrap null distribution are slower when weights are used for one and two-sample tests. Computation times when using a bootstrap null distribution also are slower for the tests lmXvsZ, lmYvsXZ, coxph.YvsXZ.
To execute the bootstrap on a computer cluster, a cluster object generated with makeCluster in the package snow may be used as the argument for cluster. Alternatively, the number of nodes to use in the computer cluster can be used as the argument to cluster. In this case, type must be specified and a cluster will be created. In both cases, Biobase and multtest will be loaded onto each cluster node if these libraries are located in a directory in the standard search path. If these libraries are in a non-standard location, it is necessary to first create the cluster, load Biobase and multtest on each node and then to use the cluster object as the argument to cluster. See documentation for snow package for additional information on creating and using a cluster.
Finally, note that the old argument csnull is now DEPRECATED as of multtest v. 2.0.0 given the expanded null distribution options described above. Previously, this argument was an indicator of whether the bootstrap estimated test statistics distribution should be centered and scaled (to produce a null distribution) or not. If csnull=FALSE, the (raw) non-null bootstrap estimated test statistics distribution was returned. If the non-null bootstrap distribution should be returned, this object is now stored in the 'rawdist' slot when keep.rawdist=TRUE in the original MTP function call.
An object of class MTP, with the following slots:
statistic |
Object of class |
estimate |
For the test of single-parameter null hypotheses using t-statistics (i.e., not the F-tests), the numeric vector of estimated parameters corresponding to each hypothesis, e.g. means, differences in means, regression parameters. |
sampsize |
Object of class |
rawp |
Object of class |
adjp |
Object of class |
conf.reg |
For the test of single-parameter null hypotheses using t-statistics (i.e., not the F-tests), the numeric array of lower and upper simultaneous confidence limits for the parameter vector, for each value of the nominal Type I error rate |
cutoff |
The numeric matrix of cut-offs for the vector of test statistics for each value of the nominal Type I error rate |
reject |
Object of class |
rawdist |
The numeric matrix for the estimated nonparametric non-null test statistics distribution (returned only if |
nulldist |
The numeric matrix for the estimated test statistics null distribution (returned only if |
nulldist.type |
Character value describing which choice of null distribution was used to generate the MTP results. Takes on one of the values of the original |
marg.null |
If |
marg.par |
If |
call |
Object of class |
seed |
An integer or vector for specifying the state of the random number generator used to create the resampled datasets. The seed can be reused for reproducibility in a repeat call to |
Thank you to Peter Dimitrov for suggestions about the code.
Katherine S. Pollard and Houston N. Gilbert with design contributions from Sandra Taylor, Sandrine Dudoit and Mark J. van der Laan.
M.J. van der Laan, S. Dudoit, K.S. Pollard (2004), Augmentation Procedures for Control of the Generalized Family-Wise Error Rate and Tail Probabilities for the Proportion of False Positives, Statistical Applications in Genetics and Molecular Biology, 3(1). http://www.bepress.com/sagmb/vol3/iss1/art15/
M.J. van der Laan, S. Dudoit, K.S. Pollard (2004), Multiple Testing. Part II. Step-Down Procedures for Control of the Family-Wise Error Rate, Statistical Applications in Genetics and Molecular Biology, 3(1). http://www.bepress.com/sagmb/vol3/iss1/art14/
S. Dudoit, M.J. van der Laan, K.S. Pollard (2004), Multiple Testing. Part I. Single-Step Procedures for Control of General Type I Error Rates, Statistical Applications in Genetics and Molecular Biology, 3(1). http://www.bepress.com/sagmb/vol3/iss1/art13/
K.S. Pollard and Mark J. van der Laan, "Resampling-based Multiple Testing: Asymptotic Control of Type I Error and Applications to Gene Expression Data" (June 24, 2003). U.C. Berkeley Division of Biostatistics Working Paper Series. Working Paper 121. http://www.bepress.com/ucbbiostat/paper121
M.J. van der Laan and A.E. Hubbard (2006), Quantile-function Based Null Distributions in Resampling Based Multiple Testing, Statistical Applications in Genetics and Molecular Biology, 5(1). http://www.bepress.com/sagmb/vol5/iss1/art14/
S. Dudoit and M.J. van der Laan. Multiple Testing Procedures and Applications to Genomics. Springer Series in Statistics. Springer, New York, 2008.
EBMTP, MTP-class, MTP-methods, mt.minP, mt.maxT, ss.maxT, fwer2gfwer
#data set.seed(99) data<-matrix(rnorm(90),nr=9) group<-c(rep(1,5),rep(0,5)) #fwer control with centered and scaled bootstrap null distribution #(B=100 for speed) m1<-MTP(X=data,Y=group,alternative="less",B=100,method="sd.minP") print(m1) summary(m1) par(mfrow=c(2,2)) plot(m1,top=9) #fwer control with quantile transformed bootstrap null distribution #default settings = N(0,1) marginal null distribution m2<-MTP(X=data,Y=group,alternative="less",B=100,method="sd.minP", nulldist="boot.qt",keep.rawdist=TRUE) #fwer control with quantile transformed bootstrap null distribution #marginal null distribution and df parameters manually set, #first all equal, then varying with hypothesis m3<-update(m2,marg.null="t",marg.par=10) mps<-matrix(c(rep(9,5),rep(10,5)),nr=10,nc=1) m4<-update(m2,marg.null="t",marg.par=mps) [email protected] [email protected] [email protected] [email protected] [email protected] [email protected] [email protected] [email protected] [email protected] [email protected]#data set.seed(99) data<-matrix(rnorm(90),nr=9) group<-c(rep(1,5),rep(0,5)) #fwer control with centered and scaled bootstrap null distribution #(B=100 for speed) m1<-MTP(X=data,Y=group,alternative="less",B=100,method="sd.minP") print(m1) summary(m1) par(mfrow=c(2,2)) plot(m1,top=9) #fwer control with quantile transformed bootstrap null distribution #default settings = N(0,1) marginal null distribution m2<-MTP(X=data,Y=group,alternative="less",B=100,method="sd.minP", nulldist="boot.qt",keep.rawdist=TRUE) #fwer control with quantile transformed bootstrap null distribution #marginal null distribution and df parameters manually set, #first all equal, then varying with hypothesis m3<-update(m2,marg.null="t",marg.par=10) mps<-matrix(c(rep(9,5),rep(10,5)),nr=10,nc=1) m4<-update(m2,marg.null="t",marg.par=mps) m1@nulldist.type m2@nulldist.type m2@marg.null m2@marg.par m3@nulldist.type m3@marg.null m3@marg.par m4@nulldist.type m4@marg.null m4@marg.par
An object of class MTP is the output of a particular multiple testing procedure, for example, generated by the MTP function. It has slots for the various data used to make multiple testing decisions, such as adjusted p-values and confidence regions.
Objects can be created by calls of the form
new('MTP',
statistic = ...., object of class numeric
estimate = ...., object of class numeric
sampsize = ...., object of class numeric
rawp = ...., object of class numeric
adjp = ...., object of class numeric
conf.reg = ...., object of class array
cutoff = ...., object of class matrix
reject = ...., object of class matrix
rawdist = ...., object of class matrix
nulldist = ...., object of class matrix
nulldist.type = ...., object of class character
marg.null = ...., object of class character
marg.par = ...., object of class matrix
label = ...., object of class numeric
index = ...., object of class matrix
call = ...., object of class call
seed = ...., object of class integer
)
statisticObject of class numeric, observed test statistics for each hypothesis, specified by the values of the MTP arguments test, robust, standardize, and psi0.
estimateFor the test of single-parameter null hypotheses using t-statistics (i.e., not the F-tests), the numeric vector of estimated parameters corresponding to each hypothesis, e.g. means, differences in means, regression parameters.
sampsizeObject of class numeric, number of columns (i.e. observations) in the input data set.
rawpObject of class numeric, unadjusted, marginal p-values for each hypothesis.
adjpObject of class numeric, adjusted (for multiple testing) p-values for each hypothesis (computed only if the get.adjp argument is TRUE).
conf.regFor the test of single-parameter null hypotheses using t-statistics (i.e., not the F-tests), the numeric array of lower and upper simultaneous confidence limits for the parameter vector, for each value of the nominal Type I error rate alpha (computed only if the get.cr argument is TRUE).
cutoffThe numeric matrix of cut-offs for the vector of test statistics for each value of the nominal Type I error rate alpha (computed only if the get.cutoff argument is TRUE).
rejectObject of class 'matrix', rejection indicators (TRUE for a rejected null hypothesis), for each value of the nominal Type I error rate alpha.
rawdistThe numeric matrix for the estimated nonparametric non-null test statistics distribution (returned only if keep.rawdist=TRUE and if nulldist is one of 'boot.ctr', 'boot.cs', or 'boot.qt'). This slot must not be empty if one wishes to call update to change choice of bootstrap-based null distribution.
nulldistThe numeric matrix for the estimated test statistics null distribution (returned only if keep.nulldist=TRUE); option not currently available for permutation null distribution, i.e., nulldist='perm'). By default (i.e., for nulldist='boot.cs'), the entries of nulldist are the null value shifted and scaled bootstrap test statistics, with one null test statistic value for each hypothesis (rows) and bootstrap iteration (columns).
nulldist.typeCharacter value describing which choice of null distribution was used to generate the MTP results. Takes on one of the values of the original nulldist argument in the call to MTP, i.e., 'boot.cs', 'boot.ctr', 'boot.qt', 'ic', or 'perm'.
marg.nullIf nulldist='boot.qt', a character value returning which choice of marginal null distribution was used by the MTP. Can be used to check default values or to ensure manual settings were correctly applied.
marg.parIf nulldist='boot.qt', a numeric matrix returning the parameters of the marginal null distribution(s) used by the MTP. Can be used to check default values or to ensure manual settings were correctly applied.
labelIf keep.label=TRUE, a vector storing the values used in the argument Y. Storing this object is particularly important when one wishes to update MTP objects with F-statistics using default marg.null and marg.par settings when nulldist='boot.qt'.
indexFor tests of correlation parameters a matrix corresponding to t(combn(p,2)), where p is the number of variables in X. This matrix gives the indices of the variables considered in each pairwise correlation. For all other tests, this slot is empty, as the indices are in the same order as the rows of X.
callObject of class call, the call to the MTP function.
seedAn integer or vector for specifying the state of the random number generator used to create the resampled datasets. The seed can be reused for reproducibility in a repeat call to MTP. This argument is currently used only for the bootstrap null distribution (i.e., for nulldist="boot.xx"). See ?set.seed for details.
signature(x = "MTP")
: Subsetting method for MTP class, which operates selectively on each slot of an MTP instance to retain only the data related to the specified hypotheses.
: Converts an object of class MTP to an object of class list, with an entry for each slot.
: plot methods for MTP class, produces the following graphical summaries of the results of a MTP. The type of display may be specified via the which argument.
1. Scatterplot of number of rejected hypotheses vs. nominal Type I error rate.
2. Plot of ordered adjusted p-values; can be viewed as a plot of Type I error rate vs. number of rejected hypotheses.
3. Scatterplot of adjusted p-values vs. test statistics (also known as "volcano plot").
4. Plot of unordered adjusted p-values.
5. Plot of confidence regions for user-specified parameters, by default the 10 parameters corresponding to the smallest adjusted p-values (argument top).
6. Plot of test statistics and corresponding cut-offs (for each value of alpha) for user-specified hypotheses, by default the 10 hypotheses corresponding to the smallest adjusted p-values (argument top).
The argument logscale (by default equal to FALSE) allows one to use the negative decimal logarithms of the adjusted p-values in the second, third, and fourth graphical displays. The arguments caption and sub.caption allow one to change the titles and subtitles for each of the plots (default subtitle is the MTP function call). Note that some of these plots are implemented in the older function mt.plot.
: print method for MTP class, returns a description of an object of class MTP, including sample size, number of tested hypotheses, type of test performed (value of argument test), Type I error rate (value of argument typeone), nominal level of the test (value of argument alpha), name of the MTP (value of argument method), call to the function MTP.
In addition, this method produces a table with the class, mode, length, and dimension of each slot of the MTP instance.
: summary method for MTP class, provides numerical summaries of the results of a MTP and returns a list with the following three components.
1. rejections: A data.frame with the number(s) of rejected hypotheses for the nominal Type I error rate(s) specified by the alpha argument of the function MTP. (NULL values are returned if all three arguments get.cr, get.cutoff, and get.adjp are FALSE).
2. index: A numeric vector of indices for ordering the hypotheses according to first adjp, then rawp, and finally the absolute value of statistic (not printed in the summary).
3. summaries: When applicable (i.e., when the corresponding quantities are returned by MTP), a table with six number summaries of the distributions of the adjusted p-values, unadjusted p-values, test statistics, and parameter estimates.
: update method for MTP class, provides a mechanism to re-run the MTP with different choices of the following arguments - nulldist, alternative, typeone, k, q, fdr.method, alpha, smooth.null, method, get.cr, get.cutoff, get.adjp, keep.nulldist, keep.rawdist, keep.margpar. When evaluate is 'TRUE', a new object of class MTP is returned. Else, the updated call is returned. The MTP object passed to the update method must have either a non-empty rawdist slot or a non-empty nulldist slot (i.e., must have been called with either 'keep.rawdist=TRUE' or 'keep.nulldist=TRUE').
To save on memory, if one knows ahead of time that one will want to compare different choices of bootstrap-based null distribution, then it is both necessary and sufficient to specify 'keep.rawdist=TRUE', as there is no other means of moving between null distributions other than through the non-transformed non-parametric bootstrap distribution. In this case, 'keep.nulldist=FALSE' may be used. Specifically, if an object of class MTP contains a non-empty rawdist slot and an empty nulldist slot, then a new null distribution will be generated according to the values of the nulldist= argument in the original call to MTP or any additional specifications in the call to update. On the other hand, if one knows that one wishes to only update an MTP object in ways which do not involve choice of null distribution, then 'keep.nulldist=TRUE' will suffice and 'keep.rawdist' can be set to FALSE (default settings). The original null distribution object will then be used for all subsequent calls to update.
N.B.: Note that keep.rawdist=TRUE is only available for the bootstrap-based resampling methods. The non-null distribution does not exist for the permutation or influence curve multivariate normal null distributions.
: coersion method for converting objects of class MTP to objects of class EBMTP. Slots common to both objects are taken from the object of class MTP and used to create a new object of class EBMTP. Once an object of class EBMTP is created, one may use the method EBupdate to perform resampling-based empirical Bayes multiple testing without the need for repeated resampling.
Katherine S. Pollard and Houston N. Gilbert with design contributions from Sandrine Dudoit and Mark J. van der Laan.
M.J. van der Laan, S. Dudoit, K.S. Pollard (2004), Augmentation Procedures for Control of the Generalized Family-Wise Error Rate and Tail Probabilities for the Proportion of False Positives, Statistical Applications in Genetics and Molecular Biology, 3(1). http://www.bepress.com/sagmb/vol3/iss1/art15/
M.J. van der Laan, S. Dudoit, K.S. Pollard (2004), Multiple Testing. Part II. Step-Down Procedures for Control of the Family-Wise Error Rate, Statistical Applications in Genetics and Molecular Biology, 3(1). http://www.bepress.com/sagmb/vol3/iss1/art14/
S. Dudoit, M.J. van der Laan, K.S. Pollard (2004), Multiple Testing. Part I. Single-Step Procedures for Control of General Type I Error Rates, Statistical Applications in Genetics and Molecular Biology, 3(1). http://www.bepress.com/sagmb/vol3/iss1/art13/
Katherine S. Pollard and Mark J. van der Laan, "Resampling-based Multiple Testing: Asymptotic Control of Type I Error and Applications to Gene Expression Data" (June 24, 2003). U.C. Berkeley Division of Biostatistics Working Paper Series. Working Paper 121. http://www.bepress.com/ucbbiostat/paper121
M.J. van der Laan and A.E. Hubbard (2006), Quantile-function Based Null Distributions in Resampling Based Multiple Testing, Statistical Applications in Genetics and Molecular Biology, 5(1). http://www.bepress.com/sagmb/vol5/iss1/art14/
S. Dudoit and M.J. van der Laan. Multiple Testing Procedures and Applications to Genomics. Springer Series in Statistics. Springer, New York, 2008.
MTP, MTP-methods,
EBMTP, EBMTP-methods,
[-methods, as.list-methods, print-methods, plot-methods, summary-methods, mtp2ebmtp,
ebmtp2mtp
## See MTP function: ? MTP## See MTP function: ? MTP
Summary, printing, plotting, subsetting, updating, as.list and class conversion methods were defined for the MTP and EBMTP classes. These methods provide visual and numeric summaries of the results of a multiple testing procedure (MTP) and allow one to perform some basic manipulations of objects class MTP or EBMTP.
Several of the methods with the same name will work on objects of their respective class. One exception to this rule is the difference between update and EBupdate (described below). Because of the differences in the testing procedures, separately named methods were chosen to clearly delineate which method was being applied to which type of object.
: Subsetting method for MTP and EBMTP classes, which operates selectively on each slot of an MTP or EBMTP instance to retain only the data related to the specified hypotheses.
: Converts an object of class MTP or EBMTP to an object of class list, with an entry for each slot.
: plot methods for MTP and EBMTP classes, produces the following graphical summaries of the results of a MTP. The type of display may be specified via the which argument.
1. Scatterplot of number of rejected hypotheses vs. nominal Type I error rate.
2. Plot of ordered adjusted p-values; can be viewed as a plot of Type I error rate vs. number of rejected hypotheses.
3. Scatterplot of adjusted p-values vs. test statistics (also known as volcano plot).
4. Plot of unordered adjusted p-values.
Only for objects of class MTP:
5. Plot of confidence regions for user-specified parameters, by default the 10 parameters corresponding to the smallest adjusted p-values (argument top).
6. Plot of test statistics and corresponding cut-offs (for each value of alpha) for user-specified hypotheses, by default the 10 hypotheses corresponding to the smallest adjusted p-values (argument top).
Plots (5) and (6) are not available for objects of class EBMTP because the function EBMTP returns only adjusted p-values and not confidence regions of cut-offs. The argument logscale (by default equal to FALSE) allows one to use the negative decimal logarithms of the adjusted p-values in the second, third, and fourth graphical displays. The arguments caption and sub.caption allow one to change the titles and subtitles for each of the plots (default subtitle is the MTP function call). Note that some of these plots are implemented in the older function mt.plot.
: print method for MTP and EBMTP classes, returns a description of an object of either class, including sample size, number of tested hypotheses, type of test performed (value of argument test), Type I error rate (value of argument typeone), nominal level of the test (value of argument alpha), name of the MTP (value of argument method), call to the function MTP or EBMTP.
In addition, this method produces a table with the class, mode, length, and dimension of each slot of the MTP or EBMTP instance.
: summary method for MTP and EBMTP classes, provides numerical summaries of the results of a MTP and returns a list with the following three components.
1. rejections: A data.frame with the number(s) of rejected hypotheses for the nominal Type I error rate(s) specified by the alpha argument of the function MTP or EBMTP. (For objects of class MTP, NULL values are returned if all three arguments get.cr, get.cutoff, and get.adjp are FALSE).
2. index: A numeric vector of indices for ordering the hypotheses according to first adjp, then rawp, and finally the absolute value of statistic (not printed in the summary).
3. summaries: When applicable (i.e., when the corresponding quantities are returned by MTP or EBMTP), a table with six number summaries of the distributions of the adjusted p-values, unadjusted p-values, test statistics, and parameter estimates.
: update methods for MTP class, respectively, provides a mechanism to re-run the MTP with different choices of the following arguments - nulldist, alternative, typeone, k, q, fdr.method, alpha, smooth.null, method, get.cr, get.cutoff, get.adjp, keep.nulldist, keep.rawdist, keep.margpar. When evaluate is 'TRUE', a new object of class MTP is returned. Else, the updated call is returned. The MTP object passed to the update method must have either a non-empty rawdist slot or a non-empty nulldist slot (i.e., must have been called with either 'keep.rawdist=TRUE' or 'keep.nulldist=TRUE').
: update method for EBMTP class, provides a mechanism to re-run the MTP with different choices of the following arguments - nulldist, alternative, typeone, k, q, alpha, smooth.null, bw, kernel, prior, keep.nulldist, keep.rawdist, keep.falsepos, keep.truepos, keep.errormat, keep.margpar. When evaluate is 'TRUE', a new object of class EBMTP is returned. Else, the updated call is returned. The EBMTP object passed to the update method must have either a non-empty rawdist slot or a non-empty nulldist slot (i.e., must have been called with either 'keep.rawdist=TRUE' or 'keep.nulldist=TRUE').
Additionally, when calling EBupdate for any Type I error rate other than FWER, the typeone argument must be specified (even if the original object did not control FWER). For example,
typeone="fdr", would always have to be specified, even if the original object also controlled the FDR. In other words, for all function arguments, it is safest to always assume that you
are updating from the EBMTP default function settings, regardless of the original call to the EBMTP function. Currently, the main advantage of the EBupdate method is that it prevents the need for repeated estimation of the test statistics null distribution.
To save on memory, if one knows ahead of time that one will want to compare different choices of bootstrap-based null distribution, then it is both necessary and sufficient to specify 'keep.rawdist=TRUE', as there is no other means of moving between null distributions other than through the non-transformed non-parametric bootstrap distribution. In this case, 'keep.nulldist=FALSE' may be used. Specifically, if an object of class MTP or EBMTP contains a non-empty rawdist slot and an empty nulldist slot, then a new null distribution will be generated according to the values of the nulldist= argument in the original call to (EB)MTP or any additional specifications in the call to (EB)update. On the other hand, if one knows that one wishes to only update an (EB)MTP object in ways which do not involve choice of bootstrap null distribution, then 'keep.nulldist=TRUE' will suffice and 'keep.rawdist' can be set to FALSE (default settings). The original null distribution object will then be used for all subsequent calls to update.
N.B.: Note that keep.rawdist=TRUE is only available for the bootstrap-based resampling methods. The non-null distribution does not exist for the permutation or influence curve multivariate normal null distributions.
: coersion method for converting objects of class MTP to objects of class EBMTP. Slots common to both objects are taken from the object of class MTP and used to create a new object of class EBMTP. Once an object of class EBMTP is created, one may use the method EBupdate to perform resampling-based empirical Bayes multiple testing without the need for repeated resampling.
: coersion method for converting objects of class EBMTP to objects of class MTP. Slots common to both objects are taken from the object of class EBMTP and used to create a new object of class MTP. Once an object of class MTP is created, one may use the method update to perform resampling-based multiple testing (as would have been done with calls to MTP) without the need for repeated resampling.
Katherine S. Pollard and Houston N. Gilbert with design contributions from Sandrine Dudoit and Mark J. van der Laan.