Effectively are you able to consider it? I caught Covid. There’s a pressure going round Boston and I drew the quick straw. Forgot how nice this was — pounding complications, eyes that really feel like they’re burning from the within out, congestion and scratchy throat, and fever. So what higher factor to do than write a substack about Lalonde and diff in diff and covariates. So right here’s the take a look at.
Lalonde famously confirmed that many econometric evaluators had been biased, of their software not of their properties and estimators. An estimator could be unbiased and but be misapplied and thus biased. And he confirmed cleverly that that was the case in his 1986 job market paper.
He took an RCT, first, and in contrast 1978 earnings for a handled and management group. That distinction was $800. Which means that the therapy — a job coaching program aimed toward poor staff — raised the true earnings of staff by $800.
Then to point out the issue going through the present econometric evaluators, he dropped the management group — which was recall exchangeable with the therapy teams counterfactual earnings because of the randomization — and changed it with that of survey information of Individuals, notably the CPS and the PSID. He then used a gamut of econometric estimators to estimate the impact of this system which was nonetheless, expressed as an ATT, $800. And his estimates didn’t get well it. In reality estimates ranged from -16,000 to +3,000 and the few that acquired shut didn’t have any apparent ex ante cause to favor them.
And thus the credibility revolution was born.
I wish to do the identical factor, solely I wish to make three modifications.
-
I wish to use the Dehejia and Wahba (2002) pattern which is a subset of the employees who had two years pre therapy earnings accessible. The impact there may be increased — $1,794. That’ll be our goal.
-
We are going to use distinction in variations with the CPS pattern for our non random comparability group. The goal parameter there may be the ATT and the impact remains to be $1,794 as a result of the ATT is the impact on the handled staff, and the handled staff remained within the pattern (simply not the unique management staff). It is a delicate level possibly however it’s price saying — the ATE and the ATT had been underneath randomization $1,794 as a result of underneath randomizing, they’re the identical. However the ATE modifications once we exchange the experimental management with the non experimental management group — it’s simply the ATT doesn’t.
-
We are going to use covariates with the diff in diff estimators however underneath totally different specs, together with the 2 approach fastened results mannequin with additive covariates which is the most typical most likely of all utilized panel estimators.
I’m going to evaluation six estimators, all diff-in-diff. And bear in mind — the impact is $1,794. It’s recognized. It’s the bottom fact. So we decide efficiency in opposition to it. And my competition in that is to point out you {that a}) covariates can matter loads, even with diff in diff and b) that doesn’t subsequently imply that each one methods of introducing Covariates right into a diff in diff mannequin are equal.
Recall from our JEL that there are 4 separate regressions which are numerically the identical for calculating the estimated OLS coefficient. They’re saturated dummies the place you could have a therapy dummy, a submit dummy and an interplay. The interplay is the numerical equal of calculating 4 averages and three subtractions. Then there may be the inclusion of employee and 12 months fastened results with an interplay of submit and therapy. Similar precise calculation on the interplay. Then there’s a first distinction regression onto a therapy dummy. Similar factor. After which there may be an across-comparison group regression onto a submit dummy.
-
Interplay OLS with out covariates
Since they’re all the identical we use the one that may let me later embody time invariant covariates: the saturated regression. And once we do that, we get $3,621. That is constructive however it’s two occasions too massive.
-
Interplay OLS with additive (however not saturated) covariates
Subsequent we run the everyday diff in diff of merely together with covariates additively. This additionally yields $3,621. Word that I didn’t use TWFE, which might have eradicated the time invariant variables totally by way of the demeaning course of, however because the TWFE and the interplay of all dummies is similar numerically, it doesn’t matter. Once we do that, we get this 3,621.
However discover that it’s the similar quantity. Why? It is because they don’t change over time. For the reason that TWFE and saturated regression are numerically equivalent, then even when TWFE deletes the time invariant covariates through employee fastened results, the precise did coefficient is unchanged as a result of the estimate by no means trusted them anyway.
-
Regression adjustment I: work together the covariates with time
The primary adjustment we’ll contemplate will handle the issue of a misspecified regression wherein covariate tendencies don’t change over time. What is that this although?
Let’s say in our pattern there’s two forms of staff for simplicity. These with a highschool diploma, and people with out. Covariates tendencies asks if the highschool educated staff are on totally different Y(0) tendencies than those that don’t. If they’re, then the above regression I did is inaccurate insofar because the therapy group has extra of 1 than the opposite. And it does — our management group is from the CPS, a consultant pattern of America. The common American had extra education than the people within the MDRC job coaching program. So, if we predict their tendencies are totally different, then the imbalance in education will violate the parallel tendencies bias time period.
Conditional parallel tendencies will fail because of the imbalance however we are able to appropriate at the least a part of it with partial saturation, which is to work together the 12 months dummies with all of the covariates.
Once we work together all of the covariates that Dehejia and Wahba use, curiously our estimate drops to $1,711, which is nearly precisely the true ATT of $1,794. This means that covariate imbalance is a nontrivial drawback with these information and that if the teams with totally different covariate mixtures are on totally different earnings tendencies, Delta Y(0), then the imbalance breaks conditional parallel tendencies mechanically. It’s a must to work together with X to repair it, not merely embody them additively, as a result of it is advisable mannequin these tendencies explicitly first.
-
Regression adjustment II: heterogenous therapy results and saturated regression
However imbalance is simply probably a part of the issue with specification #2. It additionally imposes fixed therapy results with respect to the covariates themselves. Utilizing our education instance, however permitting education to enter in additively, the therapy results are the identical for highschool and no highschool educated staff.
We appropriate this by way of saturation which is to dummy out each worth of the covariates, work together it with the therapy, and run that regression. Then you definitely get well the ATT by summing over all of the therapy indicator coefficients (each it alone and the interactions) multiplied the identical imply of the covariate within the therapy group.
Curiously, this did nothing to our estimates. We discover utilizing this type of interplay the identical as what we discovered utilizing the additive technique — $3,621. Why? As a result of saturating on a regression involving ranges did nothing to handle the conditional parallel tendencies assumption, which is in tendencies, not ranges. Till we try this, we’re going to proceed to endure from a violation of conditional parallel tendencies.
-
Regression adjustment III: saturated in time and covariates
So, now we’ll do each. That’s, we’ll work together time with the covariates, in addition to work together covariates with the therapy. After which we’ll get well the ATT by summing the therapy coefficients (together with the interactions) occasions the pattern imply of the covariates within the therapy group solely. In diff-in-diff, that is usually referred to as final result regression, primarily based on Heckman, Ichimura and Todd (1997, Restud). However right here I do it manually, and to do it manually, I first take first variations in earnings, then I regress that onto the therapy dummy, every covariate by itself, and the therapy interacted with every covariate. There’s no submit dummy right here although. The primary-differencing eliminated it.
Right here’s the delicate half: as a result of the end result is now a change, a covariate’s personal coefficient is not a degree impact. Moderately, it’s the covariate’s impact on the development. That most important impact is the X-by-time interplay from #3, free of charge. So the first-difference regression does two jobs without delay: the covariate most important results bend the management group’s counterfactual development (what #3 did), and the therapy interactions loosen up fixed therapy results (what #4 tried to do). That’s why it lands at $1,770 whereas the pure-levels model in #4 sat at $3,621. As a result of the degrees saturation by no means touched the development and thus the imbalance turned problematic.
-
Heckman, Ichimura and Todd (1997)
However guess what. The Heckman, Ichimura and Todd (1997) estimator which strikes in a couple of levels is similar because the regression adjustment “each” we simply did. The levels are:
-
Regress the primary distinction, Delta Y, onto X for management group solely
-
Acquire the coefficients on X, transport them to a regression mannequin aimed on the therapy group solely, and predict Delta Y for the therapy teams solely
-
Calculate particular person therapy results because the therapy group first distinction minus the imputed counterfactual first distinction for the therapy group solely after which take its common over the therapy group.
That’s what I discover so fascinating and never terribly intuitive. That whenever you try this totally interacted X on a primary distinction regression for all models (therapy and management), it’s numerically the identical as should you had gone by way of a collection of steps beginning with the management group solely, then predicting on the therapy group solely. That is, although, what the Oaxaca-Blinder estimator had accomplished which is itself too a therapy impact estimator for the ATT, the ATE and the ATU relying on the way you wish to go about utilizing these interactions. Nonetheless, doing this we get $1,770.
-
Propensity scores and double strong
I may also introduce the covariates through an inverted propensity scores reweighting process developed by Alberto Abadie (2005, Restud), or I can introduce each regression adjustment and the propensity scores reweighting process in what’s referred to as double strong developed by Sant’Anna and Zhao (2020). I’ll clarify these two later, however for now I’ll simply present it to you, and the code shall be under on the finish of this substack.
So let’s summarize this. Initially, the Lalonde information is the reward that retains on giving. We all know the ATT is $1,794, and we subsequently can convincingly present (for my part anyway) that the exclusion of covariates in a diff-in-diff could be the incorrect spec. Moreover, simply because together with covariates in diff in diff modified the coefficients doesn’t imply the unique was incorrect, as a result of covariates are usually not a robustness train. In the event that they had been robustness train then sure — we’d anticipate them to not matter. But when covariates carry out an precise perform, which is satisfying conditional parallel tendencies and enjoyable fixed therapy results assumption, then naturally they each will matter underneath imbalance and heterogenous therapy results and must be included.
However that’s the second level I wish to go away you with. Simply because they’re wanted doesn’t subsequently imply that any specification is equal. As we noticed right here, essentially the most versatile ones recovered the true ATT or acquired very shut, however the TWFE one even with covariates didn’t.
There’s extra to this than I cowl right here — like why use regression adjustment? Why not use the propensity rating technique that Abadie developed? What precisely makes them totally different from each other? What precisely makes another interesting than one other? That’s for a distinct day. For now, I simply needed to get this accomplished.
Stata code
********************************************************************************
* identify: lalonde_all_specs.do
* objective: Run each estimator from the deck on the LaLonde-DW non-experimental
* panel and evaluate to the RCT benchmark (~$1,794).
*
* Estimators:
* Spec 0: Naive TWFE (no covariates) -> 3,621
* Spec A: Additive X (X within the degree; time-invariant) -> 3,621 (inert)
* Spec BT: X × therapy (ranges; X within the impact) -> 3,621 (inert)
* Spec B: Put up × X (X × time; X within the development) -> 1,711 (corrects)
* Spec C: Absolutely Saturated TWFE (FD with D × X) = HIT -> 1,770 (each)
* OR by hand (HIT 1997): FD regression on controls solely, impute
* OR (DRDID regadj)
* IPW (DRDID, Abadie 2005)
* DR (DRDID dripw, Sant'Anna-Zhao 2020)
*
* Information: lalonde_nonexp_panel.dta — NSW handled + CPS controls, years 74/75/78.
* We use years 75 (pre) and 78 (submit) for the 2-period DiD.
* Covariates: time-invariant baseline traits.
********************************************************************************
* Run this do-file from the folder that holds lalonde_nonexp_panel.dta
* (all paths under are relative to the working listing).
clear all
seize log shut
set extra off
log utilizing "lalonde_all_specs.log", exchange textual content
use lalonde_nonexp_panel.dta, clear
* Prohibit to 75 (pre) and 78 (submit) for a 2-period DiD
preserve if 12 months == 75 | 12 months == 78
gen submit = (12 months == 78)
* "deal with" within the panel = ever_treated × submit; ever_treated is the group indicator
* Sanity test
show _n "===== Sanity test: cell counts ====="
tab ever_treated submit
* X set matches Scott's canonical LaLonde spec (has agecube; no re75/u75)
world X "age agesq agecube educ educsq marr nodegree black hisp re74 u74"
show _n "================================================================"
show " LaLonde DiD specs — non-experimental panel (NSW + CPS controls)"
show " Pre = 1975, Put up = 1978. RCT benchmark ~$1,794"
show "================================================================"
********************************************************************************
* Spec 0: Naive TWFE (no covariates)
********************************************************************************
show _n "===== Spec 0: Naive TWFE ====="
reg re i.submit##i.ever_treated, strong
native spec0 = _b[1.post#1.ever_treated]
native spec0_se = _se[1.post#1.ever_treated]
show "Spec 0 estimate = " %9.0fc `spec0' " (SE " %6.0fc `spec0_se' ")"
********************************************************************************
* Spec A: Additive X (X enters as linear controls)
********************************************************************************
show _n "===== Spec A: Additive X ====="
reg re i.submit##i.ever_treated $X, strong
native specA = _b[1.post#1.ever_treated]
native specA_se = _se[1.post#1.ever_treated]
show "Spec A estimate = " %9.0fc `specA' " (SE " %6.0fc `specA_se' ")"
********************************************************************************
* Spec B: Put up × X (double-## so X will get each pre and submit results)
********************************************************************************
show _n "===== Spec B: Put up × X ====="
reg re i.submit##i.ever_treated ///
i.submit##c.age i.submit##c.agesq i.submit##c.agecube ///
i.submit##c.educ i.submit##c.educsq ///
i.submit##i.marr i.submit##i.nodegree i.submit##i.black i.submit##i.hisp ///
i.submit##c.re74 i.submit##i.u74, strong
native specB = _b[1.post#1.ever_treated]
native specB_se = _se[1.post#1.ever_treated]
show "Spec B estimate = " %9.0fc `specB' " (SE " %6.0fc `specB_se' ")"
********************************************************************************
* Spec BT: X × TREATMENT, in ranges (heterogeneous therapy results).
* Work together each covariate with the treatment-on change T = submit×D, then
* get well the ATT the textbook approach: coef on T plus every interplay coef occasions
* the treated-group imply of that covariate (= predict at T=1 minus T=0, averaged
* over the treated-post cell).
*
* KEY RESULT: that is INERT. It returns the naive quantity to the greenback
* (3,621), as a result of the covariates enter the therapy EFFECT, by no means the management
* group's counterfactual TREND. Solely X × submit touches the development. SE is
* bootstrapped (margins cannot take a discrete impact right here; extra params than the
* naive spec, so its variance differs).
********************************************************************************
show _n "===== Spec BT: X × therapy (ranges, inert) ====="
seize drop T
gen T = submit*ever_treated
seize program drop _attXT
program outline _attXT, rclass
seize drop _rh1 _rh0 _tauXT
reg re i.submit i.ever_treated i.T ///
c.age c.agesq c.agecube c.educ c.educsq c.re74 ///
i.marr i.nodegree i.black i.hisp i.u74 ///
i.T#c.age i.T#c.agesq i.T#c.agecube i.T#c.educ i.T#c.educsq i.T#c.re74 ///
i.T#i.marr i.T#i.nodegree i.T#i.black i.T#i.hisp i.T#i.u74
tempvar e
gen `e' = T
exchange T = 1
predict _rh1, xb
exchange T = 0
predict _rh0, xb
exchange T = `e'
gen _tauXT = _rh1 - _rh0
sum _tauXT if ever_treated==1 & submit==1, meanonly
return scalar att = r(imply)
drop _rh1 _rh0 _tauXT
finish
quietly _attXT
native specBT = r(att)
set seed 90210
quietly bootstrap att=r(att), reps(199) nodots: _attXT
native specBT_se = _se[att]
show "Spec BT estimate = " %9.0fc `specBT' " (SE " %6.0fc `specBT_se' ")"
seize drop T
********************************************************************************
* Spec C: Absolutely Saturated TWFE (FD with D × X interactions)
* In 2-period panel with unit FE, that is equal to first-differencing,
* then regressing dy on X with full D × X interactions and predicting.
********************************************************************************
show _n "===== Spec C: Absolutely Saturated TWFE (FD) ====="
protect
preserve id ever_treated re age agesq agecube educ educsq marr nodegree black hisp re74 u74 submit
reshape large re, i(id) j(submit)
gen dy = re1 - re0
* Saturated regression on FULL pattern with D × X interactions
reg dy $X ///
i.ever_treated#c.age i.ever_treated#c.agesq i.ever_treated#c.agecube ///
i.ever_treated#c.educ i.ever_treated#c.educsq ///
i.ever_treated#i.marr i.ever_treated#i.nodegree ///
i.ever_treated#i.black i.ever_treated#i.hisp ///
i.ever_treated#c.re74 i.ever_treated#i.u74 ///
i.ever_treated, strong
* ATT: predict dy underneath ever_treated=1 minus dy underneath ever_treated=0
gen evt_orig = ever_treated
exchange ever_treated = 1
quietly predict dy_hat_1, xb
exchange ever_treated = 0
quietly predict dy_hat_0, xb
exchange ever_treated = evt_orig
drop evt_orig
gen tau_sat = dy_hat_1 - dy_hat_0
quietly sum tau_sat if ever_treated == 1
native specC = r(imply)
show "Spec C (Saturated TWFE) estimate = " %9.0fc `specC'
restore
********************************************************************************
* OR by hand (HIT 1997): FD on controls, impute counterfactual for handled
* + cross-check with drdid regadjust
********************************************************************************
show _n "===== OR by hand (HIT 1997, FD on controls solely) ====="
protect
preserve id ever_treated re age agesq agecube educ educsq marr nodegree black hisp re74 u74 submit
reshape large re, i(id) j(submit)
gen dy = re1 - re0
reg dy $X if ever_treated == 0, strong
predict dy_hat
gen delta_hat = dy - dy_hat if ever_treated == 1
quietly sum delta_hat if ever_treated == 1
native or_hand = r(imply)
show "OR by hand estimate = " %9.0fc `or_hand'
restore
show _n "===== OR cross-check: drdid regadjust ====="
seize noisily drdid re $X, time(12 months) ivar(id) tr(ever_treated) reg
if _rc == 0 {
matrix Bor = e(b)
matrix Vor = e(V)
native or_drdid = Bor[1,1]
native or_se = sqrt(Vor[1,1])
show "drdid regadjust = " %9.0fc `or_drdid' " (SE " %6.0fc `or_se' ")"
show "distinction = " %9.0fc `or_hand' - `or_drdid'
}
else {
show as error "drdid regadjust failed (rc=" _rc ")."
native or_drdid = .
native or_se = .
}
********************************************************************************
* IPW (Abadie 2005) — by hand, + cross-check with drdid ipw
********************************************************************************
show _n "===== IPW by hand (Abadie 2005) ====="
protect
preserve id ever_treated re age agesq agecube educ educsq marr nodegree black hisp re74 u74 submit
reshape large re, i(id) j(submit)
gen dy = re1 - re0
* Step 1: propensity rating
logit ever_treated $X
predict phat, pr
* Step 2: Abadie 2005 ATT weights
quietly sum ever_treated
native pr_treat = r(imply)
gen w_ipw = (ever_treated - phat) / (1 - phat) / `pr_treat'
gen ipw_dy = w_ipw * dy
quietly sum ipw_dy
native ipw_byhand = r(imply)
show "IPW by hand estimate = " %9.0fc `ipw_byhand'
* Step 3: DR (Sant'Anna-Zhao 2020)
reg dy $X if ever_treated == 0
predict dy_hat_dr
gen dr_t = ever_treated * (dy - dy_hat_dr) / `pr_treat'
gen dr_c = (1 - ever_treated) * (phat / (1 - phat)) * (dy - dy_hat_dr) / `pr_treat'
quietly sum dr_t
native dr_treated = r(imply)
quietly sum dr_c
native dr_byhand = `dr_treated' - r(imply)
show "DR by hand estimate = " %9.0fc `dr_byhand'
restore
show _n "===== IPW cross-check: drdid ipw ====="
seize noisily drdid re $X, time(12 months) ivar(id) tr(ever_treated) ipw
if _rc == 0 {
matrix Bipw = e(b)
matrix Vipw = e(V)
native ipw_drdid = Bipw[1,1]
native ipw_se = sqrt(Vipw[1,1])
show "drdid ipw = " %9.0fc `ipw_drdid' " (SE " %6.0fc `ipw_se' ")"
show "distinction = " %9.0fc `ipw_byhand' - `ipw_drdid'
}
else {
show as error "drdid ipw failed (rc=" _rc ")."
native ipw_drdid = .
native ipw_se = .
}
show _n "===== DR cross-check: drdid dripw ====="
seize noisily drdid re $X, time(12 months) ivar(id) tr(ever_treated) dripw
if _rc == 0 {
matrix Bdr = e(b)
matrix Vdr = e(V)
native dr_drdid = Bdr[1,1]
native dr_se = sqrt(Vdr[1,1])
show "drdid dripw = " %9.0fc `dr_drdid' " (SE " %6.0fc `dr_se' ")"
show "distinction = " %9.0fc `dr_byhand' - `dr_drdid'
}
else {
show as error "drdid dripw failed (rc=" _rc ")."
native dr_drdid = .
native dr_se = .
}
********************************************************************************
* Abstract
********************************************************************************
show _n
show "================================================================"
show " LaLonde — all estimators side-by-side"
show "================================================================"
show "RCT benchmark (goal) ~ $1,794"
show "Spec 0 Naive TWFE = " %9.0fc `spec0'
show "Spec A Additive X = " %9.0fc `specA'
show "Spec BT X × therapy (ranges) = " %9.0fc `specBT'
show "Spec B Put up × X (X × time) = " %9.0fc `specB'
show "Spec C Saturated TWFE (FD) = HIT = " %9.0fc `specC'
show "OR by hand (HIT 1997, FD) = " %9.0fc `or_hand'
show "OR drdid (regadjust) = " %9.0fc `or_drdid'
show "IPW by hand (Abadie 2005) = " %9.0fc `ipw_byhand'
show "IPW drdid (ipw) = " %9.0fc `ipw_drdid'
show "DR by hand (Sant'Anna-Zhao 2020) = " %9.0fc `dr_byhand'
show "DR drdid (dripw) = " %9.0fc `dr_drdid'
show "================================================================"
********************************************************************************
* Forest plot: the covariate arc, ordered as an argument.
* Three specs that go away covariates OUT of the development all return the naive
* 3,621 (gray). The second covariates enter the development (X × submit, then full
* saturation) the estimate snaps to the $1,794 experimental benchmark (inexperienced).
* Propensity-based estimators (blue) land close by. Grouping = the thesis.
********************************************************************************
protect
clear
set obs 7
gen str32 estimator = ""
gen double att = .
gen double se = .
gen byte ord = .
gen byte grp = . // 1 inert, 2 trend-corrected, 3 propensity
* --- inert: covariates within the degree or the therapy impact (gray) ---
exchange estimator="No covariates" in 1
exchange att=`spec0' in 1
exchange se=`spec0_se' in 1
exchange ord=7 in 1
exchange grp=1 in 1
exchange estimator="Additive X" in 2
exchange att=`specA' in 2
exchange se=`specA_se' in 2
exchange ord=6 in 2
exchange grp=1 in 2
exchange estimator="X × therapy" in 3
exchange att=`specBT' in 3
exchange se=`specBT_se' in 3
exchange ord=5 in 3
exchange grp=1 in 3
* --- trend-corrected: covariates enter the counterfactual development (inexperienced) ---
exchange estimator="X × submit (development)" in 4
exchange att=`specB' in 4
exchange se=`specB_se' in 4
exchange ord=4 in 4
exchange grp=2 in 4
exchange estimator="Saturated = HIT" in 5
exchange att=`specC' in 5
exchange se=`or_se' in 5
exchange ord=3 in 5
exchange grp=2 in 5
* --- propensity-based (blue) ---
* Use the BY-HAND Abadie/SZ factors so Stata and R reproduce identically; the
* drdid/DRDID comfort instructions differ throughout languages on default internals.
* SEs are drdid's analytic values (IPW level equals drdid's precisely; DR inside $14).
exchange estimator="IPW (Abadie)" in 6
exchange att=`ipw_byhand' in 6
exchange se=`ipw_se' in 6
exchange ord=2 in 6
exchange grp=3 in 6
exchange estimator="Doubly strong" in 7
exchange att=`dr_byhand' in 7
exchange se=`dr_se' in 7
exchange ord=1 in 7
exchange grp=3 in 7
gen decrease = att - 1.96*se
gen higher = att + 1.96*se
record estimator att se decrease higher grp, sep(0) noobs
* Gov 2001 palette
native gray "113 128 150"
native inexperienced "39 103 73"
native blue "43 108 176"
native crim "165 28 48"
twoway ///
(rcap higher decrease ord if grp==1, horizontal lcolor("`gray'") lwidth(medthick)) ///
(scatter ord att if grp==1, mcolor("`gray'") msize(massive) msymbol(O)) ///
(rcap higher decrease ord if grp==2, horizontal lcolor("`inexperienced'") lwidth(medthick)) ///
(scatter ord att if grp==2, mcolor("`inexperienced'") msize(massive) msymbol(D)) ///
(rcap higher decrease ord if grp==3, horizontal lcolor("`blue'") lwidth(medthick)) ///
(scatter ord att if grp==3, mcolor("`blue'") msize(massive) msymbol(S)) ///
, ///
xline(1794, lcolor("`crim'") lpattern(sprint) lwidth(medthick)) ///
ylabel(7 "No covariates" 6 "Additive X" 5 "X × therapy" ///
4 "X × submit (development)" 3 "Saturated = HIT" 2 "IPW (Abadie)" ///
1 "Doubly strong", angle(0) labsize(medsmall) nogrid) ///
yscale(vary(0.5 7.6)) ///
xlabel(0 1000 1794 "{bf:1,794}" 3000 4000, format(%9.0fc) labsize(small)) ///
xtitle("ATT estimate ($), 95% CI", measurement(small)) ytitle("") ///
title("Covariates rescue LaLonde solely after they enter the {it:development}", ///
measurement(medium) colour("45 55 72")) ///
legend(order(2 "Within the degree / impact (inert)" 4 "Within the development (corrected)" ///
6 "Propensity-based") rows(1) measurement(small) area(lstyle(none)) ///
place(6)) ///
graphregion(colour(white) margin(medium)) plotregion(colour(white)) bgcolor(white) ///
ysize(4.4) xsize(7.6)
graph export lalonde_forest.png, as(png) exchange width(2000)
show _n "Forest plot saved to lalonde_forest.png"
restore
seize log shut
exit
R code
# ==============================================================================
# lalonde_all_specs.R
# R replication of lalonde_all_specs.do — each DiD estimator on the LaLonde-DW
# non-experimental panel, in comparison with the RCT benchmark (~$1,794).
#
# Spec 0 : Naive TWFE (no covariates) -> 3,621
# Spec A : Additive X (X within the degree) -> 3,621 (inert)
# Spec BT: X x therapy (ranges; X within the impact) -> 3,621 (inert)
# Spec B : Put up x X (X x time; X within the development) -> 1,711 (corrects)
# Spec C : Absolutely saturated FD (D x X) = HIT -> 1,770 (each)
# HIT by hand (control-only FD regression, impute) -> 1,770
# IPW (Abadie 2005) / DR (Sant'Anna-Zhao 2020), by hand -> 1,861 / 1,993
#
# IPW/DR are written out by hand (not delegated to a bundle) so Stata and R
# reproduce identically; the DRDID/drdid comfort instructions differ throughout
# languages on default internals and are printed solely as a cross-check.
# Strong SEs are HC1 (= Stata's ", strong"). Spec BT / IPW / DR SEs are bootstrapped.
# ==============================================================================
suppressMessages({
library(haven); library(sandwich); library(lmtest)
library(DRDID); library(ggplot2); library(dplyr)
})
## ---- run this from the folder holding lalonde_nonexp_panel.dta --------------
## (paths under are relative to the working listing)
xvars <- c("age","agesq","agecube","educ","educsq",
"marr","nodegree","black","hisp","re74","u74")
xf <- reformulate(xvars) # ~ age + agesq + ... + u74
d <- read_dta("lalonde_nonexp_panel.dta") |>
filter(12 months %in% c(75, 78)) |>
mutate(submit = as.integer(12 months == 78),
D = as.integer(ever_treated))
# strong (HC1) ATT + SE for the submit:D interplay in a ranges DiD
did_int <- perform(type, information) D:submit$", rownames(ct))
c(att = ct[ix, 1], se = ct[ix, 2])
## ---- Spec 0: naive -----------------------------------------------------------
s0 <- did_int(re ~ submit * D, d)
## ---- Spec A: additive X ------------------------------------------------------
sA <- did_int(reformulate(c("submit*D", xvars), "re"), d)
## ---- Spec B: X x submit (covariate most important results PLUS every covariate x time) ----
sB <- did_int(reformulate(c("submit*D", xvars, paste0("submit:", xvars)), "re"), d)
## ---- Spec BT: X x therapy in ranges (g-computation over treated-post) ------
# T = submit*D ; work together each covariate with T ; ATT = imply over treated-post of
# [ yhat(T=1) - yhat(T=0) ] = coef(T) + sum(coef(T:X) * Xbar_treated). INERT.
att_XT <- perform(information) {
information$T <- information$submit * information$D
m <- lm(reformulate(c("submit","D","T", xvars, paste0("T:", xvars)), "re"), information)
d1 <- remodel(information, T = 1); d0 <- remodel(information, T = 0)
tau <- predict(m, d1) - predict(m, d0)
imply(tau[data$D == 1 & data$post == 1])
}
sBT_att <- att_XT(d)
set.seed(90210)
ids <- distinctive(d$id)
idx_map <- cut up(seq_len(nrow(d)), d$id) # id -> row indices (as soon as)
bsBT <- replicate(199, { # cluster bootstrap on id
samp <- pattern(ids, exchange = TRUE)
att_XT(d[unlist(idx_map[as.character(samp)], use.names = FALSE), ])
})
sBT <- c(att = sBT_att, se = sd(bsBT))
## ---- Spec C: totally saturated FD (D x X), g-computation -----------------------
w <- d |>
choose(id, D, re, submit, all_of(xvars)) |>
tidyr::pivot_wider(names_from = submit, values_from = re,
names_prefix = "re") |>
mutate(dy = re1 - re0)
mC <- lm(reformulate(c(xvars, paste0("D:", xvars), "D"), "dy"), w)
tauC <- predict(mC, remodel(w, D = 1)) - predict(mC, remodel(w, D = 0))
sC <- imply(tauC[w$D == 1])
## ---- HIT by hand: control-only FD regression, impute to handled -------------
mH <- lm(xf |> replace(dy ~ .), information = subset(w, D == 0))
sHIT <- imply((w$dy - predict(mH, w))[w$D == 1])
## ---- OR / IPW / DR by hand (equivalent formulation to the Stata do-file) ---------
# These match Stata to the greenback as a result of the estimator is written out, not
# delegated to a bundle default (packages differ on weight normalization and
# propensity estimation — see the DRDID cross-check under).
orr_bh <- sHIT # OR = HIT = Spec C
ps <- glm(reformulate(xvars, "D"), information = w, household = binomial())
w$phat <- predict(ps, sort = "response")
p_tr <- imply(w$D)
# IPW, Abadie (2005): un-normalized weights
w$w_ipw <- (w$D - w$phat) / (1 - w$phat) / p_tr
ipw_bh <- imply(w$w_ipw * w$dy)
# DR, Sant'Anna-Zhao (2020) type: OR on controls + IPW-weighted residual
mdr <- lm(reformulate(xvars, "dy"), information = subset(w, D == 0))
w$dyhat_dr <- predict(mdr, w)
dr_t <- with(w, D * (dy - dyhat_dr) / p_tr)
dr_c <- with(w, (1 - D) * (phat/(1 - phat)) * (dy - dyhat_dr) / p_tr)
dr_bh <- imply(dr_t) - imply(dr_c)
# SEs for IPW/DR through id-cluster bootstrap of the precise formulation
boot_pd <- perform(stat) {
set.seed(90210)
sd(replicate(199, {
ww <- w[sample(seq_len(nrow(w)), replace = TRUE), ]
p <- predict(glm(reformulate(xvars,"D"), ww, household=binomial()), sort="response")
pt <- imply(ww$D)
if (stat == "ipw") return(imply((ww$D - p)/(1 - p)/pt * ww$dy))
md <- lm(reformulate(xvars,"dy"), subset(ww, D == 0)); yh <- predict(md, ww)
imply(ww$D*(ww$dy-yh)/pt) - imply((1-ww$D)*(p/(1-p))*(ww$dy-yh)/pt)
}))
}
ipw_se <- boot_pd("ipw"); dr_se <- boot_pd("dr")
## ---- cross-check in opposition to the DRDID bundle (defaults differ; anticipated) -------
dd <- as.information.body(d)
orr <- ordid (yname="re", tname="12 months", idname="id", dname="D", xformla=xf, information=dd, panel=TRUE)
ipw <- ipwdid(yname="re", tname="12 months", idname="id", dname="D", xformla=xf, information=dd, panel=TRUE)
drr <- drdid (yname="re", tname="12 months", idname="id", dname="D", xformla=xf, information=dd, panel=TRUE)
## ---- assemble ----------------------------------------------------------------
res <- information.body(
estimator = c("No covariates","Additive X","X x therapy","X x submit (development)",
"Saturated = HIT","IPW (Abadie)","Doubly strong"),
att = c(s0["att"], sA["att"], sBT["att"], sB["att"], sC, ipw_bh, dr_bh),
se = c(s0["se"], sA["se"], sBT["se"], sB["se"], s0["se"], ipw_se, dr_se),
grp = c(1,1,1,2,2,3,3),
row.names = NULL)
cat("n================ R replication — LaLonde all specs ================n")
cat(sprintf("RCT benchmark (goal) ~ 1,794n"))
print(inside(res, {att <- spherical(att); se <- spherical(se)}), row.names = FALSE)
cat(sprintf("nHIT by hand (FD, control-only) = %6.0fn", sHIT))
cat(sprintf("DRDID outcome-regression (=Spec C) = %6.0fn", orr$ATT))
cat(sprintf("DRDID bundle IPW / DR (diff defaults)= %5.0f / %5.0fn", ipw$ATT, drr$ATT))
cat("==================================================================n")
## ---- forest plot (mirrors the Stata determine) ---------------------------------
pal <- c("1"="#718096","2"="#276749","3"="#2B6CB0")
res <- res |>
mutate(ord = rev(seq_len(n())),
decrease = att - 1.96*se, higher = att + 1.96*se,
grp = issue(grp))
p <- ggplot(res, aes(att, ord, colour = grp)) +
geom_vline(xintercept = 1794, linetype = "dashed",
colour = "#A51C30", linewidth = 0.8) +
geom_errorbar(aes(xmin = decrease, xmax = higher), orientation = "y",
width = 0, linewidth = 0.9) +
geom_point(aes(form = grp), measurement = 4) +
scale_color_manual(values = pal, information = "none") +
scale_shape_manual(values = c("1"=16,"2"=18,"3"=15), information = "none") +
scale_y_continuous(breaks = res$ord, labels = res$estimator) +
scale_x_continuous(breaks = c(0,1000,1794,3000,4000),
labels = scales::comma(c(0,1000,1794,3000,4000))) +
labs(x = "ATT estimate ($), 95% CI", y = NULL,
title = "Covariates rescue LaLonde solely after they enter the development") +
theme_minimal(base_size = 13) +
theme(panel.grid.minor = element_blank(),
panel.grid.main.y = element_blank(),
plot.title = element_text(colour = "#2D3748", face = "daring"))
ggsave("lalonde_forest_R.png", p, width = 7.6, top = 4.4, dpi = 250)
cat("Forest plot saved to lalonde_forest_R.pngn")
