dBH_mvgauss.RddBH_mvgauss computes the rejection set of \(dBH_\gamma(\alpha)\) and \(dBH^2_\gamma(\alpha)\), as well as
\(dSU_{\gamma, \Delta}(\alpha)\) and \(dSU^2_{\gamma, \Delta}(\alpha)\) a broad class of step-up procedures for
multivariate z-statistics \(z\sim N(\mu, \Sigma)\). dBH_mvgauss can handle both one-sided tests
\(H_i: \mu_i \le 0\) or \(H_i: \mu_i \ge 0\), and two-sided tests \(H_i: \mu_i = 0\).
dBH_mvgauss( zvals, Sigma = NULL, Sigmafun = NULL, vars = NULL, side = c("right", "left", "two"), alpha = 0.05, gamma = NULL, niter = 1, tautype = "QC", avals = NULL, avals_type = c("BH", "geom", "bonf", "manual"), geom_fac = 2, eps = 0.05, qcap = 2, gridsize = 20, exptcap = 0.9, is_safe = NULL, verbose = FALSE )
| zvals | a vector of z-values. The marginal variances do not have to be 1. |
|---|---|
| Sigma | the covariance matrix. For very large m, it is preferable to use |
| Sigmafun | the function that takes the row index as the input and outputs the i-th row of the covariance matrix
divided by \(\Sigma_{i, i}\). This is an alternative for |
| vars | the vector of marginal variances \((\sigma_1^2, \ldots, \sigma_m^2)\). Used only with |
| side | a string that takes values in {"right", "left", "two"}, with "right" for one-sided tests \(H_i: \mu_i \le 0\), "left" for one-sided tests, \(H_i: \mu_i \ge 0\), and "two" for two-sided tests \(H_i: \mu_i = 0\). |
| alpha | the target FDR level. |
| gamma | the parameter for the dBH and dSU procedures. The default is |
| niter | the number of iterations. In the current version it can only be 1, for dBH/dSU, or 2, for \(dBH^2\)/\(dSU^2\). |
| tautype | the type of tau function. In the current version, only "QC" (q-value-calibration) is supported with \(\tau(c; X) = cR_{BH}(c) / m\). |
| avals | the a-values in the step-up procedures defined in Appendix C.1. The default is NULL, in which case
|
| avals_type | a string that takes values in {"BH", "geom", "bonf", "manual"}, which determines the type of thresholds in the step-up procedures. See Details. |
| geom_fac | a real number that is larger than 1. This is the growing rate of thresholds when
|
| eps | a real number in [0, 1], which is used to determine \(t_{hi}\) discussed in Section 4.2. |
| qcap | a real number that is larger than 1. It is used to filter out hypotheses with q-values above
|
| gridsize | an integer for the size of the grid used when |
| exptcap | a real number in [0, 1]. It is used by \(dBH^2\)/\(dSU^2\) to filter out hypotheses with
\(g_i*(q_i | S_i)\) below |
| is_safe | a logical or NULL indicating whether the procedure is taken as safe. DON'T set |
| verbose | a logical indicating whether a progress bar is shown. |
a list with the following attributes
rejs: the indices of rejected hypotheses (after the randomized pruning step if any);
initrejs: the indices of rejected hypotheses (before the randomized pruning step if any);
cand: the set of candidate hypotheses for which \(g_i^{*}(q_i | S_i)\) is evaluated;
expt: \(g_i^{*}(q_i | S_i)\) for each hypothesis in cand;
safe: TRUE iff the procedure is safe;
secBH: TRUE iff the randomzied pruning step (a.k.a. the secondary BH procedure) is invoked;
secBH_fac: a vector that gives \(\hat{R}_i / R_{+}\) that is defined in Section 2.2. It only shows up
in the output if secBH = TRUE.
dBH_mvgauss supports two types of inputs for the covariance matrix.
When the covariance matrix fits into the memory, it can be inputted through Sigma. The diagonals do not
have to be 1. In this case, Sigmafun and vars should be left as their default;
When the covariance matrix does not fit into the memory, it can be inputted through Sigmafun, for the
correlation matrix, and vars, for marginal variances. Sigmafun should be a function with a single input
i that gives the row index and a vector output \(\Sigma_{i, }/\sqrt{\Sigma_{i, i}\Sigma_{j, j}}\) where
\(\Sigma_{i, }\) is the i-th row of \(\Sigma\). The marginal variances \((\Sigma_{1, 1}, \ldots, \Sigma_{m, m})\)
are inputted through vars. If all marginal variances are 1, vars can be left as its default.
dBH_mvgauss can handle all dSU procedures that are defined in Appendix C.1. with
$$\Delta_{\alpha}(r) = \frac{\alpha a_{\ell}}{m}, r\in [a_{\ell}, a_{\ell + 1}), \ell = 0, 1, \ldots, L$$
for any set of integer a-values \(1\le a_1 < \ldots < a_L\le m\) (with \(a_0 = 0, a_{L+1} = m+1\)). There are two-ways to input the a-values.
(Recommended) use avals_type while leave avals as its default. Three types of avals_type
are supported.
When avals_type = "BH", the BH procedure is used (i.e. \(a_\ell = \ell, \ell = 1, \ldots, m\));
When avals_type = "geom", the geometrically increasing a-values that are defined in Appendix C.1.
are used with the growing rate specified by geom_fac, whose default is 2;
When avals_type = "bonf", the Bonferroni procedure is used (i.e. \(a_1 = 1, L = 1\)).
(Not recommended) use avals while set avals_type = "manual". This option allows any set of
increasing a-values. However, it can be much slower than the recommended approach.
# Generate mu and Sigma for an AR process n <- 100 rho <- 0.8 Sigma <- rho^(abs(outer(1:n, 1:n, "-"))) mu1 <- 2.5 nalt <- 10 mu <- c(rep(mu1, nalt), rep(0, n - nalt)) # Generate the z-values set.seed(1) zvals <- rep(NA, n) zvals[1] <- rnorm(1) for (i in 2:n){ zvals[i] <- zvals[i - 1] * rho + rnorm(1) * sqrt(1 - rho^2) } zvals <- zvals + mu # Run dBH_1(\alpha) for one-sided tests alpha <- 0.05 res <- dBH_mvgauss(zvals = zvals, Sigma = Sigma, side = "right", alpha = alpha, gamma = 1, niter = 1, avals_type = "BH") # Run dBH_1(\alpha) for two-sided tests res <- dBH_mvgauss(zvals = zvals, Sigma = Sigma, side = "right", alpha = alpha, gamma = 1, niter = 1, avals_type = "BH") # Run dBH^2_1(\alpha) for one-sided tests alpha <- 0.05 res <- dBH_mvgauss(zvals = zvals, Sigma = Sigma, side = "right", alpha = alpha, gamma = 1, niter = 2, avals_type = "BH") # Run dBH^2_1(\alpha) for one-sided tests res <- dBH_mvgauss(zvals = zvals, Sigma = Sigma, side = "right", alpha = alpha, gamma = 1, niter = 2, avals_type = "BH") # Run dSU_1(\alpha) with the geometrically increasing a-values for one-sided tests res <- dBH_mvgauss(zvals = zvals, Sigma = Sigma, side = "right", alpha = alpha, gamma = 1, niter = 1, avals_type = "geom", geom_fac = 2) # Run dBY(\alpha) for one-sided tests res <- dBH_mvgauss(zvals = zvals, Sigma = Sigma, side = "right", alpha = alpha, gamma = NULL, niter = 1, avals_type = "BH") # Run dBY^2(\alpha) for one-sided tests res <- dBH_mvgauss(zvals = zvals, Sigma = Sigma, side = "right", alpha = alpha, gamma = NULL, niter = 2, avals_type = "BH") # Input the covariance matrix through \code{Sigmafun} Sigmafun <- function(i){ rho^(abs(1:n - i)) } vars <- rep(1, n) res1 <- dBH_mvgauss(zvals = zvals, Sigmafun = Sigmafun, vars = vars, side = "right", alpha = alpha, gamma = NULL, niter = 1, avals_type = "BH") res2 <- dBH_mvgauss(zvals = zvals, Sigma = Sigma, side = "right", alpha = alpha, gamma = NULL, niter = 1, avals_type = "BH") identical(res1, res2)#> [1] TRUE