clear
set seed 2
// Create data set with experimental conditions
input str4 educ str6 sex str9 treatmentstr str4 health1 number
Low Male Exposed Poor 333
Low Male Exposed Good 917
Low Male Unexposed Poor 1000
Low Male Unexposed Good 2750
Low Female Exposed Poor 500
Low Female Exposed Good 3250
Low Female Unexposed Poor 167
Low Female Unexposed Good 1083
High Male Exposed Poor 83
High Male Exposed Good 542
High Male Unexposed Poor 584
High Male Unexposed Good 3791
High Female Exposed Poor 125
High Female Exposed Good 1750
High Female Unexposed Poor 208
High Female Unexposed Good 2917
end
expand number // Create full number of cases
// Transform strings to numerical variables
generate loeduc = (educ == "Low")
generate female = (sex == "Female")
generate treatment = (treatmentstr == "Exposed")
generate good1 = (health1 == "Good")
label define loeduc 0 "High" 1 "Low"
label val loeduc loeduc
label define female 0 "Male" 1 "Female"
label val female female
label define treatment 0 "Unexposed" 1 "Exposed"
label val treatment treatment
label define health 0 "Poor" 1 "Good"
label val good1 health
// Simulate outcome variable
generate good2 = good1
replace good2 = 1 if good1 == 0 & loeduc == 1 & (runiform() <= .05)
replace good2 = 1 if good1 == 0 & loeduc == 0 & (runiform() <= .20)
replace good2 = 1 if good1 == 0 & treatment == 1 & (runiform() <= .30)
label val good2 health
// Transform some more and clean up
generate poor2 = (good2 == 0)
generate poor1 = (good1 == 0)
drop educ sex treatmentstr health1 number good1 good2
// Table 1-ish
table poor2, by(loeduc female treatment) contents(freq)
// 1) Regression adjustment
logit poor2 treatment female if loeduc == 0, or
logit poor2 treatment female if loeduc == 1, or
logit poor2 treatment##loeduc female##loeduc, or
// 2) Propensity score matching
teffects nnmatch (poor2 female) (treatment) if loeduc == 0
teffects nnmatch (poor2 female) (treatment) if loeduc == 1
// 3) Difference in difference
// Transform to long format
gen id = _n
reshape long poor, i(id) j(year)
logit poor treatment##c.year if loeduc == 0, or
logit poor treatment##c.year if loeduc == 1, or
logit poor loeduc##treatment##c.year, or
// 4) Fixed effects model
replace treatment = 0 if year == 1
xtset id year
xtreg poor treatment year if loeduc == 0, fe
xtreg poor treatment year if loeduc == 1, fe
xtreg poor treatment##loeduc year##loeduc, fe
Showing posts with label xtset. Show all posts
Showing posts with label xtset. Show all posts
Aug 2, 2017
Analyzing natural policy experiments
Hu et al. 2017 simulate data of a natural policy experiment and show how to analyze the data with regression adjustment, propensity score matching, difference-in-differences, and fixed effects regression. (The paper also includes IV, regression discontinuity, and interrupted time series, but does not describe the data created for these analyses.)
Feb 26, 2017
Chapter 2 of Singer and Willett's (2003) book on longitudinal data analysis
There is a different take on the chapter here, but I like mine better.
// Open data set use http://www.ats.ucla.edu/stat/stata/examples/alda/data/tolerance, clear // Figure 2.1 list, sep(0) noobs reshape long tol, i(id) j(age) list if inlist(id, 9, 45, 1653), sepby(id) noobs // Table 2.1 reshape wide quietly estpost cor tol*, matrix listwise esttab, unstack not noobs compress nostar b(%6.2f) // Figure 2.2
reshape long twoway (scatter tol age), by(id) ytitle(Tolerance) xtitle(Age) name(Figure22, replace) // Figure 2.3
twoway (scatter tol age) ///
(lowess tol age) ///
, by(id, note("Graphs by id, lowess curves") legend(off)) ///
ytitle(Tolerance) xtitle(Age) name(Figure23, replace)
// Table 2.2
preserve
tempfile table22
generate time = age - 11
statsby initial = _b[_cons] ///
initial_se = _se[_cons] ///
change = _b[time] ///
change_se = _se[time] ///
residual = (e(rmse)^2) ///
explained = e(r2) ///
, by(id) saving(`table22', replace): regress tol time
drop time
reshape wide
merge 1:1 id using `table22', nogenerate
format initial initial_se change change_se residual explained exposure %6.2f
keep id initial initial_se change change_se residual explained male exposure
order id initial initial_se change change_se residual explained male exposure
list, sep(0) noobs abbrev(15)
tempfile figure28
save `figure28', replace
// Figure 2.4
stem initial, round(.01)
stem change, round(.01)
stem residual, round(.01)
stem explained, round(.01)
restore
// Figure 2.5
twoway (scatter tol age) ///
(lfit tol age) ///
, by(id, note("Graphs by id, OLS curves") legend(off)) ///
ytitle(Tolerance) xtitle(Age) name(Figure25, replace)
// Figure 2.6
xtset id age
xtline tol, overlay t(age) i(id) legend(off) ///
ytitle(Tolerance) ylabel(1 (1) 4) ///
addplot(lowess tol age, lwidth(thick) lpattern(solid)) ///
xtitle(Age) title("Observed data and lowess smoother") ///
xsize(3) name(Figure26A, replace) nodraw
quietly regress tol i.id##c.age
predict tolhat
xtline tolhat, overlay t(age) i(id) legend(off) ///
ytitle(Tolerance) ylabel(1 (1) 4) ///
addplot(lfit tol age, lwidth(thick) lpattern(solid)) ///
xtitle(Age) title("OLS trajectories") ///
xsize(3) name(Figure26B, replace) nodraw
graph combine Figure26A Figure26B, col(2) name(Figure26, replace)
// Table 2.2
preserve
use `table22', clear
tabstat initial change, stat(mean sd) format(%6.2f)
cor initial change
restore
// Figure 2.7
quietly sum exposure, detail
generate highexposure = (exposure >= r(p50))
generate lowexposure = (exposure < r(p50))
generate female = !male
label var male "Males"
label var female "Females"
label var lowexposure "Low exposure"
label var highexposure "High exposure"
foreach z of varlist female male lowexposure highexposure {
capture drop tolhat
quietly regress tol i.id##c.age if `z'
predict tolhat if `z'
xtline tolhat if `z', ///
overlay t(age) i(id) legend(off) ///
ytitle("Predicted tolerance") ylabel(1 (1) 4) ///
addplot(lfit tol age if `z', lwidth(thick) lpattern(solid)) ///
xtitle(Age) title(`: variable label `z'') ///
xsize(3) ysize(2) name(Figure27_`z', replace) nodraw
}
graph combine Figure27_female Figure27_male Figure27_lowexposure Figure27_highexposure, col(2) name(Figure27, replace)
// Figure 2.8
use `figure28', clear
quietly correlate initial male
local pmca = round(r(rho), .01)
quietly correlate initial exposure
local pmcb = round(r(rho), .01)
quietly correlate change male
local pmcc = round(r(rho), .01)
quietly correlate change exposure
local pmcd = round(r(rho), .01)
label define male -1 " " 1 "Male" 0 "Female" 2 " ", modify
label val male male
twoway dot initial male, ///
ytitle("Predicted intercept") ///
xtitle("Gender") ///
xlabel(-1 0 1 2, val) ///
ysize(2) xsize(3) ///
note("{it:r} = `pmca'") ///
name(Figure28A, replace) nodraw
twoway scatter initial exposure, ///
ytitle("Predicted intercept") ///
xtitle("Exposure") ///
ysize(2) xsize(3) ///
note("{it:r} = `pmcb'") ///
name(Figure28B, replace) nodraw
twoway dot change male, ///
ytitle("Predicted change") ///
xtitle("Gender") ///
ysize(2) xsize(3) ///
xlabel(-1 0 1 2, val) ///
note("{it:r} = `pmcc'") ///
name(Figure28C, replace) nodraw
twoway scatter change exposure, ///
ytitle("Predicted change") ///
xtitle("Exposure") ///
ysize(2) xsize(3) ///
note("{it:r} = `pmcd'") ///
name(Figure28D, replace) nodraw
graph combine Figure28A Figure28B Figure28C Figure28D, col(2) row(2) name(Figure28, replace)
Reference
Singer, Judith D., and John B. Willett. 2003. Applied Longitudinal Data Analysis. Modeling Change and Event Occurrence. Oxford University Press. doi: 10.1093/acprof:oso/9780195152968.001.0001
Feb 22, 2016
Random graphs (60): Line plots
import delimited C:\internet\isoc_bdek_di_1_Data.csv, clear
// Generate country variable
kountryadd "Germany (until 1990 former territory of the FRG)" to "Germany" add
kountry geo, from(other) stuck marker
ren _ISO3N_ country
kountry country, from(iso3n) to(iso2c)
ren _ISO2C_ country_str
replace country_str = "UK" if country_str == "GB"
replace country_str = "EU25" if geo == "European Union (25 countries)"
replace country_str = "EU27" if geo == "European Union (27 countries)"
replace country_str = "EU28" if geo == "European Union (28 countries)"
replace country_str = "EU15" if geo == "European Union (15 countries)"
drop if geo == "Euro area (EA11-2000, EA12-2006, EA13-2007, EA15-2008, EA16-2010, EA17-2013, EA18-2014, EA19)"
drop MARKER country
*list country_str geo
// Fix variables
drop ind_type // Drop constant
replace value = "" if value == ":" // Fix missing data indicator
destring value, replace // Convert to numeric
encode indic_is, gen(indic) // Convert from string
label var time "Time"
// Keep relevant cases
keep if unit == "Percentage of individuals"
drop if inlist(country_str, "CH", "EU25", "EU28", "ME", "RS")
twoway (line value time if indic == 1, by(country_str, note("{it:Source:} Eurostat, isoc_bdek_di", span))), ///
ytitle("% who access internet at least once a week") ///
xlabel(2003 2005 2010 2015, ang(h) alternate) xtitle("") ///
name(graph1, replace)
// All in one plot
encode country_str, gen(country)
keep if indic == 1
xtset country time, yearly
keep if inlist(country, 6, 7, 8, 10, 16, 19, 25, 29, 30, 34) // Reduce number of countries
xtline value, overlay xlabel(2003 2005 2010 2015, ang(h)) xtitle("") ///
ytitle("% who access internet at least once a week") ///
legend(pos(2)) name(graph2, replace) ///
note("{it:Source:} Eurostat, isoc_bdek_di", span)xlabel(2003 2005 2010 2015, ang(h) alternate) xtitle("") ///
name(graph2, replace)
Nov 2, 2014
Random graphs (35): Line plots
clear
import excel "art births across countries.xlsx", sheet("Data") firstrow case(lower) clear
ren a geo // Rename country ID
// Rename variables: perc_YEAR format
foreach v of varlist b-o {
local x : variable label `v'
rename `v' perc`x'
}
// Reshape data set to long format
reshape long perc, i(geo) j(year)
// Cross-sectional bar graph at latest point in time
by geo: egen lastobs = max(cond(!missing(perc), year, .)) // Identify year of latest observation
gen byte last = year == lastobs // ID variable of latest observation
// According to http://www.stata.com/support/faqs/data-management/first-and-last-occurrences/
graph hbar perc if year >= 2007 & last == 1, over(geo, sort(perc) descending) nofill ///
ytitle("% ART live births among all births") ylabel(0(1)6) ///
note("{it:Note:} Data refer to 2010 or latest available (no earlier than 2007)", span) ///
/*title("Country differences in ART birth rates")*/ ///
name(crosssec, replace)
// Create classical twoway line plot
twoway (line perc year if geo == "DK", cmissing(no)) ///
(line perc year if geo == "SI", cmissing(no)) ///
(line perc year if geo == "BE", cmissing(no)) ///
(line perc year if geo == "SE", cmissing(no)) ///
(line perc year if geo == "DE", cmissing(no)) ///
(line perc year if geo == "FR", cmissing(no)) ///
(line perc year if geo == "IT", cmissing(no)) ///
(line perc year if geo == "UK", cmissing(no)) ///
(line perc year if geo == "US", cmissing(no)) ///
, xlabel(1997(1)2010, ang(h) labsize(small)) xtitle("") ///
ylabel(0(1)6) ytitle("% ART live births among all births") ///
legend(label(1 "DK") label(2 "SI") label(3 "BE") ///
label(4 "SE") label(5 "DE") label(6 "FR") ///
label(7 "IT") label(8 "UK") label(9 "US") ///
pos(2))
//cmissing(no) makes sure that gaps are not filled
// Use panel line plot functionality keep if inlist(geo, "DK", "SI", "BE", "SE", "DE", "FR", "IT", "UK", "US") encode geo, gen(country) xtset country year, yearly xtline perc, overlay xlabel(1997(1)2010, ang(h) labsize(small)) xtitle("") /// ylabel(0(1)6) ytitle("% ART live births among all births") ///
legend(pos(2)) // Note that the gaps are much harder to spot in this plot
// Small multiples version
twoway (line perc year, by(geo, note("") /*title("Development of ART birth rate, selected countries")*/) cmissing(no)), ///
xlabel(1998 (2) 2010, ang(h) labsize(small) alternate) xtitle("") ///
ylabel(0(1)6) yscale(r(0)) ytitle("% ART live births among all births") ///
name(longit, replace)
// yscale(r(0)) adds missing gridline according to
// http://www.stata.com/statalist/archive/2013-04/msg00904.html
Jul 2, 2014
Random graphs (24): Small-multiples and overlayed line plots
clear
// SG.GEN.PARL.ZS - Proportion of seats held by women in national parliaments (%)
wbopendata, indicator(SG.GEN.PARL.ZS) clear long
// Keep European countries
keep if regioncode == "ECS" // "Europe & Central Asia (all income levels)"
// Keep EU-28
keep if inlist(iso2code, "AT", "BE", "BG", "CY", "CZ", "DE", "DK", "EE", "ES") | ///
inlist(iso2code, "FI", "FR", "GB", "GR", "HR", "HU", "IE", "IT", "LI") | ///
inlist(iso2code, "LT", "LU", "LV", "NL", "PL", "PT", "RO", "SE", "SI") | ///
inlist(iso2code, "SK") // EU-28
// Identify and rename salient variables
ren sg_gen_parl_zs women
label var women "% of national parliament seats held by women"
ren iso2code cntry
drop countrycode region regioncode
// Reshape to drop empty observations
reshape wide women, i(cntry) j(year)
dropmiss women1960-women2014, force
reshape long
// Calculate EU-28 mean
preserve
collapse (mean) women, by(year)
gen cntry = "{bf:EU-28}"
tempfile eu28
save `eu28', replace
restore
append using `eu28'
// Plot small-multiples figure
twoway scatter women year ///
, by(cntry ///
, note(" ") ///
caption("{it:Source:} World Development Indicators, date of extraction: 2015-05-25", span size(small))) ///
cmissing(no) ///
connect(l) ///
msymbol(o) ///
xtitle("") xlabel(1990 1997 2000 2005 2010 2014, ang(45)) ///
ylabel(0 (10) 50) ///
name(small_multiples, replace)
// Plot overlayed figure
encode cntry, gen(country) // No strings for -xtline-
xtset country year, yearly // Declare panel
xtline women if cntry != "{bf:EU-28}" ///
, overlay ///
xlabel(1990 1997 2000 2005 2010 2014) xtitle("Year") ///
plot1opts(lpattern(dash)) plot2opts(lpattern(dash)) ///
plot3opts(lpattern(dash)) plot4opts(lpattern(dash)) ///
plot5opts(lpattern(dash)) plot6opts(lpattern(dash)) ///
plot7opts(lpattern(dash)) plot8opts(lpattern(dash)) ///
plot9opts(lpattern(dash)) plot10opts(lpattern(dash)) ///
plot11opts(lpattern(dash)) plot12opts(lpattern(dash)) ///
plot13opts(lpattern(dash)) plot14opts(lpattern(dash)) ///
plot15opts(lpattern(dash)) plot16opts(lpattern(dash)) ///
plot17opts(lpattern(dash)) plot18opts(lpattern(dash)) ///
plot19opts(lpattern(dash)) plot20opts(lpattern(dash)) ///
plot21opts(lpattern(dash)) plot22opts(lpattern(dash)) ///
plot23opts(lpattern(dash)) plot24opts(lpattern(dash)) ///
plot25opts(lpattern(dash)) plot26opts(lpattern(dash)) ///
plot27opts(lpattern(dash)) plot28opts(lpattern(dash)) ///
name(overlay, replace) ///
addplot(line women year if cntry == "{bf:EU-28}", ///
lwidth(thick) lpattern(solid) sort(cntry year)) ///
legend(order(29 "EU-28 (unweighted average)") pos(6) ring(0)) ///
caption("{it:Source:} World Development Indicators, date of extraction: 2015-05-25", span size(small)) ///
ytitle("% of national parliament seats" "held by women, EU-28 countries")
Subscribe to:
Posts (Atom)