A function to perform resampling-based multiple hypothesis testing
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)
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:
|
Object of class |
|
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. |
|
Object of class |
|
Object of class |
|
Object of class |
|
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 |
|
The numeric matrix of cut-offs for the vector of test statistics for each value of the nominal Type I error rate |
|
Object of class |
|
The numeric matrix for the estimated nonparametric non-null test statistics distribution (returned only if |
|
The numeric matrix for the estimated test statistics null distribution (returned only if |
|
Character value describing which choice of null distribution was used to generate the MTP results. Takes on one of the values of the original |
|
If |
|
If |
|
Object of class |
|
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.
#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
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.