1.1.2 ATE: DIM in practice

df <- fabricatr::fabricate(N = 100, Z = rep(0:1, N/2), Y = rnorm(N) + Z)

# by hand
df |>
  summarize(Y1 = mean(Y[Z==1]), 
            Y0 = mean(Y[Z==0]), 
            diff = Y1 - Y0) |> kable(digits = 2)

# with estimatr
estimatr::difference_in_means(Y ~ Z, data = df) |>
  tidy() |> kable(digits = 2)

1.1.3 ATE: DIM in practice

We can also do this with regression:

estimatr::lm_robust(Y ~ Z, data = df) |>
  tidy() |> kable(digits = 2)

See Freedman (2008) on why regression is fine here

1.1.4 ATE: Blocks

Say now different strata or blocks \(\mathcal{S}\) had different assignment probabilities. Then you could estimate:

\[ \widehat{ATE} = \sum_{S\in \mathcal{S}}\frac{n_{S}}{n} \left(\frac{1}{n_{S1}}\sum_{S\cap T}y_i - \frac{1}{n_{S0}}\sum_{S\cap C}y_i \right) \]

Note: you cannot just ignore the blocks because assignment is no longer independent of potential outcomes: you might be sampling units with different potential outcomes with different probabilities.

However, the formula above works fine because selecting is random conditional on blocks.

1.1.5 ATE: Blocks

As a DAG this is just classic confounding:

make_model("Block -> Z ->Y <- Block") |> 
  plot(x_coord = c(2,1,3), y_coord = c(2, 1, 1))

1.1.6 ATE: Blocks in practice

Data with heterogeneous assignments:

df <- fabricatr::fabricate(
  N = 500, X = rep(0:1, N/2), 
  prob = .2 + .3*X,
  Z = rbinom(N, 1, prob),
  ip = 1/(Z*prob + (1-Z)*(1-prob)), # discuss
  Y = rnorm(N) + Z*X)

True effect is 0.5, but:

estimatr::difference_in_means(Y ~ Z, data = df) |>
  tidy() |> kable(digits = 2)

1.1.7 ATE: Blocks in practice

Averaging over effects in blocks

# by hand
estimates <- 
  df |>
  group_by(X) |>
  summarize(Y1 = mean(Y[Z==1]), 
            Y0 = mean(Y[Z==0]), 
            diff = Y1 - Y0,
            W = n())

estimates$diff |> weighted.mean(estimates$W)

[1] 0.7236939

# with estimatr
estimatr::difference_in_means(Y ~ Z, blocks = X, data = df) |>
  tidy() |> kable(digits = 2)

1.1.8 ATE with IPW

This also corresponds to the difference in the weighted average of treatment outcomes (with weights given by the inverse of the probability that each unit is assigned to treatment) and control outcomes (with weights given by the inverse of the probability that each unit is assigned to control).

• The average difference in means estimator is the same as what you would get if you weighted inversely by shares of units in different conditions inside blocks.

1.1.9 ATE with IPW in practice

# by hand
df |>
  summarize(Y1 = weighted.mean(Y[Z==1], ip[Z==1]), 
            Y0 = weighted.mean(Y[Z==0],  ip[Z==0]), # note !
            diff = Y1 - Y0)|> 
  kable(digits = 2)

# with estimatr
estimatr::difference_in_means(Y ~ Z, weights = ip, data = df) |>
  tidy() |> kable(digits = 2)

1.1.10 ATE with IPW

• But inverse propensity weighting is a more general principle, which can be used even if you do not have blocks.

• The intuition for it comes straight from sampling weights — you weight up in order to recover an unbiased estimate of the potential outcomes for all units, whether or not they are assigned to treatment.

• With sampling weights however you can include units even if their weight was 1. Why can you not include these units when doing inverse propensity weighting?

1.1.11 Illustration: Estimating treatment effects with terrible treatment assignments: Fixer

Say you made a mess and used a randomization that was correlated with some variable, \(U\). For example:

• The randomization is done in a way that introduces a correlation between Treatment Assignment and Potential Outcomes
• Then possibly, even though there is no true causal effect, we naively estimate a large one — enormous bias
• However since we know the assignment procedure we can fully correct for the bias

1.1.12 Illustration: Estimating treatment effects with terrible treatment assignments: Fixer

• In the next example, we do this using “inverse propensity score weighting.”
• This is exactly analogous to standard survey weighting — since we selected different units for treatment with different probabilities, we weight them differently to recover the average outcome among treated units (same for control).

1.1.13 Basic randomization: Fixer

Bad assignment, some randomization process you can’t understand (but can replicate) that results in unequal probabilities.

N <- 400
U <- runif(N, .1, .9)

 design <- 
   declare_model(N = N,
                 Y_Z_0 = U + rnorm(N, 0, .1),
                 Y_Z_1 = U + rnorm(N, 0, .1),
                 Z = rbinom(N, 1, U)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_estimator(Y ~ Z, label = "naive") 

1.1.14 Basic randomization: Fixer

Results is a sampling distribution not centered on the true effect (0)

diagnosis <- diagnose_design(design)

diagnosis$simulations_df |>
  ggplot(aes(estimate)) + geom_histogram() + facet_grid(~estimator) +
  geom_vline(xintercept = 0, color = "red")

1.1.15 A fix

To fix you can estimate the assignment probabilities by replicating the assignment many times:

probs <- replicate(1000, design |> draw_data() |> pull(Z)) |> apply(1, mean)

and then use these assignment probabilities in your estimator

design_2 <-
  design +
  declare_measurement(weights = Z/probs + (1-Z)/(1-probs)) +
  declare_estimator(Y ~ Z, weights = weights, label = "smart") 

1.1.16 Basic randomization: Fixer

Implied weights

draw_data(design_2) |> 
  ggplot(aes(probs, weights, color = factor(Z))) + 
  geom_point()

1.1.18 IPW with one unit!

This example is surprising but it helps you see the logic of why inverse weighting gets unbiased estimates (and why that might not guarantee a reasonable answer)

Imagine there is one unit with potential outcomes \(Y(1) = 2, Y(0) = 1\). So the unit level treatment effect is 1.

You toss a coin.

• If you assign to treatment you estimate: \(\hat\tau = \frac{2}{0.5} = 4\)
• If you assign to control you estimate: \(\hat\tau = -\frac{1}{0.5} = -2\)

So your expected estimate is: \[0.5 \times 4 - 0.5 \times (-2) = 1\]

Great on average but always lousy

1.1.19 Generalization: why IPW works

• Say a given unit is assigned to treatment with probability \(\pi_i\)
• We estimate the average \(Y(1)\) using

\[\hat{\overline{Y_1}} = \frac{1}n\left(\sum_i \frac{Z_iY_i(1)}{\pi_i}\right)\] With independent assignment the expected value of \(\hat{\overline{Y_1}}\) is just:

\[\mathbb{E}[\hat{\overline{Y_1}}] =\frac1n\left( \left(\pi_1 \frac{1\times Y_1(1)}{\pi_1} + (1-\pi_1) \frac{0\times Y_1(1)}{\pi_1}\right) + \left(\pi_2 \frac{1\times Y_2(1)}{\pi_2} + (1-\pi_2) \frac{0\times Y_1(1)}{\pi_2}\right) + \dots\right)\]

\[\mathbb{E}[\hat{\overline{Y_1}}] =\frac1n\left( Y_1(1) + Y_2(1) + \dots\right) = \overline{Y_1}\]

and similarly for \(\mathbb{E}[\hat{\overline{Y_0}}]\) and so using linearity of expectations:

\[\mathbb{E}[\widehat{\overline{Y_1 - Y_0}}] = \overline{Y_1 - Y_0}\]

1.1.20 Generalization: why IPW works

• Note we needed \(\pi_i >0\) and also \(\pi_i <1\) everywhere. Why?
• We used independence here; sampling theory is used to show similar results for e.g. complete randomization
• For blocked randomization this is easy to see

1.2.3 Var(ATE)

• In classical statistics we characterize our uncertainty over an estimate using an estimate of variance of the sampling distribution of the estimator.

• Key idea is we want to be able to say: how likely are we to have gotten such an estimate if the distribution of estimates associated with our design looked a given way.

• More specifically we want to estimate “standard error” or the “standard deviation of the sampling distribution”

(See Woolridge (2023) where the standard error is understood as the “estimate of the standard deviation of the sampling distribution”)

1.2.4 Variance and standard errors

Given:

• \(\hat\tau\) is an estimate for \(\tau\)
• \(\overline{x}\) is the average values of \(x\)

The variance of the estimator of \(n\) repeated ‘runs’ of a design is: \(Var(\hat{\tau}) = \frac{1}n\sum_i(\hat\tau_i - \overline{\hat\tau_i})^2\)

And the standard error is:

\(se(\hat{\tau}) = \sqrt{\frac{1}n\sum_i(\hat\tau_i - \overline{\hat\tau_i})^2}\)

1.2.5 Variance and standard errors

If we have a good measure for the shape of the sampling distribution we can start to make statements of the form:

• What are the chances that an estimate would be this large or larger?

If the sampling distribution is roughly normal, as it may be with large samples, then we can use procedures such as: “there is a 5% probability that an estimate would be more than 1.96 standard errors away from the mean of the sampling distribution”

1.2.6 Var(ATE)

• Key idea: You can estimate variance straight from the data, given knowledge of the assignment process and assuming well defined potential outcomes?

• Recall in general \(Var(x) = \frac{1}n\sum_i(x_i - \overline{x})^2\). here the \(x_i\)s are the treatment effect estimates we might get under different random assignments, the \(n\) is number of different assignments (assumed here all equally likely, but otherwise we can weight) and \(\overline{x}\) is the truth.

• For intuition imagine we have just two units \(A\), \(B\), with potential outcomes \(A_1\), \(A_0\), \(B_1\), \(B_0\).

• When there are two units with outcomes \(x_1, x_2\), the variance simplifies like this:

\[Var(x) = \frac{1}2\left(x_1 - \frac{x_1 + x_2}{2}\right)^2 + \frac{1}2\left(x_2 - \frac{x_1 + x_2}{2}\right)^2 = \left(\frac{x_1 - x_2}{2}\right)^2\]

1.2.7 Var(ATE)

In the two unit case the two possible treatment estimates are: \(\hat{\tau}_1=A_1 - B_0\) and \(\hat{\tau}_2=B_1 - A_0\), depending on what gets put into treatment. So the variance is:

\[Var(\hat{\tau}) = \left(\frac{\hat{\tau}_1 - \hat{\tau}_2}{2}\right)^2 = \left(\frac{(A_1 - B_0) - (B_1 - A_0)}{2}\right)^2 =\left(\frac{(A_1 - B_1) + (A_0 - B_0)}{2}\right)^2 \] which we can re-write as:

\[Var(\hat{\tau}) = \left(\frac{A_1 - B_1}{2}\right)^2 + \left(\frac{A_0 - B_0}{2}\right)^2+ 2\frac{(A_1 - B_1)(A_0-B_0)}{2}\] The first two terms correspond to the variance of \(Y(1)\) and the variance of \(Y(0)\). The last term is a bit pesky though, it corresponds to twice the covariance of \(Y(1)\) and \(Y(0)\).

1.2.8 Var(ATE)

How can we go about estimating this?

\[Var(\hat{\tau}) = \left(\frac{A_1 - B_1}{2}\right)^2 + \left(\frac{A_0 - B_0}{2}\right)^2+ 2\frac{(A_1 - B_1)(A_0-B_0)}{2}\]

In the two unit case it is quite challenging because we do not have an estimate for any of the three terms: we do not have an estimate for the variance in the treatment group or in the control group because we have only one observation in each case; and we do not have an estimate for the covariance because we don’t observe both potential outcomes for any case.

Things do look a bit better however with more units…

1.2.9 Var(ATE): Generalizing

From Freedman Prop 1 / Example 1 (using combinatorics!) we have:

\(V(\widehat{ATE}) = \frac{1}{n-1}\left[\frac{n_C}{n_T}V_1 + \frac{n_T}{n_C}V_0 + 2C_{01}\right]\)

… where \(V_0, V_1\) denote variances and \(C_{01}\) covariance

This is usefully rewritten as:

\[ \begin{split} V(\widehat{ATE}) & = \frac{1}{n-1}\left[\frac{n - n_T}{n_T}V_1 + \frac{n - n_C}{n_C}V_0 + 2C_{01}\right] \\ & = \frac{n}{n-1}\left[\frac{V_1}{n_T} + \frac{V_0}{n_C}\right] - \frac{1}{n-1}\left[V_1 + V_0 - 2C_{01}\right] \end{split} \]

where the final term is positive

1.2.12 Illustration of Neyman Conservative Estimator

An illustration of how conservative the conservative estimator of variance really is (numbers in plot are correlations between \(Y(1)\) and \(Y(0)\).

We confirm that:

1. the estimator is conservative
2. the estimator is more conservative for negative correlations between \(Y(0)\) and \(Y(1)\) — eg if those cases that do particularly badly in control are the ones that do particularly well in treatment, and
3. with \(\tau\) and \(V(Y(0))\) fixed, high positive correlations are associated with highest variance.

1.2.13 Illustration of Neyman Conservative Estimator

Here \(\rho\) is the unobserved correlation between \(Y(1)\) and \(Y(0)\); and \(\Delta\) is the final term in the sample variance equation that we cannot estimate.

1.2.14 Illustration of Neyman Conservative Estimator

1.2.16 Tighter Bounds On Variance Estimate

Example:

• Take a million-observation dataset, with treatment randomly assigned
• Assume \(Y(0)=0\) for everyone and \(Y(1)\) distributed normally with mean 0 and standard deviation of 1000.
• Note here the covariance of \(Y(1)\) and \(Y(0)\) is 0.
• Note the true variance of the estimated sample average treatment effect should be (approx) \(\frac{Var(Y(1))}{{1000000}} + \frac{Var(Y(0))}{{1000000}} = 1+0=1\), for an se of \(1\).
• But using the Neyman estimator (or OLS!) we estimate (approx) \(\frac{Var(Y(1))}{({1000000/2})} + \frac{Var(Y(0))}{({1000000/2})} = 2\), for an se of \(\sqrt{2}\).
• But we can recover the truth knowing the covariance between \(Y(1)\) and \(Y(0)\) is 0.

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

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

1.3.7 Randomization Inference

• Randomization inference can get more complicated when you want to test a null other than the sharp null of no effect.
• Say you wanted to test the null that the effect is 2 for all units. How do you do it?
• Say you wanted to test the null that an interaction effect is zero. How do you do it?
• In both cases by filling in a potential outcomes schedule given the hypothesis in question and then generating a test statistic

1.3.8 ri and CIs

It is possible to use this procedure to generate confidence intervals with a natural interpretation.

• The key idea is that we can use the same procedure to assess the probability of the data given a sharp null of no effect, but also a sharp null of any other **constant* effect.
• We can then see what set of effects we reject and what set we accept
• We are left with a set of values that we cannot reject at the 0.05 level.

1.3.9 ri and CIs in practice

candidates <- seq(-.05, .3, length = 50)

get_p <- function(j)
  (conduct_ri(
      formula = Y ~ Z,
      declaration = assignment,
      sharp_hypothesis = j,
      data = df,
      sims = 5000,
    ) |> summary()
  )$two_tailed_p_value

# Implement
ps <- candidates |> sapply(get_p)

1.3.10 ri and CIs in practice

Warning: calculating confidence intervals this way can be computationally intensive

1.3.11 ri with DeclareDesign

• DeclareDesign can do randomization inference natively
• The trick is to ensure that when calculating the \(p\) values the only stochastic component is the assignment to treatment when calculating the \(p\) values

1.3.12 ri with DeclareDesign (advanced)

Here we get minimal detectable effects by using a design that has two stage simulations so we can estimate the sampling distribution of summaries of the sampling distribution generated from reassignments.

test_stat <- function(data)
  with(data, data.frame(estimate = mean(Y[Z==1]) - mean(Y[Z==0])))

b <- 0

design <- 
  declare_model(N = 100,  Z = complete_ra(N), Y = b*Z + rnorm(N)) +
  declare_estimator(handler = label_estimator(test_stat), label = "actual")+
  declare_measurement(Z = sample(Z))  + # this is the permutation step
  declare_estimator(handler = label_estimator(test_stat), label = "null")

1.3.13 ri with DeclareDesign (advanced)

Simulations data frame from two step simulation. Note computational intensity as number of runs is the product of the sims vector. I speed things up by using a simple estimation function and also using parallelization.

library(future)
options(parallelly.fork.enable = TRUE) 
plan(multicore) 

simulations <- 
  design |> redesign(b = c(0, .25, .5, .75, 1)) |> 
  simulate_design(sims = c(200, 1, 1000, 1)) 

1.3.14 ri with DeclareDesign (advanced)

A snap shot of the simulations dataframe: we have multiple step 3 draws for each design and step 1 draw.

simulations |> head() |> kable(digits = 2)

1.3.15 ri with DeclareDesign (advanced)

Power for each value of b.

simulations |> group_by(b, step_1_draw) |> 
  summarize(p = mean(abs(estimate[estimator == "null"]) >= abs(estimate[estimator == "actual"]))) |>
  group_by(b) |> summarize(power = mean(p <= .05)) |>
  ggplot(aes(b, power)) + geom_line() + geom_hline(yintercept = .8, color = "red")

If you want to figure out more precisely what b gives 80% or 90% power you can narrow down the b range.

1.3.16 ri interactions

Lets now imagine a world with two treatments and we are interested in using ri for assessing the interaction. (Code from Coppock, ri2)

set.seed(1)

N <- 100
declaration <- randomizr::declare_ra(N = N, m = 50)

data <- fabricate(
  N = N,
  Z = conduct_ra(declaration),
  X = rnorm(N),
  Y = .9 * X + .2 * Z + .1 * X * Z + rnorm(N))

1.3.17 ri interactions

The approach is to declare a null model that is nested by the full model. Then \(F\) test statistic from the model comparisons is taken as the test statistic and distribution of this is built up under re-randomizations.

conduct_ri(
    model_1 = Y ~ Z + X,
    model_2 = Y ~ Z + X + Z * X,
    declaration = declaration,
    assignment = "Z",
    sharp_hypothesis = coef(lm(Y ~ Z, data = data))[2],
    data = data, 
    sims = 1000
  )  |> summary() |> kable()

1.3.18 ri interactions with DeclareDesign

Let’s imagine a true model with interactions. We take an estimate. We then ask how likely that estimate is from a null model with constant effects

Note: this is quite a sharp hypothesis

df <- fabricate(N = 1000, Z1 = rep(c(0,1), N/2), Z2 = sample(Z1), Y = Z1 + Z2 - .15*Z1*Z2 + rnorm(N))

my_estimate <- (lm(Y ~ Z1*Z2, data = df) |> coef())[4]

null_model <-  function(df) {
  M0 <- lm(Y ~ Z1 + Z2, data = df) 
  d1 <- coef(M0)[2]
  d2 <- coef(M0)[3]
  df |> mutate(
    Y_Z1_0_Z2_0 = Y - Z1*d1 - Z2*d2,
    Y_Z1_1_Z2_0 = Y + (1-Z1)*d1 - Z2*d2,
    Y_Z1_0_Z2_1 = Y - Z1*d1 + (1-Z2)*d2,
    Y_Z1_1_Z2_1 = Y + (1-Z1)*d1 + (1-Z2)*d2)
  }

1.3.19 ri interactions with DeclareDesign

Let’s imagine a true model with interactions. We take an estimate. We then ask how likely that estimate is from a null model with constant effects

Imputed potential outcomes look like this:

df <- df |> null_model()

df |>  head() |> kable(digits = 2, align = "c")

1.3.20 ri interactions with DeclareDesign

design <- 
  declare_model(data = df) +
  declare_measurement(Z1 = sample(Z1), Z2 = sample(Z2),
                      Y = reveal_outcomes(Y ~ Z1 + Z2)) +
  declare_estimator(Y ~ Z1*Z2, term = "Z1:Z2")

diagnose_design(design, sims = 1000, diagnosands = ri_ps(my_estimate))

1.3.24 ri Applications

1.3.25 ri Applications

If you think hard about assignment you might come up with an allocation like this.

1.3.26 ri Applications

That allocation balances as much as possible. Given the allocation you might randomly assign individuals to different days as well as randomly assigning them to treatments within days. If you then figure out assignment propensities, this is what you would get:

1.3.27 ri Applications

Even under the assumption that the day of measurement does not matter, these assignment probabilities have big implications for analysis.

• Only the type \(A\) subjects could have received any of the three treatments.

• There are no two treatments for which it is possible to compare outcomes for subpopulations \(B\) and \(C\)

• A comparison of \(T1\) versus \(T2\) can only be made for population \(A \cup B\)

• However subpopulation \(A\) is assigned to \(A\) (versus \(B\)) with probability 4/5; while population \(B\) is assigned with probability 3/8

1.3.28 ri Applications

• Implications for design: need to uncluster treatment delivery

• Implications for analysis: need to take account of propensities

Idea: Wacky assignments happen but if you know the propensities you can do the analysis.

1.3.29 ri Applications: Indirect assignments

A particularly interesting application is where a random assignment combines with existing features to determine an assignment to an “indirect” treatment.

• For instance: \(n\) of \(N\) are assigned to a treatment.
• You are interested in whether “having a friend assigned to treatment” makes a difference to a subject. Or maybe “a friend of a friend”
• That means the subject has a complex clustered assignment that depends on how many friends they have
• A bit mind-boggling, but:
  • Rerun your assignment many times and each time figure out whether a subject is assigned to an indirect treatment or not
  • Calculate the implied quantity of interest for each assignment
  • Assess the place of the actual quantity in the sampling distribution

1.4.1 Example

Consider for example this data.

• You randomly pair offerers and receivers in a dictator game (in which offerers decide how much of $1 to give to receivers).
• Your population comes from two groups (80% Baganda and 20% Banyankole) so in randomly assigning partners you are randomly determining whether a partner is a coethnic or not.
• You find that in non-coethnic pairings 35% is offered, in coethnic pairings 48% is offered.

Should you believe it?

1.4.2 Covariate Adjustment

• Population: randomly matched Baganda (80% of pop) and Banyankole (20% of pop)
• You find: in non-coethnic pairings 35% is offered, in coethnic pairings 48% is offered.
• But a closer look at the data reveals…

So that’s a problem

1.4.5 Covariate adjustment via saturated regression

• Alternatively you can use propensity weights.
• Alternatively you can use a regression that includes both the treatment and the treatment interacted with the covariates.
  • In practice this is best done by demeaning the covariates; doing this lets you read off the average effect from the main term. Key resource: Lin (2012)

You should have noticed that the logic for controlling for a covariate here is equivalent to the logic we saw for heterogeneous assignment propensities. These are really the same thing.

1.4.6 Covariate adjustment via saturated regression

Returning to prior example:

df <- fabricatr::fabricate(
  N = 500, 
  X = rep(0:1, N/2), 
  Z = rbinom(N, 1, .2 + .3*X),
  Y = rnorm(N) + Z*X)

lm_robust(Y ~ Z*X_c, data = df |> mutate(X_c = X - mean(X))) |>
  tidy() |> kable(digits = 2)

1.4.7 Covariate adjustment via saturated regression

Returning to prior example:

lm_lin(Y ~ Z, ~ X, data = df) |>
  tidy() |> kable(digits = 2)

1.4.9 Demeaning and saturating

Let’s:

• Declare a factorial design in which Y is generated according to

f_Y <- function(X1, X2, u) .1 + .2*X1 + .3*X2 + u*X1*X2

where u is distributed \(U[0,1]\).

• Specify estimands carefully
• Run analyses in which we do and do not demean the treatments; compare and explain results

1.4.10 Demeaning interactions

f_Y <- function(X1, X2, u) .1 + .2*X1 + .3*X2 + u*X1*X2

design <-
  declare_model(N = 1000, u = runif(N),
                X1 = complete_ra(N), X2 = block_ra(blocks = X1),
                X1_demeaned = X1 - mean(X1),
                X2_demeaned = X2 - mean(X2),
                Y = f_Y(X1, X2, u)) +
  declare_inquiry(
    base = mean(f_Y(0, 0, u)),
    average = mean(f_Y(0, 0, u) + f_Y(0, 1, u)  + f_Y(1, 0, u)  + f_Y(1, 1, u))/4,
    CATE_X1_given_0 = mean(f_Y(1, 0, u) - f_Y(0, 0, u)),
    CATE_X2_given_0 = mean(f_Y(0, 1, u) - f_Y(0, 0, u)),
    ATE_X1 = mean(f_Y(1, X2, u) - f_Y(0, X2, u)),
    ATE_X2 = mean(f_Y(X1, 1, u) - f_Y(X1, 0, u)),
    I_X1_X2 = mean((f_Y(1, 1, u) - f_Y(0, 1, u)) - (f_Y(1, 0, u) - f_Y(0, 0, u)))
  ) +
  declare_estimator(Y ~ X1*X2, 
                    inquiry = c("base", "CATE_X1_given_0", "CATE_X2_given_0", "I_X1_X2"), 
                    term = c("(Intercept)", "X1", "X2", "X1:X2"),
                    label = "natural") +
  declare_estimator(Y ~ X1_demeaned*X2_demeaned, 
                    inquiry = c("average", "ATE_X1", "ATE_X2", "I_X1_X2"), 
                    term = c("(Intercept)", "X1_demeaned", "X2_demeaned", "X1_demeaned:X2_demeaned"),
                    label = "demeaned")

1.4.11 Demeaning interactions: Solution

It’s all good. But you need to match the estimator to the inquiry: demean for average marginal effects; do not demean for conditional marginal effects.

1.4.12 Recap

If you have different groups with different assignment propensities you can do any or all of these:

1. Blocked differences in means
2. Inverse propensity weighting
3. Saturated regression (Lin)
4. More… (coming)

You cannot (reliably):

1. Ignore the groups
2. Include them in a regression (without interactions)

1.5.1 Considerations

• Even though randomization ensures no bias, you may sometimes want to “control” for covariates in order to improve efficiency (see the discussion of blocking above).
• Or you may have to take account of the fact that the assignment to treatment is correlated with a covariate (as above).
• In observational work you might also figure out you have to control for a covariate to justify inferences (Refer to our discussion of the backdoor criteria)

1.5.2 Observational work

• Observational motivation: Controls can provide grounds for identification
• But recall – they can also destroy identification

For a great walk through of what you can draw from graphical models for the decision to control see:

A Crash Course in Good and Bad Controls by Cinelli, Forney, and Pearl (2022)

Aside: these implications generally refer to use controls as covariates – e.g. by implementing blocked differences in means or similar. For a Bayesian model of the form used in CausalQueries the information from “bad controls” is used wisely.

1.5.4 Precision Gains from Controls

However: Introducing controls can create complications

• As argued by Freedman (summary from Lin (2012)), we can get: “worsened asymptotic precision, invalid measures of precision, and small-sample bias”\(^*\)

• These adverse effects are essentially removed with an interacted model

• See discussions in Imbens and Rubin (2015) (7.6, 7.7) and especially Theorem 7.2 for the asymptotic variance of the estimator

\(^*\) though note that the precision concern does not hold when treatment and control groups are equally sized

1.5.5 To control or not

We will illustrate by comparing:

• DIM
• OLS (linear controls)
• Lin (saturated)

With:

• Varying fraction assigned to treatment
• Varying relation between \(Y(0)\) and \(Y(1)\)

1.5.6 Declaration (from RDSS)

# https://book.declaredesign.org/library/experimental-causal.html

prob <- 0.5
control_slope <- -1

declaration_18.3 <-
  declare_model(N = 100, X = runif(N, 0, 1),
                U = rnorm(N, sd = 0.1),
                Y_Z_1 = 1*X + U, Y_Z_0 = control_slope*X + U
  ) +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_assignment(Z = complete_ra(N = N, prob = prob)) + 
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z, inquiry = "ATE", label = "DIM") +
  declare_estimator(Y ~ Z + X, .method = lm_robust, inquiry = "ATE", label = "OLS") +
  declare_estimator(Y ~ Z, covariates = ~X, .method = lm_lin, 
                    inquiry = "ATE", label = "Lin")

1.5.8 To control or not

Diagnosis:

declaration_18.3 |> 
  redesign(
    control_slope = seq(-1, 1, 0.5), 
    prob = seq(0.1, 0.9, 0.1)) |> 
  diagnose_designs()

1.5.9 To control or not

• Lin always does well
• OLS does fine with equal probability assignments
• Otherwise the ranking of DIM and OLS is design dependent

1.5.13 Standard errors

We see the standard errors are larger when you control in cases in which the control is not predictive of the outcome and it is correlated with the treatment. Otherwise they can be smaller.

See Mutz et al

1.6.1 Doubly robust estimation

Doubly robust estimation combines:

1. A model for how the covariates predict the potential outcomes
2. A model for how the covariates predict assignment propensities

Using both together to estimate potential outcomes using propensity weighting lets you do well even if either model is wrong.

Each part can be done using nonparameteric methods resulting in an overall semi-parametric procedure.

• \(\pi(Z) = \Pr(Z=1|X)\): Estimate \(\hat\pi\)
• \(Y_z = \mathbb{E}[Y|Z=z, X]\): Estimate \(\hat{Y}_z\)
• Estimate of causal effect: \(\frac{1}{n}\sum_{i=1}^n\left(\left(\frac{Z_i}{\hat{\pi}_i}(Y_i - \hat{Y}_{i1}\right) - \left(\frac{1-Z_i}{1-\hat{\pi}_i}(Y_i - \hat{Y}_{i0}\right) + \left(\hat{Y}_{i1} - \hat{Y}_{i0}\right) \right)\)

1.6.3 Doubly robust estimation

• More subtly say the \(\hat{\pi}\)s are correct, but your imputations are wrong; then we again have an unbiased estimator.

To see this imagine with probability \(\pi\) we assign unit 1 to treatment and 2 to control (otherwise 1 to control and 2 to treatment).

Then our expected estimate is:

\(\frac12\pi\left(\left(\frac{1}{\pi}(Y_{11} - \hat{Y}_{11}\right) - \left(\frac{1}{\pi}(Y_{20} - \hat{Y}_{20}\right) \right) + (1-\pi)\left(\left(\frac{1}{1-\pi}(Y_{21} - \hat{Y}_{21}\right) - \left(\frac{1}{1-\pi}(Y_{10} - \hat{Y}_{10}\right) \right) + \left(\hat{Y}_{11} - \hat{Y}_{20}\right) + \left(\hat{Y}_{21} - \hat{Y}_{10}\right)\)

\(\frac12\left(Y_{11} - Y_{10} + Y_{21}- Y_{20} +\pi\left(\left(\frac{1}{\pi}( - \hat{Y}_{11}\right) - \left(\frac{1}{\pi}( - \hat{Y}_{20}\right) \right) + (1-\pi)\left(\left(\frac{1}{1-\pi}( - \hat{Y}_{21}\right) - \left(\frac{1}{1-\pi}(- \hat{Y}_{10}\right) \right)\right) + \left(\hat{Y}_{11} - \hat{Y}_{20}\right) + \left(\hat{Y}_{21} - \hat{Y}_{10}\right)\)

\(\frac12\left(Y_{11} - Y_{10} + Y_{21}- Y_{20}\right)\)

Robins, Rotnitzky, and Zhao (1994)