Chapter 7 Code call-out 7.1: Power curves

In this code call-out we will explore calculations on power, and in particular will focus on considerations of the minimum detectable effect and power curves. This will be based on the study of Alan and Kubilay (2025) who conduct a clustered randomisation, and study a behavioural intervention in schools and its impact on a range of educational, psycho-social, and well-being outcomes.

In particular, the authors document their power calculations in their pre-registration, which can be found on the AEA Social Science Repository. As they are considering a clustered randomisation, we need to know various elements to determine power including the sample size in treated and control clusters, the number of clusters themselves, and expectations about the standard deviation of outcomes of interest as well as the intra-cluster correlation of these outcomes.

The key details which the authors use to calculate power are laid out in their pre-analysis plan, and reproduced below for ease of access. Of note, they propose to randomly assign 32 clusters to receive treatment and 33 to act as controls, with 350 observations per cluster. Additionally, the note the assumed standard deviations of of each variable as SD, and the intra-cluster correlation as ICC.

Table 1: Power and Minimum Detectable Effects with Clustering

Variables	Clusters in Treatment	Clusters in Control	Cluster Size	N	Mean	MDE	SD	ICC	Percent Change
Turkish Score	32	33	350	22750	4.294	0.571	2.456	0.109	0.133
Math Score	32	33	350	22750	3.914	0.535	2.372	0.103	0.137
Bullying in Class	32	33	350	22750	0.470	0.047	0.499	0.016	0.100
Bullying in School	32	33	350	22750	0.490	0.046	0.500	0.015	0.094
Sensitivity to World Issues	32	33	350	22750	3.413	0.095	0.881	0.021	0.028
Locus of Control	32	33	350	22750	3.200	0.073	0.542	0.035	0.023
Impulsivity	32	33	350	22750	2.116	0.035	0.570	0.005	0.017
Perspective Taking	32	33	350	22750	3.166	0.105	0.659	0.049	0.033
Mental Wellbeing	32	33	350	22750	3.040	0.049	0.524	0.015	0.016
Belonging	32	33	350	22750	2.858	0.068	0.660	0.019	0.024
Autonomy	32	33	350	22750	2.593	0.052	0.581	0.014	0.020
In-Degree Friendship	32	33	350	22750	2.767	0.373	2.613	0.039	0.135

The details above can be used to calculate minimum detectable effect sizes (MDE in the table), as well as more generally used to determine the power to detect specific effects. To see a specific calculation, remember that as laid in (7.14) and related discussion in the book, to calculating power in a clustered randomisation we can apply the following: \[ (1-\beta)=\Phi\left(\frac{\tau}{\sqrt{\sigma^2\left[\frac{1+(N_{c0}-1)\rho}{N_{c0}c_{0}}+\frac{1+(N_{c1}-1)\rho}{N_{c1}c_{1}}\right]}}-\Phi^{-1}(1-\alpha/2)\right), \tag{1}\] where \(N_{c0}\) and \(N_{c1}\) refer to the number of units per control and treated cluster respectively, \(c_0\) and \(c_1\) refer to the number of control and treated clusters, and \(\rho\) refers to the intra-cluster correlation.

Let’s confirm that this all makes sense, by conindering a specific row from the table above. Below we can confirm that based on their assumptions, the authors do indeed have 80% power to detect the reported MDE in row 1 (0.571). We do this below, seeing that indeed, the reported effect of 0.571 can be detected with the stated level of power.

s1 <- (1 + (350 - 1) * 0.109) / (32 * 350)
s0 <- (1 + (350 - 1) * 0.109) / (33 * 350)

# Combined standard error
SE <- sqrt(2.456^2 * (s0 + s1))

power <- pnorm(0.571 / SE - qnorm(0.975))
power

[1] 0.8011777

It is also easy enough for us to see how we can arrive at this minimum detectable effect directly. In this case, rather than having an effect of interest \(\tau\) in (Equation 1) and wishing to find the power, we wish to find the minimum detectable effect is power is set to 0.8. It is easy enough to see that if we apply an inverse normal transformation to each side of (Equation 1) and re-arrange, we can calculate our MDE (\(\tau\)) as: \[ \tau = \left[\Phi^{-1}\left(1-\beta\right) + \Phi^{-1}\left(1-\alpha/2\right)\right]\times SE, \] where SE is simply the standard error term in the denominator of (Equation 1). We do this below, seeing how we can arrive to the value of 0.571 (approximately) below:

threshold <- (qnorm(0.8) + qnorm(0.975)) * SE

threshold

[1] 0.5701424

More generally, instead of working with a specific effect to achieve a desired power, we could generate power curves to plot the power to detect an effect if the true effect were indeed some specific value. This in essence maps the formula in (Equation 1) to a series of effects \(\tau\), plotting the resulting power. Below we will generate and plot such a power curve. We do so for effects ranging from 0 to 0.7:

s1 <- (1 + (350 - 1) * 0.109) / (32 * 350)
s0 <- (1 + (350 - 1) * 0.109) / (33 * 350)
SE <- sqrt(2.456^2 * (s0 + s1))

effect <- seq(0, 1, length.out = 51)
power  <- pnorm(effect / SE - qnorm(0.975))
df     <- data.frame(effect, power)


library(ggplot2)
library(scales)

ggplot(df, aes(x = effect, y = power)) +
  geom_line() +
  scale_x_continuous(
    name = expression("True effect " (tau)),
    breaks = seq(0, 1, by = 0.2),
    labels = number_format(accuracy = 0.1)
  ) +
  scale_y_continuous(
    name = "Power",
    limits = c(0, 1),
    breaks = seq(0, 1, by = 0.2),
    labels = number_format(accuracy = 0.1)
  ) +
  theme_minimal()

Power curve for effect on test scores with clustered randomisation

Upon doing this, we see the standard features of a power curve with an ‘upside-down’ normal shape, and see that this reaches the commonly employed 80% power threshold at the value indicated above (0.571), approximately reaching 100% power at effect sizes of 1.

Of course, the power in this setting depends upon a number of quantities, and not the true effect. We can see this by considering the power curve corresponding to each outcome laid out in Table Table 1. Below we loop through each of the outcomes from the table taking their assumed standard deviation (SDs) and intra-cluster correlation (ICC), calculating the resulting standard error, and then calculating the power curve at the values for effect defined above. Here, simultaneous iteration over SDs and ICC is handled naturally since both are parallel vectors of the same length, so SEs <- sqrt(SDs^2 * (n0 + n1)) simply works element-wise with all 12 standard errors. We then use sapply to apply the power formula across each of these 12 standard errors in turn: for each value se in SEs, it computes the full power curve over our effect grid and returns the result as a column, producing a matrix power_mat we can use below.

SDs <- c(2.456, 2.372, 0.499, 0.500, 0.881, 0.542, 
         0.570, 0.659, 0.524, 0.660, 0.581, 2.613)
ICC <- c(0.109, 0.103, 0.016, 0.015, 0.021, 0.035, 
         0.005, 0.049, 0.015, 0.019, 0.014, 0.033)

n1 <- (1 + (350 - 1) * ICC) / (32 * 350)
n0 <- (1 + (350 - 1) * ICC) / (33 * 350)

SEs <- sqrt(SDs^2 * (n0 + n1))
effect <- seq(0, 1, length.out = 51)


power_mat <- sapply(SEs, function(se) {
    pnorm(effect / se - qnorm(0.975))
})

df <- data.frame(effect = effect, power_mat)
names(df)[-1] <- paste0("power_y", seq_along(SEs))

We can now plot these power curves for each of the outcomes, where we of course see that the power to detect specific effects varies considerably by outcomes, as laid out in Table Table 1.

library(tidyr)
library(dplyr)


Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

df_long <- df %>%
    pivot_longer(
        cols      = starts_with("power_y"),
        names_to  = "outcome",
        values_to = "power"
    ) %>%
    filter(effect < 0.6) %>%
    mutate(
        outcome = factor(
            outcome,
            levels = paste0("power_y", 1:12),
            labels = c(
                "Turkish score", "Math score",
                "Bullying (class)", "Bullying (school)",
                "Sensitivity", "Locus of control",
                "Impulsivity", "Perspective",
                "Wellbeing", "Belonging",
                "Autonomy", "Friendship"
            )
        )
    )

ggplot(df_long, aes(x = effect, y = power, color = outcome)) +
    geom_line() +
    scale_x_continuous(
        name   = expression("True effect " (tau)),
        breaks = seq(0, 0.6, by = 0.1),
        labels = number_format(accuracy = 0.1)
    ) +
    scale_y_continuous(
        name   = "Power",
        labels = number_format(accuracy = 0.1)
    ) +
    theme_minimal() +
    theme(
        legend.title = element_blank()
    )

The above assumes clustered randomisation. Of course if we were working with unit-level randomisation we could proceed similarly, but in this case our standard error calculation is far simpler, and we would be very likely be more powered given the larger effective sample size.

# individual‐randomisation standard error
SEind <- sqrt(2.456^2 / 22750)
effect <- seq(0, 1, length.out = 51)

powerNoCluster <- pnorm(effect / SEind - qnorm(0.975))
df_ind <- data.frame(effect, powerNoCluster)


ggplot(df_ind, aes(x = effect, y = powerNoCluster)) +
    geom_line() +
    scale_x_continuous(
        name   = expression("True effect " (tau)),
        breaks = seq(0, 1, by = 0.2),
        labels = number_format(accuracy = 0.1)
    ) +
    scale_y_continuous(
        name   = "Power",
        limits = c(0, 1),
        breaks = seq(0, 1, by = 0.2),
        labels = number_format(accuracy = 0.1)
    ) +
    theme_minimal()

Power curve for effect on test scores with individual randomization

Finally note that if instead of asking what the MDE is for a given case we would rather determine the minimum required sample size to determine an effect, we can do so by applying (7.12) in the book. Let’s see this for a power of 0.8:

z80  <- qnorm(0.8)
z975 <- qnorm(0.975)

4*(z80+z975)^2/(0.4^2/2.456^2)

[1] 1183.599

4*(z80+z975)^2/(0.3^2/2.456^2)

[1] 2104.175

4*(z80+z975)^2/(0.2^2/2.456^2)

[1] 4734.394

4*(z80+z975)^2/(0.1^2/2.456^2)

[1] 18937.58

Let’s confirm that if we use those sample sizes we do indeed see such minimum detectable effect sizes.

effect <- seq(0, 1, length.out = 51)


SE_04 <- sqrt(2.456^2 / 1184)
SE_03 <- sqrt(2.456^2 / 2104)
SE_02 <- sqrt(2.456^2 / 4734)
SE_01 <- sqrt(2.456^2 / 18938)

z975 <- qnorm(0.975)

power_04 <- pnorm(effect / SE_04 - z975)
power_03 <- pnorm(effect / SE_03 - z975)
power_02 <- pnorm(effect / SE_02 - z975)
power_01 <- pnorm(effect / SE_01 - z975)


df <- data.frame(
    effect,
    `N = 592`  = power_04,
    `N = 1,052` = power_03,
    `N = 2,367` = power_02,
    `N = 9,469` = power_01
)


df_long <- df %>%
    pivot_longer(-effect, names_to = "Sample", values_to = "Power") %>%
    filter(effect < 0.6)


ggplot(df_long, aes(x = effect, y = Power, color = Sample)) +
    geom_line() +
    scale_x_continuous(
        name   = expression("True effect " (tau)),
        breaks = seq(0, 0.6, by = 0.1),
        labels = number_format(accuracy = 0.1)
    ) +
    scale_y_continuous(
        name   = "Power",
        labels = number_format(accuracy = 0.1)
    ) +
    theme_minimal() +
    theme(legend.title = element_blank())

Minimum sample size to detect specific effects

While there are various other calculations we can make (or indeed, we could generate power curves by simulation in cases where settings are more complex), above we see that it is relatively straightforward to both generate power curves as well as arrive at minimum detectable effect sizes or minimum sample sizes required to detect effects. This is a key step in experimental design where we will wish to ensure that our study is actually sufficiently powerful to detect any real effects that may exist with sufficient precision.

Code call-out 7.2: Bonferronni, Holm and Romano-Wolf Procedures

We are interested in examining a number of corrections for multiple hypothesis tests which provide control of the family-wise error rate of tests. To do so, we will use data from Charness and Gneezy (2009). In their paper, they examine whether a group of individuals who receive an incentive to exercise are observed to have greater improvements in a number of markers of health compared with individuals who are not incentivised to go to the gym. Below we load the data from their “Experiment 2”, which has be most complete set of health and biometric measures, make a small number of changes to correct for a missing value, and remove a variable which we will regenerate ourselves below:

library(haven)
library(dplyr)

df <- read_dta("data/Charness_Gneezy_2009.dta") %>%
    mutate(bmi2 = na_if(bmi2, 0)) %>%
    select(-bmi_diff) %>%
    filter(one != 1)

In the last step, we have dropped all observations flagged as having one==1. Charness and Gneezy (2009) actually consider two treatments: one in which individuals are given one reminder and one in which they are given eight reminders, as well as a control condition. In the interests of simplicity here we consider only those treated more intensively with eight reminders (eight==1), though if we wished, we could reproduce the entire analysis below by comparing those treated with one reminder to the control condition.

Here we will consider the tests documented in Table III of Charness and Gneezy (2009) in which changes in biometric measures of treated individuals (eight==1) are compared with measures among untreated individuals (eight==0). Below we will define a local which contains each of the variables studied, and then we generate an outcome which captures the change in these measures over time for each individual. This is based on differences between the third measure taken (a follow-up wave), and a baseline measure.

vars <- c("bodyfat","pulse","weight","bmi","waist","sbp","dbp")
for (v in vars) {
    df[[paste0(v, "_diff")]] <- df[[paste0(v, "1")]] - df[[paste0(v, "3")]]
}

We will seek to calculate 4 p-values for each test. Specifically, we will calculate p-values not considering the fact that we are conducting multiple tests, and we will generate p-values correcting for multiple testing by applying the procedures of Bonferroni, Holm, and Romano and Wolf. We will start by simply generating an empty matrix which we will fill out with each procedure below:

pvalues <- matrix(
    NA_real_,
    nrow = length(vars),
    ncol = 4,
    dimnames = list(
        vars,
        c("uncorrected", "Bonferroni", "Holm", "Romano-Wolf")
    )
)
pvalues

        uncorrected Bonferroni Holm Romano-Wolf
bodyfat          NA         NA   NA          NA
pulse            NA         NA   NA          NA
weight           NA         NA   NA          NA
bmi              NA         NA   NA          NA
waist            NA         NA   NA          NA
sbp              NA         NA   NA          NA
dbp              NA         NA   NA          NA

Let’s start by simply running the regression which essentially replicates the results of Charness and Gneezy (2009). Specifically, we will regress each differenced measure on treatment receipt, which, given random assignment, should allow us to calculate an ATT. In practice, Charness and Gneezy (2009) document these differences in a more extended manner in their Table III, but if we compare the below to this Table, we can see that estimated effects line up very well.

library(sandwich)
library(lmtest)

Loading required package: zoo


Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

for (v in vars) {
  # Fit model
  formula <- as.formula(paste0(v, "_diff ~ eight"))
  model <- lm(formula, data = df)
  
  # Compute robust (HC1) SEs and display
  robust <- coeftest(model, vcov = vcovHC(model, type = "HC1"))
  print(robust)
}


t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -1.41026    0.41470 -3.4006 0.0009775 ***
eight        2.18859    0.46479  4.7088   8.3e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -3.8974     2.0771 -1.8764  0.06361 .
eight         5.1474     2.7401  1.8786  0.06331 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.57179    0.54379 -1.0515   0.2956
eight        0.91179    0.59976  1.5203   0.1317


t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.23137    0.19078 -1.2128   0.2282  
eight        0.35497    0.21189  1.6753   0.0971 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.071795   0.359328 -0.1998   0.8421  
eight        0.796795   0.425436  1.8729   0.0641 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -5.2308     1.7024 -3.0726 0.002754 **
eight         3.4474     2.6190  1.3163 0.191173   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -2.87179    1.21621 -2.3613  0.02021 *
eight        0.28846    1.80280  0.1600  0.87321  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Because in certain multiple-hypothesis correction procedures below we will wish to work from the most to least significant hypothesis, below we we-order our variable list in terms of their statistical significance, ie in ordered from lowest p-value to highest.

vars <- c("bodyfat", "pulse", "waist", "bmi", "weight", "sbp", "dbp")

Now, let’s actually start to examine p-values corresponding to each test. Below we will consider the most simple case where no correction is implemented at all. Below, we simply iterate through variables, running the regression with each outcome as a dependent variable, and saving the resulting p-value in the pvalues matrix.

for (i in seq_along(vars)) {
    var <- vars[i]
    # Fit model
    model <- lm(as.formula(paste0(var, "_diff ~ eight")), data = df)
    # Robust SEs
    robust <- coeftest(model, vcov = vcovHC(model, type = "HC1"))
    # Extract p‐value for 'eight'
    p_val <- robust["eight", "Pr(>|t|)"]
    # Store in pvalues matrix (column 1 = "uncorrected")
    pvalues[i, "uncorrected"] <- p_val
}

pvalues

         uncorrected Bonferroni Holm Romano-Wolf
bodyfat 8.299681e-06         NA   NA          NA
pulse   6.330761e-02         NA   NA          NA
weight  6.409514e-02         NA   NA          NA
bmi     9.709828e-02         NA   NA          NA
waist   1.316954e-01         NA   NA          NA
sbp     1.911727e-01         NA   NA          NA
dbp     8.732082e-01         NA   NA          NA

Here we see that we have simply filled in the first column of the matrix, observing an ascending set of values corresponding to each outcome.

Bonferroni’s correction

Let’s now implement the Bonferroni correction. As laid out in Section 7.3.2.1 of the book, this is a relatively simple procedure in which we simply multiply all p-values by the total number of tests (the logic of this is discussed in the chapter itself). We do this below, noting just one specific point of care, which is that we want to ensure that p-values are capped at 1, and so we use the min() function below.

# Bonferroni adjustment: p * 7, capped at 1
pvalues[, "Bonferroni"] <- pmin(pvalues[, "uncorrected"] * 7, 1)
pvalues

         uncorrected   Bonferroni Holm Romano-Wolf
bodyfat 8.299681e-06 5.809776e-05   NA          NA
pulse   6.330761e-02 4.431533e-01   NA          NA
weight  6.409514e-02 4.486660e-01   NA          NA
bmi     9.709828e-02 6.796880e-01   NA          NA
waist   1.316954e-01 9.218680e-01   NA          NA
sbp     1.911727e-01 1.000000e+00   NA          NA
dbp     8.732082e-01 1.000000e+00   NA          NA

Here we see that–as expected–our Bonferroni p-values simply inflate all baseline p-values constantly, though suggest that even with this most demanding criteria, impacts of treatment on bodyfat are still statistically significant.

Holm’s correction

As we saw in the book, we can gain power while still maintaining the familywise error rate by following step-wise techniques. A first of these is Holm’s correction (see section 7.3.2). We implement this below, where the key point to note is that here instead of multiplying all p-values by 7, we mutliply the first (lowest) p-valuy by 7, then the second lowest by 6, and so on and so forth until the final p-value is simply returned as is.

Perhaps the only important thing to note in this procedure is that now we must ensure that we enforce monotonicity. Monotonicity simply implies that if we have some original p-value ordering from lowest to highest, our resulting corrected p-values should still maintain this ordering. For this reason, if some variable with a higher uncorrected p-value returns a lower corrected p-value, we simply replace the corrected p-value with the value of the p-value observed in the preceding variable. This is implemented below on the lines following the comment # enforce monotonicity.

prev <- 0
for (i in seq_len(nrow(pvalues))) {
    # 8 refer to number of tests + 1, so this goes from 7, 6, ..., 1
    m <- 8 - i
    pval <- min(pvalues[i, "uncorrected"] * m, 1)
  
    # enforce monotonicity
    if (pval < prev) {
        pvalues[i, "Holm"] <- prev
    } else {
        pvalues[i, "Holm"] <- pval
        prev <- pval
    }
}
pvalues

         uncorrected   Bonferroni         Holm Romano-Wolf
bodyfat 8.299681e-06 5.809776e-05 5.809776e-05          NA
pulse   6.330761e-02 4.431533e-01 3.798457e-01          NA
weight  6.409514e-02 4.486660e-01 3.798457e-01          NA
bmi     9.709828e-02 6.796880e-01 3.883931e-01          NA
waist   1.316954e-01 9.218680e-01 3.950863e-01          NA
sbp     1.911727e-01 1.000000e+00 3.950863e-01          NA
dbp     8.732082e-01 1.000000e+00 8.732082e-01          NA

Here we see the additional power of Holm’s correction over Bonferroni. With the exception of the lowest p-value, p-values are all strictly lower than their counterpart in Bonferroni. We also observe a second regularity which is that the lowest p-value corresponds to that of the Bonferroni procedure, while the highest corresponds to the uncorrected p-value.

Romano-Wolf

While Holm offers benefits over Bonferroni in terms of power, this is still based on the “worst-case scenario” that all test statistics are entirely independent. As laid out in Romano and Wolf, We can improve on this using bootstrap-based procedures which directly estimate the dependence among tests. A full description of this test and its logic can be found at the end of Section 7.3.2.1 of the book.

Here we will implement this procedure “by hand”. To do so, we must estimate many bootstrap replicates of each model, and in each case save the resulting test-statistic of each test once the null hypothesis (in this case of a 0 effect) has been imposed. We will begin this procedure by setting up a new frame with a large number of observations to store the test statistic from each bootstrap replicate, and in this frame we will generate a variable to store the statistic corresponding to each variable:

# Number of bootstrap replications
B <- 5000

# Create a data frame with B rows and the specified columns, initialized to NA
bstraps <- data.frame(matrix(nrow = B, ncol = length(vars)))
colnames(bstraps) <- vars

Now, before we actually go about implementing the bootstrap procedure, we will just re-estimate each of our true regressions, saving the estimated point estimate and t-statistic for use below.

beta  <- setNames(numeric(length(vars)), vars)
tstat <- setNames(numeric(length(vars)), vars)

for (v in vars) {
    model <- lm(as.formula(paste0(v, "_diff ~ eight")), data = df)
    coef  <- coef(model)["eight"]
    se    <- sqrt(vcovHC(model, type = "HC1")["eight","eight"])
    beta[v]  <- coef
    tstat[v] <- coef / se
}

beta

  bodyfat     pulse     waist       bmi    weight       sbp       dbp 
2.1885897 5.1474359 0.7967949 0.3549698 0.9117949 3.4474359 0.2884615

tstat

  bodyfat     pulse     waist       bmi    weight       sbp       dbp 
4.7087603 1.8785602 1.8728893 1.6752900 1.5202730 1.3163080 0.1600071

With this in hand, we will do the bootstrap procedure itself. Specifically, we will run a series of B resamples, and in each loop we will calculate the t-statistic where we impose the null hypothesis of a zero effect. This is calculated (in each bootstrap replicate \(b\)), as: \[ t = \frac{\widehat\tau^{*b}-\widehat\tau}{\widehat\sigma^{*b}}, \] where \(\widehat\tau^{*b}\) and \(\widehat\sigma^{*b}\) refer to the bootstrap estimate and standard error respectively, and \(\widehat\tau\) the original estimate. The idea of this is that on average this quantity will be centred around 0 across all bootstraps, but it will correctly capture the variation inherent in our data.

We do this below, iterating through each our B boostrap replicates, and within each replicate, estimate each of our outcomes and the resulting set of t-statistics. We take the absolute value of these as below we will be interested in considering whether our original t-statistic is extreme (either high or low) compared to our bootstrap t-statistics, and this can easily be done by just taking absolute values of each, and comparing the absolute values. We then store each t-statistic in our data frame called bstraps which we will return to use below.

# Bootstrap loop
for (b in seq_len(B)) {
    df_boot <- df[sample(nrow(df), replace = TRUE), ]
    for (v in vars) {
        model_boot <- lm(as.formula(paste0(v, "_diff ~ eight")), data = df_boot)
        se_boot    <- sqrt(vcovHC(model_boot, type = "HC1")[ "eight","eight" ])
        t_val      <- abs((coef(model_boot)["eight"] - beta[v]) / se_boot)
        bstraps[b, v] <- t_val
    }
}

To see what we have done so far, let’s examine the distribution of the bootstrap replicates for one of the variables:

library(ggplot2)
ggplot(bstraps, aes(x = bodyfat)) +
    geom_histogram(binwidth = 0.1) +
    labs(
        x       = "Bootstrap null |t-statistic|",
        y       = "Density"
    ) +
    theme_minimal()

Here we see that this is effectively a half normal distribution, as we would expect given that we have centred this around 0 and standardised. This in essence provides us a null distribution for a single test: if we observe that our true t-value is extreme when compared to this null distribution, we can view this as strong evidence against a null of a zero effect.

However, in Romano-Wolf, what we want is not to consider the single test, but to correct for multiple testing. And here we follow a step-wise procedure in which for our most significant variable, we compare our test-statistic with the maximum t-value across all hypotheses tested, for our second most significant, consider only the next 6 variables, and so on and so forth. We will generate these comparison distributions for each test below:

library(tidyr)

bstraps <- bstraps %>%
  mutate(
    null_bodyfat = pmax(bodyfat, pulse, waist, bmi, weight, sbp, dbp),
    null_pulse   = pmax(pulse, waist, bmi, weight, sbp, dbp),
    null_waist   = pmax(waist, bmi, weight, sbp, dbp),
    null_bmi     = pmax(bmi, weight, sbp, dbp),
    null_weight  = pmax(weight, sbp, dbp),
    null_sbp     = pmax(sbp, dbp),
    null_dbp     = dbp
  )


bstraps_long <- bstraps %>%
  select(starts_with("null_")) %>%
  pivot_longer(
    cols      = everything(),
    names_to  = "variable",
    values_to = "value"
  ) %>%
  mutate(variable = factor(
    variable,
    levels = c("null_bodyfat","null_pulse","null_waist","null_bmi",
               "null_weight","null_sbp","null_dbp"),
    labels = c("Body fat","Pulse","Waist","BMI",
               "Weight","Sistolic BP","Diastolic BP")
  ))

This very much follows the logic of Holm, and indeed, if each variable is entirely independent, this will essentially revert to the same logic as Holm’s correction. However, if the variables are not independent among themselves, our null distributions here will account for this dependence, and generally offer a less demanding null distribution (and more powerful test). We can see how our null distributions get progressively less demanding as the p-value of the original test rises by plotting these null distributions below:

ggplot(bstraps_long, aes(x = value, color = variable)) +
  geom_density() +
  labs(
    x     = "Bootstrap |t-statistic| under the null",
    y     = "Density",
    color = NULL
  ) +
  theme_minimal() +
  theme(legend.position = c(0.85, 0.85))

Now, finally, let’s just calculate our Romano-Wolf p-values. Here–and as laid out in the book–we simply wish to compare our original t-values with those in the null distribution in each case. If our original t-values are extreme compared to the null distribution, this results in a low p-value (and strong evidence against the null), while if our original t-values appear likely to occur when the null is true, we will have little evidence to reject the null. We calculate our p-values below, where as above, we enforce monotonocity to avoid unordered p-values in final tests:

prev <- 0
for (i in seq_along(vars)) {
  null_vals <- do.call(pmax, c(bstraps[vars[i:length(vars)]], na.rm = TRUE))
  p <- mean(null_vals > abs(tstat[vars[i]]))
  if (p < prev) p <- prev
  pvalues[i, "Romano-Wolf"] <- p
  prev <- p
}
pvalues

         uncorrected   Bonferroni         Holm Romano-Wolf
bodyfat 8.299681e-06 5.809776e-05 5.809776e-05      0.0002
pulse   6.330761e-02 4.431533e-01 3.798457e-01      0.2854
weight  6.409514e-02 4.486660e-01 3.798457e-01      0.2854
bmi     9.709828e-02 6.796880e-01 3.883931e-01      0.2854
waist   1.316954e-01 9.218680e-01 3.950863e-01      0.3628
sbp     1.911727e-01 1.000000e+00 3.950863e-01      0.3628
dbp     8.732082e-01 1.000000e+00 8.732082e-01      0.8686

Here we see that as expected, Romano-Wolf p-values are lower, and at times substantially lower, to the values from Holm. We also see that for the final test, these p-values are the same (subject to some minor bootstrap variance), as the uncorrected p-values. In general, we see here the benefits in terms of power which we may expect in this case.

Finally note, that while it is very useful to see that we can calculate such p-values by hand, there are also libraries available which implement these automatically. Below we use the wildrwolf package which automates the procedure we have conducted above, returning to us the Romano-Wolf p-values alongside uncorrected p-values. These values are essentially identical (subject to bootstrap variation) to what we calculated by hand ourselves. Note that wildrwolf can be installed from GitHub via remotes::install_github("s3alfisc/wildrwolf") provided that the remotes package is installed.

library(fixest)


Attaching package: 'fixest'

The following object is masked from 'package:scales':

    pvalue

library(wildrwolf)

set.seed(121316)

fits <- lapply(vars, function(v) {
  feols(as.formula(paste0(v, "_diff ~ eight")), data = df, se = "hetero")
})

rw_tbl <- rwolf(models = fits, param = "eight", B = 5000)


  |                                                                            
  |                                                                      |   0%

Warning: Please note that the seeding behavior for random number generation for
`boottest()` has changed with `fwildclusterboot` version 0.13.

It will no longer be possible to exactly reproduce results produced by versions
lower than 0.13.

If your prior results were produced under sufficiently many bootstrap
iterations, none of your conclusions will change. For more details about this
change, please read the notes in
[news.md](https://cran.r-project.org/web/packages/fwildclusterboot/news/news.html).
This warning is displayed once per session.


  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |======================================================================| 100%

print(rw_tbl)

  model  Estimate Std. Error   t value     Pr(>|t|) RW Pr(>|t|)
1     1   2.18859  0.4647911   4.70876 8.299681e-06  0.00019996
2     2  5.147436   2.740096   1.87856   0.06330761  0.27174565
3     3 0.7967949  0.4254362  1.872889   0.06409514  0.27174565
4     4 0.3549698  0.2118856   1.67529   0.09709828  0.27174565
5     5 0.9117949  0.5997573  1.520273    0.1316954  0.34093181
6     6  3.447436   2.619019  1.316308    0.1911727  0.34093181
7     7 0.2884615   1.802804 0.1600071    0.8732082  0.87562488

Code call-out 7.3: Principal Components, Summary Indexes and Anderson’s Index

In this code call-out, we will consider multiple hypothesis testing with some simulated data and use these simulations to understand why the decision of which type of index to use is not innocuous.

A Simple Illustration of Over-Rejection with Multiple Hypotheses

To begin, let’s just consider how we may explore a controlled simulation with multiple hypothesis testing. Specifically, let’s imagine that we wish to run the following regression: \[ y_i^k=\alpha + \tau^k Treat_i + \varepsilon_i^k \forall k\in\{1,\ldots,K\}. \] Here we use the super-index \(k\) to indicate that these are different outcome variables (in particular \(K\) different outcome variables), and these are regressed on a single treatment value. Let’s conduct some simulations where the true treatment effect for each of these models \(\tau^k\) is indeed 0. Let’s also simulate a treatment variable for 100 observations such that half of the population is randomly assigned treatment, and the other half is not assigned treatment. We will simulate each \(\varepsilon_i\sim\mathcal{N}(0,1)\). Then, let’s estimate the effect \(\tau^k\) in each case, and see if we reject the (true) null hypothesis that \(\tau^k=0\).

We can do this quite simply below, where we first set 100 observations and then generate a single treatment variable as treat. We will then simulate for \(K=10\) outcomes where the treatment effect is indeed 0, and run the regression of each outcome on treatment one at a time:

rm(list = ls())
set.seed(121316)

K <- 10
n <- 100
treat <- as.numeric(runif(n) > 0.5)

for (k in seq_len(K)) {
    y <- 2 + 0*treat + rnorm(n)
    model <- lm(y ~ treat)
  
    beta <- coef(model)["treat"]
    se   <- summary(model)$coefficients["treat", "Std. Error"]
    z    <- beta / se
    p_val <- pnorm(abs(z))
    cat(sprintf("p-value of test %d is: %g\n", k, p_val))
}

p-value of test 1 is: 0.948836
p-value of test 2 is: 0.633301
p-value of test 3 is: 0.991094
p-value of test 4 is: 0.932929
p-value of test 5 is: 0.76734
p-value of test 6 is: 0.799159
p-value of test 7 is: 0.939336
p-value of test 8 is: 0.508681
p-value of test 9 is: 0.576478
p-value of test 10 is: 0.94407

In this particular example, we can see that of these 10 simulated outcomes where the treatment has no effect on each outcome \(y^k\), we correctly fail to reject the null hypothesis that the effect is 0 based on the 10 p-values resulting from the regression above if a critical value of \(\alpha=0.05\) is used. Of course, this is entirely dependent upon the particular simulation which we define with the random seed using in set seed. For this reason, to see the concerns about multiple testing more generally, we will repeat this process 5000 times, observing in each of the 5000 cases whether we (falsely) reject at least one hypothesis within our \(K\) hypotheses, or falsely reject at least two hypotheses within our \(K\) tests.

We do this below, where we have essentially replicated the above block of code, but we now repeat the test \(S=5000\) times and in each case we determine whether any hypotheses are rejected, saving these in a series of variables. Finally, we can calculate the quantity of rejected hypotheses of the null of a zero effect (which, from our DGP, we know to be true), as well as the proportion of times we reject at least 1 null hypothesis and at least 2 null hypotheses within our class of \(K\) tests.

set.seed(1213)
S <- 5000
K <- 10

rejectK <- propK <- prop2K <- 0

for(s in 1:S) {
    # binary 0/1 treatment
    treat <- rbinom(100, 1, 0.5)
    rejectAny <- 0
    for(k in 1:K) {
        y <- 2 + 0*treat + rnorm(100)
        m <- lm(y ~ treat)
        z <- coef(m)["treat"] / summary(m)$coefficients["treat", "Std. Error"]
        if (abs(z) > qnorm(0.975)) {
            rejectK   <- rejectK + 1
            rejectAny <- rejectAny + 1
        }
    }
    if (rejectAny > 0) propK  <- propK + 1
    if (rejectAny > 1) prop2K <- prop2K + 1
}

mtr <- rejectK / (K * S)
ptr <- propK   / S
p2r <- prop2K  / S

cat("Total Tests Rejected: ", rejectK, "\n")

Total Tests Rejected:  2715

cat("Mean Tests Rejected:  ", sprintf("%04.3f", mtr), "\n")

Mean Tests Rejected:   0.054

cat("Proportion ≥1 rejections: ", sprintf("%04.3f", ptr), "\n")

Proportion ≥1 rejections:  0.428

cat("Proportion ≥2 rejections: ", sprintf("%04.3f", p2r), "\n")

Proportion ≥2 rejections:  0.098

In this case, we see that (as expected), while only 5% of all hypotheses are rejected (or, to be exact, 5.4%), when we consider the number of times that at least 1 hypothesis is rejected, this is much higher, at 41.8%. What’s more, if we consider cases with at least 2 rejections, we see that this is also far higher that 5%, at around 9.4%.

Of course, the degree to which such overrejections occurs will depende upon the number of false null hypotheses considered. Below we can see this where we simply replicate the above analysis, however here instead of doing this just for \(K=10\) variables, we do it for a range of values. This simply adds an outer loop around the procedure above, considering \(K=1, 2, \ldots, 20\) variables within a family of tests, where in each case there is no effect of treatment on any of the variables considered.

set.seed(1213)
S <- 5000

results <- data.frame(
    numvars    = 1:20,
    propreject = NA_real_,
    prejectAny = NA_real_
)

for (K in seq_len(20)) {
    rejectK    <- 0
    propAnyK   <- 0
  
    for (s in seq_len(S)) {
          treat     <- rbinom(100, 1, 0.5)
          rejectAny <- 0
      
          for (k in seq_len(K)) {
                y     <- 2 + 0*treat + rnorm(100)
                m     <- lm(y ~ treat)
                z     <- coef(m)["treat"] / summary(m)$coefficients["treat", "Std. Error"]
                if (abs(z) > qnorm(0.975)) {
                    rejectK    <- rejectK + 1
                    rejectAny  <- rejectAny + 1
                }
          }
      
          if (rejectAny > 0) propAnyK <- propAnyK + 1
    }
  
    results$propreject[K] <- rejectK / (K * S)
    results$prejectAny[K] <- propAnyK / S
}

results

   numvars propreject prejectAny
1        1 0.05120000     0.0512
2        2 0.05680000     0.1114
3        3 0.05486667     0.1546
4        4 0.05245000     0.1940
5        5 0.05520000     0.2496
6        6 0.05236667     0.2732
7        7 0.05171429     0.3090
8        8 0.05230000     0.3548
9        9 0.05333333     0.3896
10      10 0.05250000     0.4174
11      11 0.05300000     0.4530
12      12 0.05306667     0.4822
13      13 0.05293846     0.5034
14      14 0.05432857     0.5478
15      15 0.05333333     0.5600
16      16 0.05316250     0.5846
17      17 0.05396471     0.6144
18      18 0.05333333     0.6296
19      19 0.05370526     0.6554
20      20 0.05215000     0.6558

Finally, to view the output, we can simply plot the rate of rejections of all hypotheses considered by the number of hypotheses tested, as well as the proportion of times that at least 1 hypothesis is rejected. We do this in the graph below:

library(ggplot2)
library(tidyr)
library(dplyr)

df_long <- results %>%
    pivot_longer(
        cols = c(propreject, prejectAny),
        names_to = "Metric",
        values_to = "Proportion"
    ) %>%
    mutate(
        Metric = factor(
          Metric,
          levels = c("propreject", "prejectAny"),
          labels = c("Proportion of all hypotheses",
                     "Proportion at least 1 rejection"
          )
      )
  )

ggplot(df_long, aes(x = numvars, y = Proportion, color = Metric)) +
    geom_line() +
    scale_y_continuous(
        name = "Proportion rejected hypotheses",
        labels = scales::number_format(accuracy = 0.1)
    ) +
    scale_x_continuous(
        name = "Number of variables in class",
        breaks = seq(1, 20, by = 1)
    ) +
    labs(color = NULL) +
    theme_minimal()

As expected, we can see that while the proportion of all rejected hypotheses is stable at around 5%, the number of times at least 1 hypothesis is rejected accumulates as the number of hypotheses tested grows, reaching very high levels (around 65%) when 20 hypotheses are in a class.

Indexes and Dimension Reduction

As laid out in the book, there are multiple ways to deal with this. One way is to seek to reduce the number of outcomes (and hence hypotheses) to a single aggregate measure. We will explore this here. If all variables are indeed entirely independent as in the simulations above, the way an index is formed is a relatively innocuous choice. However, when variables are correlated among themselves, as is likely the case in real empirical settings, we will see here that the choice of indexes is very much not innocuous. In order to explore this below, we will set up some simulations where unobserved errors corresponding to each outcome of interest are correlated among themselves. Specifically, we will seek to simulate a series of variables which depend upon unobservables drawn from a multivariate normal distribution such that: \((\varepsilon_i^1, \varepsilon_i^2, \ldots \varepsilon_i^K)\sim\mathcal{N}([0,0,\ldots,0], \Omega)\), where \(\Omega\) is a correlation matrix: \[ \Omega = \begin{bmatrix} 1 & \rho & \ldots & \rho \\ \rho & 1 & \ldots & \rho \\ \vdots & \vdots & \ddots & \vdots \\ \rho & \rho & \ldots & 1 \\ \end{bmatrix} \] Thus, if \(\rho=0\) unobserved error terms will be completely independent, while if \(0<\rho\leq 1\), variables will be positively correlated. We will then simulate a single treatment variable \(Treat_i\) which takes 1 for around half of the observations and 0 for all others, and finally, simulate the outcomes themselves as: \[ y_i^k = \tau^k Treat_i + \varepsilon^k_i \]

Below we will set up such a simulation based on 1,000 observations. In this simulation we will consider \(K=10\) strongly correlated variables (setting \(\rho=0.9\)), and additionally simulate one variable which is entirely uncorrelated with other variables. For 9 of the 10 uncorrelated variables we will set \(\tau^{k}=0\), while for the tenth variable, as well as the uncorrelated variable, we will set \(\tau=0.5\):

set.seed(1213)
n <- 1000
r <- 0.9

corr <- matrix(r, nrow = 10, ncol = 10)
diag(corr) <- 1

library(MASS)


Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select

u <- mvrnorm(n, mu = rep(0, 10), Sigma = corr)
df <- as.data.frame(u)
names(df) <- paste0("u", 1:10)

df$treat <- as.numeric(runif(n) > 0.5)

for (i in 1:9) {
    df[[paste0("y", i)]] <- df[[paste0("u", i)]]
}

df$y10 <- 0.5 * df$treat + df$u10
df$y11 <- 0.5 * df$treat + rnorm(n)

Now with these variables in hand, let’s consider three of the potential indexes discussed in the book. These are Anderson’s index which overweights variables which bring independent variation, a principal-components based index, and a simple summary index. Be we generate these below, we define our own function to implement the index described by Anderson (2008), and which takes as arguments a dataframe, a series of variables to be used in the index, and an indicator of which observations are control units.

anderson_index <- function(df, y_vars, control_indicator) {
    # Subset to control group
    df_control <- df[control_indicator, y_vars]
    
    # Standardise using control group mean and SD
    means <- colMeans(df_control, na.rm = TRUE)
    sds   <- apply(df_control, 2, sd, na.rm = TRUE)
    
    y_std <- sweep(sweep(df[y_vars], 2, means, "-"), 2, sds, "/")
    
    # Covariance matrix from standardised control group outcomes
    y_std_control <- y_std[control_indicator, ]
    cov_matrix    <- cov(y_std_control)
    
    # Invert the covariance matrix
    cov_inv <- solve(cov_matrix)
    
    # Weights are row sums of the inverse covariance matrix
    weights <- rowSums(cov_inv)
    
    # Weighted average across outcomes for each observation
    index <- as.matrix(y_std) %*% weights / sum(weights)
    
    return(as.numeric(index))
}

Finally, we generate our three indexes, and examine whether we detect any effect of treatment on each of the aggregate indexes using a simple OLS regression:

library(sandwich)
library(lmtest)

y_vars <- paste0("y", 1:10)

# Anderson's index
control_idx <- df$treat == 0
df$AndersonIndex <- anderson_index(df, y_vars, control_idx)

# Principal‐component index
pca_model <- prcomp(df[y_vars], center = TRUE, scale. = FALSE)
df$PCIndex <- pca_model$x[, 1]

# Summary index
df$SumIndex <- rowSums(df[y_vars])

# Regressions with robust SEs
for (idx in c("AndersonIndex", "SumIndex", "PCIndex")) {
    m <- lm(as.formula(paste(idx, "~ treat")), data = df)
    print(coeftest(m, vcov = vcovHC(m, type = "HC1")))
}


t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)
(Intercept)  4.3281e-17  4.2925e-02  0.0000   1.0000
treat       -1.9721e-02  5.9449e-02 -0.3317   0.7402


t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.20248    0.45773  0.4424   0.6583
treat       -0.21643    0.63265 -0.3421   0.7323


t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.032417   0.144757 -0.2239   0.8228
treat        0.064575   0.200073  0.3228   0.7469

In this case, where there is a single true effect, we observe that if we consult the summary index or the principal components index we find some evidence of a significant effect of treatment on aggregate outcomes. On the other hand, if we consult the Anderson index we observe very little evidence of a significant effect. However, below we consider precisely the same setting, but now instead of examining an index based on the 10 correlated variables (with one true effect), we consider 9 correlated variables (with no effects), and 1 uncorrelated variable with a true effect.

vars2 <- c(paste0("y", 1:9), "y11")
control_idx <- df$treat == 0

pca2 <- prcomp(df[vars2], center = TRUE, scale. = FALSE)
df$PCIndex2       <- pca2$x[, 1]
df$SumIndex2      <- rowSums(df[vars2])
df$AndersonIndex2 <- anderson_index(df, vars2, control_idx)

summary(lm(AndersonIndex2 ~ treat, data = df))


Call:
lm(formula = AndersonIndex2 ~ treat, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.63020 -0.50798  0.00036  0.47773  2.29713 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.577e-18  3.183e-02   0.000        1    
treat       1.915e-01  4.492e-02   4.263  2.2e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7102 on 998 degrees of freedom
Multiple R-squared:  0.01789,   Adjusted R-squared:  0.0169 
F-statistic: 18.18 on 1 and 998 DF,  p-value: 2.205e-05

summary(lm(SumIndex2      ~ treat, data = df))


Call:
lm(formula = SumIndex2 ~ treat, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-26.177  -6.171  -0.217   5.993  32.724 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.2042     0.4066   0.502    0.616
treat        -0.1742     0.5738  -0.304    0.761

Residual standard error: 9.073 on 998 degrees of freedom
Multiple R-squared:  9.235e-05, Adjusted R-squared:  -0.0009096 
F-statistic: 0.09217 on 1 and 998 DF,  p-value: 0.7615

summary(lm(PCIndex2       ~ treat, data = df))


Call:
lm(formula = PCIndex2 ~ treat, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9970  -2.0022   0.0525   2.0935   9.0686 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.1032     0.1341  -0.769    0.442
treat         0.2055     0.1893   1.086    0.278

Residual standard error: 2.993 on 998 degrees of freedom
Multiple R-squared:  0.00118,   Adjusted R-squared:  0.0001788 
F-statistic: 1.179 on 1 and 998 DF,  p-value: 0.2779

In this case we see that our conclusion is very different. Now the Anderson index strongly suggests some statistically significant treatment effect, while the principal components based index suggests little evidence of a statistically significant effect. Of course this should not surprise us if we consider what the indexes do “under the hood”: Anderson’s index overweights variables which bring independent variation, while principal components seeks to pick up underlying signals, and hence will downweight variables if they are quite unrelated to others within the class.

To see that this is not a statisical curio with this specific simulation, below we can generalise this, running \(S=500\) simulations. In each simulation we will consider these three types of indexes: the Anderson index, principal components, and a summary index. We will also consider 4 specific cases. A first case where there are no significant treatment effects for any of the variables, a second case where there is one significant treatment effect from a correlated variable, a third case where there is one significant treatment effect from an uncorrelated variable, and a fourth case where there is a significant treatment effect on all of the 10 variables considered. We will then examine the distribution of t-statistics from each of the simulations and each of the index types. We set up these simulations below:

set.seed(1213)
S <- 500
n <- 1000
r <- 0.9
library(MASS)
corr <- matrix(r, 10, 10)
diag(corr) <- 1
results <- data.frame(matrix(NA_real_, nrow=S, ncol=12))
colnames(results) <- c(
  "AndersonA_t","AndersonB_t","AndersonC_t","AndersonD_t",
  "PCA_t","PCB_t","PCC_t","PCD_t",
  "SumA_t","SumB_t","SumC_t","SumD_t"
)
for(s in 1:S){
    u <- mvrnorm(n, rep(0,10), corr)
    df <- as.data.frame(u)
    names(df) <- paste0("u",1:10)
    df$treat <- rbinom(n,1,0.5)
    for(i in 1:10){
        df[[paste0("y",i)]]  <- 1 + 0*  df$treat + df[[paste0("u",i)]]
        df[[paste0("yr",i)]] <- 1 + 0.5*df$treat + df[[paste0("u",i)]]
    }
    df$yr11 <- 1 + 0.5*df$treat + rnorm(n)
    idxC <- df$treat==0

    #define variable groups
    ys   <- paste0("y",1:10)
    Bvars <- c(paste0("y",1:9),"yr10")
    Cvars <- c(paste0("y",1:9),"yr11")
    Dvars <- paste0("yr",1:10)

    #define Anderson Index
    df$AndersonIndexA <- anderson_index(df, ys, idxC)
    df$AndersonIndexB <- anderson_index(df, Bvars, idxC)
    df$AndersonIndexC <- anderson_index(df, Cvars, idxC)
    df$AndersonIndexD <- anderson_index(df, Dvars, idxC)

    #define Principal Components
    df$PCIndexA  <- prcomp(df[ys], center=TRUE)$x[,1]
    df$PCIndexB  <- prcomp(df[Bvars], center=TRUE)$x[,1]
    df$PCIndexC  <- prcomp(df[Cvars], center=TRUE)$x[,1]
    df$PCIndexD  <- prcomp(df[Dvars], center=TRUE)$x[,1]

    #define Summary index
    df$SumIndexA <- rowSums(df[ys])
    df$SumIndexB <- rowSums(df[Bvars])
    df$SumIndexC <- rowSums(df[Cvars])
    df$SumIndexD <- rowSums(df[Dvars])

    tvals <- numeric(12)
    m <- lm(AndersonIndexA ~ treat, df);  tvals[1] <- summary(m)$coefficients["treat","t value"]
    m <- lm(AndersonIndexB ~ treat, df);  tvals[2] <- summary(m)$coefficients["treat","t value"]
    m <- lm(AndersonIndexC ~ treat, df);  tvals[3] <- summary(m)$coefficients["treat","t value"]
    m <- lm(AndersonIndexD ~ treat, df);  tvals[4] <- summary(m)$coefficients["treat","t value"]
    m <- lm(PCIndexA       ~ treat, df);  tvals[5] <- summary(m)$coefficients["treat","t value"]
    m <- lm(PCIndexB       ~ treat, df);  tvals[6] <- summary(m)$coefficients["treat","t value"]
    m <- lm(PCIndexC       ~ treat, df);  tvals[7] <- summary(m)$coefficients["treat","t value"]
    m <- lm(PCIndexD       ~ treat, df);  tvals[8] <- summary(m)$coefficients["treat","t value"]
    m <- lm(SumIndexA      ~ treat, df);  tvals[9] <- summary(m)$coefficients["treat","t value"]
    m <- lm(SumIndexB      ~ treat, df); tvals[10] <- summary(m)$coefficients["treat","t value"]
    m <- lm(SumIndexC      ~ treat, df); tvals[11] <- summary(m)$coefficients["treat","t value"]
    m <- lm(SumIndexD      ~ treat, df); tvals[12] <- summary(m)$coefficients["treat","t value"]
    results[s,] <- tvals
}

This results in a series of variables which contain the t-statistics of resulting tests, and we can plot these to see how each index performs. We do this below, starting with the Anderson Index:

xlims <- c(-3, 11)
xbreaks <- seq(-3, 11, by = 1)

# Variables to plot
vars_Anderson <- c("AndersonA_t","AndersonB_t","AndersonC_t","AndersonD_t")

for (var in vars_Anderson) {
    ggplot(results, aes_string(x = var)) +
        geom_histogram(aes(y = ..density..), bins = 30, fill = "grey", color = "black") +
        stat_function(fun = dnorm, args = list(mean = 0, sd = 1), color = "red") +
        scale_x_continuous(limits = xlims, breaks = xbreaks) +
        labs(
            x = expression(italic(t) ~ "statistic"),
            y = "Density"
        ) +
        theme_minimal() +
        theme(legend.position = "none")
}

In this case we see that – as we may hope – the distribution of t-statistics we observe in panel (a), based on pure null effects, lines up nearly perfectly with a normal distribution centred around 0 (ie consistent with zero effects). This is clearly different in panel (b) where all non-zero effects are considered, with t-statistics based on an Anderson Index all suggestive of non-zero effects. Panel (c), based on a single true effect which is highly correlated with all zero effects suggests relatively little evidence to reject null hypotheses of zero effects. However, panel (d), which is based on the same structure as panel (c) but now with an uncorrelated effect shows clear evidence against the null of zero effects.

Below we compare this to the case where we consider a principal components analysis. We build identical graphs based on principal components instead of the Anderson Index:

# Variables to plot
vars_pc <- c("PCA_t", "PCB_t", "PCC_t", "PCD_t")

for (var in vars_pc) {
    ggplot(results, aes_string(x = var)) +
        geom_histogram(aes(y = ..density..), bins = 30, fill = "grey", color = "black") +
        stat_function(fun = dnorm, args = list(mean = 0, sd = 1), color = "red") +
        scale_x_continuous(limits = xlims, breaks = xbreaks) +
        labs(
            x = expression(italic(t) ~ "statistic"),
            y = "Density"
        ) +
        theme_minimal() +
        theme(legend.position = "none")
}

Now we see quite different patterns in panel (c) and (d). While panels (a) and (b) are virtually identical, in cases where only certain effects are observed, our conclusions appear to be quite different. While there is some evidence that t-statistics are slightly shifted in the case of a single correlated effect, we see no evidence at all to reject the null where our single effect is independent of other variables, as documented in panel (d).

If we consider the simply summary index we observe a case between the Anderson Index and the Principal Components index. In this case the t-statistics in both panels (c) and (d) appear to be slightly shifted, though this is not as extreme as the cases above which (respectively) mechanically give more weight to the uncorrelated non-zero effect and correlated non-zero effect.

# Variables to plot
vars_sum <- c("SumA_t", "SumB_t", "SumC_t", "SumD_t")

for (var in vars_sum) {
    ggplot(results, aes_string(x = var)) +
        geom_histogram(aes(y = ..density..), bins = 30, fill = "grey", color = "black") +
        stat_function(fun = dnorm, args = list(mean = 0, sd = 1), color = "red") +
        scale_x_continuous(limits = c(-3, 11), breaks = seq(-3, 11, by = 1)) +
        labs(x = expression(italic(t) - " statistic"), y = "Density") +
        theme_minimal() +
        theme(legend.position = "none")
}

These patterns illustrate the importance of considering what type of index should be used if seeking to reduce multiple hypothesis tests to a single dimension. These are of course extreme cases based on a single non-zero effect among a class of outcomes (or cases where all or no significant treatment effects exist). You may wish to further explore these simulations by varying correlation structures, the number of significant effects, effect sizes and so forth. What is clear even from these simple structures however is that the choice of indexes is not innocuous, and we should take into account the nature of the effects being tested, whether we believe it is correct to overweight certain types of variables in our analysis, and where variation in effects comes from when determining how to aggregate underlying variables.

References

Alan, Sule, and Elif Kubilay. 2025. “Empowering Adolescents to Transform Schools: Lessons from a Behavioral Targeting.” American Economic Review 115 (2): 365–407. https://doi.org/10.1257/aer.20240374.

Anderson, Michael L. 2008. “Multiple Inference and Gender Differences in the Effects of Early Intervention: A Reevaluation of the Abecedarian, Perry Preschool, and Early Training Projects.” Journal of the American Statistical Association 103 (484): 1481–95.

Charness, Gary, and Uri Gneezy. 2009. “Incentives to Exercise.” Econometrica 77 (3): 909–31. https://doi.org/https://doi.org/10.3982/ECTA7416.