bonferroni.RdFunction to carry out the Bonferroni method.
bonferroni(p, adjust = "none", R, m,
size = 10000, threshold, side = 2, batchsize, nearpd = TRUE, ...)vector of length \(k\) with the (one- or two-sided) p-values to be combined.
character string to specify an adjustment method to account for dependence. The default is "none", in which case no adjustment is applied. Methods "nyholt", "liji", "gao", or "galwey" are adjustments based on an estimate of the effective number of tests (see meff). Adjustment method "empirical" uses an empirically-derived null distribution using pseudo replicates. See ‘Details’.
a \(k \times k\) symmetric matrix that reflects the dependence structure among the tests. Must be specified if adjust is set to something other than "none". See ‘Details’.
optional scalar (between 1 and \(k\)) to manually specify the effective number of tests (instead of estimating it via one of the methods described above).
size of the empirically-derived null distribution. Can also be a numeric vector of sizes, in which case a stepwise algorithm is used. This (and the following arguments) are only relevant when adjust = "empirical".
numeric vector to specify the significance thresholds for the stepwise algorithm (only relevant when size is a vector).
scalar to specify the sidedness of the \(p\)-values that are used to simulate the null distribution (2, by default, for two-sided tests; 1 for one-sided tests).
optional scalar to specify the batch size for generating the null distribution. When unspecified (the default), this is done in a single batch.
logical indicating if a negative definite R matrix should be turned into the nearest positive definite matrix (only relevant when adjust = "empirical").
other arguments.
Bonferroni Method
By default (i.e., when adjust = "none"), the function applies the Bonferroni method to the \(p\)-values. Letting \(p_1, p_2, \ldots, p_k\) denote the individual (one- or two-sided) \(p\)-values of the \(k\) hypothesis tests to be combined, the combined \(p\)-value is then computed with \[p_c = \min(1, \min(p_1, p_2, \ldots, p_k) \times k).\]
The Bonferroni method does not assume that the \(p\)-values to be combined are independent. However, if the \(p\)-values are not independent, the method can become quite conservative (not reject often enough), depending on the dependence structure among the tests. In this case, one can adjust the method to account for such dependence (to bring the Type I error rate closer to some desired nominal significance level).
Adjustment Based on the Effective Number of Tests
When adjust is set to "nyholt", "liji", "gao" or "galwey", the Bonferroni method is adjusted based on an estimate of the effective number of tests (see meff for details on these methods for estimating the effective number of tests). In this case, argument R needs to be set to a matrix that reflects the dependence structure among the tests.
There is no general solution for constructing such a matrix, as this depends on the type of test that generated the \(p\)-values and the sidedness of these tests. If the \(p\)-values are obtained from tests whose test statistics can be assumed to follow a multivariate normal distribution and a matrix is available that reflects the correlations among the test statistics, then the mvnconv function can be used to convert this correlation matrix into the correlations among the (one- or two-sided) \(p\)-values, which can then be passed to the R argument. See ‘Examples’.
Once the effective number of tests, \(m\), is estimated based on R using one of the four methods described above, the combined \(p\)-value is then computed with \[p_c = \min(1, \min(p_1, p_2, \ldots, p_k) \times m).\]
Alternatively, one can also directly specify the effective number of tests via the m argument (e.g., if some other method not implemented in the poolr package is used to estimate the effective number of tests). Argument R is then irrelevant and doesn't need to be specified.
Adjustment Based on an Empirically-Derived Null Distribution
When adjust = "empirical", the combined \(p\)-value is computed based on an empirically-derived null distribution using pseudo replicates (using the empirical function). This is appropriate if the test statistics that generated the \(p\)-values to be combined can be assumed to follow a multivariate normal distribution and a matrix is available that reflects the correlations among the test statistics (which is specified via the R argument). In this case, test statistics are repeatedly simulated from a multivariate normal distribution under the joint null hypothesis, converted into one- or two-sided \(p\)-values (depending on the side argument), and the Bonferroni method is applied. Repeating this process size times yields a null distribution based on which the combined \(p\)-value can be computed, or more precisely, estimated, since repeated applications of this method will yield (slightly) different results. To obtain a stable estimate of the combined \(p\)-value, size should be set to a large value (the default is 10000, but this can be increased for a more precise estimate). If we consider the combined \(p\)-value an estimate of the ‘true’ combined \(p\)-value that would be obtained for a null distribution of infinite size, we can also construct a 95% (pseudo) confidence interval based on a binomial distribution.
If batchsize is unspecified, the null distribution is simulated in a single batch, which requires temporarily storing a matrix with dimensions [size,k]. When size*k is large, allocating the memory for this matrix might not be possible. Instead, one can specify a batchsize value, in which case a matrix with dimensions [batchsize,k] is repeatedly simulated until the desired size of the null distribution has been obtained.
One can also specify a vector for the size argument, in which case one must also specify a corresponding vector for the threshold argument. In that case, a stepwise algorithm is used that proceeds as follows. For j = 1, ..., length(size),
estimate the combined \(p\)-value based on size[j]
if the combined \(p\)-value is \(\ge\) than threshold[j], stop (and report the combined \(p\)-value), otherwise go back to 1.
By setting size to increasing values (e.g., size = c(1000, 10000, 100000)) and threshold to decreasing values (e.g., threshold = c(.10, .01, 0)), one can quickly obtain a fairly accurate estimate of the combined \(p\)-value if it is far from significant (e.g., \(\ge\) .10), but hone in on a more accurate estimate for a combined \(p\)-value that is closer to 0. Note that the last value of threshold should be 0 (and is forced to be inside of the function), so that the algorithm is guaranteed to terminate (hence, one can also leave out the last value of threshold, so threshold = c(.10, .01) would also work in the example above). One can also specify a single threshold (which is replicated as often as necessary depending on the length of size).
An object of class "poolr". The object is a list containing the following components:
combined \(p\)-value.
confidence interval for the combined \(p\)-value (only when adjust = "empirical"; otherwise NULL).
number of \(p\)-values that were combined.
estimate of the effective number of tests (only when adjust is one of "nyholt", "liji", "gao" or "galwey"; otherwise NULL).
chosen adjustment method.
value of the (adjusted) test statistic.
name of calling function.
The method underlying adjust = "empirical" assumes that the test statistics that generated the \(p\)-values to be combined follow a multivariate normal distribution. Hence, the matrix specified via R must be positive definite. If it is not and nearpd = TRUE, it will be turned into one (based on Higham, 2002, and a slightly simplified version of nearPD from the Matrix package).
Bland, J. M., & Altman, D. G. (1995). Multiple significance tests: The Bonferroni method. British Medical Journal, 310(6973), 170.
Cinar, O. & Viechtbauer, W. (2022). The poolr package for combining independent and dependent p values. Journal of Statistical Software, 101(1), 1–42. https://doi.org/10.18637/jss.v101.i01
Higham, N. J. (2002). Computing the nearest correlation matrix: A problem from finance. IMA Journal of Numerical Analysis, 22(3), 329–343.
# copy p-values and LD correlation matrix into p and r
# (see help(grid2ip) for details on these data)
p <- grid2ip.p
r <- grid2ip.ld
# apply the Bonferroni method
bonferroni(p)
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: none
#> combined p-value: 0.03881585
# use mvnconv() to convert the LD correlation matrix into a matrix with the
# correlations among the (two-sided) p-values assuming that the test
# statistics follow a multivariate normal distribution with correlation
# matrix r (note: 'side = 2' by default in mvnconv())
mvnconv(r, target = "p", cov2cor = TRUE)[1:5,1:5] # show only rows/columns 1-5
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.00000000 0.02160864 0.022809124 0.010804322 0.097238896
#> [2,] 0.02160864 1.00000000 0.013205282 0.000000000 0.042016807
#> [3,] 0.02280912 0.01320528 1.000000000 0.006002401 0.006002401
#> [4,] 0.01080432 0.00000000 0.006002401 1.000000000 0.000000000
#> [5,] 0.09723890 0.04201681 0.006002401 0.000000000 1.000000000
# adjustment based on estimates of the effective number of tests
bonferroni(p, adjust = "nyholt", R = mvnconv(r, target = "p", cov2cor = TRUE))
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: effective number of tests (m = 22; Nyholt, 2004)
#> combined p-value: 0.0371
bonferroni(p, adjust = "liji", R = mvnconv(r, target = "p", cov2cor = TRUE))
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: effective number of tests (m = 21; Li & Ji, 2005)
#> combined p-value: 0.0354
bonferroni(p, adjust = "gao", R = mvnconv(r, target = "p", cov2cor = TRUE))
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: effective number of tests (m = 23; Gao, 2008)
#> combined p-value: 0.0388
bonferroni(p, adjust = "galwey", R = mvnconv(r, target = "p", cov2cor = TRUE))
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: effective number of tests (m = 20; Galwey, 2009)
#> combined p-value: 0.0338
# setting argument 'm' manually
bonferroni(p, m = 12)
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: effective number of tests (m = 12; user-defined)
#> combined p-value: 0.0203
# adjustment based on an empirically-derived null distribution (setting the
# seed for reproducibility)
set.seed(1234)
bonferroni(p, adjust = "empirical", R = r)
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: empirical distribution (size = 10000)
#> combined p-value: 0.0323 (95% CI: 0.0289, 0.036)
# generate the empirical distribution in batches of size 100
bonferroni(p, adjust = "empirical", R = r, batchsize = 100)
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: empirical distribution (size = 10000)
#> combined p-value: 0.0322 (95% CI: 0.0288, 0.0358)
# using the stepwise algorithm
bonferroni(p, adjust = "empirical", R = r, size = c(1000, 10000, 100000), threshold = c(.10, .01))
#> combined p-values with: Bonferroni method
#> number of p-values combined: 23
#> minimum p-value: 0.002
#> adjustment: empirical distribution (size = 10000)
#> combined p-value: 0.0305 (95% CI: 0.0272, 0.0341)