# Ensemble power or sample size under selected control of the FDP This is a function for calculating two differing notions of power, or deriving sample sizes for specified requisite power in multiple testing experiments under a variety of methods for control of the distribution of the False Discovery Proportion (FDP). More specifically, one can choose to control the FDP distribution according to control of its (i) mean, e.g. the usual BH-FDR procedure, or via the probability that it exceeds a given value, delta, via (ii) the Romano procedure, or via (iii) my procedure based upon asymptotic approximation. Likewise, we can think of the power in multiple testing experiments in terms of a summary of the distribution of the True Positive Proportion (TPP). The package will compute power, sample size or any other missing parameter required for power based upon (i) the mean of the TPP which is the average power (ii) the probability that the TPP exceeds a given value, lambda, via my asymptotic approximation procedure. The theoretical results are described in Izmirlian, G. (2020), and an applied paper describing the methodology with a simulation study is in preparation. ```r pwrFDR(effect.size, n.sample, r.1, alpha, delta = NULL, groups = 2, N.tests, average.power, TPX.power, lambda, type = c("paired","balanced","unbalanced"), grpj.per.grp1=NULL, corr.struct=list(type=c("CS-Blocks","Toeplitz-Blocks"), block.size=NA, rho = NA), FDP.control.method = c("BHFDR","BHFDX","Romano","Auto","both","Holm","Hochberg", "Bonferroni"), distopt=1, method = c("Theoretical","simulation"), n.sim = 1000, temp.file, control=list(tol=1e-08, max.iter=c(1000, 20), sim.level=2, low.power.stop=TRUE, FDP.meth.thresh=FDP.cntl.mth.thrsh.def, ast.le.a=TRUE, verb=FALSE, show.footer=TRUE)) ``` ## Arguments - `effect.size`: The effect size (mean over standard deviation) for test statistics having non-zero means. Assumed to be a constant (in magnitude) over non-zero mean test statistics. - `n.sample`: The number of experimental replicates. Required for calculation of power - `r.1`: The proportion of simultaneous tests that are non-centrally located - `alpha`: The false discovery rate (in the BH case) or the upper bound on the probability that the FDP exceeds delta (BHFDX and Romano case) - `delta`: If the "FDP.control.method" is set to 'Romano' or 'BHFDX', then the user can set the exceedance thresh-hold for the FDP tail probability control $P\{ FDP > \delta \} < \alpha$. The default value is $\alpha$. - `groups`: The number of experimental groups to compare. Must be integral and >=1. The default value is 2. - `N.tests`: The number of simultaneous hypothesis tests. - `average.power`: The desired average power. Sample size calculation requires specification of either 'average.power' or 'TPX.power'. - `TPX.power`: The desired tp-power (see details for explanation). Sample size calculation requires specification of either 'average.power' or 'TPX.power'. - `lambda`: The tp-power threshold, required when calculating the tp-power (see details for explanation) or when calculating the sample size required for tp-power. - `type`: A character string specifying, in the groups=2 case, whether the test is 'paired', 'balanced', or 'unbalanced' and in the case when groups >=3, whether the test is 'balanced' or 'unbalanced'. The default in all cases is 'balanced'. Left unspecified in the one sample (groups=1) case. - `grpj.per.grp1`: Required when `type`="unbalanced", specifies the group 0 to group 1 ratio in the two group case, and in the case of 3 or more groups, the group j to group 1 ratio, where group 1 is the group with the largest effect under the alternative hypothesis. - `corr.struct`: Specifies a block correlation structure between test statistics which is used in both the simulation routine, and in the computations based upon asymptotic approximation, e.g. the AFDX control and the ATPP method. Its form is specified via the following named elements. "type": a string taking the value "CS-Blocks" for compound compound symmetry within blocks or "Toeplitz-Blocks" for toeplitz within blocks, respectively. "block.size": the size of the correlated blocks "rho": When type="CS-Blocks", then 'rho' is a correlation of length 1 when type="Toeplitz-Blocks", then 'rho' is a vector of correlations of length 'block.size' - 1 - `FDP.control.method`: A character string specifying how the false discovery proportion (FDP) is to be controlled. You may specify the whole word or any shortened uniquely identifying truncation. "BHFDR": the usual BH-FDR "BHFDX": guarantees control of the FDX e.g. probability that the FDP exceeds delta will be less than alpha. Uses asymptotic approximationto the distribution of the FDP to find alpha* such that the BHFDR procedure controls the FDX at alpha. Only slightly more conservative than BHFDR control at alpha. Allows dependence which must be specified. Not guaranteed to have a solution in all cases, in which case Romano's procedure is used as a fall-back. "Romano": use Romano's method which guarantees control of the FDX, e.g. the probability that the FDP exceeds delta will be less than alpha. "Auto": in 'FixedPoint' mode, the program will use its own wisdom to determine which choice above to make. The order of conservatism is Bonferroni > Holm > Hochberg > Romano > BHFDX > BHFDR, but BHFDR offers only expected control, the most conservative 3 control the FWER and the other two guarantee bounds on the excedance probabilty. If the distribution of the FDP is nearly degenerate, then BHFDR is the best option. Otherwise, if it can be reliably used, BHFDX would be the best choice. The 'effective' denominator, gamma*N.tests, in the CLT determines when the approximation is good enough and the asymptotic standard error of the FDP determines when the distribution is dispersed enough to matter. Use "Auto" to run through these checks and determine the best. A return argument, 'Auto', displays the choice made. See output components and details. "both": in 'simulation' mode, compute statistics R and T under BHFDX and Romano (in addition to BHFDR). Corresponding statistics are denoted R.st, T.st corresponding to BHFDX control of the FDP, and R.R and T.R corresponding to Romano control of the FDP. If sim.level is set to 2, (default) the statistics R.st.ht and T.st.ht, which are the number rejected and number true positives under BHFDX where r_0 = 1-r_1, gamma, and alpha.star have been estimated from the P-value data and then alpha.star computed from these. "Holm": Use Holm's step-down procedure which guarantees control of the family-wise error rate (FWER) "Hochberg": Use Hochberg's step-up procedure which guarantees control of the FWER if the tests are independent. Less conservative than Holm and Bonferroni but more conservative than all other procedures here. "Bonferroni": Use the Bonferroni procedure. Guarantees control of the FWER. the most conservative of procedures offered here. - `distopt`: Test statistic distribution in among null and alternatively distributed sub-populations. distopt=0 gives normal (2 groups), distop=1 gives t- (2 groups) and distopt=2 gives F- (2+ groups) - `method`: Specify the method whereby the average power is calculated. You may specify the whole word or any unqiuely indentifying truncation. "Theoretical": for all procedures but Bonferroni, use the fixed point method, e.g., first find the solution to the equation u = G(alpha u) where G is the CDF of the pooled P-values. This solution gives 'gamma', the positive proportion. The average power and TPX power can then be determined e.g. average.power = G_1( gamma alpha), where G_1 is the CDF of the P-values corresponding to statistics drawn under H_A. TPX power is determined via asymptotic approximation. "simulation": uses brute force simulation to determine the average power and TPX power. - `control`: Optionally, a list with components with the following components: 'tol' is a convergence criterion used in iterative methods which is set to 1e-8 by default. 'max.iter' is an iteration limit, set to 20 for the iterated function limit and 1000 for all others by default. 'sim.level' sim level 2 (default) stipulates, when FDP.control.method is set to "BHFDX", or "both", R.st.ht and T.st.ht are computed in addition to R.st and T.st (see above). 'low.power.stop' in simulation option, will result in an error message if the power computed via FixedPoint method is too low, which result in no solution for the BHFDX option. Default setting is TRUE. Set to FALSE to over-ride this behavior. 'FDP.meth.thresh' fine-tunes the 'Auto' voodoo (see above). Leave this alone. 'ast.le.a' leaving this at the default value TRUE forces 'alpha.star', the solution under FDP.method.control="BHFDX", to be less than the specified 'alpha'. 'verb' vebosity level. 'show.footer' optionally, show a textual output - `n.sim`: If 'simulation' method is chosen you may specify number of simulations. Default is 1000. - `temp.file`: If 'simulation' method is chosen you may specify a tempfile where the current simulation replicate is updated. Very usefull for batch runs. You can use the included utility 'gentempfilenm' ## Details This function will compute one of a variety of ensemble powers under a given choice of FDP control methods. The underlying model assumes that the m simultaneous test statistics are i.i.d., each being formed from k samples which can be paired (k=2), balanced or unbalanced (k>=2), k=1,2,..., and distributed according to one of the available relevent distribution types (see above). The location parameter for each of the statistical tests is either 0 (null hypothesis) or a specified constant effect size (alternative hypothesis), with the identity of these two possibilities in each of the m cases being an i.i.d. unmeasured latent bernouli variable with density r.1, the mixing proportion. The m simultaneous statistical tests partition into those which are distributed according to the alternative, numbering M_m, and those distributed according to the NULL, numbering m-M_m. Once a selected thresholding method is applied, the m statistics can also be partitioned into those which are called significant, numbering R_m, and those which are not, numbering m-R_m. Each of the test statistics is thus given two labels, alternative hypothesis membership and whether a significant call was made. Of the R_m significant calls, T_m are true positives and V_m are false positives. This results in the following table. |||||| |--:|--:|:-:|:-:|:-:| |1.||rej H0|acc H0|row Total| |2.|H0 is FALSE|T_m|M.1-T|M_m| |3.|H0 is TRUE|R_m-T_m|(m-M_m)-(R_m-T_m)|m - M_m| |4.|col Total|R|m-R_m|m| The ratio of the false positive count to the significant call count, V_m/R_m, is called the False Discovery Proportion (FDP). Thresholding methods which result in the most reproducibility seek to control the FDP distribution. The most well known is the Benjamini-Hochberg False Discovery Rate (BH-FDR) procedure. It guarantees that the FDR, which is the expected FDP, will be less than a stipulated alpha $E[ V_m / R_m ] < \alpha$ While it is true that for large m, the distribution of the FDP, V_m/M_m will become spiked at its mean, $(1-r_1)\alpha$, in many commonly occuring situations, there will still be non-negligible dispersion in the distribution of the FDP. For this reason, any validity promised by the BH-FDR procedure does not actually apply on a case to case basis, and individual FDP's may differ non-negligibly from the FDR. For this reason, the function supplies two other methods of FDP control in addition to FDP.control.method="BHFDR". These two alternate methods, FDP.control.method="Romano" and FDP.control.method="BHFDX" guarantee control of the tail probability of the FDP distribution: $P\{ V_m/R_m > \delta \} < \alpha$ The lower bound $\delta$ is left arbitrary for greater flexibility, $\delta=\alpha$ being the default. There is also an automatic option, FDP.control.method="Auto", which lets the function decide which of the three FDP control methods is the most advisable in a given situation. The two tail probability control options are preferred when the standard error of the FDP exceeds a cutoff given in the default 'control' settings: se[V_m/R_m] / alpha > FDP.cntl.mth.thrsh.def[1] The default is 10%. When the standard error to alpha ratio is 10% or less then the BHFDR, being the least conservative, is preferred. When the se to alpha ratio is 10% or more, then Romano and BHFDX are decided between, with the BHFDX (asymptotic approximation) being less conservative than Romano and therefor preferred if the CLT approximation is adequate. This will be the case provided $m$ is large enough, $m \geq$ FDP.cntl.mth.thrsh.def[2]. The default is 50. The concept of ensemble power for the purposes of this function, concern the distribution of the true positive proportion (TPP), $T_m/M_m$. The most well known is the average power, which is the expected value of the TPP, which is called the true positve rate (TPR): $E[ T_m/M_m ] = average power$ For large m, the distribution of the TPP will be spiked at its mean, which is the asymptotic average power. This is used in the function in the average power computation. As was the case for the FDP, there are many commonly occuring situations when the distribution of the TPP will still be non-negligibly dispersed. For this reason, we provide an alternate notion of power which is based upon the tail probability of the TPP distribution: $P\{ T_m/M_m > \lambda \} = tp-power$ This is computed via asymptotic approximation and also requires that $m$ be large enough: $m > 50$. The user decides when the tp-power is to be preferred. A good check is to look at the ratio of the se[TPP] to average power ratio which is the sigma.rtm.TPP/average.power/N.tests^0.5 If this ratio is unacceptibly large (10% or so) than the tp-power is preferred. For the "FixedPoint" method (default) and for any specified choice of FDP.control method, the function can be used in the following ways: 1. Specify 'n.sample', 'effect.size', 'r.1' and 'alpha'. Calculates 'average power' 2. Specify 'n.sample', 'effect.size', 'r.1', 'alpha' and 'lambda'. 'N.tests' is also required. The function wil calculate the 'TPX.power' in addition to the 'average power'. 3. Specify the 'average.power' or the pair 'TPX.power' and 'lambda'. Specify all but one of the parameters, 'n.sample', 'effect.size', 'r.1' and 'alpha'. The function will calculate the value of the missing parameter required for the specified 'average power' or 'tp-power'. Note: a solution is guaranteed for missing 'n.sample' and missing 'effect.size', but not necessarily for missing 'r.1' or 'alpha'. The Holm, Hochberg and Bonferroni procedures for controlling the FWER under the assumption of independent test statistics are also provided for sake of completeness. ## Returns An object of class "pwr" with with components including: - **call**: The call which produced the result - **average.power**: Resulting average power. - **TPX.power**: When 'lambda' is specified, the tp-power is also computed - **L.eq**: The lambda at which the tp-power and average-power are equal. - **n.sample**: If 'n.sample' is missing from the argument list, then the sample size required for the specified average- or lambda- power. - **alpha.star**: If 'FDP.control.method' was set to "BHFDX" or it resulted from the "Auto" setting, the alpha at which the probability that the FDP exceeds alpha.star is less than or equal to the originally specified alpha. - **c.g**: The FDP control method threshold on the scale of the test statistics. - **gamma**: The proportion of all 'm' tests declared significant. - **objective**: Result of optimization yielding the average or tp- power. - **err.III**: Mass on the wrong side of the threshold. - **sigma.rtm.ToM**: Asymptotic standard deviation of the true positive fraction. - **Auto**: If 'FDP.control.method' was set to "Auto", this returns the resulting choice (a string) which was made internally. - **se.by.a**: The ratio of the standard error of the FDP to alpha, the nominal FDR, which gives an indication of the dispersion of its distribution relative to the nominal FDP. Used by the "Auto" specification. - **gma.Ntsts**: the effective denominator in the CLT asymptotic approximation to the distribution of the FDP, which equals the positive proportion, 'gamma', times the number of simultaneous tests, 'm'. - **detail**: The extractor function, `detail`, will return simulation replicates. See the linked documentation ## References Izmirlian G. (2020) Strong consistency and asymptotic normality for quantities related to the Benjamini-Hochberg false discovery rate procedure. Statistics and Probability Letters; 108713, Izmirlian G. (2017) Average Power and $\lambda$-power in Multiple Testing Scenarios when the Benjamini-Hochberg False Discovery Rate Procedure is Used. Jung S-H. (2005) Sample size for FDR-control in microarray data analysis. Bioinformatics; 21:3097-3104. Kluger D. M., Owen A. B. (2023) A central limit theorem for the Benjamini-Hochberg false discovery proportion under a factor model. Bernoulli; xx:xxx-xxx. Liu P. and Hwang J-T. G. (2007) Quick calculation for sample size while controlling false discovery rate with application to microarray analysis. Bioinformatics; 23:739-746. Lehmann E. L., Romano J. P.. Generalizations of the familywise error rate. Ann. Stat.. 2005;33(3):1138-1154. Romano Joseph P., Shaikh Azeem M.. Stepup procedures for control of generalizations of the familywise error rate. Ann. Stat.. 2006;34(4):1850-1873. ## Author(s) Grant Izmirlian ## See Also `pwrFDR.grid` `controlFDP` ## Examples ```r ## Example 1a: average power rslt.avgp <- pwrFDR(effect.size=0.79, n.sample=46, r.1=2000/54675, alpha=0.15) rslt.avgp ## Example 1b: average power, FDP.control.method set to "Auto", N.tests=1000 rslt.avgp.auto <- pwrFDR(effect.size = 0.79, n.sample = 46, r.1 = 2000/54675, alpha = 0.15, N.tests = 1000, FDP.control.method = "Auto") rslt.avgp.auto ## Example 1c: average power, FDP.control.method set to "Auto", N.tests=2000 rslt.avgp.auto <- update(rslt.avgp.auto, N.tests = 2000) rslt.avgp.auto ## Example 1d: tp-power rslt.lpwr <- pwrFDR(effect.size=0.79, n.sample=46, r.1=2000/54675, alpha=0.15, lambda=0.80, N.tests=54675) rslt.lpwr ## Example 1e: sample size required for given average power rslt.ss.avgp <- pwrFDR(effect.size=0.79, average.power=0.82, r.1=2000/54675, alpha=0.15) rslt.ss.avgp ## Example 1f: sample size required for given tp-power rslt.ss.lpwr <- pwrFDR(effect.size=0.79, TPX.power=0.82, lambda=0.80, r.1=2000/54675, alpha=0.15, N.tests=54675) rslt.ss.lpwr ## Example 1g: simulation rslt.sim <- update(rslt.avgp, method="sim", n.sim=500, N.tests=1000) rslt.sim ## Example 1h: simulation rslt.sim <- update(rslt.avgp, method="sim", FDP.control.method="both", n.sim=500, N.tests=1000) rslt.sim ## Example 2: methods for adding, subtracting, multiplying, dividing, exp, log, ## logit and inverse logit rslt.avgp - rslt.sim logit(rslt.avgp) ## etc ## Example 3: Compare the asymptotic distribution of T/M with kernel ## density estimate from simulated data pdf <- with(detail(rslt.sim)$reps, density(T/M1)) med <- with(detail(rslt.sim)$reps, median(T/M1)) avg <- rslt.sim$average.power sd <- rslt.sim$se.ToM rng.x <- range(pdf$x) rng.y <- range(c(pdf$y, dnorm(pdf$x, mean=avg, sd=sd))) plot(rng.x, rng.y, xlab="u", ylab="PDF for T/M", type="n") with(pdf, lines(x, y)) lines(rep(rslt.sim$average.power, 2), rng.y, lty=2) lines(pdf$x, dnorm(pdf$x, mean=avg, sd=sd), lty=3) ```

Ch01-pwrFDR function