Package illustration
In this section, we illustrate how to use the PUMP package, using our example motivated by the Diplomas Now study. Given the study's design, we ask a natural initial question: What size of impact could we reasonably detect after using a MTP to adjust \(p\)-values to account for our multiple outcomes?
We mimic the planning process one might use for planning a study similar to Diplomas Now (e.g., if we were planning a replication trial in a slightly different context). To answer this question we therefore first have to decide on our experimental design and modeling approach. We also have to determine values for the associated design parameters that accompany these choices. In the following sections we walk through selecting these parameters (sample size, control variables, intraclass correlation coefficients, impact variation, and correlation of outcomes). We calculate MDES for the resulting context and determine how necessary sample sizes change depending on what kind of power we desire. We finally illustrate some sensitivity checks, looking at how MDES changes as a function of rho, the correlation of the test statistics.
Establishing needed design parameters
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 H. S. Bloom (2006) and Dong and Maynard (2013). For discussion in the multiple testing context, especially with regards to the overall power measures such as \(1-\)minimal or complete power, see Porter (2018); 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.
Calculating MDES
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:
m <- pump_mdes(
d_m = "d3.2_m3fc2rc", # choice of design and analysis strategy
MTP = "HO", # multiple testing procedure
target.power = 0.80, # desired power
power.definition = "D1indiv", # power type
M = 5, # number of outcomes
J = 3, # number of schools per block
K = 21, # number districts
nbar = 258, # average number of students per school
Tbar = 0.50, # prop treated
alpha = 0.05, # significance level
numCovar.1 = 5, # number of covariates at level 1
numCovar.2 = 3, # number of covariates at level 2
R2.1 = 0.1, R2.2 = 0.7, # explanatory power of covariates for each level
ICC.2 = 0.05, ICC.3 = 0.4, # intraclass correlation coefficients
rho = 0.4 ) # how correlated outcomes areThe results are easily made into a nice table via the knitr kable() command:
knitr::kable( m, digits = 3 ) %>%
kableExtra::kable_styling( position = "center" )| MTP | Adjusted.MDES | D1indiv.power |
|---|---|---|
| HO | 0.106 | 0.797 |
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:
m2 <- update( m, power.definition = "min1" )#> 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.2 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:
m3 <- update( m2, numZero = 2 )#> 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.
Determining necessary sample size
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):
smp <- pump_sample(
d_m = "d3.2_m3fc2rc",
MTP = "HO",
typesample = "K",
target.power = 0.80, power.definition = "min1", tol = 0.01,
MDES = 0.10, M = 5, nbar = 258, J = 3,
Tbar = 0.50, alpha = 0.05, numCovar.1 = 5, numCovar.2 = 3,
R2.1 = 0.1, R2.2 = 0.7, ICC.2 = 0.05, ICC.3 = 0.40, rho = 0.4 )
print( smp )#> 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:
p_check <- update( smp, type = "power", tnum = 20000,
long.table = TRUE )| power | None | HO |
|---|---|---|
| individual outcome 1 | 0.7 | 0.52 |
| individual outcome 2 | 0.7 | 0.52 |
| individual outcome 3 | 0.7 | 0.53 |
| individual outcome 4 | 0.7 | 0.52 |
| individual outcome 5 | 0.7 | 0.53 |
| mean individual | 0.7 | 0.53 |
| 1-minimum | 0.81 | |
| 2-minimum | 0.65 | |
| 3-minimum | 0.51 | |
| 4-minimum | 0.39 | |
| complete | 0.32 |
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.3
We can look at a power curve of our pump_sample() call to assess how sensitive power is to our level two sample size:4
plot( smp )
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.
pow <- update( p_check, tnum = 10000 )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.
Comparing adjustment procedures
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:
p2 <- update( pow,
MTP = c( "BF", "HO", "WY-SD" ) )
plot( p2 )
To speed up computation, we could also 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 Porter (2018) 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.
Exploring sensitivity to design parameters
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.
Exploring power with update()
Update allows for a quick change of some of the set of parameters used in a prior call; we saw update() used several times above. As a further example, here we examine what happens if the ICCs are more equally split across levels two and three:
p_b <- update( pow, ICC.2 = 0.20, ICC.3 = 0.25 )
print( p_b )#> power result: d3.2_m3fc2rc d_m with 5 outcomes
#> MTP D1indiv D2indiv D3indiv D4indiv D5indiv indiv.mean min1 min2 min3
#> None 0.2381 0.2435 0.2408 0.2463 0.2437 0.24248 NA NA NA
#> HO 0.0943 0.0963 0.0945 0.0995 0.0943 0.09578 0.2671 0.1119 0.0552
#> min4 complete
#> NA NA
#> 0.0285 0.0241
#> 0.000 <= SE <= 0.002We immediately see that our assumption of substantial variation in level three matters a great deal for power.
When calculating power for a given scenario, it is also easy to vary many of our design parameters by outcome. For example, if we thought we had better predictive covariates for our second outcome, we might try:
p_d <- update( pow,
R2.1 = c( 0.1, 0.3, 0.1, 0.2, 0.2 ),
R2.2 = c( 0.4, 0.8, 0.3, 0.2, 0.2 ) )
print( p_d )#> power result: d3.2_m3fc2rc d_m with 5 outcomes
#> MTP D1indiv D2indiv D3indiv D4indiv D5indiv indiv.mean min1 min2 min3
#> None 0.4345 0.8525 0.3912 0.3416 0.3380 0.47156 NA NA NA
#> HO 0.2472 0.6542 0.2147 0.1902 0.1884 0.29894 0.7139 0.3776 0.2145
#> min4 complete
#> NA NA
#> 0.1201 0.0874
#> 0.001 <= SE <= 0.002Notice how the individual powers are heavily impacted. The \(d\)-minimal powers naturally take the varying outcomes into account as we are calculating a joint distribution of test statistics that will have the correct marginal distributions based on these different design parameter values.
After several update()s, we may lose track of where we are; to find out, we can always check details with print_context() or summary():
summary( p_d )
#> power result: d3.2_m3fc2rc d_m with 5 outcomes
#> MDES vector: 0.1, 0.1, 0.1, 0.1, 0.1
#> nbar: 258 J: 3 K: 15 Tbar: 0.5
#> alpha: 0.05
#> Level:
#> 1: R2: 0.1 / 0.3 / 0.1 / 0.2 / 0.2 (5 covariates)
#> 2: R2: 0.4 / 0.8 / 0.3 / 0.2 / 0.2 (3 covariates) ICC: 0.05 omega: 0
#> 3: fixed effects ICC: 0.4 omega: 0
#> rho = 0.4
#> MTP D1indiv D2indiv D3indiv D4indiv D5indiv indiv.mean min1 min2 min3
#> None 0.4345 0.8525 0.3912 0.3416 0.3380 0.47156 NA NA NA
#> HO 0.2472 0.6542 0.2147 0.1902 0.1884 0.29894 0.7139 0.3776 0.2145
#> min4 complete
#> NA NA
#> 0.1201 0.0874
#> 0.001 <= SE <= 0.002
#> (tnum = 10000)Using update allows for targeted comparison of major choices, but if we are interested in how power changes across a range of options, we can do this more systematically with the grid() functions, as we do next.
Exploring the impact of the ICC
We above saw that the ICC does impact power considerably. We next extend this evaluation by exploring a range of options for both level two and three ICCs, so we can assess whether our power is sufficient across a set of plausible values. The update_grid() call makes this straightforward: we pass our baseline scenario along with lists of parameters to additionally explore. To decrease computation time and achieve less precise estimates in this exploratory phase, we decrease tnum. We can then easily visualize the variation in min1 power by calling plot() on the object.
grid <- update_grid( pow,
ICC.2 = seq( 0, 0.3, 0.05 ),
ICC.3 = seq( 0, 0.60, 0.20 ) )
plot( grid, power.definition = "min1" )
Our plots show the impact of one varying factor averaged over the other factors, like you would get with a plot for a main effect in a fully interacted experiment. Note that in addition to update_grid(), there are also base functions pump_power_grid(), pump_mdes_grid(), and pump_sample_grid().
We see that higher ICC.2 radically reduces power to detect anything and ICC.3 does little. To understand why, we turn to our standard error formula for this design and model: \[ \begin{aligned} SE( \hat{\tau} ) = \sqrt{ \frac{\text{ICC}_{2}(1 - R^2_{2})}{\bar{T}(1 - \bar{T}) JK} + \frac{(1-\text{ICC}_{2} - \text{ICC}_{3})(1-R^2_{1})}{\bar{T}(1 - \bar{T}) J K\bar{n}} } . \end{aligned} \] In the above, the \(\bar{n} = 258\) students per group makes the second term very small compared to the first, regardless of the ICC.3 value. The first term, however, is a direct scaling of ICC.2; changing it will change the standard error, and therefore power, a lot. All provided designs and models implemented in the package are discussed, along with corresponding formula such as these, in our technical supplement accompanying this paper and package.
For grid searches we recommend reducing the number of permutations, via tnum, to speed up computation. As tnum shrinks, we will get increasingly rough estimates of power, but even these rough estimates can help us determine trends.
The grid() functions provide easy and direct ways of exploring how power changes as a function of the design parameters. We note, however, that in order to keep syntax simple, they do not allow different design parameters, including MDES, by outcome. This is to keep package syntax simpler. When faced with contexts where it is believed that these parameters do vary, we recommend using average values for the broader searches, and then double-checking a small set of potential final designs with the pump_power() method.
Exploring the impact of rho
The correlation of test statistics, \(\rho\), is a critical parameter for how power will play out across the multiple tests. For example, with Westfall-Young, we saw that the correlation can improve our individual power, as compared to Bonferroni. We might not know what will happen to \(2-\)minimal power, however: on one hand, correlated statistics make individual adjustment less severe, and on the other correlation means we succeed or fail all together. We can explore this question relatively easily by letting rho vary as so:
gridRho <- update_grid( pow,
MTP = c( "BF", "WY-SD" ),
rho = seq( 0, 0.9, by = 0.15 ),
tnum = 1000,
B = 3000 )We then plot our results.
plot( gridRho )
First, we see the benefit of the Westfall-Young single-step procedure is minimal, as compared to Bonferroni. Second, the impact on individual adjustment is flat, as anticipated. Third, across a very broad range of rho, we maintain good \(1-\)minimal power. Complete power climbs as correlation increases, and \(2-\)minimal power is generally unchanged.
Exploring the impact of null outcomes
We finally explore varying the number of outcomes with no effects. This exploration is an important way to hedge a design against the possibility that some number of the identified outcomes are measured poorly, or are simply not impacted by treatment. We use a grid search, varying the number of outcomes that have no treatment impact via the numZero design parameter:
gridZero <- update_grid( pow,
numZero = 0:4,
M = 5 )
plot( gridZero, nrow = 1 )
There are other ways of exploring the impact of weak or null effects on some outcomes. In particular, the pump_power() and pump_sample() methods allow the researcher to provide an MDES vector with different values for each outcome, including 0s for some outcomes. The grid() functions, by contrast, take a single MDES value for the non-null outcomes, with a separate specification of how many of the outcomes are 0. (This single value plus numZero parameter also works with pump_power() if desired.)