Chapter 2

Code Call-out 2.1: Regression, Comparison of Means and Covariates

In this code call-out we will explore some basic elements of regression as a means to estimate treatment effects. First we will, very simply, confirm the equivalance between a regression and difference of means estimator for the pooint estimate of binary treatment receipt on outcomes. And secondly, we will explore effect heterogeneity by a single covariate. We will do this with data from Bari et al. (2024), who experimentally assess the impact of access to large business assets in microenteprises in Pakistan, where a randomly treated sample was offered a larger microfinance route to purchase such large assets.

Equivalence between regression and comparison of means

To understand the equivalence between regression analysis and the comparison of means in a binary regression set-up, we will first open the data from Bari et al. (2024) and run a basic regression such as that documented in Table 4 of their paper:

library(dplyr)
#library(readxl)
#library(plyr)
library(ggplot2)
#library(xtable)
#library(tidyr)
library(haven)
library(sandwich)
library(lmtest)

data <- read_dta("data/Bari_et_al_2024.dta")

# Keep only observations where wave == 1
data_wave1 <- subset(data, wave == 1)

# Linear regression model with robust errors
model <- lm(biz_ta ~ A, data = data_wave1)
robust_se <- vcovHC(model, type = "HC2")

coeftest(model, vcov = robust_se)

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  1149.01     110.06 10.4399 < 2.2e-16 ***
A             507.04     134.21  3.7779 0.0001709 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef_treatment <- coef(model)["A"]

print(coef_treatment)
       A 
507.0368

Here biz_ta is an outcome measuring total business assets, and this is measured at wave 1 (ie the wave in data immediately after treatment was first provided).

To confirm that this is–by definition–simply equivalent to a difference in means we can consider the quantities which would be observed if we substituted \(A=1\) and \(A=0\) into the above equation. In this case we can see that the regression consists of a constant term: \(E[Y_i|Treat_i=0]\), and a second term which capture the difference between units receiving treatment and those not receiving treatment: \(E[Y_i|Treat_i=1]-E[Y_i|Treat_i=0]\). This latter is just the treatment effect of interest. Below we confirm that the coefficient on \(A\) above is equivalent to the difference in means.

# mean for the treatment group
mean_treatment <- mean(data_wave1$biz_ta[data_wave1$A == 1], na.rm = TRUE)

# mean for the control group
mean_control <- mean(data_wave1$biz_ta[data_wave1$A == 0], na.rm = TRUE)

# difference in means
diff_means <- mean_treatment - mean_control


cat("Coefficient for treatment in regression:", coef_treatment, "\n")
Coefficient for treatment in regression: 507.0368 
cat("Difference in means (Treatment - Control):", diff_means, "\n")
Difference in means (Treatment - Control): 507.0368 
cat("The coefficient from the regression should be equal to the difference in means to demonstrate equivalence.\n")
The coefficient from the regression should be equal to the difference in means to demonstrate equivalence.

Here we see that, as expected, our mean for the control group is equivalent to the constant in the regression, and the difference in means is equivalent to the coefficient on the treatment indicator in the regression.

Heterogeneity by a baseline covariate \(X\)

Let’s now consider the incorporation of a covariate into this regression and confirm that we can also quite simply break down (linear) heterogeneity, as discussed in (2.10) of the book. Here we will consider a variable which measures the management capacity of the small businesses in the sample. We will consider this management capacity at baseline, which below we generate, standardise so that it is mean zero, and then interact with treatment:

data$temp    <- ifelse(data$wave == 0, data$mpall, NA)
data$manage  <- ave(data$temp, data$id, FUN = function(x) mean(x, na.rm = TRUE))
data$manage  <- data$manage - mean(data$manage, na.rm = TRUE)
data$Amanage <- data$A * data$manage

summary(data$Amanage)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-1.0770 -0.4426  0.0000  0.0122  0.1438  2.3057

We can now consider whether there are differential treatment effects by a firm’s management capacity. To do so, we estimate the interacted regression below which allows for (a) a differential intercept for treated and untreated groups by the inclusion of the treatment dummy and a constant, and (b) a differential slope by management capacity, by including both a control for management capacity, as well as an interaction with treatment:

model <- lm(biz_ta ~ A + manage + Amanage, data = data, subset = (wave == 1))
coeftest(model, vcov = vcovHC(model, type = "HC2"))

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  1167.59     110.75 10.5423 < 2.2e-16 ***
A             485.21     134.66  3.6034 0.0003353 ***
manage        442.00     164.48  2.6873 0.0073669 ** 
Amanage      -229.22     194.02 -1.1814 0.2378188    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Above we eastimate this regression, and observe that while in general business with higher management capacity at baseline have much higher business assets, there is not (at least in wave 1) a statistically significantly different return to management among treated and untreated individuals. There is, however, a quite large negative effect on the interaction, weakly suggestive of smaller treatment effects among firms with higher management capacity. We can visualise this as in Figure 2.1 of the book. There are many ways we can do this, and you may wish to explore alternative ways of setting this up by considering marginal effects, or by confirming that the four quantities estimated in the regression above correspond to the two intercepts and two slopes of the lines plotted below which can be similarly plotted “by hand” based on these values. However, below we can see that it is sufficient to plot the predicted effect among treated and control units across the range of the manage variable. We do this in the code block below, first generating this predicted effect from the regression, and then plotting resulting curves for treated and control units.

df_pred <- data %>%
  filter(wave == 1) %>%
  mutate(yhat = predict(model, newdata = .))

ggplot(df_pred, aes(x = manage, y = yhat, color = factor(A), linetype = factor(A))) +
  geom_line(size = 1) +
  scale_color_manual(values = c("blue", "red"),
                     labels = c("Untreated", "Treated")) +
  scale_linetype_manual(values = c("dashed", "solid"),
                        labels = c("Untreated", "Treated")) +
  labs(x = "manage",
       y = "Linear prediction",
       color = "",
       linetype = "") +
  theme_minimal() +
  theme(legend.position = "bottom")

While this is just a simple illustration, it points to the flexibility of regression for capturing treatment effects when treatment is randomly assigned, and assumptions of unconditional unconfoundedness are reasonable.

Code Call-out 2.2: Randomization Inference

Randomization inference is perhaps best-illustrated with practical examples. A particularly illuminating approach to understand how randomisation inference works is visualization through tabular permutation. In this code call-out we will first consider a made-up example based on 3 treated units and 3 control units, before working with a larger number of control and treated units in data from an experiment implemented by Banerjee, Duflo, and Sharma (2021), which we discuss below.

An Exact p-value

It is perhaps useful to see a simple example. Consider the case of 6 units, with 3 observations randomly assigned treatment. Imagine that after having been exposed to treatment, the observed outcomes were, in the treatment group: \((34,27,29)\) and in the control group: \((14,18,24)\). A simple comparison of means estimator suggests that the treatment effect is 11.33. To calculate a p-value, we can permute all the possible combinations, and ask what proportion of these are greater than or equal to this treatment effect. If we consider random orderings of 6 units, this suggests that there are \(6!\) possible combinations, but in reality, as we are randomly choosing 3 units from these 6 to assign a permuted treatment status, the actual value of different combinations is \(6\choose 3\) \(=\frac{6!}{3!*(6-3)!}=20\). We document each of these possible permutations, as well as their permuted treatment effect in the Table below. In this case, we can see that only 1 of the 20 different permutations is greater than or equal to 11.33 (the original treatment assignment). Suggesting an exact p-value of \(1/20=0.05\) if a one-sided test is considered. If, however, we wish to consider a two-sided p-value, there are two values as extreme as 11.33, which is permutation 1, and permutation 20 below, suggesting a two-sided p-value of \(2/10=0.10\).

A Simple Illustration of Randomization Inference
Permutation T1 T2 T3 C1 C2 C3 Estimate
Original (1) 34 27 29 14 18 24 11.33
2 34 27 14 29 18 24 1.33
3 34 27 18 14 29 24 4
4 34 27 24 14 18 29 8
5 34 14 29 27 18 24 2.67
6 34 18 29 14 27 24 5.33
7 34 24 29 14 18 27 9.33
8 14 27 29 34 18 24 -2
9 18 27 29 14 34 24 0.67
10 24 27 29 14 18 34 4.67
11 34 14 18 27 29 24 -4.67
12 34 14 24 27 18 29 -0.67
13 34 18 24 14 27 29 2
14 14 27 18 34 29 24 -9.33
15 14 27 24 34 18 29 -5.33
16 18 27 24 14 34 29 -2.67
17 14 18 29 34 27 24 -8
18 14 24 29 34 18 27 -4
19 18 24 29 14 34 27 -1.33
20 14 18 24 34 27 29 -11.33

While this is so simple that we can set it up by hand, it is also useful to see how we can compute this in R. First, we will import our “data”:

data <- data.frame(
  Y = c(34, 27, 29, 14, 18, 24),
  W = c(1, 1, 1, 0, 0, 0)
)

and then calculate a comparison of means estimator of the ATT:

mean_Y1 <- mean(data$Y[data$W == 1])
mean_Y0 <- mean(data$Y[data$W == 0])

tau_hat <- mean_Y1 - mean_Y0

cat("Treatment effect is:", tau_hat, "\n")
Treatment effect is: 11.33333

Now let’s generate our p-value by permuting all possible treatment combinations. Rather than iterating manually through nested loops, we can use R’s combn(1:6, 3) function, which generates all possible combinations of 3 treated units drawn from 6 directly, returning these as a matrix of 20 columns. Below we begin by preparing a data frame in which to store our resulting treatment effects:

permutations <- data.frame(
  permutation = 1:20,
  effect = NA_real_
)
print(permutations)
   permutation effect
1            1     NA
2            2     NA
3            3     NA
4            4     NA
5            5     NA
6            6     NA
7            7     NA
8            8     NA
9            9     NA
10          10     NA
11          11     NA
12          12     NA
13          13     NA
14          14     NA
15          15     NA
16          16     NA
17          17     NA
18          18     NA
19          19     NA
20          20     NA

This is simply a frame containing 20 lines where permutation increases from 1 to 20, with an empty variable effect that we can fill in with our effects. Now, let’s begin permuting our effect, where we will generate a variable Wperm which indicates treatment assignment for each permutation. We loop over each column of the combinations matrix, assign treatment to the selected units, estimate the permuted treatment effect via regression, and store it in our permutations frame. In each iteration we also print out which units are receiving treatment so that we can convince ourselves that we are indeed seeing the 20 treatment combinations mentioned in the table above.

combs <- combn(1:6, 3)

effect <- numeric(ncol(combs))

# Loop through permutations
for (j in seq_len(ncol(combs))) {
  
  treated_indices <- combs[, j]
  
  
  Wperm <- rep(0, 6)
  Wperm[treated_indices] <- 1
  
  # Regression
  model <- lm(Y ~ Wperm, data = data.frame(Y = data$Y, Wperm = Wperm))
  effect[j] <- coef(model)["Wperm"]

  cat("Treated units are:", paste(treated_indices, collapse = ", "), "\n")
}
Treated units are: 1, 2, 3 
Treated units are: 1, 2, 4 
Treated units are: 1, 2, 5 
Treated units are: 1, 2, 6 
Treated units are: 1, 3, 4 
Treated units are: 1, 3, 5 
Treated units are: 1, 3, 6 
Treated units are: 1, 4, 5 
Treated units are: 1, 4, 6 
Treated units are: 1, 5, 6 
Treated units are: 2, 3, 4 
Treated units are: 2, 3, 5 
Treated units are: 2, 3, 6 
Treated units are: 2, 4, 5 
Treated units are: 2, 4, 6 
Treated units are: 2, 5, 6 
Treated units are: 3, 4, 5 
Treated units are: 3, 4, 6 
Treated units are: 3, 5, 6 
Treated units are: 4, 5, 6 
permutations <- data.frame(
  permutation = seq_along(effect),
  effect = effect
)

In the code above, within each loop we have estimated our permuted treatment effect, and placed it in the permutations frame. We can now open up this frame, and compare these permutation-based effects to the original effect (which is the first effect stored in our data), and calculate the one-sided and two-sided p-values. We do this below.

effect1 <- permutations$effect[1]

# one-sided test
p1side <- mean(permutations$effect >= effect1)

# two sided test
p2side <- mean(abs(permutations$effect) >= abs(effect1))

# Show p-values
cat("One sided p-value:", p1side, "\n")
One sided p-value: 0.05 
cat("Two sided p-value:", p2side, "\n")
Two sided p-value: 0.1 
library(ggplot2)

ggplot(permutations, aes(x = effect)) +
  geom_histogram(bins = 10, fill = "lightblue", color = "black") +
  geom_vline(xintercept = 11.333, color = "red", linewidth = 1, linetype = "solid") +
  labs(
    title = "Permutation Distribution of Treatment Effects",
    x = "Effect",
    y = "Count"
  ) +
  theme_minimal()

Here we can confirm that we find the p-values we calculated by hand above, and view these permutations graphically in the resulting histogram above.

A Real Example

It turns out, that to translate this to a case with real data, we do not need to add many additional elements. In this section, we work with data from a randomized control trial that examines asset transfers to poor households in India, as discussed in the paper by Banerjee, Duflo, and Sharma (2021). Here we will consider the impact of receipt of asset transfers on two outcomes: a financial index ind_fin_el1 and an asset index asset_ind_tot_el1.

Before we explore randomisation inference, let’s open these data and estimate a simple regression, considering the effect of treatment receipt on the asset index:

library(readr)
library(lmtest)
library(sandwich)

data <- read_csv("data/Banerjee_et_al_2021.csv")
Rows: 991 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (8): treatment, el1, index_ctotal_el1, index_foodsecurity_el1, ind_fin_e...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
model <- lm(asset_ind_tot_el1 ~ treatment, data = data)

coeftest(model, vcov = vcovHC(model, type = "HC1"))

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -0.196104   0.075692 -2.5908 0.0097802 ** 
treatment    0.407753   0.117979  3.4561 0.0005821 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
assetEffect <- coef(model)["treatment"]

print(assetEffect)
treatment 
0.4077528

From the summary we can note that \(\widehat{\tau}_{ATE} = 0.408\) and its standard error equals 0.118. This suggests a p-value of around 0.0005. We have saved the original estimated effect in a local (assetEffect) which we will use below.

We can similarly consider the financial index:

model_fin <- lm(ind_fin_el1 ~ treatment, data = data)

coeftest(model_fin, vcov = vcovHC(model_fin, type = "HC1"))

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  0.1354387  0.0295983  4.5759 5.482e-06 ***
treatment   -0.0095893  0.0415171 -0.2310    0.8174    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
finEffect <- coef(model_fin)["treatment"]

print(finEffect)
   treatment 
-0.009589323

where we observe a point estimate quite close to zero, and p-value of 0.817.

Let’s now consider how we conduct permutation inference in this particular case where there are too many permutations of treatment to feasibly consider them all. Here, rather than considering all possible permutations, we will consider randomly assigned permutations of treatment which respect the original treatment assignment structure. Namely, if we inspect these data, we can see that there are 991 units, 525 of whom receive treatment. We can thus conduct a single permutation by randomly re-assigning treatment to 525 units, and control to the remaning units. We do this below, estimating the correspnding effect.

table(data$treatment)

  0   1 
466 525 
set.seed(121316)
data$epsilon <- rnorm(nrow(data))
data <- data[order(data$epsilon), ]

data$Wperm <- c(rep(1, 525), rep(0, nrow(data) - 525))

model_perm <- lm(asset_ind_tot_el1 ~ Wperm, data = data)
coeftest(model_perm, vcov = vcovHC(model_perm, type = "HC1"))

t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.067511   0.089229  0.7566   0.4495
Wperm       -0.130803   0.118640 -1.1025   0.2706
data$Wperm <- NULL
data$epsilon <- NULL

Here we see that, perhaps unsurprisingly, our randomly assigned treatment results in a small estimated “effect”, of 0.082, which is substantially smaller than the original treatment effect itself.

To estimate a permutation-based p-value all we need to do is repeat this procedure a large number of times (below 10,000 times), and count the proportion of permuted effects which are larger than our estimated treatment effect. We set this up below, incrementing pval by one any time the absolute value of the permuted effect exceeds the absolute value of our observed treatment effect:

Npermut <- 10000
pval <- 0

set.seed(121316)

for (i in 1:Npermut) {
  epsilon <- rnorm(nrow(data))
  idx <- order(epsilon)
  Wperm <- numeric(nrow(data))
  # N = 525 treated units hardcoded in.  To generalise, calculate from original data
  Wperm[idx[1:525]] <- 1
  
  model_i <- lm(asset_ind_tot_el1 ~ Wperm, data = data.frame(asset_ind_tot_el1 = data$asset_ind_tot_el1, Wperm = Wperm))
  
  if (abs(coef(model_i)["Wperm"]) > abs(assetEffect)) {
    pval <- pval + 1
  }
}

cat(pval, "\n")
5 
cat(pval / Npermut, "\n")
5e-04

Here we see that among all of our 10,000 placebos, only 5 effects were larger than the effect we actually observed in our data, suggesting a p-value of 5/10000, or 0.0005.

We can generalise this idea below (for this specific data set-up), and below we do this defining a program which we call randomisation_inference. This program takes a single argument, which is the outcome we wish to consider, and packages the procedure we have explored above, finally printing the p-value and a number of other pieces of information. While we could further optimise this program, for example to automatically infer how to permute the treatment assignment and to return resulting p-values to users, the below is sufficient for our interests here of estimating a randomisation-based p-value.

randomisation_inference <- function(varname, data, Npermut = 10000) {
  formula_str <- as.formula(paste(varname, "~ treatment"))
  model_obs <- lm(formula_str, data = data)
  obs_effect <- coef(model_obs)["treatment"]

  pval <- 0
  set.seed(121316)

  for (i in 1:Npermut) {
    epsilon <- rnorm(nrow(data))
    idx <- order(epsilon)
    Wperm <- numeric(nrow(data))
    Wperm[idx[1:525]] <- 1

    model_perm <- lm(data[[varname]] ~ Wperm)
    if (abs(coef(model_perm)["Wperm"]) > abs(obs_effect)) {
      pval <- pval + 1
    }
  }

  p_value <- pval / Npermut

  cat("The observed effect of", varname, "is", obs_effect, "and its p-value is", p_value, "\n")
}

Finally, we can run this program, simply passing as an argument the two variables we have examined above. We do this below:

randomisation_inference("ind_fin_el1", data)
The observed effect of ind_fin_el1 is -0.009589323 and its p-value is 0.8171 
randomisation_inference("asset_ind_tot_el1", data)
The observed effect of asset_ind_tot_el1 is 0.4077528 and its p-value is 5e-04

If we compare these p-values to regression-based estimates, we will see that, perhaps unsurpringly, they are very similar.

Code Call-out 2.3: Bootstrap

For an introduction to the bootstrap, refer to Section 2.4.4 of the book. Here we will examine a computational implementation of a bootstrap standard error and confidence intervals, using data from Chong et al. (2016). First we will consider the main regression \[Y_i = \mu + \tau_{ATE}W_i + \varepsilon_i.\] Here we will work with the data Chong_et_al_2016, and examine the total cognitive score (wii_total) as our outcome measure \(Y_i\). We can load the data and examine selected variables on the first few lines of this dataset, before estimating the treatment effect via regression below.

library(haven)
library(dplyr)

data <- read_dta("data/Chong_et_al_2016.dta")

head(select(data, student_id, hhid, treat), 5)
# A tibble: 5 × 3
  student_id  hhid treat
  <chr>      <dbl> <dbl>
1 001001         1     0
2 001002         1     0
3 001003         1     0
4 002001         2     0
5 002002         2     1
data_clean <- data %>%
  filter(!is.na(wii_total) & !is.na(treat))

model <- lm(wii_total ~ treat, data = data_clean)

summary(model)

Call:
lm(formula = wii_total ~ treat, data = data_clean)

Residuals:
    Min      1Q  Median      3Q     Max 
-432.66 -103.91  -12.66  134.99  389.34 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   455.66      14.76  30.874   <2e-16 ***
treat          57.46      25.26   2.275    0.024 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 172.8 on 206 degrees of freedom
Multiple R-squared:  0.0245,    Adjusted R-squared:  0.01977 
F-statistic: 5.174 on 1 and 206 DF,  p-value: 0.02395

From the summary we can note that \(\widehat{\tau}_{ATE} = 57.46\) and its standard error equals 25.26. This suggests a p-value of 0.024, and 95% confidence intervals of [7.66;107.267]. In this case the standard errors are computed via the standard OLS variance estimator, assuming homoscedastic errors, although heteroscedasticity-robust standard errors can be requested quite simply using coeftest, and requesting a HC1 or HC2 variance-coviariance matrix. This is shown below, resulting in slight variations to the reported standard error (and correspondingly, p-value and resulting confidence intervals).

library(lmtest)
library(sandwich)

model <- lm(wii_total ~ treat, data = data_clean)

coeftest(model, vcov = vcovHC(model, type = "HC2"))

t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  455.664     15.397 29.5948  < 2e-16 ***
treat         57.463     24.193  2.3752  0.01846 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

An alternative to these closed-form methods for variance estimation is to use bootstrap resampling methods. There are multiple ways a bootstrap can be implemented, but the simplest here will be to simply conduct a paired bootstrap, which is also robust to heteroscedasticity. The pairs bootstrap consists of, first, sampling with replacement the ordered pairs \(\left\{(y_i,w_i)\right\}_{i=1}^{N} = \left\{(y_1,w_1),\ldots,(y_N,w_N)\right\}\) from original data, obtaining a “new” dataset of \(N\) resampled pairs \(\left\{(y_i^*,w_i^*)\right\}_{i=1}^{N} = \left\{(y_1^*,w_1^*),\ldots,(y_N^*,w_N^*)\right\}\). The new dataset thus simply consists of (potentially repeated) randomly selected rows of the original data. Let’s see what this looks like with a single bootstrap replicate. First, we will generate a vector of size \(N\), drawn with replacement. We do this below, having a look at the first 5 rows of the vector.

set.seed(121316)

data_bootstrap <- data_clean %>%
  slice_sample(n = nrow(data_clean), replace = TRUE)

head(select(data_bootstrap, student_id, hhid, treat), 5)
# A tibble: 5 × 3
  student_id  hhid treat
  <chr>      <dbl> <dbl>
1 162001       162     1
2 088001        88     0
3 110003       110     1
4 144001       144     0
5 086001        86     1

With this bootstrap sample we estimate again the coefficient of interest with standard OLS procedure:

model_bootstrap <- lm(wii_total ~ treat, data = data_bootstrap)
summary(model_bootstrap)

Call:
lm(formula = wii_total ~ treat, data = data_bootstrap)

Residuals:
    Min      1Q  Median      3Q     Max 
-425.56  -97.97   -2.64  131.57  388.44 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   456.56      13.90  32.845  < 2e-16 ***
treat          65.16      23.15   2.815  0.00535 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 160.3 on 206 degrees of freedom
Multiple R-squared:  0.03704,   Adjusted R-squared:  0.03237 
F-statistic: 7.924 on 1 and 206 DF,  p-value: 0.005351

With this sample, we have obtained an estimate of \(\widehat{\tau}_{ATE}^* =23.41\); a different value to the original value, in line with variation in the sample used for estimation. If we want some idea of how much estimates vary as the sample changes (ie the sampling variation of our estimate), the next step simply consists of repeating this bootstrap process \(B\) times, resulting in \(B\) estimates of \(\widehat{\tau}_{ATE,(b)}^*\) with \(b\in\{1,\ldots,B\}\). Let’s set \(B=1,000\), and do this:

set.seed(121316)

taus <- numeric(1000)

for (i in 1:1000) {
  data_bootstrap <- data_clean %>%
    slice_sample(n = nrow(data_clean), replace = TRUE)
  
  model <- lm(wii_total ~ treat, data = data_bootstrap)
  
  taus[i] <- coef(model)["treat"]
}

You may wish to step through each of the lines in this loop to ensure that you can see what is going on, but the end product of this code is a column matrix called taus which has been filled in with \(B\) estimates of \(\tau\), \(\left\{\widehat{\tau}^{(b)}\right\}_{b=1}^{B}\). The next step is to convert this column matrix into a variable in our dataset

taus_df <- data.frame(taus = taus)

We can have a look at the entire empirical distribution of these estimates as follows:

library(ggplot2)

p25 <- quantile(taus_df$taus, 0.025)
p975 <- quantile(taus_df$taus, 0.975)

ggplot(taus_df, aes(x = taus)) +
  geom_density(fill = "lightblue", alpha = 0.5) +
  geom_vline(xintercept = c(p25, p975), color = "red", linetype = "dashed") +
  geom_vline(xintercept = 57.46253, color = "grey40", linetype = "solid", linewidth = 1) +
  labs(
    title = expression(tau ~ " Empirical Distribution"),
    x = expression(tau),
    y = "Density"
  ) +
  scale_x_continuous(labels = scales::number_format(accuracy = 0.01)) +
  theme_minimal()

As we can see, this empirical distribution is close to centred on the original estimate (in the limit, they will be exactly the same, which is something you may wish to confirm by using a larger value for \(B\)), while also giving us some idea of the variation of the estimate over alternative resamples. We additionally plot the empirical 2.5th and 97.5th quantiles of the distribution, which provides a 95% confidence interval for \(\widehat\tau\). Clearly, and in line with the regression results observed above, we can reject the null that \(\tau=0\) with some certainty. Using these values \(\left\{\widehat{\tau}^{(b)}\right\}_{b=1}^{B}\), finally we display the estimate’s standard error as the standard deviation of the bootstrap estimates, and the confidence intervals can be generated from empirical quantiles. We illustrate this below:

# Var
var_estimate <- var(taus_df$taus)
sd_estimate <- sd(taus_df$taus)

# Percentils 2.5% y 97.5%
ci <- quantile(taus_df$taus, c(0.025, 0.975))

cat(
  "The variance estimate is:", round(var_estimate, 4), 
  " as such, the standard error estimate is:", round(sd_estimate, 4), "\n"
)
The variance estimate is: 567.7202  as such, the standard error estimate is: 23.8269 
cat(
  "The empirical 95% confidence interval is: [",
  round(ci[1], 4), ",", round(ci[2], 4), "].\n"
)
The empirical 95% confidence interval is: [ 12.4403 , 104.096 ].

We can see that the standard error, and hence 95% confidence interval, is very closely alligned to that from regression documented above, given both are valid under broadly similar assumptions.

References

Banerjee, Abhijit, Esther Duflo, and Garima Sharma. 2021. “Long-Term Effects of the Targeting the Ultra Poor Program.” American Economic Review: Insights 3 (4): 471–86. https://doi.org/10.1257/aeri.20200667.
Bari, Faisal, Kashif Malik, Muhammad Meki, and Simon Quinn. 2024. “Asset-Based Microfinance for Microenterprises: Evidence from Pakistan.” American Economic Review 114 (2): 534–74. https://doi.org/10.1257/aer.20210169.
Chong, Alberto, Isabelle Cohen, Erica Field, Eduardo Nakasone, and Maximo Torero. 2016. “Iron Deficiency and Schooling Attainment in Peru.” American Economic Journal: Applied Economics 8 (4): 222–55. https://doi.org/10.1257/app.20140494.