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 equivalence 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:
clearallsetmoreoffuse"data/Bari_et_al_2024.dta", clearreg biz_ta A if wave==1, robustscalar coef_treatment = _b[A]
Linear regression Number of obs = 739
F(1, 737) = 14.28
Prob > F = 0.0002
R-squared = 0.0192
Root MSE = 1712.5
------------------------------------------------------------------------------
| Robust
biz_ta | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
A | 507.0368 134.1654 3.78 0.000 243.645 770.4287
_cons | 1149.011 109.9867 10.45 0.000 933.0865 1364.936
------------------------------------------------------------------------------
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.
// Calculate mean for the treatment groupsummarize biz_ta if A==1 & wave==1scalar mean_treatment = r(mean)// Calculate mean for the control groupsummarize biz_ta if A==0 & wave==1scalar mean_control = r(mean)// Calculate and display the difference in meansscalar diff_means = mean_treatment - mean_controldisplay"Coefficient for treatment in regression: " coef_treatmentdisplay"Difference in means (Treatment - Control): " diff_means// Display the comparisondi"The coefficient from the regression should be equal to the difference in means to demonstrate equivalence."
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
biz_ta | 491 1656.048 1701.928 0 7666.667
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
biz_ta | 248 1149.011 1733.225 0 7666.667
Coefficient for treatment in regression: 507.03681
Difference in means (Treatment - Control): 507.03681
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:
//Generate baseline valuegen temp = mpall if wave==0bys id: egenmanage = mean(temp)//Standardise so mean-zerosummanagereplacemanage = manage-r(mean)//Interact variable with treatmentgen Amanage = A*manage
(3,608 missing values generated)
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
manage | 4,365 .0263182 .7326993 -1.050634 2.332034
(4,365 real changes made)
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:
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.
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 Stata. First, we will load our “data”:
clearinput Y W34 127 129 114 018 024 0end
and then calculate a comparison of means estimator of the ATT:
sum Y if W == 1local mean_Y1 = r(mean)sum Y if W == 0local mean_Y0 = r(mean)// Treatment effectlocal tau_hat = `mean_Y1' - `mean_Y0'dis "Treatment effect is:"`tau_hat'
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
Y | 3 30 3.605551 27 34
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
Y | 3 18.66667 5.033223 14 24
Treatment effect is:11.333333
Now let’s generate our p-value by permuting all possible treatment combinations. While there are many ways we could consider doing this, one is to note that we could start our permutation by assigning units 1, 2 and 3 to treatment and the last three units to control, and subsequently increase the units assigned to treatment until we are assigning units 4, 5 and 6 to treatment, and the remaining units to control. Below we begin by preparing a new frame of data in which to store our resulting treatment effects:
Number of observations (_N) was 0, now 20.
(20 missing values generated)
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 is our new treated unit. The below code will iterate through the three possible units receiving treatment (t1, t2 and t3), starting from units 1, 2 and 3, and going until these reach units 3, 4 and 5 respectively. It is worth reading through this code carefully given that this is done “manually”, but in each loop we will print out the units receiving treatment so that we can convince ourselves that we are indeed seeing the 20 treatment combinations mentioned in the table above.
gen Wperm = .local t1 = 1local t2 = 2local t3 = 3local j=1while`t1'<=4 {while`t2'<=5 {while`t3'<= 6 { dis "Treated units are `t1', `t2', `t3'"quireplace Wperm = 0quireplace Wperm = 1 in`t1'quireplace Wperm = 1 in`t2'quireplace Wperm = 1 in`t3'quireg Y Wperm// Store permuted effect in permutation framequi frame permutations: replace effect = _b[Wperm] in`j'// Increase to next permutationlocal ++j// Increase so unit 3 goes up by one unitlocal ++t3 }// We have now reached unit 6 for the third treated unit, so increase second treated unit and start againlocal ++t2local t3=`t2'+1 }// We have now reached unit 5 for 2nd unit and 6 for the third treated unit// so increase first treated unit and start again local ++t1local t2=`t1'+1local t3=`t2'+1}
(6 missing values generated)
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
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.
// Change to permutations framecwf permutationslist// Calculate a one sided p-valuecountif effect >= effect[1]local p1side = `r(N)' / 20// Calculate a two sided p-valuecountifabs(effect) >= effect[1]local p2side = `r(N)' / 20display"One sided p-value: "`p1side'display"Two sided p-value: "`p2side'histogram effect, bins(10) xline(11.333, lcolor(red) lwidth(thick))
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:
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.
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 remaining units. We do this below, estimating the corresponding effect.
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, using the ++ notation to increase our local pval by one any time the condition on line 12 is met:
local Npermut = 10000local pval = 0forvalues i=1/`Npermut' {qui {gen epsilon = rnormal()sort epsilongenerate Wperm = 1 in 1/525replace Wperm = 0 in 526/991regress asset_ind_tot_el1 Wperm }// Increase pval by 1 if permuted effect is larger than true effectifabs(_b[Wperm])>abs(`assetEffect') local ++pvaldrop epsilon Wperm}dis `pval'dis `pval'/`Npermut'
8
.0008
Here we see that among all of our 10,000 placebos, only 8 effects were larger than the effect we actually observed in our data, suggesting a p-value of 8/10000, or 0.0008.
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.
program randomisation_inferenceargs varnamequi: reg`varname' treatmentlocal obs_effect = _b[treatment]// Set number of permutationslocal Npermut = 10000local pval = 0forvalues i = 1/`Npermut' {// Generate a permuted treatment variablequi {gen epsilon = rnormal()sort epsilongenerate Wperm = 1 in 1/525replace Wperm = 0 in 526/991regress`varname' Wperm }// Increase pval by 1 if permuted effect is larger than true effectifabs(_b[Wperm])>abs(`obs_effect') local ++pvaldrop epsilon Wperm }// Calculate p-valuelocal p_value = `pval' / `Npermut'display"The observed effect of ""`varname'"" is "`obs_effect'" and its p-value is "`p_value'end
Finally, we can run this program, simply passing as an argument the two variables we have examined above. We do this below:
The observed effect of ind_fin_el1 is -.00958932 and its p-value is .8245
The observed effect of asset_ind_tot_el1 is .40775282 and its p-value is .0004
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.
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 the robust option from the command regress This is shown below, resulting in slight variations to the reported standard error (and correspondingly, p-value and resulting confidence intervals).
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 this new dataset, having a look at the first 5 rows of selected variables in the new dataset.
preservesetseed 121316bsamplelist student_id hhid treat in 1/5
With this bootstrap sample we estimate again the coefficient of interest with standard OLS procedure:
reg wii_total treatrestore
(Note: Below code run with echo to enable preserve/restore functionality.)
Source | SS df MS Number of obs = 208
-------------+---------------------------------- F(1, 206) = 0.81
Model | 25627.5269 1 25627.5269 Prob > F = 0.3680
Residual | 6485700.15 206 31483.9813 R-squared = 0.0039
-------------+---------------------------------- Adj R-squared = -0.0009
Total | 6511327.67 207 31455.6892 Root MSE = 177.44
------------------------------------------------------------------------------
wii_total | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
treat | 23.40968 25.94699 0.90 0.368 -27.74602 74.56539
_cons | 494.365 15.15949 32.61 0.000 464.4773 524.2526
------------------------------------------------------------------------------
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:
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
clearsvmat taus, names("taus")
number of observations will be reset to 1000
Press any key to continue, or Break to abort
Number of observations (_N) was 0, now 1,000.
We can have a look at the entire empirical distribution of these estimates as follows:
As we can see, this empirical distribution is close to centered 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:
quisum taus1diastext"The variance estimate is: "as result round(r(Var), 0.0001) ///astext" as such, the standard error estimates is: "as result round(r(sd), 0.0001)_pctile taus1, percentiles(2.5 97.5)diastext"The empirical 95% confidence interval is: ["///as result round(r(r1), 0.0001) astext","///as result round(r(r2), 0.0001) astext"]."
The variance estimate is: 599.7058 as such, the standard error estimates is: 24
> .4889
The empirical 95% confidence interval is: [10.7098,104.6838].
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.
