Power Under Multiplicity Project (PUMP): Estimating Power, Minimum Detectable Effect Size, and Sample Size When Adjusting for Multiple Outcomes in Multi-level Experiments

As introduced above, multiple testing procedures adjust $p$-values to counteract the multiple testing problem.⁴⁴4Alternatively, MTPs can decrease the critical values for rejecting hypothesis tests. For ease of presentation, this paper focuses only on the approach of adjusting $p$-values. We next describe how using a MTP protects against false positives.

Considering multiple outcomes presents both challenges and opportunities. First, we discuss the impact of MTPs on individual power. The power of an individual hypothesis test is the probability of correctly rejecting a null hypothesis when the effect is at least a specified size. We refer to a setting in which the true impact is at least as large as the desired effect size as a “false null” hypothesis, while a “true null” is a setting in which the true impact is zero. In the proceeding explanations, when we refer to rejecting a null hypothesis or detecting an effect, we assume that we are detecting an effect of a certain pre-specified size. If $p$-values are adjusted upward, one is less likely to reject true nulls, which reduces the probability of Type I errors, or false positive findings. At the same time, MTPs increase the probability of a Type II error, or false negative findings, when the test fails to reject a false null. Individual power is $1-Pr(\text{Type II error})$, so MTPs have the tradeoff of reducing Type I errors but also reducing individual power.

Next, we consider the impact of multiple outcomes on other definitions of power. Applying a MTP reduces power according to all definitions of power relative to the case when no MTP is applied to adjust $p$-values. However, as discussed previously, having multiple outcomes also allows for a wider variety of definitions of success. Recall that $1-$minimal power is the probability of detecting an effect on at least one outcome. Typically, $1-$minimal power, even after applying a MTP, is higher than individual power for a hypothesis test on a single, pre-specified outcome. Depending on the study, other definitions of power, such as $1/2$ or $1/3$-minimal power, may or may not have higher power than the power of a single hypothesis test.

The MTPs that are the focus of this paper have three key features that affect statistical power: (1) whether the MTP is a familywise procedure or a false discovery rate procedure; (2) whether the MTP is single-step or stepwise; and (3) whether the MTP takes the correlation between test statistics into account. Below we explain each of these features of MTPs and provide discussion of the new parameter specifications they require when estimating power.

Familywise procedures “reframe Type I error as a rate across the entire set or “family” of multiple hypothesis tests. This rate is called the familywise error rate (FWER) [15]. The FWER is typically set to the same value as the probability of a Type I error for a single test, e.g., $\alpha$. MTPs that control the FWER at $5\%$ adjust $p$-values in a way that ensures that the probability of at least one Type I error across the entire set of hypothesis tests is no more than $5\%$. The MTPs introduced by Bonferroni [16, 17], [18], and [19] all control the FWER” [7].

MTPs that control false discovery rate (FDR) take an entirely different approach to the multiple testing problem. FDR, introduced by [20], is a less stringent criteria than FWER. It is the expected proportion of all rejected hypotheses that are erroneously rejected. As laid out in [7], the two-by-two representation in Table 1 is often found in articles on multiple hypothesis testing, and helps to illustrate the difference between FWER and FDR. Let $M$ be the total number of tests. Therefore, we have $M$ unobserved truths: whether or not each null hypothesis is true or false. We also have $M$ observed decisions: whether or not the null hypotheses were rejected, because the $p$-values were less than $\alpha$. In Table 1, $U$, $V$, $W$, and $X$ are four possible scenarios: the numbers of true or false hypotheses not rejected or rejected. $M_{0}$ and $M_{1}$ are the unobservable numbers of true null and false null hypotheses. $R$ is the number of null hypotheses that were rejected, and $M-R$ is the number of null hypotheses that were not rejected.

In Table 1, $V$ is the number of erroneously rejected null hypotheses, or the number of false positive findings. Therefore, the FWER is equivalent to $Pr(V>0)$, the probability of at least one false positive finding. As noted in [7], “recall that Type I error is inflated when testing for effects when no MTPs are applied. Consider the setting when all the outcomes are independent of each other. The Type I error is almost $10\%$ when testing effects on two independent outcomes and $23\%$ when testing effects on five independent outcomes. These Type I error rates both correspond to the FWER. The goal of MTPs that control the FWER is to bring these percentages back down to $5\%$.”

Also in Table 1, the FDR is equal to $E(V/R)$ but is defined to be $0$ when $R=0$, or when no hypotheses are rejected. “As is frequently noted in the literature (e.g., [21, 14], the FWER and FDR have different objectives. Control of the FWER protects researchers from spurious findings and so may be preferred when even a single false positive could lead to the wrong conclusion about the effectiveness of an intervention. On the other hand, the FDR is more lenient with false positives” [7]. Researchers may be willing to accept a few false positives, $V$, when the total number of rejected hypotheses, $R$, is large. Note that under the complete null hypothesis that all $M$ null hypotheses are true, the FDR is equal to the FWER. Referring back to Table 1, under the complete null we have $V=R$, so

However, if any effects truly exist, then FWER $\geq$ FDR. As a result of the difference in objective between FWER and FDR, in the case where there is at least one false null hypothesis (at least one true effect), a MTP that controls the FDR at $5\%$ will have a Type I error rate that is greater than $5\%$.

A side remark is that MTPs may provide either “weak control” or “strong control” of the error rate they target. A MTP “provides weak control of the FWER or the FDR at level $\alpha$ if the control can only be guaranteed when all null hypotheses are true, e.g. when the effects on all outcomes are zero. A MTP provides strong control of the FWER or FDR at level $\alpha$ if the control is guaranteed when some null hypotheses are true and some are false, e.g. when there may be effects on at least some outcomes. Of course, strong control is preferred” [7].⁵⁵5The single-step and step-down Westfall Young MTPs (which we discuss below) always provide at least weak control of the FWER. In order for these procedures to provide strong control of the FWER, they require the assumption of subset pivotality [22].The distribution of the unadjusted test statistics or $p$-values is said to have subset pivotality if for any subset of null hypotheses, the joint distribution of the test statistics or of the $p-$values for the subset is identical to the distribution under the complete null. A consequence of this assumption is that the permutation of test statistics or $p$-values can be done under the complete null hypothesis rather than under the unknown partial hypothesis [22].

An additional feature of a MTP that affects its statistical power is whether it is a “single-step” or “stepwise” procedure. “Single-step procedures adjust each $p$-value independently of the other $p$-values. For example, the Bonferroni MTP multiplies all raw $p$-values by $M$. Therefore, one $p$-value adjustment does not depend on other $p$-value adjustments, only on the number of tests” [7].

In contrast, “stepwise procedures first order raw $p$-values (or test statistics), and then adjust according to the order of the tests. The adjustments depend on the null hypotheses already rejected in previous steps. For example, the Holm MTP — the stepwise counterpart to the Bonferroni MTP — orders raw $p$-values from smallest to largest. The procedure then multiplies the smallest $p$-value by $M$, the second smallest $p$-value by $M-1$, and so on. The Holm MTP, like most other stepwise procedures, also enforces monotonicity: each adjusted $p$-value is greater than or equal to the previous adjusted $p$-value, and enforces that any $p$-values is not greater than one. Overall, stepwise MTPs allow for less adjustment than single-step MTPs in later steps, and therefore preserve more power (for outcomes in the later steps). The Bonferroni and Westfall-Young single-step procedure are single-step; the Holm and Benjamini-Hochberg MTPs and the Westfall-Young step-down procedure are stepwise” [7]. Note that stepwise procedures may be “step-down” or “step-up,” referring to whether a procedure begins with the smallest $p$-value, and thus the largest effect size (step-down) or the largest $p$-value (step-up)“.

Due to the dependencies of adjustments in stepwise MTPs, a new assumption must be considered when estimating power under multiplicity: the proportion of outcomes on which there are truly impacts, or, equivalently, the number of false null hypotheses. “Researchers may be inclined to assume that there will be effects on all outcomes, as hypotheses of effects probably drive the selection of outcomes in the first place.(…) However, if the researchers are incorrect, the probability of detecting the extant effects can be diminished, sometimes substantially” [7].

As noted in [7], it is important to point out that under the assumption that some effects are truly null, we must change our notion of power for $d-$minimal powers (e.g., $1-$minimal power, $1/3-$minimal power, etc.) and complete power. While individual power is defined based on the probability of correctly rejecting false nulls, the definition is loosened here and includes the probability of erroneous rejections of true nulls. For example, $1/3-$minimal power is defined as the probability of detecting effects on at least $1/3$ of the total outcomes $M$, regardless of the number of outcomes with true effects. That is, $1/3-$minimal power is not defined as the probability of detecting effects among the $M$ outcomes on which the effects truly exist. This reframing of power is necessary for power to be consistent. If $d-$minimal power were defined based on false nulls, then the value and interpretation would change depending on what assumption the researcher is making about the number of false nulls, which is an unknown quantity. For example, with $M=5$ outcomes, the probability of detecting at least one effect would be very different depending on if we assume all five outcomes are false nulls, or if we assume only two of them are false nulls. Complete power, which is the probability of detecting effects on all outcomes, has similar issues. We define complete power only in the context where all effects are assumed to be false nulls — if any outcomes are assumed to be true nulls, then complete power is undefined.

There is an additional technical note about the calculation of complete power.⁶⁶6Complete power has also been referred to as “conjunctive power” [23] and “all pairs power” [24]. To calculate complete power, we do not need to adjust the $p$-values, and can instead reject each individual test based on the unadjusted $p$-values. Complete power is the power of the omnibus test constructed by whether or not we reject all the null hypotheses. This test was originally introduced as the intersection-union test because the null hypothesis is expressed as a union and the alternative hypothesis is expressed as an intersection [25, 26]. [25] showed that if all the individual tests are level $\alpha$, the intersection-union test is also a level $\alpha$ test. To provide some intuition, we do not need to adjust $p$-values for complete power because it is a special case where we must reject all the hypothesis tests. Thus, there is no way for the omnibus test to be rejected by chance because of a favorable configuration [3]. For example, consider if we have four tests, with two false nulls and two true nulls. If we consider $3-$minimal power, we just need one of the two true nulls to be rejected by chance alone, and there are two ways for this to occur. For complete power, there is only one way for us to reject all of the nulls. The downside of an intersection-union test is that it is conservative: the FWER is generally less than $\alpha$. For example, if we have two independent tests with Type I error $\alpha$, then if both of are true nulls, the probability of a Type I error for the omnibus test (the probability of rejecting both null hypotheses) is $\alpha^{2}$ [27].

The final feature of a MTP that affects its statistical power is whether or not it takes into account the correlation between test statistics. “The Bonferroni and Holm procedures strongly control the FWER in all cases, even when the test statistics are correlated, but they adjust $p$-values more than is necessary in that case. Along with the Bonferroni and Holm MTPs, the Benjamin-Hochberg MTP also does not take correlations into account.⁷⁷7The Benjamini-Hochberg procedure was originally shown to control the FDR for independent test statistics. However, [28] showed that it also controls the FDR for true null hypotheses with “positive regression dependence.” This condition is satisfied for most applications in practice. In contrast, both of the Westfall-Young MTPs rely on the estimation of the joint distribution of test statistics when the complete null hypothesis (that there are not effects on any of the outcomes) is true. This joint distribution of the test statistics is estimated from the study’s data” [7]. For example, random permutations of the treatment indicator break the association between treatment status and outcome. Repeating these permutations a large number of times results in a distribution of test statistics under the complete null. “Because the actual data are used to generate this null distribution, correlations among the test statistics are captured. Then observed test statistics can be compared with the distribution of test statistics under the complete null hypothesis” [7].⁸⁸8Instead of using test statistics, the Westfall-Young MTPs can alternatively compare raw $p$-values with the estimated joint null distribution of $p$-values.

The correlation between test statistics is a parameter a researcher must specify in order to estimate power, MDES or sample size requirements when using a MTP. When fitting a separate regression model for the impact on each outcome, the $\binom{M}{2}$ correlations between test statistics are equal to the “pairwise correlations between the residuals in the $M$ impact models” [7]. Then, “if there are no covariates in the impact models or if the $R^{2}$’s of the covariates are equivalent in all impact models, then the correlations between the test statistics are equal to the correlations between the outcomes. However, having different $R^{2}$’s across the impact models reduces the correlations between the residuals and therefore between test statistics.⁹⁹9For example, one of the multiple outcomes may have a baseline covariate with a high $R^{2}$ while another may have a baseline covariate with a smaller $R^{2}$. Also, block dummies may explain more variation in some outcomes than in others. Models of outcomes that are highly correlated are more likely to have residuals that are highly correlated because baseline covariates will tend to have similar $R^{2}$’s. The gaps between the correlations between outcomes and the correlations between residuals — and therefore the test statistics — may be wider for moderately or weakly correlated outcomes. In any case, the upper bounds of correlations between the test statistics are the correlations between the outcomes” [7].

We take an innovative simulation-based approach to estimating power, as introduced in [7]. This approach is then also applied to estimate MDES and sample size. In order to estimate power for a single outcome, we can often use closed-form algebraic expressions, which are derived from the assumed model. However, with multiple outcomes, finding such expressions can be quite difficult, or even impossible. In cases where it is possible to find a closed-form expression, we would need to find expressions for every design and model, MTP, and definition of power. Importantly, we would also need to find new expressions for any possible number of outcomes, which quickly becomes an intractable problem. Furthermore, in some cases, such as permutation-based procedures like Westfall-Young approaches, a closed-form solution does not exist. To avoid these complexities, we rely on simulation to calculate estimated power. The approach outlined below can estimate power for any scenario.

If we were to rely on a full simulation approach, we could use the following method to estimate power. We introduce this full simulation approach to provide intuition, but use a simplified and far less computationally intensive approach in the package, as discussed below.

1. 1.

   Simulate a data sample according to the joint alternative hypothesis. First, we formulate what we will refer to as the joint alternative hypothesis, which is the set of outcomes we assume to have treatment effects above the desired size. We define $\psi_{m}$ to be the treatment effect for outcome $m$, with $M$ total outcomes. If we have $M=5$ outcomes, as in the Diplomas Now study, one possible joint alternative hypothesis is that all outcomes have effects above specified sizes: $H_{A}:\psi_{1}>0.125,\psi_{2}>0.2,\psi_{3}>0.1,\psi_{4}>0.1,\psi_{5}>0.05$. Another possible joint alternative hypothesis is one where only the first two outcomes have effects above the desired sizes: $H_{A}:\psi_{1}>0.125,\psi_{2}>0.2,\psi_{3}=\psi_{4}=\psi_{5}=0$. Once our joint alternative hypothesis is specified, we would generate simulated data under this hypothesis. To simulate data, we need to specify the full set of parameters as mentioned in Section 5.3 that allow for data generation. The Technical Appendix contains more details about the assumed data-generating process. For example, for the Diplomas Now experiment, we would assume a specific data generating process to allow us to simulate synthetic students, schools, and districts, including covariates, outcomes, and treatment assignment. This process would involve specifying parameter values such as $R^{2}$, the amount of outcome variation explained by covariates at a particular level, and translating these parameter choices into data-generating parameters, such as the coefficient values for covariates in a linear model.

2. 2.

   Estimate impacts on the simulated data. Given simulated data, we could fit $M$ regression models (specified to match the experimental design and model assumptions). For the models supported by PUMP, the relevant functions would be lm(), lmer() from the lme4 library [29], and interacted_linear_estimators() from the blkvar library.¹⁰¹⁰10This package is currently under development on GitHub; see https://github.com/lmiratrix/blkvar From the model output we extract the test statistics $t_{m}$ for the estimated impacts, one statistic for each outcome, along with estimated standard errors.

3. 3.

   Calculate unadjusted $p$-values. The test statistics and standard errors would in turn give raw (unadjusted) $p$-values. We can either calculate these by hand, or use the $p$-values routinely returned by regression functions. For Diplomas Now we could run a regression model of each attendance measure on treatment status and student and school covariates, and extract $p$-values from the regression outputs.

4. 4.

   Repeat above steps (1 through 3) for a large number of iterations. Denote the number of iterations tnum. Repeating steps $1$-$3$ tnum times results in a matrix of unadjusted $p$-values which we call $\mathbf{F}$, and is of dimension $tnum\times M$. One row corresponds to one set of simulated raw $p$-values from regressions for the $5$ attendance outcomes of interest for Diplomas Now.

5. 5.

   Adjust $p$-values. For each row, corresponding to one simulated dataset, the $M$ raw $p$-values corresponding to the $M$ hypothesis tests can be adjusted according to the desired multiple testing procedure. This process generates a new matrix $\mathbf{G}$ of adjusted $p$-values. For Bonferroni, Holm, and Benjamini-Hochberg adjustments, we use the function p.adjust in R (found in the stats package). We developed our own functions for implementing adjustment using the Westfall-Young procedures. One row corresponds to one set of simulated adjusted $p$-values for the $5$ attendance outcomes of interest for Diplomas Now.

6. 6.

   Calculate hypothesis rejection indicators. For any MTP, the matrix of adjusted $p$-values $\mathbf{G}$ can then be compared with a specified value of $\alpha$ (the default is $0.05$, but the value can be changed by the user). For each row, corresponding to one iteration of simulated data, we record whether or not the null hypothesis was rejected for each outcome. This process results in a new matrix $\mathbf{H}$, which contains hypothesis rejection indicators (still of dimension $tnum\times M$). Using $\mathbf{H}$, we can compute all definitions of power.

7. 7.

   Calculate power. To compute the different definitions of power:

• •

  Individual power for outcome $m$ is the proportion of the tnum rows in which the null hypothesis $m$ was rejected (the mean of column $m$ of $\mathbf{H}$). We would have $5$ different individual power values for Diplomas Now, corresponding to each outcome of interest.

• •

  $d$-minimal power is the proportion of the $tnum$ rows in which at least $d$ of the $M$ hypotheses were rejected.¹¹¹¹11Note that others refer to $1-$minimal power simply as “minimal power” [30, 3, 6], “disjunctive power” [23], or “any pair” power [24]. [3] use the terminology of “r-power” for what is referred to here as $d-$minimal power for $d>1$. For Diplomas Now, we could consider $1$-minimal power through $4$-minimal power.

• •

  Complete power is the proportion of the $tnum$ rows in which all of the null hypotheses were rejected based on the raw $p$-values rather than adjusted $p$-values (based on the matrix $\mathbf{G}$ rather than $\mathbf{H}$.) We would be interested in complete power if we want to evaluate whether Diplomas Now resulted in improvement for every single attendance outcome of interest. With $5$ outcomes, this criteria is a relatively strict indicator of success.

Above, we outlined a full simulation-based approach for calculating power. This approach would be computationally intensive because of the need to generate and analyze a full simulated dataset at each iteration. We can simplify this process by skipping the simulation of data and modeling steps. Given an assumed model and correlation structure for the test statistics, we can directly sample from $f(t_{1},\dots,t_{M})$, the joint alternative distribution of the test statistics. This shortcut vastly improves both the simplicity and the speed of computation. In summary, our approach is:

1. 1.

   Generate draws of test statistics $t_{1},\dots,t_{M}$ under the joint alternative hypothesis. This step produces a $tnum\times M$ matrix $\mathbf{E}$.

2. 2.

   Calculate unadjusted $p$-values. This produces the matrix $\mathbf{F}$, as in the procedure above.

3. 3.

   Adjust $p$-values. This produces the matrix $\mathbf{G}$, as in the procedure above.

4. 4.

   Calculate hypothesis rejection indicators. This produces the matrix $\mathbf{H}$, as in the procedure above.

5. 5.

   Calculate power.

We now describe how to sample from $f(t_{1},\dots,t_{M})$ directly. First, we assume a particular research design and model. In our example based on the Diplomas Now study, the research design is a $3$-level experiment, with randomization at level $2$. We plan for analyzing our data with a linear regression model with fixed intercepts at the district level, random intercepts at the school level, and a constant treatment effect across schools and districts. As previously, denote $\psi_{m}$ as the treatment effect for outcome $m$. We express treatment effects in terms of effect sizes:

where $\sigma_{m}$ is the standard deviation of outcome $Y_{m}$ in the control group. In order to calculate power, we also need the standard error of the impact in effect size units, which we denote as

The quantity $Q_{m}$ is a consequence of the assumed model, the number of units at different levels, the percent of units treated, the assumed $R^{2}$, and other parameters; our technical appendix shows formulae for $Q_{m}$ for all the designs and models our package supports. In our Diplomas Now example, $Q_{m}$ will be a function of the number of students, schools, and districts; the proportion of treated units; the number of student and school covariates; the explanatory power of the student and school covariates; the proportion of variation in the outcome explained by schools and districts; and the amount of impact variation relative to the amount of mean variation. Some parameters, such as the percent of units treated, will generally be known, while others, such as the $R^{2}$ at different levels, would need to be supplied by the user through either estimation on pilot data or assumptions based on prior knowledge.

Given the effect sizes $ES_{m}$ and the standard errors $Q_{m}$, we can determine the distribution of the vector of test statistics. When testing the hypothesis for outcome $m$, the test statistic for a $t$-test is:

with degrees of freedom $df$, also defined by the assumed model. Under the alternative hypothesis for outcome $m$, $t_{m}$ has a $t$ distribution with degrees of freedom $df$ and mean $ES_{m}/Q_{m}$. Finally, in addition to the parameters above, we also need to choose the correlation matrix between test statistics $\rho$ to sample from the joint distribution of $f(t_{1},\dots,t_{M})$. With these distributions specified, we can calculate $p$-values.

Note that this approach of simulating test statistics builds on work by [31], who use simulated test statistics to identify critical values based on the distribution of the maximum test statistics. Their approach produces the same estimates as the approach described here for the single-step Westfall-Young MTP. As an alternative to a simulation-based approach, [3] derived explicit formulae for $d-$minimal powers of stepwise procedures and for complete power of single-step procedures, but only for $1$, $2$, or $3$ tests. The approach presented here is more generally applicable, as it can be used for all MTPs, for any number of tests, and for all definitions of power discussed in the present paper.

Remark. The $p$-value adjustment using Westfall-Young procedures is the most complex correction procedure, so we briefly outline it here. Similar to above, we first explain a full simulation approach, and then discuss our simplification. Under a full simulation approach, we would first generate a single dataset under the joint alternative hypothesis and calculate a set of $M$ observed test statistics. Then, we would permute the single simulated dataset, say $B=3,000$ times, assuming the joint null hypothesis, and calculate test statistics on each of these permuted datasets. This process generates an empirical distribution of $B$ test statistics under the joint null distribution. Next, we compare the distribution of observed test statistics to the generated distribution of test statistics under the joint null distribution to calculate $p$-values. We would then re-generate a new simulated dataset, and repeat the process. If we were to generate $tnum=10,000$ datasets under the joint alternative hypothesis, for each of these datasets we generate $B=3,000$ permuted datasets under the joint null, so we would have to generate $10,000\times 3,000$ datasets!

When we skip the step of simulating data, then for each iteration $t$ in $1,\dots,tnum$ we first generate a set of $M$ observed test statistics from the joint alternative distribution. Then, we draw $B$ samples of test statistics under the joint null rather than permuting the data $B$ times. Under the null hypothesis, $t_{m}$ has a $t$ distribution with degrees of freedom $df$ and mean $0$. As before, we then compare the distribution of observed test statistics to the distribution of test statistics under the joint null distribution to calculate $p$-values. Westfall-Young procedures are computationally intensive, so the approach of skipping the simulated data step is particularly helpful here. This approach substantially reduces computational time by drawing test statistics directly rather than permutating the data.

Frequently, a researcher’s main concern with power is calculating either the MDES for each outcome in a given study, or determining the necessary sample size to achieve a target power given a specified set of MDES values. In Diplomas Now, we might want to know what sample sizes we would need to detect at least one significant effect across our outcomes if all the outcomes had a specified effect size (corresponding to $1-$minimal power) and we were planning on using the Holm procedure.

For pump_mdes() and pump_sample(), the user provides a particular target power, say $80\%$. The method then conducts a stochastic optimization problem to determine a value (of sample size or MDES) that is within a specified tolerance of the target power with high probability. We discuss the algorithm for MDES, although the approach for determining sample size is the same.

The algorithm first determines an initial range of MDES values that likely contain the target MDES. This initial range is calculated using formulae for unadjusted power based on the standard errors and degrees of freedom. In particular, from [8], in general the MDES for a single outcome can be estimated as

where $MT_{df}$, the “multiplier,” is the sum of two $t$ statistics with degrees of freedom $df$. For one-tailed tests, $MT_{df}=t_{\alpha}^{\star}+t_{1-\beta}^{\star}$ where $\alpha$ is the Type I error rate and $\beta$ is the desired power. For two-tailed tests, $MT_{df}=t_{\alpha/2}^{\star}+t_{1-\beta}^{\star}$. We do not explain the details of the derivations of the multiplier here; for more details and understanding, see [8] or [12]. These expressions can be further manipulated to obtain sample size formulae; see our technical appendix for all formulae used in the package.

We can calculate our initial bounds by manipulating the $\alpha$ and $\beta$ values in the above. First, to calculate the preliminary lower bound, we apply the formula above as given, assuming individual unadjusted power will give the smallest MDES; to calculate the preliminary upper bound, we apply the formula using $\alpha/M$ to correspond to a Bonferroni correction. We also adjust $\beta$ to account for different power types. For example, if we are interested in complete power, we need a larger upper bound than for individual power; in order to have a complete power of $80\%$, we would need each outcome to have an individual power of $\text{0.8}^{(1/M)}$, assuming independence. If we are interested in minimal power, we must have a smaller lower bound; in order to have $1-$minimal power of $80\%$, each outcome would only need to have individual power of $1-(1-\text{0.8})^{(1/M)}$. We ignore correlation in the setting of the initial bounds; the bounds do not need to be strict, given the adaptive nature of the subsequent search.

Once the initial range is established, we use pump_power() with the complete array of design parameters, including the correlation between test statistics, to obtain rough (using a small tnum, or number of simulation trials) estimates of power for five initial values across this range. We then fit a scaled logistic curve to these five points, and identify where the curve crosses the desired power level. After fitting an initial curve, we iterate, repeatedly calculating power for the targeted point and using the result to update the logistic curve model. At any point, if the current fitted curve’s range does not contain the target power, the algorithm extrapolates beyond the initial bounds for the next step. With each iteration we increase tnum to increase precision as we narrow in on the final answer; with each update to our estimated power curve, we weigh the collection of observations by their precisions (determined by corresponding tnum value). Once a test point achieves the target power to within tolerance, we conduct an additional simulation check using a high number of replicates to verify the proposed answer is within a specified tolerance of the target power; if it is not, we continue the iterative search. The default tolerance is $1\%$, so given a target power of $80\%$, we stop when we find a MDES that gives an estimated power between $79\%$ and $81\%$.

In practice, due to the monotonic nature of the logistic functional form, our algorithm generally converges fairly rapidly. However, in certain corner cases the algorithm may fail to converge on a value within tolerance. For more information on applying the search algorithm, see the sample size vignette.

We completed extensive validation checks to ensure our power calculation procedures are correct. First, we compared our power estimates in scenarios with only one outcome, $M=1$, to those from the PowerUpR package. Without a multiple testing procedure adjustment, our estimates match. Second, in order to validate our estimates under multiplicity, we followed the full simulation approach outlined above, in Section 4.1. The simulation approach involves generating many iterations of full datasets according to the assumed design and model, calculating $p$-values, and calculating an empirical estimate of power. Using a binomial distribution we constructed Monte Carlo confidence intervals for the power estimates from the full simulation approach. Then, we validated that the PUMP estimates fall within these confidence intervals.

A more detailed explanation of the validation procedure can be found in the Appendix, and full validation code and results are in a supplementary GitHub repository pump_validate. For some scenarios, we have some apparent discrepancies from PowerUp, but these result from different modeling choices. For example, for certain models PowerUp assumes the intraclass correlation is zero, while we allow for nonzero values. These discrepancies are noted in the appendix.

When planning a study, the researcher first has to identify the design of the experiment, including the number of levels, and the level at which randomization occurs. These decisions can be a mix of the realities of the context (e.g., the treatment must be applied at the school level, and students are naturally nested in schools, making for a cluster randomization), or deliberate (e.g., the researcher groups similar schools to block their experiment in an attempt to improve power). Second, based on the design and the inferential goals of the study, the researchers chooses an assumed model, including whether intercepts and treatment effects should be treated as constant, fixed, or random. For the same experimental design, the analyst can sometimes choose from a variety of possible models, and these two decisions should be kept conceptually separated from each other.

The design. The PUMP package supports designs with one, two, or three levels, with randomization occurring at any level. For example, a design with two levels and randomization at level one is a blocked design (or equivalently a multisite experiment), where level two forms the blocks (blocks being groups of units, some of which are treated and some not). Ideally, the blocks in a trial will be groups of relatively homogenous units, but frequently they are a consequence of the units being studied (e.g., evaluations of college supports, with students, the units, nested in colleges, the blocks). A design with two levels and randomization at level two is commonly called a cluster design (e.g., a collection of schools, with treatment applied to a subset of the schools, with outcomes at the student level); here the schools are the clusters, with a cluster being a collection of units which is entirely treated or entirely not. We can also have both blocking and clustering: randomizing schools within districts, creating a series of cluster-randomized experiments, would be a blocked (by district), cluster-randomized experiment, with randomization at level two.

The model. Given a design, the researcher can select a model via a few modeling choices. In particular the researcher has to decide, for each level beyond the first, about the intercepts and the treatment impacts:

• •

  Whether level two and level three intercepts are:

  • –

    fixed: we have a separate intercept for each unit.

  • –

    random: we have a separate intercept for each unit as above, but model the collection of intercepts as Normally distributed, allowing for partial pooling.

• •

  Whether level two and level three treatment effects are:

  • –

    constant: we model all units within a group as having the same single average impact.

  • –

    fixed: we allow each block or cluster within a level to have its own individual estimated impact (we can only do this if we have treated and control units within said block or cluster).

  • –

    random: we allow variation as with fixed, but model the collection of treatment impacts as Normally distributed around a grand mean mean impact. This is implicitly allowing for the sample as being representative of a larger super-population, in terms of treatment impact estimation.

We denote the research design by $d$, followed by the number of levels and randomization level, so d3.1 is a three level design with randomization at level one. The model is denoted by $m$, followed by the level and the assumption for the intercepts, either $f$ or $r$ and then the assumption for the treatment impacts, $c$, $f$, or $r$. For example, m3ff2rc means at level $3$, we assume fixed intercepts and fixed treatment impacts, and at level two we assume random intercepts and constant treatment impacts. The full design and model are specified by concatenating these together, e.g. d2.1_m3fc. The Diplomas Now model, for example, is d3.2_m3fc2rc.

The full list of supported design and model combinations is below. The user can see the list by calling pump_info(), which provides the designs and models, MTPs, power definitions, and model parameters. We also include the corresponding names from the PowerUP! package where appropriate. For more details about each combination of design and model, see the Technical Appendix.

The supported multiple testing procedures were covered in more detail in Section 3. Here we provide a review of the multiple testing procedures supported by the PUMP package:

• •

  Bonferroni: adjusts $p$-values by multiplying them by $M$ to ensure strong control of the FWER. Bonferroni is a simple procedure, but the most conservative.

• •

  Holm: a step-down version of Bonferroni. Starting from smallest to largest, $p$-values are sequentially adjusted by different multipliers. Holm is less conservative than Bonferroni for larger $p$-values.

• •

  Benjamini-Hochberg: A sequential, step-up procedure that controls the FDR. Using the BH method, only null hypotheses with $p$-values below a certain threshold are rejected, where the threshold is determined by the number of tests and the level $\alpha$.

• •

  Single-step Westfall-Young: A permutation-based procedure for controlling the FWER, which directly takes into account the joint correlation structure of the outcomes. In the single-step approach, all outcomes are adjusted by using the permuted distribution of the minimum $p$-value. Although Westfall-Young procedures are less conservative while still protecting against false discoveries, they are computationally very intensive.

• •

  Step-down Westfall-Young: A similar approach to the single-step procedure, except that outcomes are adjusted sequentially from smallest to largest according to the permuted distributions of the corresponding sequential $p$-values.

For a more detailed explanation of each MTP, see Appendix A of [7].

The following table from [7] summarizes the important features for each of the MTPs supported by PUMP.

The table below shows the parameters that influence $Q_{m}$, the standard error, for different designs and models.

A few parameters warrant more explanation.

• •

  The quantity ICC is the unconditional Intraclass Correlation, and gives a measure of variation at different levels of the model. For each outcome, the ICC for each level is defined as the ratio of the variance at that level divided by the overall variance of the individual outcomes. The ICC includes the variation due to covariates.

• •

  For each outcome, the quantity omega ($\omega$) for each level is the ratio between impact variation at that level and variation in intercepts (including covariates) at that level. It is a measure of treatment impact heterogeneity.

• •

  The $R^{2}$ expressions are the percent of variation at a particular level predicted by covariates specific to that level. For simplicity we assume covariates at a level are group mean centered, so only covariates at a particular level explain variance at that level.

For precise formulae of these expressions, see the Technical Appendix, which outlines the assumed data-generating process, and the resulting expressions for ICC, $\omega$, and $R^{2}$.

In addition to design parameters, there are additional parameters that control the precision of the power estimates themselves:

• •

  tnum is the number of test statistics generated in order to estimate power. A larger number of test statistics results in greater computation time, but also a more precise estimate of power. Note that the pump_mdes() and pump_sample() have multiple tnum parameters controlling the precision of the search.

• •

  B is the number of Westfall-Young permutations. Again, there is a tradeoff between precision and computation time.

• •

  parallel.WY.cores specifies the number of cores to use for parallel computation of the Westfall-Young Step-Down procedure, which is the most computationally intensive. The default of 1 does not result in parallel computation. Parallelization is done using parApply from the parallel package.

To conduct power, MDES, and sample size calculations, we first specify the design, sample sizes, analytic model, and level of statistical significance. We also must specify parameters of the data generating distribution that match the selected design and model. All of these numbers have to be determined given resource limitations, or estimated using prior knowledge, pilot studies, or other sources of information.

We next discuss selection of all needed design parameters and modeling choices. For further discussion of selecting these parameters see, for example [12] and [8]. For discussion in the multiple testing context, especially with regards to the overall power measures such as $1-$minimal or complete power, see [7]; the findings there are general, as they are a function of the final distribution of test statistics. The key insight is that power is a function of only a few summarizing elements: the individual-level standard errors, the degrees of freedom, and the correlation structure of the test statistics. Once we have these elements, regardless of the design, we can proceed.

Analytic model. We first need to specify how we will analyze our data; this choice also determines which design parameters we will need to specify. Following the original Diplomas Now report, we plan on using a multi-level model with fixed effects at level three, a random intercept at level two, and a single treatment coefficient. We represent this model as “m3fc2rc.” The “3fc” means we are including block fixed effects, and not modeling any treatment impact variation at level three. The “2rc” means random intercept and no modeled variation of treatment within each block (the “c” is for “constant”). We note that the Diplomas Now report authors call their model a “two-level” model, but this is not quite aligned with the language of this package. In particular, fixed effects included at level two are actually accounting for variation at level three; we therefore identify their model as a three level model with fixed effects at level three.

Sample sizes. We assume equal size randomization blocks and schools, as is typical of most power analysis packages. For our context, this gives about three schools per randomization block; we can later do a sensitivity check where we increase and decrease this to see how power changes. The Diplomas Now report states there were 14,950 students, yielding around 258 students per school. Normally we would use the geometric means of schools per randomization block and students per school as our design parameters, but that information is not available in the report. We assume 50% of the schools are treated; our calculations will be approximate here in that we could not actually treat exactly 50% in small and odd-sized blocks.

Control variables. We next need values for the $R^{2}$ of the possible covariates. The report does not provide these quantities, but it does mention covariate adjustment in the presentation of the model. Given the types of outcomes we are working with, it is unlikely that there are highly predictive individual-level covariates, but our prior year school-average attendance measures are likely to be highly predictive of corresponding school-average outcomes. We thus set $R^{2}_{1}=0.1$ and $R^{2}_{2}=0.5$. We assume five covariates at level one and three at level two; this decision, especially for level one, usually does not matter much in practice, unless sample sizes are very small (the number of covariates along with sample size determine the degrees of freedom for our planned tests).

ICCs. We also need a measure of where variation occurs: the individual, the school, or the randomization block level. We capture this with Intraclass Correlation Coefficients (ICCs), one for level two and one for level three. ICC measures specify overall variation in outcome across levels: e.g., do we see relatively homogeneous students within schools that are quite different, or are the schools generally the same with substantial variation within them? We typically would obtain ICCs from pilot data or external reports on similar data. We here specify a level-two ICC of 0.05, and a level-three ICC of 0.40. We set a relatively high level three ICC as we expect our school type by district blocks to isolate variation; in particular we might believe middle and high school attendance rates would be markedly different.

Correlation of outcomes. We finally need to specify the number and relationship among our outcomes and associated test-statistics. For illustration, we select attendance as our outcome group. We assume we have five different attendance measures. The main decision regarding outcomes is the correlation of our test statistics. As a rough proxy, we use the correlation of the outcomes at the level of randomization; in our case this would be the correlation of school-average attendance within block. We believe the attendance measures would be fairly related, so we select rho = 0.40 for all pairs of outcomes. This value is an estimate, and we strongly encourage exploration of different values of this correlation choice as a sensitivity check for any conducted analysis. Selecting a candidate rho is difficult, and will be new for those only familiar with power analyses of single outcomes; we need to more research in the field, both empirical and theoretical, to further guide this choice.

If the information were available, we could specify different values for the design parameters such as the $R^{2}$s and $ICC$s for each outcome, if we thought they had different characteristics; for simplicity we do not do this here. The PUMP package also allows specifying different pairwise correlations between the test statistics of the different outcomes via a matrix of $\rho$s rather than a single $\rho$; also for simplicity, we do not do that here.

Once we have established initial values for all needed parameters, we first conduct a baseline calculation, and then explore how MDES or other quantities change as these parameters change.

We now have an initial planned design, with a set number of schools and students. But is this a large enough experiment to reliably detect reasonably sized effects? To answer this question we calculate the minimal detectable effect size (MDES), given our planned analytic strategy, for our outcomes.

To identify the MDES of a given setting we use the pump_mdes method, which conducts a search for a MDES that achieves a target level of power. The MDES depends on all the design and model parameters discussed above, but also depends on the type of power and target level of power we are interested in. For example, we could determine what size effect we can reliably detect on our first outcome, after multiplicity adjustment. Or, we could determine what size effects we would need across our five outcomes to reliably detect an impact on at least one of them. We set our goal by specifying the type (power.definition) and desired power (target.power).

Here, for example, we find the MDES if we want an 80% chance of detecting an impact on our first outcome when using the Holm procedure:

The results are easily made into a nice table via the knitr kable() command:

The answers pump_mdes() gives are approximate as we are calculating them via monte carlo simulation. To control accuracy, we can specify a tolerance (tol) of how close the estimated power needs to be to the desired target along with the number of iterations in the search sequence (via start.tnum, tnum, and final.tnum). The search will stop when the estimated power is within tol of target.power, as estimated by final.tnum iterations. Lower tolerance and higher tnum values will give more exact results (and take more computational time).

Changing the type of power is straightforward: for example, to identify the MDES for $1-$minimal power (i.e., what effect do we have to assume across all observations such that we will find at least one significant result with 80% power?), we simply update our result with our new power definition:

#> mdes result: d3.2_m3fc2rc d_m with 5 outcomes
#>   target min1 power: 0.80
#>  MTP Adjusted.MDES min1.power   SE
#>   HO    0.08144353    0.79075 0.01
#>  (10 steps in search)

The update() method can replace any number of arguments of the prior call with new ones, making exploration of different scenarios very straightforward.¹²¹²12The update() method re-runs the underlying call of pump_mdes(), pump_sample(), or pump_power() with the revised set of design parameters. You can even change which call to use via the type parameter. Our results show that if we just want to detect at least one outcome with 80% power, we can reliably detect an effect of size $0.08$ (assuming all five outcomes have effects of at least that size).

When estimating power for multiple outcomes, it is important to consider cases where some of the outcomes in fact have null, or very small, effects, to hedge against circumstances such as one of the outcomes not being well measured. One way to do this is to set two of our outcomes to no effect with the numZero parameter:

#> mdes result: d3.2_m3fc2rc d_m with 5 outcomes
#>   target min1 power: 0.80
#>  MTP Adjusted.MDES min1.power SE
#>   HO    0.09050718     0.8065  0
#>  (7 steps in search)

The MDES goes up, as expected: when there are not effects on some outcomes, there are fewer good chances for detecting an effect. Therefore, an increased MDES (for the nonzero outcomes) is required to achieve the same level of desired power (80%). Below we provide a deeper dive into the extent to which numZero can effect power estimates.

The MDES calculator tells us what we can detect given a specific design. We might instead want to ask how much larger our design would need to be in order to achieve a desired MDES. In particular, we might want to determine the needed number of students per school, the number of schools, or the number of blocks to detect an effect of a given size. The pump_sample method will search over any one of these.

Assuming we have three schools per block, we first calculate how many blocks we would need to achieve a MDES of 0.10 for $1-$minimal power (this answers the question of how big of an experiment do we need in order to have an 80% chance of finding at least one outcome significant, if all outcomes had a true effect size of 0.10):

#> sample result: d3.2_m3fc2rc d_m with 5 outcomes
#>   target min1 power: 0.80
#>  MTP Sample.type Sample.size min1.power   SE
#>   HO           K          15      0.798 0.01
#>  (4 steps in search)

We would need 15 blocks, rather than the originally specified 21, giving 45 total schools in our study, to achieve 80% $1-$minimal power.

We recommend checking MDES and sample-size calculators, as the estimation error combined with the stochastic search can give results a bit off the target in some cases. A check is easy to do; simply run the found design through pump_power(), which directly calculates power for a given scenario, to see if we recover our originally targeted power (we can use update() and set the type to power to pass all the design parameters automatically). When we do this, we can also increase the number of iterations to get more precise estimates of power, as well:

When calculating power directly, we get power for all the implemented definitions of power applicable to the design.

In the above, the first five rows are the powers for rejecting each of the five outcomes—they are (up to simulation error) the same since we are assuming the same MDES and other design parameters for each. The “mean individual” is the mean individual power across all outcomes. The first column is power without adjustment, and the second has our power with the listed $p$-value adjustment.

The next rows show different multi-outcome definitions of power. In particular, 1-minimum shows the chance of rejecting at least one hypotheses. The complete row shows the power to reject all hypotheses; it is only defined if all outcomes are specified to have a non-zero effect.¹³¹³13The package does not show power for these without adjustment for multiple testing, as that power would be grossly inflated and meaningless.

We can look at a power curve of our pump_sample() call to assess how sensitive power is to our level two sample size:¹⁴¹⁴14The points on the plots show the evaluated simulation trials, with larger points corresponding to more iterations and greater precision.

Though increasing tnum is useful for checking the power calculation, it also increases computation time. Thus, for future calculations we save a call with the default tnum to reduce computation time.

Remark. In certain settings, a wide range of sample sizes may result in very similar levels of power. In this case, the algorithm may return a sample size that is larger than necessary. This pattern does not occur for the sample size at the highest level of the hierarchy, and only for occurs for sample sizes at lower levels of the hierarchy; e.g. for nbar for all models, and for nbar and J for three level models. In addition, due to the nature of the search algorithm, occasionally the algorithm may not converge. For a more detailed discussion of these challenges, see the package sample size vignette.

It is easy to rerun the above using the Westfall-Young Stepdown procedure (this procedure is much more computationally intensive to run), or other procedures of interest. Alternatively, simply provide a list of procedures you wish to compare. If you provide a list, the package will re-run the power calculator for each item on the list; this can make the overall call computationally intensive. Here we obtain power for our scenario using Bonferroni, Holm and Westfall-Young adjustments, and plot the results using the default plot() method:

To speed up computation, we set an additional option, parallel.WY.cores, to the number of desired cores to parallelize the computation. We also reduce tnum to decrease computation time.

The more sophisticated (and less conservative) adjustment exploits the correlation in our outcomes (rho = 0.4) to provide higher individual power. Note, however, that we do not see elevated rates for $1-$minimal power. Accounting for the correlation of the test statistics when adjusting $p$-values can drive some power (individual power) up, but on the flip side $1-$power can be driven down as the lack of independence between tests gives fewer chances for a significant result. See [7] for further discussion; while the paper focuses on the multisite randomized trial context, the lessons learned there apply to all designs as the only substantive differences between different design and modeling choices is in how we calculate the unadjusted distribution of their test statistics.

Within the pump package we have two general ways of exploring design sensitivity. The first is with update(), which allows for quickly generating a single alternate scenario. To explore sensitivity to different design parameters more systematically, use the grid() functions, which calculate power, mdes, and sample size for all combinations of a set of passed parameter values. There are two main differences between the two approaches. First, update() allows for different values of a parameter for the different outcomes; the grid approach, by contrast, is more limited in this regard, and assumes the same parameter value across different outcomes. Second, the grid functions are a powerful tool for systematically exploring many possible combinations, while update() only allows the user to explore one value at a time.

We first illustrate the update() approach, and then turn to illustrating grid() across three common areas of exploration: Intraclass Correlation Coefficients (ICCs), the correlation of test statistics, and the assumed number of non-zero effects. The last two are particularly important for multiple outcome contexts.