*! version 4.2.2 25apr2017 E. Leuven, B. Sianesi program define pstest version 11.0 #delimit ; syntax [varlist(default=none fv)] [if] [in] [, Treated(varname) BOTH RAW SUPport(varname) MWeight(varname) DENSity BOX OUTlier NOTable DISt LABel ONLYsig GRaph HIST SCatter RUBin ATU * ]; #delimit cr marksample touse, novarlist if "`varlist'"=="" & "`r(exog)'"!="" { local varlist = "`r(exog)'" } if "`varlist'"=="" & "`r(exog)'"=="" { di as error "Error: specify covariates" exit 198 } if ("`density'"!="" & "`box'"!="") { di as error "Error: choose between " as input "density " as error "and " as input "box" exit 198 } if ("`both'"!="" & "`raw'"!="") { di as error "Error: choose between " as input "both " as error "and " as input "raw" exit 198 } if ("`graph'"!="" & "`hist'"!="") { di as error "Error: choose between " as input "dot graph " as error "and " as input "histogram" exit 198 } if ("`graph'"!="" & "`scatter'"!="") { di as error "Error: choose between " as input "dot graph " as error "and " as input "scatter" exit 198 } if ("`scatter'"!="" & "`hist'"!="") { di as error "Error: choose between " as input "scatter " as error "and " as input "histogram" exit 198 } capture confirm var _treated if (_rc & "`treated'"=="") | ("`raw'"!="" & "`treated'"=="") { di as error "Error: provide treatment indicator variable" exit 198 } else if (!_rc & "`treated'"=="") { tempvar treated qui g double `treated' = _treated } tempvar weight if ("`mweight'"=="") { if ("`raw'"!="") qui g double `weight' = 1 else { capture confirm var _weight if (!_rc) qui g double `weight' = _weight else di as error "Error: provide weight" } } else g double `weight' = `mweight' if ("`support'"=="") { tempvar support capture confirm var _support if (_rc | "`raw'"!="") qui g byte `support' = 1 else qui g byte `support'= _support } qui replace `weight' = `support' if cond("`atu'" == "", `treated'==1, `treated'==0) if ("`density'"=="" & "`box'"=="" & "`both'"!="") { breduc `varlist' , touse(`touse') mw(`weight') tr(`treated') sup(`support') `notable' `dist' `label' `graph' `hist' `scatter' options("`options'") `rubin' } if ("`density'"=="" & "`box'"=="" & "`both'"=="") { breduc1 `varlist' , touse(`touse') mw(`weight') tr(`treated') sup(`support') `notable' `dist' `label' `graph' `hist' `scatter' options("`options'") `rubin' `onlysig' `raw' } if ("`density'"!="" | "`box'"!="") { plotvar `varlist' , touse(`touse') `raw' `both' mw(`weight') tr(`treated') sup(`support') `density' `box' `outlier' options("`options'") } end program define breduc, rclass syntax varlist(min=1 fv), MWeight(varname) TReated(varname) SUPport(varname) [touse(varname) NOTable DISt LABel GRaph HIST SCatter options(string) RUBin] tempvar sumbias sumbias0 _bias0 _biasm xvar meanbiasbef medbiasbef meanbiasaft medbiasaft _vratio_bef _vratio_aft tempname Flowu Fhighu Flowm Fhighm qui count if `treated'==1 & `touse' scalar `Flowu' = invF(r(N)-1, r(N)-1, 0.025) scalar `Fhighu' = invF(r(N)-1, r(N)-1, 0.975) qui count if `treated'==1 & `support'==1 & `touse' scalar `Flowm' = invF(r(N)-1, r(N)-1, 0.025) scalar `Fhighm' = invF(r(N)-1, r(N)-1, 0.975) qui g `_bias0' = . qui g `_biasm' = . qui g str12 `xvar' = "" qui g `sumbias' = . qui g `sumbias0' = . qui g `_vratio_bef' = . qui g `_vratio_aft' = . if "`notable'"!="" { local quietly "quietly" } fvexpand `varlist' local hasfactorvars = ("`=r(fvops)'" == "true") local vnames `r(varlist)' local vlength 22 foreach v of local vnames { local vlength = max(`vlength', length("`v'")) } /* construct header */ local c = `vlength' + 4 local s = `vlength' - 22 if ("`rubin'"!="" | "`scatter'"!="") { local add "_e" } `quietly' di `quietly' di as text "{hline `c'}{c TT}{hline 34}{c TT}{hline 15}{c TT}{hline 10}" `quietly' di as text " " _s(`s') " Unmatched {c |} Mean %reduct {c |} t-test {c |} V`add'(T)/" `quietly' di as text "Variable " _s(`s') " Matched {c |} Treated Control %bias |bias| {c |} t p>|t| {c |} V`add'(C)" `quietly' di as text "{hline `c'}{c +}{hline 34}{c +}{hline 15}{c +}{hline 10}" /* get linear index and some overall stats for later*/ tempvar index0 indexm tempname r2bef r2aft chibef chiaft pchibef pchiaft qui probit `treated' `varlist' if `touse' qui predict double `index0' if e(sample), xb scalar `r2bef' = e(r2_p) scalar `chibef' = e(chi2) scalar `pchibef' = chi2tail(e(df_m), e(chi2)) return scalar r2bef = e(r2_p) return scalar chiprobbef = chi2tail(e(df_m), e(chi2)) qui probit `treated' `varlist' [iw=`mweight'] if `support'==1 & `touse' qui predict double `indexm' if e(sample), xb scalar `r2aft' = e(r2_p) scalar `chiaft' = e(chi2) scalar `pchiaft' = chi2tail(e(df_m), e(chi2)) return scalar r2aft = e(r2_p) return scalar chiprobaft = chi2tail(e(df_m), e(chi2)) /* calculate stats for varlist */ tempname m1u m0u v1u v0u m1m m0m bias biasm absreduc tbef taft pbef paft tempname v1m v0m v_ratiobef v_ratioaft v_e_1 tempvar resid0 resid1 local cnt_concbef = 0 /* counting vars with ratio of concern - rubin */ local cnt_concaft = 0 local cnt_badbef = 0 /* counting vars with bad ratio - rubin */ local cnt_badaft = 0 local cont_cnt = 0 /* counting continuous vars */ local cont_varbef = 0 /* counting continuous vars w/ excessive var ratio*/ local cont_varaft = 0 local i 0 fvrevar `varlist' foreach v in `r(varlist)' { local ++i local xlab : word `i' of `vnames' if (regexm("`xlab'", ".*b[\\.].*") == 1) continue if (`hasfactorvars'==0 & "`label'" != "") { local xlab : var label `v' if ("`xlab'" == "") local xlab `v' } qui sum `v' if `treated'==1 & `touse' scalar `m1u' = r(mean) scalar `v1u' = r(Var) qui sum `v' if `treated'==0 & `touse' scalar `m0u' = r(mean) scalar `v0u' = r(Var) qui sum `v' [iw=`mweight'] if `treated'==1 & `support'==1 & `touse' scalar `m1m' = r(mean) scalar `v1m' = r(Var) qui sum `v' [iw=`mweight'] if `treated'==0 & `support'==1 & `touse' scalar `m0m' = r(mean) scalar `v0m' = r(Var) scalar `v_ratiobef' = . scalar `v_ratioaft' = . local starbef "" local staraft "" if ("`rubin'"=="" & "`scatter'"=="") { capture assert `v'==0 | `v'==1 | `v'==., fast if (_rc) { local cont_cnt = `cont_cnt' +1 /* get Var ratio*/ scalar `v_ratiobef' = `v1u'/`v0u' if `v_ratiobef'>`Fhighu' | `v_ratiobef'<`Flowu' { local cont_varbef = `cont_varbef' +1 local starbef "*" } scalar `v_ratioaft' = `v1m'/`v0m' if `v_ratioaft'>`Fhighm' | `v_ratioaft'<`Flowm' { local cont_varaft = `cont_varaft' +1 local staraft "*" } } } if ("`rubin'"!="" | "`scatter'"!="") { cap drop `resid1' cap drop `resid0' qui regress `v' `index0' if `treated'==1 & `touse' qui predict double `resid1' if e(sample), resid qui regress `v' `index0' if `treated'==0 & `touse' qui predict double `resid0' if e(sample), resid qui sum `resid1' scalar `v_e_1' = r(Var) qui sum `resid0' scalar `v_ratiobef' = `v_e_1'/r(Var) qui replace `_vratio_bef' = `v_ratiobef' in `i' if (`v_ratiobef'>1.25 & `v_ratiobef'<=2) | (`v_ratiobef'<0.8 & `v_ratiobef'>=0.5) { local cnt_concbef = `cnt_concbef' +1 local starbef "*" } if (`v_ratiobef'>2 & `v_ratiobef'<.) | `v_ratiobef'<0.5 { local cnt_badbef = `cnt_badbef' +1 local starbef "**" } drop `resid1' drop `resid0' qui regress `v' `indexm' [iw=`mweight'] if `treated'==1 & `support'==1 & `touse' qui predict double `resid1' if e(sample), resid qui regress `v' `indexm' [iw=`mweight'] if `treated'==0 & `support'==1 & `touse' qui predict double `resid0' if e(sample), resid qui sum `resid1' [iw=`mweight'] scalar `v_e_1' = r(Var) qui sum `resid0' [iw=`mweight'] scalar `v_ratioaft' = `v_e_1'/r(Var) qui replace `_vratio_aft' = `v_ratioaft' in `i' if (`v_ratioaft'>1.25 & `v_ratioaft'<=2) | (`v_ratioaft'<0.8 & `v_ratioaft'>=0.5) { local cnt_concaft = `cnt_concaft' +1 local staraft "*" } if (`v_ratioaft'>2 & `v_ratioaft'<.) | `v_ratioaft'<0.5 { local cnt_badaft = `cnt_badaft' +1 local staraft "**" } } qui replace `xvar' = "`v'" in `i' /* standardised % bias before matching */ scalar `bias' = 100*(`m1u' - `m0u')/sqrt((`v1u' + `v0u')/2) qui replace `_bias0' = `bias' in `i' qui replace `sumbias0' = abs(`bias') in `i' /* standardised % bias after matching */ scalar `biasm' = 100*(`m1m' - `m0m')/sqrt((`v1u' + `v0u')/2) qui replace `_biasm' = `biasm' in `i' qui replace `sumbias' = abs(`biasm') in `i' /* % reduction in absolute bias */ scalar `absreduc' = -100*(abs(`biasm') - abs(`bias'))/abs(`bias') /* t-tests before matching */ qui regress `v' `treated' if `touse' scalar `tbef' = _b[`treated']/_se[`treated'] scalar `pbef' = 2*ttail(e(df_r),abs(`tbef')) /* t-tests after matching */ qui regress `v' `treated' [iw=`mweight'] if `support'==1 & `touse' scalar `taft' = _b[`treated']/_se[`treated'] scalar `paft' = 2*ttail(e(df_r),abs(`taft')) `quietly' di as text %-`vlength's abbrev("`xlab'",`vlength') _col(`=`c'-2') "U {c |}" as result %7.0g `m1u' " " %7.0g `m0u' " " %7.1f `bias' _s(8) as text " {c |}" as res %7.2f `tbef' _s(2) as res %05.3f `pbef' " {c |}" as res %6.2f `v_ratiobef' "`starbef'" `quietly' di as text _col(`=`c'-2') "M {c |}" as result %7.0g `m1m' " " %7.0g `m0m' " " %7.1f `biasm' %8.1f `absreduc' as text " {c |}" as res %7.2f `taft' _s(2) as res %05.3f `paft' " {c |}" as res %6.2f `v_ratioaft' "`staraft'" `quietly' di as text _col(`=`c'-2') " {c |}" as text _s(31) " {c |}" as text _s(12) " {c |}" } `quietly' di as text "{hline `c'}{c BT}{hline 34}{c BT}{hline 15}{c BT}{hline 10}" if ("`rubin'"=="" & "`scatter'"=="") { `quietly' di as text "* if variance ratio outside [" %4.2f `Flowu' "; " %4.2f `Fhighu' "] for U and [" %4.2f `Flowm' "; " %4.2f `Fhighm' "] for M" } if ("`rubin'"!="" | "`scatter'"!="") { `quietly' di as text "* if 'of concern', i.e. variance ratio in [0.5, 0.8) or (1.25, 2]" `quietly' di "** if 'bad', i.e. variance ratio <0.5 or >2 " } di "" local quietly "" if "`dist'"=="" { local quietly "quietly" } `quietly' di as text "{hline 61}" `quietly' di as text _col(10) "Summary of the distribution of the abs(bias)" `quietly' di as text "{hline 61}" label var `sumbias0' "BEFORE MATCHING" `quietly' sum `sumbias0', detail scalar `meanbiasbef' = r(mean) scalar `medbiasbef' = r(p50) return scalar meanbiasbef = r(mean) return scalar medbiasbef = r(p50) `quietly' di as text "{hline 61}" label var `sumbias' "AFTER MATCHING" `quietly' sum `sumbias', detail scalar `meanbiasaft' = r(mean) scalar `medbiasaft' = r(p50) return scalar meanbiasaft = r(mean) return scalar medbiasaft = r(p50) `quietly' di as text "{hline 61}" `quietly' di // Rubin's stats: // Rubin's B = absolute standardized differences of means of the linear index of the propensity score // Rubin's R = ratio of treated to non-treated variance of the propensity score index tempname mi1 vi1 mi0 vi0 iratiobef ibiasbef iratioaft ibiasaft qui sum `index0' if `treated'==1 & `touse' scalar `mi1' = r(mean) scalar `vi1' = r(Var) qui sum `index0' if `treated'==0 & `touse' scalar `mi0' = r(mean) scalar `vi0' = r(Var) scalar `ibiasbef' = 100*(`mi1' - `mi0')/sqrt((`vi1' + `vi0')/2) scalar `iratiobef' = `vi1'/`vi0' return scalar Bbef = `ibiasbef' return scalar Rbef = `iratiobef' if (`ibiasbef'>=25) local starBbef "*" if !inrange(`iratiobef', 0.5, 2) local starRbef "*" qui sum `indexm' [iw=`mweight'] if `treated'==1 & `support'==1 & `touse' scalar `mi1' = r(mean) scalar `vi1' = r(Var) qui sum `indexm' [iw=`mweight'] if `treated'==0 & `support'==1 & `touse' scalar `mi0' = r(mean) scalar `vi0' = r(Var) scalar `ibiasaft' = 100*(`mi1' - `mi0')/sqrt((`vi1' + `vi0')/2) scalar `iratioaft' = `vi1'/`vi0' return scalar Baft = `ibiasaft' return scalar Raft = `iratioaft' if (`ibiasaft'>=25) local starBaft "*" if !inrange(`iratioaft', 0.5, 2) local starRaft "*" if ("`rubin'"=="" & "`scatter'"=="") { di as text "{hline 11}{c TT}{hline 71}" di as text " Sample {c |} Ps R2 LR chi2 p>chi2 MeanBias MedBias B R %Var" di as text "{hline 11}{c +}{hline 71}" di as text " Unmatched {c | }" as res %6.3f `r2bef' _s(1) as res %9.2f `chibef' _s(1) as res %8.3f `pchibef' _s(3) as res %6.1f `meanbiasbef' _s(4) as res %6.1f `medbiasbef' _s(4) as res %6.1f `ibiasbef' "`starBbef'" _col(70) as res %5.2f `iratiobef' "`starRbef'" _col(79) as res %3.0f 100*`cont_varbef'/`cont_cnt' di as text " Matched {c | }" as res %6.3f `r2aft' _s(1) as res %9.2f `chiaft' _s(1) as res %8.3f `pchiaft' _s(3) as res %6.1f `meanbiasaft' _s(4) as res %6.1f `medbiasaft' _s(4) as res %6.1f `ibiasaft' "`starBaft'" _col(70) as res %5.2f `iratioaft' "`starRaft'" _col(79) as res %3.0f 100*`cont_varaft'/`cont_cnt' di as text "{hline 11}{c BT}{hline 71}" di as text "* if B>25%, R outside [0.5; 2]" } if ("`rubin'"!="" | "`scatter'"!="") { di as text "{hline 11}{c TT}{hline 81}" di as text " Sample {c |} Ps R2 LR chi2 p>chi2 MeanBias MedBias B R %concern %bad" di as text "{hline 11}{c +}{hline 81}" di as text " Unmatched {c | }" as res %6.3f `r2bef' _s(1) as res %9.2f `chibef' _s(1) as res %8.3f `pchibef' _s(3) as res %6.1f `meanbiasbef' _s(4) as res %6.1f `medbiasbef' _s(4) as res %6.1f `ibiasbef' "`starBbef'" _col(70) as res %5.2f `iratiobef' "`starRbef'" _col(80) as res %3.0f 100*`cnt_concbef'/`i' _s(6) as res %3.0f 100*`cnt_badbef'/`i' di as text " Matched {c | }" as res %6.3f `r2aft' _s(1) as res %9.2f `chiaft' _s(1) as res %8.3f `pchiaft' _s(3) as res %6.1f `meanbiasaft' _s(4) as res %6.1f `medbiasaft' _s(4) as res %6.1f `ibiasaft' "`starBaft'" _col(70) as res %5.2f `iratioaft' "`starRaft'" _col(80) as res %3.0f 100*`cnt_concaft'/`i' _s(6) as res %3.0f 100*`cnt_badaft'/`i' di as text "{hline 11}{c BT}{hline 81}" di as text "* if B>25%, R outside [0.5; 2]" } if ("`graph'"!="") { qui count if `xvar'!="" if r(N) > 30 { local nolabelx "label(nolabel)" } graph dot `_bias0' `_biasm', over(`xvar', sort(1) descending `nolabelx') legend(pos(5) ring(0) col(1) lab(1 "Unmatched") lab(2 "Matched")) yline(0, lcolor(gs10)) marker(1, mcolor(black) msymbol(O)) marker(2, mcolor(black) msymbol(X)) ytitle("Standardized % bias across covariates") `options' } if "`hist'"!="" { tempname grbef graft qui sum `_bias0' local bnd = round(max(-r(min), r(max)), 4) local stp = `bnd'/4 qui histogram `_bias0', xlab(-`bnd'(`stp')`bnd') xtitle("Standardized % bias across covariates") title("Unmatched") `options' saving(`grbef'.gph , replace) nodraw qui histogram `_biasm', xlab(-`bnd'(`stp')`bnd') xtitle("Standardized % bias across covariates") title("Matched") `options' saving(`graft'.gph , replace) nodraw qui graph combine `grbef'.gph `graft'.gph, xsize(6) ysize(7) col(1) scheme(s1mono) ycommon qui erase `grbef'.gph qui erase `graft'.gph } if ("`scatter'"!="") { tempname grbef graft qui sum `_bias0' local bnd = round(max(-r(min), r(max)), 4) local stp = `bnd'/4 qui scatter `_vratio_bef' `_bias0', xline(0, lw(medthick) lc(gs5)) yline(1, lw(medthick) lc(gs5)) yline(0.8 1.25, lp(dash) lw(medium) lc(gs5)) yline(0.5 2, lp(dot) lw(medium) lc(gs5)) xlab(-`bnd'(`stp')`bnd') ylab(0(0.5)2) ytitle("Variance ratio of residuals") xtitle("Standardized % bias") title("Unmatched") `options' saving(`grbef'.gph , replace) nodraw qui scatter `_vratio_aft' `_biasm', xline(0, lw(medthick) lc(gs5)) yline(1, lw(medthick) lc(gs5)) yline(0.8 1.25, lp(dash) lw(medium) lc(gs5)) yline(0.5 2, lp(dot) lw(medium) lc(gs5)) xlab(-`bnd'(`stp')`bnd') ylab(0(0.5)2) ytitle("Variance ratio of residuals") xtitle("Standardized % bias") title("Matched") `options' saving(`graft'.gph , replace) nodraw qui graph combine `grbef'.gph `graft'.gph, xsize(6) ysize(7) col(1) scheme(s1mono) ycommon qui erase `grbef'.gph qui erase `graft'.gph } return local exog = "`varlist'" end /* ************************************************************************************************************************************* */ program define breduc1, rclass syntax varlist(min=1 fv) , [ RAW MWeight(varname) TReated(varname) SUPport(varname) touse(varname) NOTable DISt LABel ONLYsig GRaph HIST SCatter options(string) RUBin] tempvar sumbias _bias xvar _vratio tempname Flow Fhigh qui count if `treated'==1 & `support'==1 & `touse' scalar `Flow' = invF(r(N)-1, r(N)-1, 0.025) scalar `Fhigh' = invF(r(N)-1, r(N)-1, 0.975) qui g `_bias' = . qui g str12 `xvar' = "" qui g `sumbias' = . qui g `_vratio' =. if "`notable'"!="" { local quietly "quietly" } fvexpand `varlist' local hasfactorvars = ("`=r(fvops)'" == "true") local vnames `r(varlist)' local vlength 22 foreach v of local vnames { local vlength = max(`vlength', length("`v'")) } /* construct header */ local c = `vlength' + 2 if ("`rubin'"!="" | "`scatter'"!="") { local add "_e" } `quietly' di `quietly' di as text "{hline `c'}{c TT}{hline 26}{c TT}{hline 15}{c TT}{hline 10}" `quietly' di as text " " _col(`c') " {c |} Mean {c |} t-test {c |} V`add'(T)/" `quietly' di as text "Variable" _col(`c') " {c |} Treated Control %bias {c |} t p>|t| {c |} V`add'(C)" `quietly' di as text "{hline `c'}{c +}{hline 26}{c +}{hline 15}{c +}{hline 10}" /* get linear index and some overall stats for later*/ tempvar index tempname r2 chi pchi qui probit `treated' `varlist' [iw=`mweight'] if `support'==1 & `touse' qui predict double `index' if e(sample), xb scalar `r2' = e(r2_p) scalar `chi' = e(chi2) scalar `pchi' = chi2tail(e(df_m), e(chi2)) return scalar r2 = e(r2_p) return scalar chiprob = chi2tail(e(df_m), e(chi2)) /* calculate stats for varlist */ tempname m1m m0m v1u v0u biasm t p meanbias medbias var1 var0 v_ratio v_e_1 tempvar resid0 resid1 local cnt_conc = 0 /* counting vars with ratio of concern - rubin */ local cnt_bad = 0 /* counting vars with bad ratio - rubin */ local cont_cnt = 0 /* counting continuous vars */ local cont_var = 0 /* counting continuous vars w/ excessive var ratio*/ local i 0 fvrevar `varlist' foreach v in `r(varlist)' { local ++i local xlab : word `i' of `vnames' if (regexm("`xlab'", ".*b[\\.].*") == 1) continue if (`hasfactorvars'==0 & "`label'" != "") { local xlab : var label `v' if ("`xlab'" == "") local xlab `v' } qui sum `v' if `treated'==1 & `touse' scalar `v1u' = r(Var) qui sum `v' if `treated'==0 & `touse' scalar `v0u' = r(Var) qui sum `v' [iw=`mweight'] if `treated'==1 & `support'==1 & `touse' scalar `m1m' = r(mean) scalar `var1' = r(Var) qui sum `v' [iw=`mweight'] if `treated'==0 & `support'==1 & `touse' scalar `m0m' = r(mean) scalar `var0' = r(Var) scalar `v_ratio' = . local star "" if ("`rubin'"=="" & "`scatter'"=="") { capture assert `v'==0 | `v'==1 | `v'==., fast if (_rc) { local cont_cnt = `cont_cnt' +1 /* get Var ratio */ scalar `v_ratio' = `var1'/`var0' if `v_ratio'>`Fhigh' | `v_ratio'<`Flow' { local cont_var = `cont_var' +1 local star "*" } } } if ("`rubin'"!="" | "`scatter'"!="") { cap drop `resid1' cap drop `resid0' qui regress `v' `index' [iw=`mweight'] if `treated'==1 & `support'==1 & `touse' qui predict double `resid1' if e(sample), resid qui regress `v' `index' [iw=`mweight'] if `treated'==0 & `support'==1 & `touse' qui predict double `resid0' if e(sample), resid qui sum `resid1' [iw=`mweight'] scalar `v_e_1' = r(Var) qui sum `resid0' [iw=`mweight'] scalar `v_ratio' = `v_e_1'/r(Var) qui replace `_vratio' = `v_ratio' in `i' if (`v_ratio'>1.25 & `v_ratio'<=2) | (`v_ratio'<0.8 & `v_ratio'>=0.5) { local cnt_conc = `cnt_conc' +1 local star "*" } if (`v_ratio'>2 & `v_ratio'<.) | `v_ratio'<0.5 { local cnt_bad = `cnt_bad' +1 local star "**" } } qui replace `xvar' = "`v'" in `i' /* standardised % bias after matching */ scalar `biasm' = 100*(`m1m' - `m0m')/sqrt((`v1u' + `v0u')/2) qui replace `_bias' = `biasm' in `i' qui replace `sumbias' = abs(`biasm') in `i' /* t-tests after matching */ qui regress `v' `treated' [iw=`mweight'] if `support'==1 & `touse' scalar `t' = _b[`treated']/_se[`treated'] scalar `p' = 2*ttail(e(df_r), abs(`t')) if ("`onlysig'" != "" & `p' >= 0.10) continue local c = `vlength' + 2 `quietly' di as text %-`vlength's substr("`xlab'",1, `vlength') _col(`c') " {c |}" as result %7.0g `m1m' " " %7.0g `m0m' " " %7.1f `biasm' as text " {c |}" as res %7.2f `t' _s(2) as res %05.3f `p' as text " {c |}" as res %6.2f `v_ratio' "`star'" } `quietly' di as text "{hline `c'}{c BT}{hline 26}{c BT}{hline 15}{c BT}{hline 10}" if ("`rubin'"=="" & "`scatter'"=="") { `quietly' di as text "* if variance ratio outside [" %4.2f `Flow' "; " %4.2f `Fhigh' "]" } if ("`rubin'"!="" | "`scatter'"!="") { `quietly' di as text "* if 'of concern', i.e. variance ratio in [0.5, 0.8) or (1.25, 2]" `quietly' di "** if 'bad', i.e. variance ratio <0.5 or >2 " } di "" local quietly "" if "`dist'"=="" { local quietly "quietly" } `quietly' di as text "{hline 61}" label var `sumbias' "Summary of the distribution of |bias|" `quietly' sum `sumbias', detail scalar `meanbias' = r(mean) scalar `medbias' = r(p50) return scalar meanbias = r(mean) return scalar medbias = r(p50) `quietly' di as text "{hline 61}" `quietly' di // Rubin's stats // Rubin's B = absolute standardized differences of means of the linear index of the propensity score tempname mi1 vi1 mi0 vi0 iratio ibias qui sum `index' [iw=`mweight'] if `treated'==1 & `support'==1 & `touse' scalar `mi1' = r(mean) scalar `vi1' = r(Var) qui sum `index' [iw=`mweight'] if `treated'==0 & `support'==1 & `touse' scalar `mi0' = r(mean) scalar `vi0' = r(Var) scalar `ibias' = 100*(`mi1' - `mi0')/sqrt((`vi1' + `vi0')/2) // Rubin's R = ratio of treated to non-treated variance of the propensity score index scalar `iratio' = `vi1'/`vi0' return scalar B = `ibias' return scalar R = `iratio' if (`ibias'>=25) local starB "*" if !inrange(`iratio', 0.5, 2) local starR "*" if ("`rubin'"=="" & "`scatter'"=="") { di as text "{hline 70}" di as text "Ps R2 LR chi2 p>chi2 MeanBias MedBias B R %Var " di as text "{hline 70}" di as text as res %5.3f `r2' _s(1) as res %9.2f `chi' _s(1) as res %8.3f `pchi' _s(3) as res %6.1f `meanbias' _s(4) as res %6.1f `medbias' _s(4) as res %6.1f `ibias' "`starB'" _s(3) as res %5.2f `iratio' "`starR'" _col(67) as res %3.0f 100*`cont_var'/`cont_cnt' di as text "{hline 70}" di as text "* if B>25%, R outside [0.5; 2]" } if ("`rubin'"!="" | "`scatter'"!="") { di as text "{hline 81}" di as text "Ps R2 LR chi2 p>chi2 MeanBias MedBias B R %concern %bad" di as text "{hline 81}" di as text as res %5.3f `r2' _s(1) as res %9.2f `chi' _s(1) as res %8.3f `pchi' _s(3) as res %6.1f `meanbias' _s(4) as res %6.1f `medbias' _s(4) as res %6.1f `ibias' "`starB'" _s(3) as res %5.2f `iratio' "`starR'" _col(68) as res %3.0f 100*`cnt_conc'/`i' _s(6) as res %3.0f 100*`cnt_bad'/`i' di as text "{hline 81}" di as text "* if B>25%, R outside [0.5; 2]" } if "`graph'"!="" { qui count if `xvar'!="" if r(N) > 30 { local nolabelx "label(nolabel)" } graph dot `_bias', over(`xvar', sort(1) descending `nolabelx') yline(0, lcolor(gs10)) marker(1, mcolor(black) msymbol(O)) ytitle("Standardized % bias across covariates") `options' } if "`hist'"!="" { qui histogram `_bias', xtitle("Standardized % bias across covariates") `options' } if ("`scatter'"!="") { qui scatter `_vratio' `_bias', xline(0, lw(medthick) lc(gs5)) yline(1, lw(medthick) lc(gs5)) yline(0.8 1.25, lp(dash) lw(medium) lc(gs5)) yline(0.5 2, lp(dot) lw(medium) lc(gs5)) ylab(0(0.5)2) ytitle("Variance ratio of residuals") xtitle("Standardized % bias") `options' } return local exog = "`varlist'" end program define plotvar syntax varname, [RAW BOTH MWeight(varname) TReated(varname) SUPport(varname) touse(varname) DENSity BOX OUTlier options(string) ] capture assert `varlist'==0 | `varlist'==1 | `varlist'==., fast if (!_rc) { di as error "Error: you can't specify a dummy variable to be plotted with options " as input "density " as error "or " as input "box" exit 198 } local Ytitle : var label `varlist' if ("`Ytitle'" == "") local Ytitle `varlist' if ("`outlier'" == "") local nooutsides "nooutsides" if ("`density'"!="") { tempname grbef graft if ("`raw'"=="" & "`both'"=="") { qui twoway (kdensity `varlist' if `touse' & `treated'==1 [aw=`mweight'], clwid(thick)) (kdensity `varlist' if `touse' & `treated'==0 [aw=`mweight'], clwid(thin) clcolor(black)), xlab(#6) xti("") yti("") title("`Ytitle'") subtitle("Matched samples") legend(order(1 "Treated" 2 "Untreated")) graphregion(color(gs16)) `options' } if ("`raw'"!="") { qui twoway (kdensity `varlist' if `touse' & `treated'==1, clwid(thick)) (kdensity `varlist' if `touse' & `treated'==0, clwid(thin) clcolor(black)), xlab(#6) xti("") yti("") title("`Ytitle'") legend(order(1 "`treated'==1" 2 "`treated'==0")) graphregion(color(gs16)) `options' } if ("`both'"!="") { qui twoway (kdensity `varlist' if `touse' & `treated'==1, clwid(thick)) (kdensity `varlist' if `touse' & `treated'==0, clwid(thin) clcolor(black)), title("Unmatched") xlab(#6) ytitle("") xtitle("") legend(off) saving(`grbef'.gph , replace) graphregion(color(gs16)) nodraw qui twoway (kdensity `varlist' if `touse' & `treated'==1 [aw=`mweight'], clwid(thick)) (kdensity `varlist' if `touse' & `treated'==0 [aw=`mweight'], clwid(thin) clcolor(black)), title("Matched") xlab(#6) ytitle("") xtitle("") legend(order(1 "Treated" 2 "Untreated")) saving(`graft'.gph , replace) graphregion(color(gs16)) nodraw qui graph combine `grbef'.gph `graft'.gph, xsize(6) ysize(7) title("`Ytitle'") scheme(s1mono) col(1) xcommon `options' qui erase `grbef'.gph qui erase `graft'.gph } } if ("`box'"!="") { tempname grbef graft if ("`raw'"=="" & "`both'"=="") { qui graph box `varlist' if `touse' [aw=`mweight'], over(`treated', sort(`treated') descending relabel(1 "Untreated" 2 "Treated")) `nooutsides' note("") yti("") title("`Ytitle'") subtitle("Matched samples") `options' } if ("`raw'"!="") { qui graph box `varlist' if `touse', over(`treated', sort(`treated') descending relabel(1 "Untreated" 2 "Treated")) `nooutsides' note("") yti("") title("`Ytitle'") `options' } if ("`both'"!="") { qui graph box `varlist' if `touse', over(`treated', sort(`treated') descending relabel(1 "Untreated" 2 "Treated")) `nooutsides' note("") yti("") title("Unmatched") saving(`grbef'.gph , replace) nodraw qui graph box `varlist' if `touse' [aw=`mweight'], over(`treated', sort(`treated') descending relabel(1 "Untreated" 2 "Treated")) `nooutsides' note("") yti("") title("Matched") saving(`graft'.gph , replace) nodraw qui graph combine `grbef'.gph `graft'.gph, xsize(6) ysize(7) title("`Ytitle'") scheme(s1mono) col(1) xcommon `options' qui erase `grbef'.gph qui erase `graft'.gph } } end