1.1.1 Experiments

• Experiments are investigations in which an intervention, in all its essential elements, is under the control of the investigator. (Cox & Reid)

• Two major types of control:

1. control over assignment to treatment – this is at the heart of many field experiments
2. control over the treatment itself – this is at the heart of many lab experiments

• Main focus today is on 1 and on the question: how does control over assignment to treatment allow you to make reasonable statements about causal effects?

1.1.3 Basic randomization

• Basic randomization is very simple. For example, say you want to assign 5 of 10 units to treatment. Here is simple code:

sample(0:1, 10, replace = TRUE)

 [1] 0 1 0 1 1 1 1 1 1 0

1.1.5 Basic randomization

Even better is to set it up so that it can reproduce lots of possible draws so that you can check the propensities for each unit.

set.seed(20111112)
P <- replicate(1000, sample(0:1, 10, replace = TRUE)) 
apply(P, 1, mean)

 [1] 0.519 0.496 0.510 0.491 0.524 0.514 0.535 0.497 0.470 0.506

Here the \(P\) matrix gives 1000 possible ways of allocating 5 of 10 units to treatment. We can then confirm that the average propensity is 0.5.

• A huge advantage of this approach is that if you make a mess of the random assignment; you can still generate the P matrix and use that for all analyses!

1.2.2 Cluster Randomization

• Often used if intervention has to function at the cluster level or if outcome defined at the cluster level.

• Disadvantage: loss of statistical power

• However: perfectly possible to assign some treatments at cluster level and then other treatments at the individual level

• Principle: (unless you are worried about spillovers) generally make clusters as small as possible

• Principle: Surprisingly, variability in cluster size makes analysis harder. Try to control assignment so that cluster sizes are similar in treatment and in control.

• Be clear about whether you believe effects are operating at the cluster level or at the individual level. This matters for power calculations.

• Be clear about whether spillover effects operate only within clusters or also across them. If within only you might be able to interpret treatment as the effect of being in a treated cluster…

1.2.3 Cluster Randomization: Block by cluster size

Surprisingly, if clusters are of different sizes the difference in means estimator is not unbiased, even if all units are assigned to treatment with the same probability.

Here’s the intuition.Say there are two clusters each with homogeneous treatment effects:

Then: What is the true average treatment effect? What do you expect to estimate from cluster random assignment?

The solution is to block by cluster size. For more see: http://gking.harvard.edu/files/cluster.pdf

1.3.1 Blocking

There are more or less efficient ways to randomize.

• Randomization helps ensure good balance on all covariates (observed and unobserved) in expectation.
• But balance may not be so great in realization
• Blocking can help ensure balance ex post on observables

1.3.2 Blocking

Consider a case with four units and two strata. There are 6 possible assignments of 2 units to treatment:

Even with a constant treatment effect and everything uniform within blocks, there is variance in the estimation of \(\widehat{\tau}\). This can be eliminated by excluding R1 and R6.

1.3.5 Blocking

• Blocking is a case of restricted randomization. Although each unit is sampled with equal probability, the profiles of possible assignments are not.
• You have to take account of this when doing analysis
• There are many other approaches.
  • Matched Pairs are a particularly fine approach to blocking
  • You could also randomize and then replace the randomization if you do not like the balance. This sounds tricky (and it is) but it is OK as long as you understand the true lottery process you are employing and incorporate that into analysis
  • It is even possible to block on covariates for which you don’t have data ex ante, by using methods in which you allocate treatment over time as a function of features of your sample (also tricky)

1.3.6 Other types of restricted randomization

• Really you can set whatever criterion you want for your set of treated units to have (eg no treated unit beside another treated unit; at least 5 from the north, 10 from the south, guaranteed balance by some continuous variable etc)
• You just have to be sure that you understand the random process that was used and that you can use it in the analysis stage
• But here be dragons
  • The more complex your design, the more complex your analysis.
  • General injunction http://www.ncbi.nlm.nih.gov/pubmed/15580598 ``as ye randomize so shall ye analyze’’)
  • In general you should make sure that a given randomization procedure coupled with a given estimation procedure will produce an unbiased estimate. DeclareDesign can help with this.

1.4.2 Factorial Designs

• Often when you set up an experiment you want to look at more than one treatment.
• Should you do this or not? How should you use your power?

Load up:

Spread out:

1.4.3 Factorial Designs

• Often when you set up an experiment you want to look at more than one treatment.
• Should you do this or not? How should you use your power?

Three arm it?:

Bunch it?:

1.4.4 Factorial Designs

• Surprisingly, adding multiple treatments does not generally eat into your power (unless you are decomposing a complex treatment – then it can. Why?)
• Especially when you use a fully crossed design like the middle one above.
• Fisher: “No aphorism is more frequently repeated in connection with field trials, than that we must ask Nature few questions, or, ideally, one question, at a time. The writer is convinced that this view is wholly mistaken.”
• However – adding multiple treatments does alter the interpretation of your average treatment effects. If T2 is an unusual treatment for example, then half the T1 effect is measured for unusual situations.

This speaks to “spreading out.” Note: the “bunching” example may not pay off and has undesireable feature of introducing a correlation between treatment assignments.

1.4.5 Factorial Designs

Two ways to do favtial assignments in DeclareDesign:

# Block the second assignment
declare_assignment(Z1 = complete_ra(N)) +
declare_assignment(Z2 = block_ra(blocks = Z1)) +
  
# Recode four arms  
declare_assignment(Z = complete_ra(N, num_arms = 4)) +
declare_measurement(Z1 = (Z == "T2" | Z == "T4"),
                      Z2 = (Z == "T3" | Z == "T4"))

1.4.6 Factorial Designs: In practice

• In practice if you have a lot of treatments it can be hard to do full factorial designs – there may be too many combinations.

• In such cases people use fractional factorial designs, like the one below (5 treatments but only 8 units!)

1.4.7 Factorial Designs: In practice

• Then randomly assign units to rows. Note columns might also be blocking covariates.

• In R, look at library(survey)

1.4.8 Factorial Designs: In practice

• But be careful: you have to be comfortable with possibly not having any simple counterfactual unit for any unit (invoke sparsity-of-effects principle).

• In R, look at library(survey)

1.4.9 Controversy?

Muralidharan, Romero, and Wüthrich (2023) write:

Factorial designs are widely used to study multiple treatments in one experiment. While t-tests using a fully-saturated “long” model provide valid inferences, “short” model t-tests (that ignore interactions) yield higher power if interactions are zero, but incorrect inferences otherwise. Of 27 factorial experiments published in top-5 journals (2007–2017), 19 use the short model. After including interactions, over half of their results lose significance. […]

1.5.1 Principle: Address external validity at the design stage

Anything to be done on randomization to address external validity concerns?

• Note 1: There is little or nothing about field experiments that makes the external validity problem greater for these than for other forms of ‘’sample based’’ research
• Note 2: Studies that use up the available universe (cross national studies) actually have a distinct external validity problem
• Two ways to think about external validity issues:
  1. Are things likely to operate in other units like they operate in these units? (even with the same intervention)
  2. Are the processes in operation in this treatment likely to operate in other treatments? (even in this population)

1.5.2 Principle: Address external validity at the design stage

• Two ways to think about external validity issues:
  1. Are things likely to operate in other units like they operate in these units? (even with the same intervention) 2.Are the processes in operation in this treatment likely to operate in other treatments? (even in this population)
• Two approaches for 1.
  • Try to sample cases and estimate population average treatment effects
  • Exploit internal variation: block on features that make the case unusal and assess importance of these (eg is unit poor? assess how effects differ in poor and wealthy components)
• 2 is harder and requires a sharp identification of context free primitives, if there are such things.

1.6.2 Sample data

Here are the first couple of rows and columns of the resulting data frame.

my_data <- draw_data(design)
kable(head(my_data), digits = 2)

1.6.3 Sample data

Here is the distribution between treatment and control:

kable(t(as.matrix(table(my_data$Z))), 
      col.names = c("control", "treatment"))

1.6.4 Complete Random Assignment using the built in function

assignment_complete <-   declare_assignment(Z = complete_ra(N))  

design_complete <- 
  replace_step(design, "assignment", assignment_complete)

1.6.5 Data from complete assignment

We can draw a new set of data and look at the number of subjects in the treatment and control groups.

set.seed(1:5)
data_complete <- draw_data(design_complete)

kable(t(as.matrix(table(data_complete$Z))))

1.6.6 Plotted

1.6.14 Illustration of efficiency gains from blocking

designs <- 
  list(
    simple = design, 
    complete = design_complete, 
    blocked = design_blocked, 
    clustered = design_clustered,  
    clustered_blocked = design_clustered_blocked) 

diagnoses <- diagnose_design(designs)

1.6.15 Illustration of efficiency gains from blocking

1.6.17 Nasty integer issues

• In many designs you seek to assign an integer number of subjects to treatment from some set.

• Sometimes however your assignment targets are not integers.

Example:

• I have 12 subjects in four blocks of 3 and I want to assign each subject to treatment with a 50% probability.

Two strategies:

1. I randomly set a target of either 1 or 2 for each block and then do complete assignment in each block. This can result in the numbers treated varying from 4 to 8
2. I randomly assign a target of 1 for two blocks and 2 for the other two blocks: Intuition–set a floor for the minimal target and then distribue the residual probability across blocks

1.6.18 Nasty integer issues

# remotes::install_github("macartan/probra")
library(probra)
set.seed(1)

blocks <- rep(1:4, each = 3)

table(blocks, prob_ra(blocks = blocks))

      
blocks 0 1
     1 1 2
     2 1 2
     3 2 1
     4 2 1

table(blocks, block_ra(blocks = blocks))

      
blocks 0 1
     1 1 2
     2 2 1
     3 1 2
     4 1 2

1.6.19 Nasty integer issues

Can also be used to set targets

# remotes::install_github("macartan/probra")
library(probra)
set.seed(1)

fabricate(N = 4,  size = c(47, 53, 87, 25), n_treated = prob_ra(.5*size)) %>%
  janitor::adorn_totals("row") |> 
  kable(caption = "Setting targets to get 50% targets with minimal variance")

1.6.20 Nasty integer issues

Can also be used to set for complete assignment with heterogeneous propensities

set.seed(1)

df <- fabricate(N = 100,  p = seq(.1, .9, length = 100), Z = prob_ra(p)) 

mean(df$Z)

[1] 0.5

df |> ggplot(aes(p, Z)) + geom_point() + theme_bw()

1.7.1 Indirect assignments

Indirect assignments are generally generated by applying a direct assignment and then figuring our an implied indirect assignment

set.seed(1)

df <-
  fabricate(
    N = 100, 
    latitude = runif(N),
    longitude = runif(N))

adjacency <- 
  sapply(1:nrow(df), function(i) 
    1*((df$latitude[i] - df$latitude)^2 + (df$longitude[i] - df$longitude)^2)^.5 < .1)

diag(adjacency) <- 0

1.7.2 Indirect assignments: Adjacency matrix

adjacency |>  
  reshape2::melt(c("x", "y"), value.name = "close") |> mutate(close = factor(close)) |>
  ggplot(aes(x=x,y=y,fill=close))+
  geom_tile() + xlab("individual") + ylab("individual") + theme_bw() +
  scale_fill_grey(start = 1, end = 0)  # 1 = white, 0 = black

1.7.3 Indirect assignments

n_assigned <- 50

design <-
  declare_model(data = df) + 
  declare_assignment(
    direct = complete_ra(N, m = n_assigned),
    indirect = 1*(as.vector(as.vector(direct) %*% adjacency) >= 1))

draw_data(design) |> with(table(direct, indirect))

      indirect
direct  0  1
     0 13 37
     1 13 37

2.2.2 Weighing Evidence

When we test a hypothesis we decide first on what sort of evidence we need to see in order to decide that the hypothesis is not reliable.

• Othello has a hypothesis that Desdemona is innocent.

• Iago confronts him with evidence:

  • See how she looks at him: would she look a him like that if she were innocent?
  • … would she defend him like that if she were innocent?
  • … would he have her handkerchief if she were innocent?
  • Othello, the chances of all of these things arising if she were innocent is surely less than 5%

2.2.3 Hypotheses are often rejected, sometimes maintained, but rarely accepted

• Note that Othello is focused on the probability of the events if she were innocent but not the probability of the events if Iago were trying to trick him.

• He is not assessing his belief in whether she is faithful, but rather how likely the data would be if she were faithful.

So:

• He assesses: \(\Pr(\text{Data} | \text{Hypothesis is TRUE})\)
• While a Bayesian would assess: \(\Pr(\text{Hypothesis is TRUE} | \text{Data})\)

2.2.4 Recap: Calculate a \(p\) value in your head

• Illustrating \(p\) values via “randomization inference”

• Say you randomized assignment to treatment and your data looked like this.

Then:

• Does the treatment improve your health?
• What’s the \(p\) value for the null that treatment had no effect on anybody?

2.4.6 Step 2: What is the power?

• Now we need to stipulate some belief about how the world really works—this is not the null hypothesis that we plan to reject, but something that we actually take to be true.

• For instance: we think that in fact sixes appear 20% of the time.

Now what’s the probability of seeing at least 190 sixes?

1 - pbinom(189, 1000, .2)

[1] 0.796066

So given I think 6s appear 20% of the time, I think it likely I’ll see at least 190 sixes and reject the hypothesis of a fair dice.

2.4.9 Subtleties

• Is a significant result from an underpowered study less credible? (only if there is a significance filter)
• What significance level should you choose for power? (Obviously the stricter the level the lower the power, so use what you will use when you actually implement tests)
• Do you really have to know the effect size to do power analysis? (No, but you should know at least what effects sizes you would want to be sure about picking up if they were present)
• Power is just one of many possible diagnosands
• What’s power for Bayesians?

2.4.10 Power analytics

Simplest intuition on power:

What is the probability of getting a significant estimate given the sampling distribution is centered on \(b\) and the standard error is 1?

• Probability below -1.96: \(F(-1.96 | \tau))\)
• Probability above -1.96: \(1-F(1.96 | \tau)\)

Add these together: probability of getting an estimate above 1.96 or below -1.96.

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

  1 - pnorm(critical, mean = abs(b)) + pnorm(-critical, mean = abs(b))

power(0)

[1] 0.05

power(1.96)

[1] 0.5000586

power(-1.96)

[1] 0.5000586

power(3)

[1] 0.8508388

2.4.15 Power analytics for a trial: by hand

• Of course the standard error will depend on the number of units and the variance of outcomes in treatment and control.

• Say \(N\) subject are divided into two groups and potential outcomes have standard deviation \(\sigma\) in treatment and control. Then the conservative variance of the treatment effect is (approx / conservatively):

\[Var(\tau)=\frac{\sigma^2}{N/2} + \frac{\sigma^2}{N/2} = 4\frac{\sigma^2}{N}\]

and so the (conservative / approx) standard error is:

\[\sigma_\tau=\frac{2\sigma}{\sqrt{N}}\]

Note here we seem to be using the actual standard error but of course the tests we actually run will use an estimate of the standard error…

2.4.16 Power analytics for a trial: by hand

se <- function(sd, N) (N/(N-1))^.5*2*sd/(N^.5)


power_2 <- function(b, alpha = .05, sd = 1, N = 100, critical = qnorm(1-alpha/2), se = 2*sd/N^.5)  

  1 - pnorm(critical, mean = abs(b)/se(sd, N)) + pnorm(-critical, mean = abs(b)/se(sd, N))

power_2(0)

[1] 0.05

power_2(.5)

[1] 0.7010827

2.4.25 Design diagnosis does it all (over multiple designs)

  diagnose_design(design)

2.4.26 Design diagnosis does it all

design |>
  redesign(b = c(0, 0.25, 0.5, 1)) |>
  diagnose_design()

2.5.6 Power when estimates are biased

bad_design <- 
  
  declare_model(N = 100, 
    U = rnorm(N),
    potential_outcomes(Y ~ 0 * X + U, conditions = list(X = 0:1)),
    X = ifelse(U > 0, 1, 0)) + 
  
  declare_measurement(Y = reveal_outcomes(Y ~ X)) + 
  
  declare_inquiry(ate = mean(Y_X_1 - Y_X_0)) + 
  
  declare_estimator(Y ~ X, inquiry = "ate", .method = lm_robust)

2.5.7 Power when estimates are biased

You can see from the null design that power is great but bias is terrible and coverage is way off.

diagnose_design(bad_design)

Power without unbiasedness corrupts, absolutely

2.5.8 Power with a more subtly biased experimental design

another_bad_design <- 
  
  declare_model(
    N = 100, 
    female = rep(0:1, N/2),
    U = rnorm(N),
    potential_outcomes(Y ~ female * Z + U)) + 
  
  declare_assignment(
    Z = block_ra(blocks = female, block_prob = c(.1, .5)),
    Y = reveal_outcomes(Y ~ Z)) + 

  declare_inquiry(ate = mean(Y_Z_1 - Y_Z_0)) + 
  
  declare_estimator(Y ~ Z + female, inquiry = "ate", 
                    .method = lm_robust)

  diagnose_design(another_bad_design)

2.5.9 Power with a more subtly biased experimental design

You can see from the null design that power is great but bias is terrible and coverage is way off.

2.5.10 Power with the wrong standard errors

clustered_design <-
  declare_model(
    cluster = add_level(N = 10, cluster_shock = rnorm(N)),
    individual = add_level(
        N = 100,
        Y_Z_0 = rnorm(N) + cluster_shock,
        Y_Z_1 = rnorm(N) + cluster_shock)) +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_assignment(Z = cluster_ra(clusters = cluster)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z, inquiry = "ATE")

What alerts you to a problem?

2.5.13 Multiple comparisons correction (complex code)

Say I run two tests and want to correct for multiple comparisons.

Two approaches. First, by hand:

b = .2

design_mc <-
  declare_model(N = 1000, 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 = b) +
  declare_estimator(Y ~ factor(Z), term = TRUE)

2.5.14 Multiple comparisons correction (complex code)

design_mc |>
  simulate_designs(sims = 1000) |>
  filter(term != "(Intercept)") |>
  group_by(sim_ID) |>
  mutate(p_bonferroni = p.adjust(p = p.value, method = "bonferroni"),
         p_holm = p.adjust(p = p.value, method = "holm"),
         p_fdr = p.adjust(p = p.value, method = "fdr")) |>
  ungroup() |>
  summarize(
    "Power using naive p-values" = mean(p.value <= 0.05),
    "Power using Bonferroni correction" = mean(p_bonferroni <= 0.05),
    "Power using Holm correction" = mean(p_holm <= 0.05),
    "Power using FDR correction" = mean(p_fdr <= 0.05)
    ) 

2.5.15 Multiple comparisons correction (approach 2)

The alternative approach (generally better!) is to design with a custom estimator that includes your corrections.

my_estimator <- function(data) 
  lm_robust(Y ~ factor(Z), data = data) |> 
  tidy() |>
  filter(term != "(Intercept)") |>
  mutate(p.naive = p.value,
         p.value = p.adjust(p = p.naive, method = "bonferroni"))
  

design_mc_2 <- design_mc |>
  replace_step(5, declare_estimator(handler = label_estimator(my_estimator))) 

run_design(design_mc_2) |> 
  select(term, estimate, p.value, p.naive) |> kable()

2.5.16 Multiple comparisons correction (Null model case)

Lets try same thing for a null model (using redesign(design_mc_2, b = 0))

design_mc_3 <- 
  design_mc_2 |> 
  redesign(b = 0) 

run_design(design_mc_3) |> select(estimate, p.value, p.naive) |> kable(digits = 3)

2.5.17 Multiple comparisons correction (Null model case)

…and power:

bothered?

2.5.19 Big takeaways

• Power is affected not just by sample size, variability and effect size but also by you data and analysis strategies.
• Try to estimate power under multiple scenarios
• Try to use the same code for calculating power as you will use in your ultimate analysis
• Basically the same procedure can be used for any design. If you can declare a design and have a test, you can calculate power
• Your power might be right but misleading. For confidence:
  • Don’t just check power, check bias and coverage also
  • Check power especially under the null
• Don’t let a focus on power distract you from more substantive diagnosands