1.2.1 Review

In the classical approach to testing a hypothesis we ask:

How likely are we to see data like this if indeed the hypothesis is true?

• If the answer is “not very likely” then we treat the hypothesis as suspect.
• If the answer is not “not very likely” then the hypothesis is maintained (some say “accepted” but this is tricky as you may want to “maintain” multiple incompatible hypotheses)

How unlikely is “not very likely”?

1.4.1 By hand

I want to test the hypothesis that a six never comes up on this dice.

Here’s my test:

• I will roll the dice once.
• If a six comes up I will reject the hypothesis.

What is the power of this test?

1.4.2 By hand

I want to test the hypothesis that a six never comes up on this dice.

Here’s my test:

• I will roll the dice twice.
• If a six comes up either time I will reject the hypothesis.

What is the power of this test?

1.4.3 Two probabilities

Power sometimes seems more complicated because hypothesis rejection involves a calculated probability and so you need the probability of a probability.

I want to test the hypothesis that this dice is fair.

Here’s my test:

• I will roll the dice 1000 times and if I see fewer than x 6s or more than y 6s I will reject the hypothesis.

Now:

• What should x and y be?
• What is the power of this test?

1.4.4 Step 1: When do you reject?

For this we need to figure a rule for rejection. This is based on identifying events that should be unlikely under the hypothesis.

Here is how many 6’s I would expect if the dice is fair:

fabricate(N = 1001, sixes = 0:1000, p = dbinom(sixes, 1000, 1/6)) |>
  ggplot(aes(sixes, p)) + geom_line()

1.4.5 Step 1: When do you reject?

I can figure out from this that 143 or fewer is really very few and 190 or more is really very many:

c(lower = pbinom(143, 1000, 1/6), upper = 1 - pbinom(189, 1000, 1/6))

     lower      upper 
0.02302647 0.02785689 

1.4.7 Rule of thumb

• 80% or 90% is a common rule of thumb for “sufficient” power
• but really, how much power you need depends on the purpose

1.4.8 Think about

• Are there other tests I could have implemented?
• Are there other ways to improve this test?

1.4.11 Power analytics: graphed

This is essentially what is done by pwrss::power.z.test – and it produces nice graphs!

See:

pwrss::power.z.test(ncp = 1.96, alpha = 0.05, alternative = "not equal", plot = TRUE)

     power ncp.alt ncp.null alpha  z.crit.1 z.crit.2
 0.5000586    1.96        0  0.05 -1.959964 1.959964

1.4.12 Power analytics: graphed

Substantively: if in expectation an estimate will be just significant, then your power is 50%

1.4.13 Equivalent

power <- function(b, alpha = 0.05, critical = qnorm(1-alpha/2))  

  1 - pnorm(critical - b) + pnorm(-critical - b)

power(1.96)

[1] 0.5000586

Intuition:

x <- seq(-3, 3, .01)
plot(x, dnorm(x), main = "power associated with effect of 1 se")
abline(v = 1.96 - 1)
abline(v = -1.96 - 1)

1.4.16 Power analytics for a trial: flexible

This can be done e.g. with pwrss like this:

pwrss::pwrss.t.2means(mu1 = .2, mu2 = .1, sd1 = 1, sd2 = 1, 
               n2 = 50, alpha = 0.05,
               alternative = "not equal")

 Difference between Two means 
 (Independent Samples t Test) 
 H0: mu1 = mu2 
 HA: mu1 != mu2 
 ------------------------------ 
  Statistical power = 0.079 
  n1 = 50 
  n2 = 50 
 ------------------------------ 
 Alternative = "not equal" 
 Degrees of freedom = 98 
 Non-centrality parameter = 0.5 
 Type I error rate = 0.05 
 Type II error rate = 0.921 

power_2(.50, N = 100)

[1] 0.7010827

1.4.17 Power for more complex trials: Analytics

Mostly involve figuring out the standard error.

Consider a cluster randomized trial, with each unit having a cluster level shock \(\epsilon_k\) and an individual shock \(\nu_i\). Say these have variances \(\sigma^2_k\), \(\sigma^2_i\).

The standard error is:

\[\sqrt{\frac{4\sigma^2_k}{K} + \frac{4\sigma^2_i}{nK}}\]

Define \(\rho = \frac{\sigma^2_k}{\sigma^2_k + \sigma^2_i}\)

\[\sqrt{\rho \frac{4\sigma^2}{K} + (1- \rho)\frac{4\sigma^2}{nK}}\]

\[\sqrt{((n - 1)\rho + 1)\frac{4\sigma^2}{nK}}\]

where

• \(((n - 1)\rho + 1)\) is the “design effect”
• \(\frac{nK}{((n - 1)\rho + 1)}\) is the “effective sample size”

Plug in these standard errors and proceed as before

1.4.18 Power via design diagnosis

Is arbitrarily flexible

N <- 100
b <- .5

design <- 
  declare_model(N = N, 
    U = rnorm(N),
    potential_outcomes(Y ~ b * Z + U)) + 
  declare_assignment(Z = simple_ra(N),
                     Y = reveal_outcomes(Y ~ Z)) + 
  declare_inquiry(ate = mean(Y_Z_1 - Y_Z_0)) + 
  declare_estimator(Y ~ Z, inquiry = "ate", .method = lm_robust)

1.4.19 Run it many times

sims_1 <- simulate_design(design) 

sims_1 |> select(sim_ID, estimate, p.value)

1.4.20 Power is mass of the sampling distribution of decisions under the model

sims_1 |>
  ggplot(aes(p.value)) + 
  geom_histogram() +
  geom_vline(xintercept = .05, color = "red")

1.4.21 Power is mass of the sampling distribution of decisions under the model

Obviously related to the estimates you might get

sims_1 |>
  mutate(significant = p.value <= .05) |>
  ggplot(aes(estimate, p.value, color = significant)) + 
  geom_point()

1.4.22 Check coverage is correct

sims_1 |>
  mutate(within = (b > sims_1$conf.low) & (b < sims_1$conf.high)) |> 
  pull(within) |> mean()

[1] 0.9573333

1.4.23 Check validity of \(p\) value

A valid \(p\)-value satisfies \(\Pr(p≤x)≤x\) for every \(x \in[0,1]\) (under the null)

sims_2 <- 
  
  redesign(design, b = 0) |>
  
  simulate_design() 

1.4.26 Diagnose over multiple moving parts (and ggplot)

design |>
  ## Redesign
  redesign(b = c(0.1, 0.3, 0.5), N = 100, 200, 300) |>
  ## Diagnosis
  diagnose_design() |>
  ## Prep
  tidy() |>
  filter(diagnosand == "power") |>
  ## Plot
  ggplot(aes(N, estimate, color = factor(b))) +
  geom_line()

1.4.27 Diagnose over multiple moving parts (and ggplot)

1.4.28 Diagnose over multiple moving parts and multiple diagnosands (and ggplot)

design |>

  ## Redesign
  redesign(b = c(0.1, 0.3, 0.5), N = 100, 200, 300) |>
  
  ## Diagnosis
  diagnose_design() |>
  
  ## Prep
  tidy() |>
  
  ## Plot
  ggplot(aes(N, estimate, color = factor(b))) +
  geom_line()+
  facet_wrap(~diagnosand)

1.4.29 Diagnose over multiple moving parts and multiple diagnosands (and ggplot)

1.5.1 Power tips

coming up:

• power everywhere
• power with bias
• power with the wrong standard errors
• power with uncertainty over effect sizes
• power and multiple comparisons

1.5.2 Power depends on all parts of MIDA

We often focus on sample sizes

But

Power also depends on

• the model – obviously signal to noise
• the assignments and specifics of sampling strategies
• estimation procedures

1.5.3 Power from a lag?

Say we have access to a “pre” measure of outcome Y_now; call it Y_base. Y_base is informative about potential outcomes. We are considering using Y_now - Y_base as the outcome instead of Y_now.

N <- 100
rho <- .5

design <- 
  declare_model(N,
                 Y_base = rnorm(N),
                 Y_Z_0 = 1 + correlate(rnorm, given = Y_base, rho = rho),
                 Y_Z_1 = correlate(rnorm, given = Y_base, rho = rho),
                 Z = complete_ra(N),
                 Y_now = Z*Y_Z_1 + (1-Z)*Y_Z_0,
                 Y_change = Y_now - Y_base) +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_estimator(Y_now ~ Z, label = "level") +
  declare_estimator(Y_change ~ Z, label = "change")+
  declare_estimator(Y_now ~ Z + Y_base, label = "RHS")

1.5.4 Power from a lag?

design |> redesign(N = c(10, 100, 1000, 10000), rho = c(.1, .5, .9)) |>
  diagnose_design() 

1.5.5 Power from a lag?

Punchline:

• if you difference: the lag has to be sufficiently information to pay its way (the \(\rho = .5\) equivalent between level and change follows from Gerber and Green (2012) equation 4.6)
• The right hand side is your friend, at least for experiments (Ding and Li (2019))
• As \(N\) grows the stakes fall

1.5.11 Let’s fix that one

clustered_design_2  <-
  clustered_design |> replace_step(5, 
  declare_estimator(Y ~ Z, clusters = cluster))

1.5.12 Power when you are not sure about effect sizes (always!)

• you can do power analysis for multiple stipulations
• or you can design with a distribution of effect sizes

design_uncertain <-
  declare_model(N = 1000, b = 1+rnorm(1), Y_Z_1 = rnorm(N), Y_Z_2 = rnorm(N) + b, Y_Z_3 = rnorm(N) + b) +
  declare_assignment(Z = complete_ra(N = N, num_arms = 3, conditions = 1:3)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_inquiry(ate = mean(b)) +
  declare_estimator(Y ~ factor(Z), term = TRUE)

draw_estimands(design_uncertain)

  inquiry   estimand
1     ate -0.3967765

draw_estimands(design_uncertain)

  inquiry  estimand
1     ate 0.7887188

1.5.18 You might try

• Power for an interaction (in a factorial design)
• Power for a binary variable (versus a continuous variable?)
• Power gains from blocked randomization
• Power losses from clustering at different levels
• Controlling the ICC directly? (see book cluster designs section)