2.2.1 Four elements of any research design

• Model: set of models of what causes what and how
• Inquiry: a question stated in terms of the model
• Data strategy: the set of procedures we use to gather information from the world (sampling, assignment, measurement)
• Answer strategy: how we summarize the data produced by the data strategy

2.2.2 Four elements of any research design

2.2.3 Declaration

Design declaration is telling the computer (and readers) what M, I, D, and A are.

2.2.4 Diagnosis

• Design diagnosis is figuring out how the design will perform under imagined conditions.

• Estimating “diagnosands” like power, bias, rmse, expected error rates, expected ethical harm, expected “amount learned”.

• Diagnosis takes account of model uncertainty: it aims to identify models for which the design works well and models for which it does not

2.2.5 Redesign

Redesign is the fine-tuning of features of the data- and answer strategies to understand how changing them affects the diagnosands

• Different sample sizes
• Different randomization procedures
• Different estimation strategies
• Implementation (e.g. gains from effort into compliance versus more effort into sample size)

2.2.6 Very often you have to simulate

• Doing all this is often too hard to work out from rules of thumb or power calculators
• Specialized formulas exist for some diagnosands, but not all

2.3.1 Key commands for design declaration

• declare_model()

• declare_inquiry()

• declare_sampling()

• declare_assignment()

• declare_measurement()

• declare_estimator()

and there are more declare_ functions!

2.3.2 Key commands for using a design

• draw_data(design)
• draw_estimands(design)
• draw_estimates(design)
• get_estimates(design, data)
• run_design(design), simulate_design(design)
• diagnose_design(design)
• redesign(design, N = 200)
• compare_designs(), compare_diagnoses()

2.3.3 Pipeable commands

design |> 
  redesign(N = c(200, 400)) |>
  diagnose_designs() |> 
  tidy() |> 
  ggplot(...) 

2.3.4 Cheat sheet

https://raw.githubusercontent.com/rstudio/cheatsheets/master/declaredesign.pdf

2.3.5 Other resources

• The website: https://declaredesign.org/
• The book: https://book.declaredesign.org
• The console: ?DeclareDesign

2.4.1 The simplest possible (diagnosable) design?

mean <- 0

simplest_design <- 
  declare_model(N = 100, Y = rnorm(N, mean)) +
  declare_inquiry(Q = mean) +
  declare_estimator(Y ~ 1)

• we draw 100 units from a standard normal distribution
• we define our inquiry as the population expectation
• we estimate the average using a regression with a constant term

2.4.2 The simplest possible design?

simplest_design <- 
  declare_model(N = 100, Y = rnorm(N, mean)) +
  declare_inquiry(Q = mean) +
  declare_estimator(Y ~ 1)

• This design has three steps, with steps connected by a +
• The design itself is just a list of steps and has class design

str(simplest_design)

List of 3
 $ model    :design_step:    declare_model(N = 100, Y = rnorm(N, mean)) 
 $ Q        :design_step:    declare_inquiry(Q = mean) 
 $ estimator:design_step:    declare_estimator(Y ~ 1) 
 - attr(*, "call")= language construct_design(steps = steps)
 - attr(*, "class")= chr [1:2] "design" "dd"

2.4.3 The design is a pipe

Each step is a function (or rather: a function that generates functions) and each function presupposes what is created by previous functions.

• The ordering of steps is quite important
• Most steps take the main data frame in and push the main dataframe out; this data frame normally builds up as you move along the pipe.

2.4.4 The design is a pipe

Each step is a function (or rather: a function that generates functions) and each function presupposes what is created by previous functions.

• The ordering of steps is quite important
• declare_estimator steps take the main data frame in and send out an estimator_df dataframe
• declare_inquiry steps take the main data frame in and send out an estimand_df dataframe.

2.4.5 The design is a pipe

• You can run these functions one at a time if you like.
• For instance the third step presupposes the data from the first step:

df <- simplest_design[[1]]()
A  <- simplest_design[[3]](df)

A |> kable(digits = 2) |> kable_styling(font_size = 20)

Estimand  <- simplest_design[[2]](df)

Estimand |> kable(digits = 2) |> kable_styling(font_size = 20)

2.4.6 Run it once

You can also just run through the whole design once by typing the name of the design:

simplest_design

Research design declaration summary

Step 1 (model): declare_model(N = 100, Y = rnorm(N, mean)) ---------------------

Step 2 (inquiry): declare_inquiry(Q = mean) ------------------------------------

Step 3 (estimator): declare_estimator(Y ~ 1) -----------------------------------

Run of the design:

 inquiry estimand estimator        term estimate std.error statistic p.value
       Q        0 estimator (Intercept)    0.155     0.108      1.43   0.155
 conf.low conf.high df outcome
  -0.0597     0.371 99       Y

2.4.7 Run it again

Or by asking for a run of the design

one_run <- simplest_design |> run_design()
one_run |> kable(digits = 2) |> kable_styling(font_size = 18)

A single run creates data, calculates estimands (the answer to inquiries) and calculates estimates plus ancillary statistics.

2.4.8 Simulation

Or by asking for many runs of the design

some_runs <- simplest_design |> simulate_design(sims = 1000)

some_runs |> head() |> kable(digits = 2) |> kable_styling(font_size = 16)

2.4.9 Diagnosis

Once you have simulated many times you can “diagnose”.

This is the next topic

2.5.1 Diagnosis by hand

Once you have simulated many times you can “diagnose”.

For instance we can ask about bias: the average difference between the estimand and the estimate:

some_runs |> 
  summarize(mean_estimate = mean(estimate), 
            mean_estimand = mean(estimand), 
            bias = mean(estimate - estimand)) 

2.5.2 diagnose_design()

diagnose_design() does this in one step for a set of common “diagnosands”:

diagnosis <-
  simplest_design |>
  diagnose_design()

2.5.3 What is the diagnosis object?

The diagnosis object is also a list; of class diagnosis

names(diagnosis)

[1] "simulations_df"       "diagnosands_df"       "diagnosand_names"    
[4] "group_by_set"         "parameters_df"        "bootstrap_replicates"
[7] "bootstrap_sims"       "duration"            

class(diagnosis)

[1] "diagnosis"

2.5.4 What is the diagnosis object?

diagnosis$simulations_df |> 
  head() 

2.5.5 What is the diagnosis object?

diagnosis$diagnosands_df |> 
  head() 

2.5.6 What is the diagnosis object?

diagnosis$bootstrap_replicates |> 
  head()

2.5.7 Diagnosis: Bootstraps

• The bootstraps dataframe is produced by resampling from the simulations dataframe and producing a diagnosis dataframe from each resampling.

• This lets us generate estimates of uncertainty around our diagnosands.

• It can be controlled thus:

diagnose_design(
  ...,
  bootstrap_sims = 100
)

2.5.8 After Diagnosis

It’s reshapeable: as a tidy dataframe, ready for graphing

diagnosis |> 
  tidy() 

2.5.9 After Diagnosis

It’s reshapeable: as a tidy dataframe, ready for graphing

diagnosis |> 
  tidy() |> 
  ggplot(aes(estimate, diagnosand)) + geom_point() + 
  geom_errorbarh(aes(xmax = conf.high, xmin = conf.low, height = .2))

2.5.10 After Diagnosis: Tables

Or turn into a formatted table:

diagnosis |> 
  reshape_diagnosis() 

2.5.11 Spotting design problems with diagnosis

Diagnosis alerts to problems in a design. Consider the following simple alternative design.

simplest_design_2 <- 
  
  declare_model(N = 100, Y = rnorm(N)) +
  declare_inquiry(Q = mean(Y)) +
  declare_estimator(Y ~ 1)

Here we define the inquiry as the sample average \(Y\) (instead of the population mean). But otherwise things stay the same.

What do we think of this design?

2.5.12 Spotting design problems with diagnosis

Here is the diagnosis

• Why is coverage so high? is that OK?
• Why is the RMSE 0? But in each run the std.error > 0? is that OK?
  • Is it because the RMSE is too low?
  • Or the standard error is too large?

2.5.13 It depends on the inquiry

• If we are really interested in the sample average then our standard error is off: we should have no error at all!
• If we are really interested in the population average then our inquiry is badly defined: it should not be redefined on each run!

2.5.14 Diagnosing multiple designs

You can diagnose multiple designs or a list of designs

list(dum = simplest_design, dee = simplest_design) |>
  diagnose_design(sims = 5) |>
  reshape_diagnosis() |> 
  kable() |> 
  kable_styling(font_size = 20)

2.6.1 What is redesign?

Redesign is the process of taking a design and modifying it in some way.

There are a few ways to do this:

1. Just make a new design using modified code
2. Take a design and alter some steps using replace_step, insert_step or delete_step
3. Modify a design parameter using redesign

we will focus on the third approach

2.6.3 Redesign parameter definition

N <- 100

simplest_design_N <- 
  
  declare_model(N = N, Y = rnorm(N)) +
  declare_inquiry(Q = 0) +
  declare_estimator(Y ~ 1)

Note that N is defined in memory; and it gets called in one of the steps. It has now become a parameter of the design and it can be modified using redesign.

2.6.4 Redesign illustration

Here is a version of the design with N = 200:

design_200 <- simplest_design_N |> redesign(N = 200)
  
design_200 |> draw_data() |> nrow()

[1] 200

2.6.5 Redesigning to a list

Here is a list of three different designs with different Ns.

design_Ns <- simplest_design_N |> redesign(N = c(200, 400, 800))

design_Ns |> lapply(draw_data) |> lapply(nrow)

$design_1
[1] 200

$design_2
[1] 400

$design_3
[1] 800

2.6.6 Redesigning to a list

The good thing here is that it is now easy to diagnose over multiple designs and compare diagnoses. The parameter names then end up in the diagnosis_df

Consider this:

N <- 100
m <- 0

design <- 
  declare_model(N = N, Y = rnorm(N, m)) +
  declare_inquiry(Q = m) +
  declare_estimator(Y ~ 1) 

Then:

designs <-  redesign(design, N = c(100, 200, 300), m = c(0, .1, .2))
  
designs |> diagnose_design() |> tidy() 

2.6.7 Redesigning to a list

Output:

designs |> diagnose_design() |> tidy() 

2.6.8 Graphing after redesigning to a list

Graphing after redesign is easy:

designs |> diagnose_design() |> 
  tidy() |>
  filter(diagnosand %in% c("power", "rmse")) |> 
  ggplot(aes(N, estimate, color = factor(m))) + 
  geom_line() + 
  facet_wrap(~diagnosand)

Power depends on N and m, rmse depends on N only

2.7.1 Using a design

We motivate with a slightly more complex experimental design (more on the components of this later)

b <-1
N <- 100
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")

2.7.2 Make data from the design

data <- draw_data(design)

data |> head () |> kable() |> kable_styling(font_size = 20)

2.7.3 Make data from the design

Play with the data:

lm_robust(Y ~ Z, data = data)

2.7.4 Draw estimands

draw_estimands(design) |>
  kable(digits = 2) |> 
  kable_styling(font_size = 20)

2.7.5 Draw estimates

draw_estimates(design) |> 
  kable(digits = 2) |> 
  kable_styling(font_size = 20)

2.7.6 Get estimates

Using your actual data:

get_estimates(design, data) |>
  kable(digits = 2) |> 
  kable_styling(font_size = 20)

2.7.7 Simulate design

simulate_design(design, sims = 3) |>
  kable(digits = 2) |> 
  kable_styling(font_size = 16)

2.7.8 Diagnose design

design |> 
  diagnose_design(sims = 100) 

2.7.9 Redesign

new_design <-
  
  design |> redesign(b = 0)

• Modify any arguments that are explicitly called on by design steps.
• Or add, remove, or replace steps

2.7.10 Compare designs

compare_diagnoses(design1 = design,
                  design2 = redesign(design, N = 50))

3.1.2 Steps: Order matters

e.g. If you sample before defining the inquiry you get a different inquiry to if you sample after you define the inquiry

design_1 <- 
  declare_model(N = 1000, X = rep(0:1, N/2), Y = X + rnorm(N)) + 
  declare_sampling(S= strata_rs(strata = X, strata_prob = c(.2, .8))) +
  declare_inquiry(m = mean(Y))

design_1 |> draw_estimands()

  inquiry  estimand
1       m 0.7606704

3.1.3 Steps: Order matters

e.g. If you sample before defining the inquiry you get a different inquiry to if you sample after you define the inquiry

design_2 <- 
  declare_model(N = 1000, X = rep(0:1, N/2), Y = X + rnorm(N)) + 
  declare_inquiry(m = mean(Y)) +
  declare_sampling(S= strata_rs(strata = X, strata_prob = c(.2, .8))) 

design_2 |> draw_estimands()

  inquiry  estimand
1       m 0.4957891

3.2.1 Hierarchical data

You can generate hierarchical data like this:

M <- 
  declare_model(
    households = add_level(
      N = 100, 
      N_members = sample(c(1, 2, 3, 4), N, 
                         prob = c(0.2, 0.3, 0.25, 0.25), replace = TRUE)
    ),
    individuals = add_level(
      N = N_members, 
      age = sample(18:90, N, replace = TRUE)
    )
  )

3.2.2 Hierarchical data

You can generate hierarchical data like this:

M() |> head() |> kable(digits = 2) |> kable_styling(font_size = 20)

3.2.4 Panel data

You can generate panel data like this:

M() |> head() |> kable(digits = 2) |> kable_styling(font_size = 20)

3.2.5 Preexisting data

M <- 
  declare_model(
    data = baseline_data,
    attitudes = sample(1:5, N, replace = TRUE)
  )

3.2.6 Steps

You can repeat steps and play with the order, always conscious of the direction of the pipe

design <- 
  declare_model(N = N, X = rep(0:1, N/2)) +
  declare_model(U = rnorm(N), potential_outcomes(Y ~ b * Z * X + U)) + 
  declare_assignment(Z = block_ra(blocks = X), Y = reveal_outcomes(Y ~ Z)) + 
  declare_inquiry(ate = mean(Y_Z_1 - Y_Z_0)) + 
  declare_inquiry(cate = mean(Y_Z_1[X==0] - Y_Z_0[X==0])) + 
  declare_estimator(Y ~ Z, inquiry = "ate", label = "ols") + 
  declare_estimator(Y ~ Z*X, inquiry = "cate", label = "fe")

3.2.7 You can generate multiple columns together

M2 <-
  declare_model(
    draw_multivariate(c(X1, X2) ~ MASS::mvrnorm(
      n = 1000,
      mu = c(0, 0),
      Sigma = matrix(c(1, 0.3, 0.3, 1), nrow = 2)
    )))

3.2.8 You can generate multiple columns together

M2() |> head() |> kable(digits = 2) |> kable_styling(font_size = 28) 

3.2.9 Cluster structures with cluster correlations

M <-
  declare_model(households = add_level(N = 1000),
                individuals = add_level(
                  N = 4,
                  X = draw_normal_icc(
                    mean = 0,
                    clusters = households,
                    ICC = 0.65
                  )
                ))

3.2.10 Cluster structures with cluster correlations

model <- lm_robust(X ~ households, data = M())
model$adj.r.squared

[1] 0.6215795

3.3.1 Inquiries using predefined potential outcomes

Many causal inquiries are simple summaries of potential outcomes:

Generating potential outcomes columns gets you far

3.3.2 Inquiries using functions

Often though we need to define inquiries as a function of continuous variables. For this generating a potential outcomes function can make life easier. This helps for:

• Continuous quantities
• Spillover quantities
• Complex counterfactuals

3.3.3 Complex counterfactuals

Here is an example of using functions to define complex counterfactuals:

f_M <- function(X, UM) 1*(UM < X)
f_Y <- function(X, M, UY) X + M - .4*X*M + UY

design <- 
  declare_model(N = 100,
                X = simple_rs(N),
                UM = runif(N),
                UY = rnorm(N),
                M = f_M(X, UM),
                Y = f_Y(X, M, UY)) +
  declare_inquiry(Q1 = mean(f_Y(1, f_M(0, UM), UY) - f_Y(0, f_M(0, UM), UY)))

design |> draw_estimands() |> kable() |> kable_styling(font_size = 20)

3.4.1 Assignment schemes

The randomizr package has a set of functions for different types of block and cluster assignments.

• Simple random assignment: “Coin flip” or Bernoulli random assignment. All units have the same probability of assignment: simple_ra(N = 100, prob = 0.25)
• Complete random assignment: Exactly m of N units are assigned to treatment, and all units have the same probability of assignment m/N complete_ra(N = 100, m = 40)

3.4.2 Assignment schemes

• Block random assignment: Complete random assignment within pre-defined blocks. Units within the same block have the same probability of assignment \(m_b / N_b\): block_ra(blocks = regions)
• Cluster random assignment: Whole groups of units are assigned to the same treatment condition. cluster_ra(clusters = households) * Block-and-cluster assignment: Cluster random assignment within blocks of clusters block_and_cluster_ra(blocks = regions, clusters = villages)

3.4.3 Assignment schemes

You can combine these in various ways. For examples with saturation random assignment first clusters are assigned to a saturation level, then units within clusters are assigned to treatment conditions according to the saturation level:

saturation = cluster_ra(clusters = villages, 
                        conditions = c(0, 0.25, 0.5, 0.75))

block_ra(blocks = villages, prob_unit = saturation)

3.5.4 A: Answers: other packages

diagnose_design(design)

if(run)
  diagnose_design(design) |> write_rds("saved/probit.rds")

read_rds("saved/probit.rds") |> 
  reshape_diagnosis() |>
  kable() |> 
  kable_styling(font_size = 20)

Why is it so terrible?

3.5.6 A: Answers: other packages

if(run)
  diagnose_design(design) |> write_rds("saved/probit_2.rds")

read_rds("saved/probit_2.rds") |> reshape_diagnosis() |> kable() |> 
  kable_styling(font_size = 20)

Much better

3.6.1 diagnosands code

DeclareDesign:::default_diagnosands

    mean_estimand <- mean(estimand)
    mean_estimate <- mean(estimate)
    bias <- mean(estimate - estimand)
    sd_estimate <- sd(estimate)
    rmse <- sqrt(mean((estimate - estimand)^2))
    power <- mean(p.value <= alpha)
    coverage <- mean(estimand <= conf.high & estimand >= conf.low)

3.6.2 More diagnosands

    mean_se = mean(std.error)
    type_s_rate = mean((sign(estimate) != sign(estimand))[p.value <= alpha])
    exaggeration_ratio = mean((estimate/estimand)[p.value <= alpha])
    var_estimate = pop.var(estimate)
    mean_var_hat = mean(std.error^2)
    prop_pos_sig = estimate > 0 & p.value <= alpha
    mean_ci_length = mean(conf.high - conf.low)

3.6.3 Custom diagnosands

my_diagnosands <-
  declare_diagnosands(median_bias = median(estimate - estimand))

diagnose_design(simplest_design, diagnosands = my_diagnosands, sims = 10) |>
  reshape_diagnosis() |> kable() |> kable_styling(font_size = 20)

3.6.4 Adding diagnosands to a design

simplest_design <- 
  set_diagnosands(simplest_design, my_diagnosands)

simplest_design |> diagnose_design(sims = 10)|>
  reshape_diagnosis() |> kable() |> kable_styling(font_size = 20)

3.6.6 Significance filter

grouped_diagnosis$simulations_df |>
  ggplot(aes(estimate, p.value, color = significant)) + geom_point()

3.6.7 Multistage simulation

• Usually a design simulation simulates “from the top”: going from the beginning to the end of the design in each run and repeating
• But sometimes you might want to follow a tree like structure and simulate different steps a different number of times

3.6.8 Multistage simulation illustration

Consider for instance this sampling design:

sampling_design <- 
  
  declare_model(N = 500, Y = 1 + rnorm(N, sd = 10)) +
  declare_inquiry(Q = mean(Y)) +
  declare_sampling(S = complete_rs(N = N, n = 100)) + 
  declare_estimator(Y ~ 1)

3.6.9 Multistage simulation illustration

Compare these two diagnoses:

diagnosis_1 <- diagnose_design(sampling_design, sims = c(5000, 1, 1, 1)) 
diagnosis_2 <- diagnose_design(sampling_design, sims = c(1, 5000, 1, 1))

In the second the estimate is drawn just once. The SD of the estimate is lower. But the RMSE is not very different.

3.6.10 Illustration with tailored estimator and diagnosand

my_estimator <- function(data) 
  data.frame(outcome = "Y", estimate = mean(data$Y), std.error = 0)

design <-
  declare_model(N = 5000, Y = rnorm(N)) +
  declare_inquiry(Y_bar_pop = mean(Y, 1)) +
  declare_sampling(S = complete_rs(N = N, n = 500)) +
  declare_inquiry(Y_bar_sample = mean(Y)) +
  declare_estimator(Y ~ 1, inquiry = "Y_bar_pop",label = "ols") +
  declare_estimator(handler = label_estimator(my_estimator),
                    inquiry = "Y_bar_sample",
                    label = "mean")
my_diagnosands <-
  declare_diagnosands(
    bias = mean(estimate - estimand),
    rmse = mean((estimate - estimand)^2)^.5,
    mean_se = mean(std.error))

3.6.11 Diagnosis

diagnose_design(design, diagnosands = my_diagnosands)

3.7.1 Redesign with vector arguments

When redesigning with arguments that are vectors, use list() in redesign, with each list item representing a design you wish to create

prob_each <- c(.1, .5, .4)

design_multi  <- 
  declare_model(N = 10) +
  declare_assignment(Z = complete_ra(N = N, prob_each = prob_each))

## returns two designs

designs <- design_multi |> 
  redesign(prob_each = list(c(.2, .5, .3), c(0, .5, .5)))
  
designs |> lapply(draw_data)

3.7.2 Redesign warnings

A parameter has to be called correctly. And you get no warning if you misname.

simplest_design_N  |> redesign(n = 200) |> draw_data() |> nrow()

[1] 100

why not 200?

3.7.3 Redesign warnings

A parameter has to be called explicitly

N <- 100

my_N <- function(n = N) n

simplest_design_N2 <- 
  
  declare_model(N = my_N(), Y = rnorm(N)) +
  declare_inquiry(Q = 0) +
  declare_estimator(Y ~ 1)

simplest_design_N2 |> redesign(N = 200) |> draw_data() |> nrow()

[1] 100

why not 200?

3.7.4 Redesign warnings

A parameter has to be called explicitly

N <- 100

my_N <- function(n = N) n

simplest_design_N2 <- 
  
  declare_model(N = my_N(N), Y = rnorm(N)) +
  declare_inquiry(Q = 0) +
  declare_estimator(Y ~ 1)

simplest_design_N2 |> redesign(N = 200) |> draw_data() |> nrow()

[1] 200

OK

3.7.5 Redesign with a function

Here is an example of redesigning where the “parameter” is a function

new_N <- function(n, factor = 1.31) n*factor

simplest_design_N2 |> redesign(my_N = new_N) |> draw_data() |> nrow()

[1] 131

4.1.1 A design: Multilevel data

A design with hierarchical data and different assignment schemes.

design <- 
  declare_model(
    school = add_level(N = 16, 
                       u_school = rnorm(N, mean = 0)),     
    classroom = add_level(N = 4,    
                  u_classroom = rnorm(N, mean = 0)),
    student =  add_level(N = 20,    
                         u_student = rnorm(N, mean = 0))
    ) +
  declare_model(
    potential_outcomes(Y ~ .1*Z + u_classroom + u_student + u_school)
    ) +
  declare_assignment(Z = simple_ra(N)) + 
  declare_measurement(Y = reveal_outcomes(Y ~ Z))  +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) + 
  declare_estimator(Y ~ Z, .method = difference_in_means)    

4.3.1 Block Random Assignment

• The treatment and control group will in expectation contain the same share of students in different classrooms.
• But as we saw this does necessarily hold in realization
• We make this more obvious by sorting the students by treatment status with schools

4.3.2 Blocked design

assignment_blocked <-   
  declare_assignment(Z = block_ra(blocks = classroom))  

estimator_blocked <- 
  declare_estimator(Y ~ Z, blocks = classroom, 
                    .method = difference_in_means)  

design_blocked <- 
  design |> 
  replace_step("assignment", assignment_blocked) |>
  replace_step("estimator", estimator_blocked)

4.3.3 Illustration of blocked assignment

• Note that subjects are sorted here after the assignment to make it easier to see that in this case blocking ensures that exactly 5 students within each classroom are assigned to treatment.

4.4.1 Clustering

But what if all students in a given class have to be assigned the same treatment?

assignment_clustered <- 
  declare_assignment(Z = cluster_ra(clusters = classroom))  
estimator_clustered <- 
  declare_estimator(Y ~ Z, clusters = classroom, 
                    .method = difference_in_means)  


design_clustered <- 
  design |> 
  replace_step("assignment", assignment_clustered) |> 
  replace_step("estimator", estimator_clustered)

4.4.2 Illustration of clustered assignment

4.5.1 Clustered and Blocked

assignment_clustered_blocked <-   
  declare_assignment(Z = block_and_cluster_ra(blocks = school,
                                              clusters = classroom))  
estimator_clustered_blocked <- 
  declare_estimator(Y ~ Z, blocks = school, clusters = classroom, 
                    .method = difference_in_means)  


design_clustered_blocked <- 
  design |> 
  replace_step("assignment", assignment_clustered_blocked) |> 
  replace_step("estimator", estimator_clustered_blocked)

4.5.2 Illustration of clustered and blocked assignment

4.6.3 Sampling distributions

diagnoses$simulations_df |> 
  mutate(design = factor(design, c("blocked", "complete", "simple", "clustered_blocked", "clustered"))) |>
  ggplot(aes(estimate)) +
  geom_histogram() + facet_grid(~design)

4.7.5 Design with two units and heterogeneous probabilities

probs <- c(.8, .2)

design <- 
  declare_model(N = 2, 
                Y_Z_1 = c(1, 1),
                Y_Z_0 = c(-1, 1)) +
  declare_inquiry(ATE = 1) +
  declare_assignment(
    Z = prob_ra(prob = probs),
    condition_prs = probs,
    Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z, label  = "ht", 
                    .method = horvitz_thompson,
                    condition_prs = condition_prs)

4.7.6 Design with two units and heterogeneous probabilities

Unbiased but very very noisy (simulations also noisy)

4.8.4 Indirect assignments: Properties

indirect_propensities <- replicate(5000, draw_data(design)$indirect) |> 
  apply(1, mean) 

4.8.5 Indirect assignments: Properties

df |> ggplot(aes(latitude, longitude, label = round(indirect_propensities_1, 2))) + geom_text()

4.8.6 Indirect assignments: Redesign

replicate(5000, draw_data(design |> redesign(n_assigned = 25))$indirect) |> 
  apply(1, mean) 

4.8.7 Indirect assignments: Redesign

df |> ggplot(aes(latitude, longitude, label = round(indirect_propensities_2, 2))) + 
  geom_text()

Looks better: but there are trade offs between the direct and indirect distributions

Figuring out the optimal procedure requires full diagnosis

4.9.1 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?

5.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”?

5.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?

5.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?

5.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?

5.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()

5.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 

5.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

5.5.1 The good

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)

5.5.2 Run it many times

sims_1 <- simulate_design(design) 

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

5.5.3 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")

5.5.4 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()

5.5.5 Check coverage is correct

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

[1] 0.9573333

5.5.6 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() 

5.5.9 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()

5.5.10 Diagnose over multiple moving parts (and ggplot)

5.5.11 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)

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

5.6.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

5.6.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

5.6.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")

5.6.4 Power from a lag?

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

5.6.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

5.6.11 Let’s fix that one

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

5.6.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