sysuse auto, clear qui regress mpg i.foreign local df = e(df_r) local t = round(((_b[1.foreign]/_se[1.foreign])), .01) local p = round((2 * ttail(`df', abs(`t'))), .001) margins, at(foreign = (0 1)) marginsplot, recast(bar) ytitle("Mileage") xtitle("Car type") /// plotopts(fcolor(gs14) lcolor(black)) /// title("{bf:A} Mileage", justification(left) bexpand span) /// addplot(scatteri 28.95 0 28.95 1, recast(line) lwidth(medium) lpattern(solid) /// text(29.5 .5 "{it:t}(`df') = `t', {it:p} = `p'") /// xlabel(0 1) /// || scatteri 28.95 0 28.05 0, recast(line) lwidth(medium) lpattern(solid) /// || scatteri 28.95 1 28.05 1, recast(line) lwidth(medium) lpattern(solid)) /// legend(off) /// name(figureA, replace) qui regress weight i.foreign local df = e(df_r) local t = round(((_b[1.foreign]/_se[1.foreign])), .01) local double p = round((2 * ttail(`df', abs(`t'))), .0000001) margins, at(foreign = (0 1)) marginsplot, recast(bar) ytitle("Weight") xtitle("Car type") /// plotopts(fcolor(gs14) lcolor(black)) /// title("{bf:B} Weight", justification(left) bexpand span) /// addplot(scatteri 3700 0 3700 1, recast(line) lwidth(medium) lpattern(solid) /// text(3760 .5 "{it:t}(`df') = `t', {it:p} = `p'") /// xlabel(0 1) /// || scatteri 3700 0 3600 0, recast(line) lwidth(medium) lpattern(solid) /// || scatteri 3700 1 3600 1, recast(line) lwidth(medium) lpattern(solid)) /// legend(off) /// name(figureB, replace) graph combine figureA figureB, col(2) name(figure, replace)
Showing posts with label marginsplot. Show all posts
Showing posts with label marginsplot. Show all posts
Jul 26, 2021
Random graphs (148): Dynamite plot
Apr 23, 2020
Random graphs (144): Automatically labeling figures with letters
// Open ESS round 5
use agea stflife cntry if inlist(cntry, "DE", "GB", "NL", "SE") using "ESS5e03_4.dta", clear
// Country variable
kountry cntry, from(iso2c)
rename NAMES_STD country
// Restrict to complete cases
keep if !missing(agea, stflife, cntry)
// Center age
qui sum agea, detail
replace agea = `r(p99)' if agea > `r(p99)' & !missing(agea)
bys country: center agea, gen(age)
drop agea
label var age "Age (centered)"
// Letters
tokenize "`c(ALPHA)'"
local i 1
levelsof country, local(country)
foreach z of local country {
qui sum age if country == `"`z'"', detail
foreach x of numlist 10 25 50 75 90 {
local x`x': di %9.1f r(p`x')
di `x`x''
}
qui regress stflife c.age##c.age if country == `"`z'"'
qui margins, at(age = (`x10' `x25' `x50' `x75' `x90'))
marginsplot, title(`"{bf:``i''} `z'"') ///
ytitle("Life satisfaction" "(predicted)") ///
xtitle("Age (percentiles)") ///
xlabel(`x10' "10th" ///
`x25' "25th" ///
`x50' "50th" ///
`x75' "75th" ///
`x90' "90th") ///
ylabel(, format(%6.1f)) ///
recastci(rarea) ciopts(color(gs10)) ///
name("``i''", replace) nodraw
local ++i
}
graph combine A B C D, col(2) ycommon
Mar 23, 2018
Random graphs (129): Offsetting markers
// Open Allbus 2016
use eastwest hs01 id02 using ZA5250_v2-0-0.dta, clear
// Region
recode eastwest (1 = 0 "West Germany") (2 = 1 "East Germany"), gen(east)
label var east "East Germany"
// Subjective social class
recode id02 (-50/-7 = .) ///
( 1 = 0 "Lower class") ///
( 2 = 1 "Working class") ///
( 3 = 2 "Middle class") ///
( 4 5 = 3 `" "Upper" "(middle)" "class" "'), gen(class)
label var class "Subjective social class"
// Health
recode hs01 (-9 = .) (5 = 0 "Bad") (4 = 1 "Less good") (3 = 2 "Satisfactory") ///
(2 = 3 "Good") (1 = 4 "Very good"), gen(health)
label var health "Self-rated health"
// Model
regress health i.class##i.east //[pw = wghtpew]
// Plot using marginsplot
margins class#east, saving(myfile, replace)
marginsplot, legend(pos(4) ring(0)) name(marginsplot, replace) title("") ///
ytitle("Predicted self-rated health") ylabel(, format(%6.1f))
// Plot using twoway with overlay
use myfile, clear
// Create offset
clonevar _mx = _m1
replace _mx = cond(_m2 == 2, _mx - 0.1, _mx + 0.1)
// Plot using marginsplot
twoway (rcap _ci_lb _ci_ub _m1 if _m2==0, sort) ///
(rcap _ci_lb _ci_ub _mx if _m2==1, sort) ///
(connected _margin _m1 if _m2==0, msymbol(oh)) ///
(connected _margin _mx if _m2==1, msymbol(o)) ///
, title("") ///
ytitle("Predicted self-rated health") xla(, valuelabel) ///
legend(pos(4) ring(0) col(1) order(3 "West Germany" 4 "East Germany")) ///
name(twoway, replace) ylabel(, format(%6.1f))
Jul 27, 2017
Random graphs (106): Interaction plot with overlaid density
// Open Allbus 2008
use V154 V151 V156 V760 V754 V755 V5 V767 using ZA4602_v1-0-0.dta, clear
// Age
generate age = V154 if V154 != 999
// Sex
generate female = (V151 == 2)
// Migrant
generate migrant = (V156 == 2)
// Interviewer ID
rename V760 id
// Interviewer sex
generate femaleinterviewer = (V754 == 2)
// Interviewer age
generate interviewerage = V755
// Attractiveness rating by interviewer before and after interview
qui factor V5 V767, pcf
predict attractiveness
// Calculate Spearman-Brown for two-item scales according
// to Eisinga et al. (https://doi.org/10.1007/s00038-012-0416-3):
spearman V5 V767
// Fit model with respondents clustered in interviewers
mixed attractiveness c.age##i.female i.migrant i.femaleinterviewer##c.interviewerage || id:
// Calculate margins and plot
qui margins, at(age = (18 (10) 98) female = (0 1))
marginsplot, recast(line) recastci(rarea) ciopt(color(gs14)) ///
plot1opts(lpattern(dash)) ///
legend(ring(0) pos(2)) ///
title("Female beauty premium disappears with age", span) ///
ytitle("Predicted attractiveness", axis(1)) ///
xtitle("Age") ///
xlabel(20 (10) 100) ///
addplot(histogram age, discrete yaxis(2) ///
lcolor(white) ///
ylabel(0 0.01 0.02, format(%6.2f) axis(2)) ///
yscale(alt range(0 0.1) axis(2)) ///
ytitle("Age density", axis(2)) ///
legend(order(3 "Men" 4 "Women"))) ///
note(" " "{it:Source:} German General Social Survey Allbus 2008, doi:10.4232/1.12345", span) ///
name(figure1, replace)
Jul 4, 2017
Chapter 5 of Singer and Willett's (2003) book on longitudinal data analysis
There is another take on the chapter here, but I like mine better.
// Table 5.1 use "C:\singer willett (2003)\reading_pp.dta", clear format age %6.2f list id wave agegrp age piat if inlist(id, 4, 27, 31, 33, 41, 49, 69, 77, 87), sep(0) noobs // Figure 5.1twoway (scatter piat age) /// (lfit piat age) /// (scatter piat agegrp) /// (lfit piat agegrp) /// if inlist(id, 4, 27, 31, 33, 41, 49, 69, 77, 87) /// , by(id, note("")) xtitle("{it:AGE} or {it:AGEGRP}") /// ylabel(0 (20) 80) xlabel(6 (1) 12, format(%6.0f)) ytitle("{it:PIAT}") /// legend(order(1 "Age" 2 "Linear fit age" /// 3 "Target age" 4 "Linear fit target age") col(2)) /// name(figure51, replace) ysize(8) // Table 5.2generate agegrp_65 = agegrp - 6.5 generate age_65 = age - 6.5 capture program drop randomslopetable program define randomslopetable eststo `1' qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance exp([lns1_1_1]_b[_cons]) * /// exp([lns1_1_2]_b[_cons]) qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance end eststo clear mixed piat agegrp_65 || id: agegrp_65, var cov(uns) mle eststo agegrp randomslopetable agegrp mixed piat age_65 || id: age_65, var cov(uns) mle eststo age randomslopetable age esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(agegrp_65 age_65) /// // Align coefficients coeflabels(_cons "Intercept" age_65 "Change") /// // Label coefficients order(_cons age_65) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(2 2 2 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("AGEGRP - 6.5" "AGE - 6.5") /// keep(piat:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.3 use "C:\singer willett (2003)\wages_pp.dta", clear list id exper lnw black hgc uerate if inlist(id, 206, 332, 1028), sep(0) noobs // Table 5.4eststo clear mixed lnw exper || id: exper, var mle cov(uns) eststo modela randomslopetable modela mixed lnw c.exper##c.hgc_9 c.exper##black || id: exper, var mle cov(uns) eststo modelb randomslopetable modelb mixed lnw exper hgc_9 c.exper#black || id: exper, var mle cov(uns) eststo modelc randomslopetable modelc esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// 1.black "BLACK" exper "Change" /// c.exper#c.hgc_9 "Change x (HGC - 9)" /// 1.black#c.exper "Change x BLACK") /// // Label coefficients order(_cons hgc_9 1.black /// exper c.exper#c.hgc_9 1.black#c.exper) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.2estimates restore modelc margins, at(exper = (0 (1) 10) black = (0 1) hgc_9 = (0 3)) marginsplot, ytitle("Predicted log hourly wage") title("") /// recastci(rarea) ciopt(color(gs14)) /// xtitle("{it:EXPER}") /// addplot(scatteri 2.35 10 "White", msymbol(none) mlabpos(0) /// || scatteri 2.23 10 "White", msymbol(none) mlabpos(0) /// || scatteri 2.185 10 "Black", msymbol(none) mlabpos(0) /// || scatteri 2.07 10 "Black", msymbol(none) mlabpos(0) /// || scatteri 1.75 0 "9th grade", msymbol(none) mlabpos(0) /// || scatteri 1.88 0 "12th grade", msymbol(none) mlabpos(0)) /// legend(off) ylabel(, format(%6.1f)) /// name(figure52, replace) // Table 5.5use "C:\singer willett (2003)\wages_small_pp.dta", clear eststo clear mixed lnw hgc_9 exper c.exper#black || id: exper, var mle cov(uns) eststo modela randomslopetable modela mixed lnw hgc_9 exper c.exper#black || id: , var mle cov(uns) eststo modelc qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(agegrp_65 age_65) /// // Align coefficients coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// 1.black "BLACK" exper "Change" /// 1.black#c.exper "Change x BLACK") /// // Label coefficients order(_cons hgc_9 1.black /// exper 1.black#c.exper) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.6 use "C:\singer willett (2003)\unemployment_pp.dta", clear list id months cesd unemp if inlist(id, 7589, 55697, 67641, 65441, 53782), sepby(id) noobs // Table 5.7eststo clear mixed cesd months || id: months, var mle cov(uns) eststo modela randomslopetable modela mixed cesd months i.unemp || id: months, var mle cov(uns) eststo modelb randomslopetable modelb mixed cesd c.months##i.unemp || id: months, var mle cov(uns) eststo modelc randomslopetable modelc version 10 // This one is a challenge to fit generate unempXmonths = unemp * months xtmixed cesd unemp unempXmonths || id: unemp unempXmonths, var mle cov(uns) eststo modeld qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v4 = exp(2*[lns1_1_2]_b[_cons]) // Unemp variance qui estadd scalar v3 = exp(2*[lns1_1_3]_b[_cons]) // Unemp x months variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance version 14 esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(unemp 1.unemp unempXmonths 1.unemp#c.months) /// coeflabels(_cons "Intercept" months "Change" /// 1.unemp "UNEMP" /// 1.unemp#c.months "Change x UNEMP") /// // Label coefficients order(_cons months /// 1.unemp 1.unemp#c.months) /// // Order coefficients stats(v_e v1 v2 v3 v4 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Var(UNEMP)" /// "Var(UNEMP x TIME)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C" "Model D") /// keep(cesd:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.3estimates restore modelb margins, at(months = (0 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("{it:Predicted CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("") /// name(figure53, replace) ** All other plots of Figure 5.3 are only variants of this one ** // Figure 5.4 estimates restore modelb margins, at(months = (0 (1) 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("Model B") subtitle("Main effects of" /// "{it:UNEMP} and {it:TIME}") /// name(figure54A, replace) estimates restore modelc margins, at(months = (0 (1) 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("Model C") subtitle("Interaction between" /// "{it:UNEMP} and {it:TIME}") /// name(figure54B, replace) estimates restore modeld // not sure how to do this with marginsplot predict pd twoway (line pd months if unemp == 0, c(L)) /// (line pd months if unemp == 1, c(L)) /// (scatteri 11.5 15 "Employed", msymbol(none) mlabpos(9)) /// (scatteri 15 15 "Unemployed", msymbol(none) mlabpos(9)) /// , ylabel(5 (5) 20) legend(off) /// xlabel(0 (2) 14) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// title("Model D") subtitle("Constraining the effect of {it:TIME}" /// "among the re-employed") /// name(figure54C, replace) graph combine figure54A figure54B figure54C, row(1) xsize(11) // Table 5.8 use "C:\singer willett (2003)\wages_pp.dta", clear eststo clear mixed lnw hgc_9 ue_7 exper c.exper#black || id: exper, mle cov(un) var eststo modela randomslopetable modela mixed lnw hgc_9 ue_mean ue_person_centered exper c.exper#black || id: exper, mle cov(un) var eststo modelb randomslopetable modelb mixed lnw hgc_9 ue1 ue_centert1 exper c.exper#black || id: exper, mle cov(un) var eststo modelc randomslopetable modelc esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(ue_mean ue_7 ue1 ue_7 ue_centert1 ue_person_centered) /// coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// ue_7 "UERATE" ue_person_centered "Deviation UERATE" /// exper "Change" 1.exper#black "Change x BLACK") /// // Label coefficients order(_cons hgc_9 ue_7 ue_person_centered /// exper c.exper#black) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.9 use "C:\singer willett (2003)\medication_pp.dta", clear list wave day timeofday time time333 time667 in 1/11, noobs sep(0) // Table 5.10 eststo clear mixed pos i.treat##c.time || id: time, mle cov(uns) var eststo modela randomslopetable modela mixed pos i.treat##c.time333 || id: time333, mle cov(uns) var eststo modelb randomslopetable modelb mixed pos i.treat##c.time667 || id: time667, mle cov(uns) var eststo modelc randomslopetable modelc esttab /// , b(2) se(2) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(time333 time time667 time /// 1.treat#c.time333 1.treat#c.time /// 1.treat#c.time667 1.treat#c.time) /// coeflabels(_cons "Intercept" time "Change" /// 1.treat "TREAT" 1.treat#c.time "Change x TREAT") /// // Label coefficients order(_cons 1.treat time 1.treat#c.time) /// // Order coefficients stats(v_e v1 v2 cov dev aic bic nc N, /// Add variance components to table fmt(2 2 2 2 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Cov(Init., Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(pos:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.5estimates restore modela margins, at(time = (0 (1) 7) treat = (1 0)) marginsplot, recastci(rarea) ciopts(color(gs14)) legend(off) /// title("") xtitle("Days") ytitle("Predicted {it:POS}") /// addplot(scatteri 187 7 "Treatment", msymbol(none) mlabpos(11) /// || scatteri 153 7 "Control", msymbol(none) mlabpos(11) /// xlabel(0 (1) 7)) /// name(figure55, 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
Jun 21, 2017
Chapter 4 of Singer and Willett's (2003) book on longitudinal data analysis
There is another take on the chapter here, but I like mine better.
use alcohol1_pp.dta, clear // Figure 4.1twoway (scatter alcuse age) (lfit alcuse age) /// if inlist(id, 4, 14, 23, 32, 41, 56, 65, 82), /// by(id, row(2) legend(off) note("")) /// xlabel(13 (1) 17) ytitle({it:ALCUSE}) xtitle({it:AGE}) /// name(figure41, replace) // Figure 4.2preserve bysort id: sample 32, count qui sum peer generate hipeer = (peer >= r(mean)) qui regress alcuse i.id##c.age predict predicted_alcuse label var predicted_alcuse "Predicted {it:ALCUSE}" xtline predicted_alcuse if coa == 0, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("{it:COA} = 0") ylabel(-1 (1) 4) /// name(coa0, replace) nodraw xtline predicted_alcuse if coa == 1, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("{it:COA} = 1") ylabel(-1 (1) 4) /// name(coa1, replace) nodraw xtline predicted_alcuse if hipeer == 0, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("Low {it:PEER}") ylabel(-1 (1) 4) /// name(hipeer0, replace) nodraw xtline predicted_alcuse if hipeer == 1, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("High {it:PEER}") ylabel(-1 (1) 4) /// name(hipeer1, replace) nodraw graph combine coa0 coa1 hipeer0 hipeer1, col(2) ysize(8) name(figure42, replace) restore // Table 4.1generate time = age - 14 // Create program that gets all the random slope parameters for esttab capture program drop randomslopetable program define randomslopetable eststo `1' qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance exp([lns1_1_1]_b[_cons]) * /// exp([lns1_1_2]_b[_cons]) qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance capture drop p1 // For R-squared based on observed and predicted outcomes qui predict p1 // For R-squared based on observed and predicted outcomes qui regress alcuse p1 // For R-squared based on observed and predicted outcomes local orsq = e(r2) // For R-squared based on observed and predicted outcomes estimates drop . // For R-squared based on observed and predicted outcomes estimates restore `1' // For R-squared based on observed and predicted outcomes qui estadd scalar orsq = `orsq' // For R-squared based on observed and predicted outcomes scalar v_e = exp(2*[lnsig_e]_b[_cons]) // For R-squared based on residual variance reduction scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar rsqe = (v_eumm - v_e) / v_eumm // For R-squared based on residual variance reduction of unconditional means model if !missing(`v2_ugm') { qui estadd scalar rsqv2 = (v2_ugm - v2) / v2_ugm // For R-squared based on residual variance reduction of unconditional means model qui estadd scalar rsqv1 = (v1_ugm - v1) / v1_ugm // For R-squared based on residual variance reduction of unconditional means model } end // This was based on: https://www.statalist.org/forums/forum/general-stata-discussion/general/1309801-saving-commands-as-macro-variables-local-and-global-to-make-do-files-shorter eststo clear mixed alcuse || id: , variance mle // Model A, unconditional means model eststo modela qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance scalar v_eumm = exp(2*[lnsig_e]_b[_cons]) // Residual variance for R-squared mixed alcuse time || id: time, cov(un) variance mle // Model B, unconditional growth model randomslopetable modelb scalar v2_ugm = exp(2*[lns1_1_1]_b[_cons]) // Slope variance scalar v1_ugm = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance mixed alcuse i.coa##c.time || id: time, cov(un) variance mle // Model C randomslopetable modelc mixed alcuse i.coa##c.time c.peer##c.time || id: time, cov(un) variance mle // Model D randomslopetable modeld mixed alcuse i.coa c.peer##c.time || id: time, cov(un) variance mle // Model E randomslopetable modele mixed alcuse i.coa c.cpeer##c.time || id: time, cov(un) variance mle // Model F randomslopetable modelf mixed alcuse ccoa c.cpeer##c.time || id: time, cov(un) variance mle // Model G randomslopetable modelg esttab /// , b(3) se(3) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define starts rename(ccoa 1.coa cpeer peer c.cpeer#c.time c.peer#c.time) /// // Align coefficients coeflabels(1.coa "COA" 1.coa#c.time "COA x time" peer "PEER" /// c.peer#c.time "PEER x time" _cons "Intercept" time "Time") /// // Label coefficients order(_cons 1.coa peer ) /// // Order coefficients stats(v_e v1 v2 cov orsq rsqe rsqv1 rsqv2 dev aic bic nc N, /// Add variance components to table fmt(3 3 3 3 3 3 3 3 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Cov(Init., Change)" /// "R-squared obs./pred." /// "R-squared var(Residual)" /// "R-squared var(Initial)" /// "R-squared var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C" "Model D" "Model E" "Model F" "Model G") /// keep(alcuse:) // Drop variance components in weird shapes // ICC of p. 96: estimates restore modela estat icc // Figure 4.3estimates restore modelb qui margins, at(time = (0 1 2)) marginsplot, xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") ylabel(0 .5 1 1.5 2) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Unconditional growth model", span) name(modelb, replace) nodraw estimates restore modelc qui margins, at(time = (0 1 2) coa = (0 1)) marginsplot, xlabel(-1 "13" 0 "14" .5 " " 1 "15" 1.5 " " 2 "16" 3 "17", ) ylabel(0 (1) 2, ) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Uncontrolled effects of {it:COA}", span) legend(off) /// addplot(scatteri .9 2 "{it:COA} = 0" 1.5 2 "{it:COA} = 1", msymbol(none)) /// name(modelc, replace) nodraw estimates restore modele margins, at(time = (0 1 2) coa = (0 1) peer = (.655 1.381)) marginsplot, xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17", format(%6.0f)) ylabel(0 1 2, ) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Controlled effects of {it:COA}", span bexpand) legend(off) /// addplot(scatteri .95 2 "{it:COA} = 0" /// 1.5 2 "{it:COA} = 1" /// .14 -0.85 "Low {it:PEER}" /// .64 -0.85 "High {it:PEER}" /// .72 -0.85 "Low {it:PEER}" /// 1.21 -0.85 "High {it:PEER}", msymbol(none)) /// name(modele, replace) nodraw graph combine modelb modelc modele, ycommon xcommon col(3) xsize(12) ysize(6) altshrink name(figure43, replace) // Figure 4.4preserve statsby, by(id) saving(temp, replace): regress alcuse time merge m:1 id using temp qui cor _b_cons coa local r = round(r(rho), .01) twoway (scatter _b_cons coa) /// (scatteri 4 .5 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(-.5 " " 0 1 1.5 " ") ytitle(Intercept) /// xtitle("{it:COA}") legend(off) name(g1, replace) nodraw qui cor _b_cons peer local r = round(r(rho), .01) twoway (scatter _b_cons peer) /// (scatteri 4 1 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(0 1 2 3) ytitle(Intercept) /// xtitle("{it:PEER}") legend(off) name(g2, replace) nodraw qui cor _b_time coa local r = round(r(rho), .01) twoway (scatter _b_time coa) /// (scatteri 4 .5 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(-.5 " " 0 1 1.5 " ") ytitle(Change) /// xtitle("{it:COA}") legend(off) name(g3, replace) nodraw qui cor _b_time peer local r = round(r(rho), .01) twoway (scatter _b_time peer) /// (scatteri 4 1 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(0 1 2 3) ytitle(Change) /// xtitle("{it:PEER}") legend(off) name(g4, replace) nodraw graph combine g1 g2 g3 g4, name(figure44, replace) restore // Figure 4.5estimates restore modelf predict e, resid predict e_time e_cons, reffects relevel(id) qnorm e, yline(0) ytitle(Var(residual)) name(g1, replace) nodraw qnorm e_cons, yline(0) ytitle(Var(intercept)) name(g2, replace) nodraw qnorm e_time, yline(0) ytitle(Var(change)) name(g3, replace) nodraw egen ze = std(e) egen ze_time = std(e_time) egen ze_cons = std(e_cons) scatter ze id, yline(0) ytitle("Standardized" "var(residual)") name(g4, replace) nodraw scatter ze_cons id, yline(0) ytitle("Standardized" "var(intercept)") name(g5, replace) nodraw scatter ze_time id, yline(0) ytitle("Standardized" "var(change)") name(g6, replace) nodraw graph combine g1 g2 g3 g4 g5 g6, cols(2) colfirst name(figure45, replace) ysize(8) // Figure 4.6twoway (scatter e time), ylabel(-2(1)2) ytitle("Var(residual)") /// xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") /// xtitle("{it:TIME}") /// yline(0) name(g1, replace) nodraw twoway (scatter e_cons coa), ylabel(-2(1)2) ytitle("Var(intercept)") /// xlabel(-.5 " " 0 1 1.5 " ") xtitle("{it:COA}") /// yline(0) name(g2, replace) nodraw twoway (scatter e_cons peer), ylabel(-2(1)2) ytitle("Var(intercept)") /// xlabel(0 1 2 3) xtitle("{it:PEER}") /// yline(0) name(g3, replace) nodraw twoway (scatter e_time coa), ylabel(-2(1)2) ytitle("Var(change)") /// xlabel(-.5 " " 0 1 1.5 " ") xtitle("{it:COA}") /// yline(0) name(g4, replace) nodraw twoway (scatter e_time peer), ylabel(-2(1)2) ytitle("Var(change)") /// xlabel(0 1 2 3) xtitle("{it:PEER}") /// yline(0) name(g5, replace) nodraw graph combine g1 g2 g3 g4 g5, cols(2) hole(2) name(figure46, replace) ysize(8) // Figure 4.7qui reg alcuse time coa cpeer c.cpeer#c.time predict pa estimates restore modelf gen bayes0 = _b[_cons] + _b[1.coa]*coa + _b[cpeer]*cpeer + e_cons gen bayes1 = _b[time] + _b[c.cpeer#c.time]*cpeer + e_time gen bayes = bayes0 + bayes1 * time twoway (scatter alcuse time) /// (lfit alcuse time) /// (line pa time, sort) /// (line bayes time, sort) /// if inlist(id, 4, 14, 23, 32, 41, 56, 65, 82), by(id, col(4) note("")) /// ytitle({it:ALCUSE}) xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") /// legend(order(1 "Data points" 2 "OLS" /// 3 "Population average" 4 "Bayes estimate") col(4)) /// xtitle("{it:TIME}") name(figure47, 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
May 29, 2016
Random graphs (88): Model predictions
use gndr fcldbrn yrbrn health agea ///
fltdpr flteeff slprl wrhpp fltlnl enjlf fltsd cldgng cntry ///
using ESS3e03_5.dta, clear
// Gender variable
gen female = (gndr == 2) if gndr !=9
// Calculate age at first birth
recode fcldbrn (6666 = .a "NA") ///
(7777 = .b "Refusal") ///
(8888 = .c "Don't know") ///
(9999 = .d "No answer") ///
, gen(yrfirstbirth)
recode yrbrn (7777 = .b "Refusal") ///
(8888 = .c "Don't know") ///
(9999 = .d "No answer") ///
, gen(yrbirth)
generate afb = yrfirstbirth - yrbirth
qui centile afb, centile(.5 99) // Truncate
replace afb = r(c_1) if afb < r(c_1)
replace afb = r(c_2) if afb >= r(c_2) & !missing(afb)
// Self-rated health variable
recode health ( 1 = 4 "Very good") ///
( 2 = 3 "Good") ///
( 3 = 2 "Fair") ///
( 4 = 1 "Bad") ///
( 5 = 0 "Very bad") ///
( 7 = .a "Refusal") ///
( 8 = .b "Don't know") ///
( 9 = .c "No answer") ///
, generate(srh)
// Depression variable
recode fltdpr flteeff slprl wrhpp fltlnl ///
enjlf fltsd cldgng ///
(7 = .a) (8 = .b) (9 = .c)
recode wrhpp enjlf (1 = 4) (2 = 3) (3 = 2) (4 = 1) (. = .)
scores depression = mean(fltdpr flteeff slprl wrhpp fltlnl enjlf fltsd cldgng)
// Age variable
recode agea (999 = .a)
// Models and plots
regress srh i.afb i.agea if female == 1, cluster(cntry)
qui margins, over(afb)
marginsplot, ytitle("Predicted self-rated health") xtitle("Age at first birth") ///
name(womensrh, replace) recastci(rarea) ciopts(color(gs12)) xsize(3) ///
title("") ylabel(2 (.2) 3.2, format(%6.1f)) nodraw
regress depression i.afb i.agea if female == 1, cluster(cntry)
qui margins, over(afb)
marginsplot, ytitle("Predicted depression") xtitle("Age at first birth") ///
name(womendep, replace) recastci(rarea) ciopts(color(gs12)) xsize(3) ///
title("") ylabel(1.2 (.2) 2.4, format(%6.1f)) nodraw
regress srh i.afb i.agea if female == 0, cluster(cntry)
qui margins, over(afb)
marginsplot, ytitle("Predicted self-rated health") xtitle("Age at first birth") ///
name(mensrh, replace) recastci(rarea) ciopts(color(gs12)) xsize(3) ///
title("") ylabel(2 (.2) 3.2, format(%6.1f)) nodraw
regress depression i.afb i.agea if female == 0, cluster(cntry)
qui margins, over(afb)
marginsplot, ytitle("Predicted depression") xtitle("Age at first birth") ///
name(mendep, replace) recastci(rarea) ciopts(color(gs12)) xsize(3) ///
title("") ylabel(1.2 (.2) 2.4, format(%6.1f)) nodraw
graph combine womensrh womendep, col(1) ysize(8) name(a, replace) ///
title("Women", span) nodraw
graph combine mensrh mendep, col(1) ysize(8) name(b, replace) ///
title("Men", span) nodraw
graph combine a b, col(2) ysize(8) xsize(8) ycommon xcommon ///
note("{it:Source:} European Social Survey, round 3." ///
"{it:Note:} All models control for age as dummy variables. 95% CI's based on robust SE's", span size(*.9))
May 26, 2016
Random graphs (86): Interaction plots
version 14
use ESS1e06_4.dta, clear
// Age variable
replace agea = . if agea == 999
qui centile(agea), centile(99.9)
replace agea = r(c_1) if age > r(c_1) & !missing(agea)
rename agea age
// Happiness variable
recode happy (77 88 99 = .)
// Education variable
recode eduyrs (77 88 99 = .)
qui centile eduyrs, centile(99.9)
replace eduyrs = r(c_1) if eduyrs > r(c_1) & !missing(eduyrs)
// Fit model
regress happy c.eduyrs##c.age##c.age, robust
// Caculate education slope as function of age
qui margins, dydx(eduyrs) over(age)
marginsplot, recastci(rarea) recast(line) ciopts(color(gs12)) ///
ytitle("Effect of education on happiness") ylabel(, format(%6.2f)) ///
xtitle("Age") xlabel(20(10)90) xmtick(14(1)93) ///
title("") name(slope, replace) xsize(4) nodraw
// Predict happiness for different values of age and education
qui margins, at(eduyrs=(9 12 15) age=(14 20 (5) 90 93))
marginsplot, xdimension(age) recastci(rarea) recast(connected) ciopts(color(gs12)) ///
ytitle("Predicted happiness") ylabel(, format(%6.1f)) ///
xtitle("Age") xlabel(20(10)90) xmtick(14(1)93) ///
legend(order(6 "15 yrs" 5 "12 yrs." 4 "9 yrs.") ///
title(Education, size(*.8)) pos(11) ring(0)) ///
title("") name(prediction, replace) xsize(4) nodraw
// Plot
graph combine slope prediction, xsize(8) row(1) note("{it:Source:} European Social Survey, round 1 (2002/2003).")
May 19, 2016
Random graphs (83): Quadratic by continuous interaction
unzipfile ZA4582_v1-0-0.dta.zip
use ZA4582_v1-0-0.dta, clear
// Fix variables
// Life satisfaction
recode v545 ( 0 = 0 "Completely dissatisfied") ///
(10 = 10 "Completely satisfied") ///
(96 = .a "Not asked") ///
(99 = .b "No response") ///
, gen(lsat)
label var lsat "Life satisfaction"
// Household income (OECD)
recode v924 (99996 = .d "Can't be determined") ///
(99997 = .a "Refused") ///
( 0 = 0 "No income") ///
(99998 = .b "Don't know") ///
(99999 = .c "No response") ///
, gen(income)
label var income "Income"
// Age, truncated at 99.9th percentile
recode v729 (997 = .a "Refused") ///
(999 = .b "No Response"), gen(age)
label var age "Age"
qui centile(age), centile(99.9)
replace age = r(c_1) if age > r(c_1) & !missing(age)
// Fit model
regress lsat c.age##c.age##c.income
// Calculate income percentiles
centile(income), centile(10 25 50 75 90)
local r10 = r(c_1)
local r25 = r(c_2)
local r50 = r(c_3)
local r75 = r(c_4)
local r90 = r(c_5)
// Calculate margins
margins, at(age=(20 (5) 90) income=(`r10' `r50' `r75' `r90'))
// Plot
marginsplot, recastci(rarea) ciopts(color(gs12)) ///
ytitle("Predicted life satisfaction") ///
legend(order(8 "90th percentile" ///
7 "75th percentile" ///
6 "50th percentile" ///
5 "10th percentile") ring(0) pos(11) ///
title("Equivalized household income", size(*.7))) ///
xtitle("Age") ylab(, format(%6.1f)) title("") ///
note(" " "{it:Source:} German General Social Survey (Allbus), doi: 10.4232/1.12439" ///
"{it:Note:} Gray-shaded areas denote 95% confidence intervals.", span)
erase ZA4582_v1-0-0.dta
May 11, 2016
Random graphs (82): Quadratic by categorical interaction
unzipfile output7138637360988071942.zip
use agea gndr happy using ESS1-6e01_1_F1, clear
// Generate age variable, truncated at 99.9th percentile
qui centile(agea), centile(99.9)
generate age = agea
replace age = r(c_1) if age > r(c_1) & !missing(age)
// Generate gender variable
generate female =(gndr == 2) if !missing(gndr)
// Fit model and plot
regress happy c.age##c.age##i.female
qui margins, at(age=(15(5)90) female=(0 1))
marginsplot, recastci(rarea) ciopts(color(gs12)) ytitle("Predicted happiness") ///
legend(order(3 "Men" 4 "Women") ring(0) pos(2)) xtitle("Age") ///
ylab(, format(%6.1f)) title("") ///
note(" " "{it:Source:} European Social Survey, round 1{c 150}6." ///
"{it:Note:} Gray-shaded areas denote 95% confidence intervals.", span)
May 8, 2016
Random graphs (81): Plot average and confidence interval for each value of other variable
use pspwght happy gndr agea if inrange(agea, 15, 86) using ESS1-6e01_1_F1.dta, clear
// Create gender variable
generate female = (gndr == 2) if !missing(gndr)
label define female 0 "Men" 1 "Women"
label val female female
// Fit model
regress happy i.agea##i.female [pweight = pspwght]
// Calculate means and confidence intervals
margins, over(agea female)
// Plot means and confidence intervals
marginsplot, recastci(rarea) ciopts(color(gs12)) recast(line) ylabel(, format(%6.1f)) ///
xtitle("Age") ytitle("Average happiness") ///
title("Happiness over the life course") ///
legend(pos(2) ring(0)) xlabel(15 (10) 86) ///
plot1opts(lpattern(dash)) plot2opts(lpattern(solid)) ///
note(" " "{it:Source:} European Social Survey, rounds 1{c 150}6, weighted data" ///
"{it:Note:} Gray areas denote 95% confidence intervals", span)
Apr 4, 2016
Random graphs (70): Bar graphs for interaction plots
unzipfile "output7717034280080927434.zip", replace
use ESS1-6e01_0_F1, clear
*use cntry essround stfjbot wkhct using ESS1-6e01_0_F1, clear
// Select Round 5
keep if inlist(essround, 5)
// Fix country variable
encode cntry, gen(country)
// Satisfaction with WLB
rename stfjbot swlb
// Contracted working hours
recode wkhct ( 0/20 = 1 "Marginal part-time") ///
( 21/34.75 = 2 "Substantial part-time") ///
(35/220 = 0 "Full-time") ///
, gen(parttime)
label var parttime "Working hours"
// Professional status
generate professional = .
replace professional = (iscoco >= 1000 & iscoco <= 2999)
label var professional "Professional status"
label define professional 0 "Non-prof." 1 "Professional"
label val professional professional
// Sex
recode gndr ( 2 = 1 "Female") ///
( 1 = 0 "Male") ///
(.a = .a "No answer") ///
, gen(female) label(female)
label var female "Gender"
// Select cases: Those in paid work and with a partner
keep if partner == 1 & mnactic == 1
// Are part-time workers more satisfied with their WLB than full-time employees?
qui regress swlb i.parttime i.country, cluster(country)
estimates store m1
// Hypothesis 1b: The effect is stronger for marginal PT than substantial PT.
qui test 1.parttime == 2.parttime
local F = round(r(F), .001)
local p = round(r(p), .001)
local r = r(df)
qui margins, at(parttime=(0 1 2))
marginsplot, recast(bar) ///
title("SWLB difference between full-time and part-time workers", size(large)) ///
plotopts(fcolor(gs14) lcolor(black)) ytitle("Predicted SWLB") ylabel(6 (.5) 7.5, format(%6.1f) grid) ///
name(noint, replace) xtitle("") ysize(3) ///
note(" " ///
"Marginal and substantial part-time" ///
"differs significantly:" ///
"{it:F}(`r', `r(df_r)') = `F', {it:p} = `p'", ///
pos(11) ring(0) bmargin(small)) ///
nodraw
// Are professional part-time workers less satisfied than non-professional part-time workers?
qui regress swlb i.parttime##i.professional i.country, cluster(country)
estimates store m2
qui margins, at(parttime=(0 1 2) professional=(0 1))
marginsplot, recast(bar) xdimension(professional) ///
bydimension(parttime) byopts(row(1) noiyaxes imargin(zero) ///
title("Interaction part-time status and professional status")) ///
subtitle(, pos(6)) /// // Place label of by dimensions below plot
plotopts(fcolor(gs14) lcolor(black)) ytitle("Predicted SWLB") ylabel(6 (.5) 7.5, format(%6.1f)) ///
name(int1, replace) xtitle("") ysize(3) nodraw
// Are women working part-time more satisfied with their SWLB than part-time working men?
qui regress swlb i.parttime##i.female i.country, cluster(country)
estimates store m3
qui margins, at(parttime=(0 1 2) female=(0 1))
marginsplot, recast(bar) xdimension(female) ///
bydimension(parttime) byopts(row(1) noiyaxes imargin(zero) ///
title("Interaction part-time status and gender")) ///
subtitle(, pos(6)) /// // Place label of by dimensions below plot
plotopts(fcolor(gs14) lcolor(black)) ytitle("Predicted SWLB") ylabel(6 (.5) 7.5, format(%6.1f)) ///
name(int2, replace) xtitle("") ysize(3) nodraw
// Output table and figure
esttab m1 m2 m3 using test.tex, compress replace se label nomtitles ///
indicate(Country dummies = *country) ///
varwidth(30) interaction(" X ") ///
title(Regression table\label{tab1}) ///
booktabs
graph combine noint int1 int2, col(1) ysize(9) ///
note("95% CI's based on cluster-robust standard errors")
Feb 4, 2016
Interactions with and without main effects
set scheme s1mono
clear
set obs 10000
set seed 1
gen a = rnormal()
gen b = rnormal()
gen e = 2*rnormal()
gen y = 2 + 3*a + 4*b + 5*a*b + e
summarize y a b
estimates clear
regress y c.a##c.b
estimates store M1
regress y c.a c.b
estimates store M2
regress y c.b c.a#c.b
estimates store M3
regress y c.a c.a#c.b
estimates store M4
regress y c.a#c.b
estimates store M5
estimates restore M1
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M1, replace) noci ///
title("M1: Main terms and interaction") legend(col(1))
estimates restore M2
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M2, replace) noci ///
title("M2: Main terms, no interaction") legend(off)
estimates restore M3
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M3, replace) noci ///
title("M3: Interaction, one main term missing") legend(off)
estimates restore M4
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M4, replace) noci ///
title("M4: Interaction, other main term missing") legend(off)
estimates restore M5
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M5, replace) noci ///
title("M5: Interaction term, both main terms missing") legend(off)
grc1leg M1 M2 M3 M4 M5, col(2) ysize(11) xsize(8) ring(0) pos(5) ycommon
estimates table M1 M2 M3 M4 M5, b(%9.2f) stat(F r2_a)
Oct 6, 2015
Random graphs (54): Cross-level interactions
use deleteme.dta
// Cross-level interaction model
qui mixed lsat c.wfc##c.gdp || country: wfc , cov(uns)
estimates store cross_lvl
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Create table esttab cross_lvl /// , se /// stats(v1 v2 v_e cov dev nc N, /// Add variance components to table labels("Var(Intercept)" /// "Var(Slope)" /// "Var(Residual)" /// "Cov(Int., Slope)" /// "Deviance" /// "No. clusters" /// "No. individuals")) /// coeflabel(wfc "Work-family conflict" /// gdp "GDP" /// c.wfc#c.gdp "Work-family conflict X GDP" /// _cons "Intercept") /// mtitles("Cross-level interaction") /// nonumbers varwidth(28) modelwidth(24) /// keep(lsat:) // Drop variance components in weird shapes // Interaction plot (A) preserve estimates restore cross_lvl qui sum wfc // get standard deviation and mean local wfcminsd = r(mean) - r(sd) local wfcmean = r(mean) local wfcplusd = r(mean) + r(sd) qui sum gdp // get standard deviation local gdpmin2sd = r(mean) - 2*r(sd) local gdpminsd = r(mean) - r(sd) local gdpmean = r(mean) local gdpplusd = r(mean) + r(sd) local gdpplu2sd = r(mean) + 2*r(sd) margins, at(gdp = (`gdpmin2sd' `gdpminsd' `gdpmean' `gdpplusd' `gdpplu2sd') /// wfc = (`wfcminsd' `wfcmean' `wfcplusd')) vsquish marginsplot, x(gdp) noci /// xlabel(`gdpmin2sd' "-2 SD" /// `gdpminsd' "-1 SD" /// `gdpmean' `""Average" "GDP""' /// `gdpplusd' "+1 SD" /// `gdpplu2sd' "+2 SD") /// xtitle("") /// ylabel(, format(%6.1f)) /// ytitle("Predicted happiness") /// plotopts(msymbol(none)) /// // Turn off markers plot1opts(lpattern(longdash)) /// // Define line types here plot2opts(lpattern(solid)) /// plot3opts(lpattern(shortdash)) /// legend(subtitle("Work{c 150}family conflict" , size(small)) /// order(1 "- 1 SD" 2 "Mean" 3 "+ 1 SD") size(small) /// pos(11) ring(0)) /// title("(A) Interaction plot") /// name(xlvl_plot1, replace) restore // Interaction plot (B) preserve estimates restore cross_lvl predict u1 u0, reffect capture drop pickone egen pickone = tag(country) keep if pickone gen wfceffect = u1 + _b[wfc] + _b[c.wfc#c.gdp] * gdp qui sum gdp // get standard deviation local gdpmin2sd = r(mean) - 2*r(sd) local gdpminsd = r(mean) - r(sd) local gdpmean = r(mean) local gdpplusd = r(mean) + r(sd) local gdpplu2sd = r(mean) + 2*r(sd) twoway (scatter wfceffect gdp, mlabel(country) mlabpos(0) msymbol(none)) /// (function y = _b[wfc] + _b[c.wfc#c.gdp] * x, range(gdp)) /// , ytitle("Work{c 150}family conflict coefficient") /// xlabel(`gdpmin2sd' "-2 SD" `gdpminsd' "-1 SD" /// `gdpmean' `""Average" "GDP""' /// `gdpplusd' "+1 SD" `gdpplu2sd' "+2 SD") /// title("(B) Interaction plot") /// legend(order(2 "Work{c 150}family conflict coefficient by GDP") /// pos(7) ring(0) size(small)) /// name(xlvl_plot2, replace) restore // Combine graphs graph combine xlvl_plot1 xlvl_plot2, row(1) ysize(3) xsize(5.5) altshrink
Sep 29, 2015
Scaling predictors for models with interactions
The scaling of predictors is an important issue for the interpretation of models with interaction terms. The Table below presents a simple OLS regression model of happiness regressed on age, sex, and their interaction.
When none of the predictors are standardized and 0 is not a meaningful value on a predictor, interpretation of the regression coefficients can become awkward. For instance, in the model with none of the predictors standardized, the coefficient for female – .348 – states that women score .348 points higher in terms of happiness than men among those who are 0 years old. The coefficient for age – -.004 – indicates that with every additional year, men's happiness decreases by .004 points. The coefficient for the interaction term – -.008 – denotes that women's happiness decrease more than men's with each additional year of age. The intercept is the average happiness of men who are 0 years old.
A first step towards giving the coefficients from an interaction model an easier interpretation is to standardize continuous variables. For the second model shown in the Table, age was standardized by subtracting its mean and dividing it by its standard deviation. Thus, for the standardized variable, the value of 0 refers to the average age, and a change of 1 refers to a change of one standard deviation. The age coefficient (-.07) now indicates that for every additional standard deviation of age, men's happiness decreases by .07. For men and women of average age, the happiness differential is .02. Women's happiness decreases by .15 per standard deviation of age stronger than men's. The intercept shows the average happiness of average-age men.
While it is a rather controversial technique, it is in principle also possible to standardize the dummy variable indicating the difference between men and women. Then, the intercept is the overall mean of happiness. However, an increase of 1 in the gender coefficient becomes difficult to interpret, as there is no such thing as a one-standard deviation increase in gender. When the dummy variable is approximately evenly distributed, then the coefficient reflects half of the difference in the outcome between the two groups.
Gelman and Hill (2007, p. 56) suggest to standardize variables by two rather than one standard deviation to keep binary and continuous variables on a roughly common scale. An increase of one in age then refers to the difference between one standard deviation below the mean and one standard deviation above the mean; an increase of one in age reflects roughly the difference between men and women.
As the plot above shows, none of these different scaling decisions affect the form of the interaction.
When none of the predictors are standardized and 0 is not a meaningful value on a predictor, interpretation of the regression coefficients can become awkward. For instance, in the model with none of the predictors standardized, the coefficient for female – .348 – states that women score .348 points higher in terms of happiness than men among those who are 0 years old. The coefficient for age – -.004 – indicates that with every additional year, men's happiness decreases by .004 points. The coefficient for the interaction term – -.008 – denotes that women's happiness decrease more than men's with each additional year of age. The intercept is the average happiness of men who are 0 years old.
A first step towards giving the coefficients from an interaction model an easier interpretation is to standardize continuous variables. For the second model shown in the Table, age was standardized by subtracting its mean and dividing it by its standard deviation. Thus, for the standardized variable, the value of 0 refers to the average age, and a change of 1 refers to a change of one standard deviation. The age coefficient (-.07) now indicates that for every additional standard deviation of age, men's happiness decreases by .07. For men and women of average age, the happiness differential is .02. Women's happiness decreases by .15 per standard deviation of age stronger than men's. The intercept shows the average happiness of average-age men.
While it is a rather controversial technique, it is in principle also possible to standardize the dummy variable indicating the difference between men and women. Then, the intercept is the overall mean of happiness. However, an increase of 1 in the gender coefficient becomes difficult to interpret, as there is no such thing as a one-standard deviation increase in gender. When the dummy variable is approximately evenly distributed, then the coefficient reflects half of the difference in the outcome between the two groups.
Gelman and Hill (2007, p. 56) suggest to standardize variables by two rather than one standard deviation to keep binary and continuous variables on a roughly common scale. An increase of one in age then refers to the difference between one standard deviation below the mean and one standard deviation above the mean; an increase of one in age reflects roughly the difference between men and women.
Reference
Gelman, Andrew, and Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.use agea gndr happy using ESS1e06.3_F1.dta, clear
// Set up variables
rename agea age
recode gndr (1 = 0) (2 = 1) (.a = .), gen(female)
label var female "Female sex (ref. male)"
label var age "Age"
label var happy "Happiness"
// Listwise deletion
drop if missing(female, age, happy)
// Standardize variables
qui sum female
gen female_1z = (female - r(mean)) / r(sd)
gen female_2z = (female - r(mean)) / (2*r(sd))
qui sum age
gen age_1z = (age - r(mean)) / r(sd)
gen age_2z = (age - r(mean)) / (2*r(sd))
// Fit models
reg happy c.age##female
estimates store unstandardized
reg happy c.age_1z##female
estimates store cont_standardized
reg happy c.age_1z##c.female_1z
estimates store standardized_1z
reg happy c.age_2z##c.female_2z
estimates store standardized_2z
// Create table
esttab unstandardized cont_standardized standardized_1z standardized_2z, ///
rename(age_1z age ///
age_2z age ///
female_1z 1.female ///
female_2z 1.female ///
c.age_1z#c.female_1z 1.female#c.age ///
c.age_2z#c.female_2z 1.female#c.age ///
1.female#c.age_1z 1.female#c.age) ///
drop(0b.female 0b.female#co.age 0b.female#co.age_1z) ///
coeflabel(age "Age" ///
1.female "Female (ref. male)" ///
1.female#c.age "Age X Female" ///
_cons "Intercept") ///
mtitles("Unstandardized" "Age standardized (1 SD)" ///
"Standardized (1 SD)" "Standardized (2 SD)") ///
nonumbers varwidth(18) modelwidth(24) ///
se r2 obslast
// Plot models
estimates restore unstandardized
qui sum age
local age_minsd = r(mean) - r(sd)
local age_mean = r(mean)
local age_plusd = r(mean) + r(sd)
qui margins, at(c.age = (`age_minsd' `age_mean' `age_plusd') ///
female = (0 1)) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
xlabel(`age_minsd' "-1 SD" `age_mean' "Average age" `age_plusd' "+1 SD") ///
title("Unstandardized variables", size(*.9)) ///
ytitle("Predicted happiness") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(dash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
legend(subtitle(Gender, size(*.8)) ///
order(3 "Males" 4 "Females") ring(0) pos(7) ///
size(*.8) symysize(*.7) symxsize(*.7)) ///
name(unstandardized, replace)
estimates restore cont_standardized
qui margins, at(c.age_1z = (-1 0 1) ///
female = (0 1)) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
xlabel(-1 "-1 SD" 0 "Average age" 1 "+1 SD") ///
title("Age standardized at 1 SD", size(*.9)) ///
ytitle("Predicted happiness") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(dash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
legend(subtitle(Gender, size(*.8)) ///
order(3 "Males" 4 "Females") ring(0) pos(7) ///
size(*.8) symysize(*.7) symxsize(*.7)) ///
name(cont_standardized, replace)
estimates restore standardized_1z
qui sum female
local fem_hi = (1 - r(mean)) / r(sd)
local fem_lo = (0 - r(mean)) / r(sd)
qui margins, at(c.age_1z = (-1 0 1) ///
c.female_1z = (`fem_lo' `fem_hi')) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
xlabel(-1 "-1 SD" 0 "Average age" 1 "+1 SD") ///
title("All variables standardized at 1 SD", size(*.9)) ///
ytitle("Predicted happiness") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(dash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
legend(subtitle(Gender, size(*.8)) ///
order(3 "Males" 4 "Females") ring(0) pos(7) ///
size(*.8) symysize(*.7) symxsize(*.7)) ///
name(standardized_1z, replace)
estimates restore standardized_2z
qui sum female
local fem_hi = (1 - r(mean)) / (2*r(sd))
local fem_lo = (0 - r(mean)) / (2*r(sd))
qui margins, at(c.age_2z = (-.5 0 .5) ///
c.female_2z = (`fem_lo' `fem_hi')) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
xlabel(-.5 "-1 SD" 0 "Average age" .5 "+1 SD") ///
title("All variables standardized at 2 SD", size(*.9)) ///
ytitle("Predicted happiness") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(dash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
legend(subtitle(Gender, size(*.8)) ///
order(3 "Males" 4 "Females") ring(0) pos(7) ///
size(*.8) symysize(*.7) symxsize(*.7)) ///
name(standardized_2z, replace)
graph combine unstandardized cont_standardized standardized_1z standardized_2z, ///
row(2) col(2) ysize(7) name(combined, replace)
May 24, 2015
Random graphs (48): Interaction plot with overlaid data points
// Simulate data
clear
set seed 1
set obs 500
generate e = 0 + (500 - 0) * runiform() // To generate random variates over the
generate x1 = 0 + (800 - 0) * runiform() // interval [a,b), a+(b-a)*runiform()
generate x2 = round(0 + (1 - 0) * runiform()) // Binary variable
generate y = (x1 + x2 + (x1*x2) + e) / 1000
// Calculations for Aiken & West (1991) style plot
regress y c.x1##x2
qui sum x1
local x1_minsd = r(mean) - (2*r(sd))
local x1_mean = r(mean)
local x1_plusd = r(mean) + (2*r(sd))
// Calculate predicted values for plot
margins, at(c.x1 = (`x1_minsd' `x1_mean' `x1_plusd') ///
c.x2 = (0 1)) vsquish
// Plot
qui marginsplot, recastci(rarea) ciopts(color(gs10)) ///
title("Interaction plot with overlaid data points") ///
ytitle("y") ///
ylabel(, format(%6.1f)) ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(longdash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
addplot(scatter y x1 ///
, symbol(o) ///
xlabel(`x1_minsd' "-2 SD" /// // The addplot seems to override
`x1_mean' "Average x1" /// // the regular axis label
`x1_plusd' "+2 SD") /// // command
legend(subtitle(x2) /// // And the legend command
order(3 "x2 = 0" 4 "x2 = 1" 2 "95 % CI"))) ///
legend(pos(5) ring(0)) /// // But not completely
name(plot, replace)
May 6, 2015
Random graphs (46): Plotting the contours of interactions
// Simulate data
clear
set seed 1
set obs 100
generate e = 0 + (100000 - 0) * runiform() // To generate random variates over the
generate x1 = 0 + (600 - 0) * runiform() // interval [a,b), a+(b-a)*runiform()
generate x2 = 0 + (600 - 0) * runiform()
generate y = (x1 + x2 + (x1*x2) + e) / 1000
// Calculations for Aiken & West (1991) style plot
qui regress y c.x1##c.x2
qui sum x1
local x1_minsd = r(mean) - r(sd)
local x1_mean = r(mean)
local x1_plusd = r(mean) + r(sd)
qui sum x2
local x2_minsd = r(mean) - r(sd)
local x2_mean = r(mean)
local x2_plusd = r(mean) + r(sd)
// Calculate predicted values for plot
margins, at(c.x2 = (`x2_minsd' `x2_mean' `x2_plusd') ///
c.x1 = (`x1_minsd' `x1_mean' `x1_plusd')) vsquish
// Plot
qui marginsplot, recastci(rarea) ciopts(color(gs10)) ///
xlabel(`x2_minsd' "-1 SD" `x2_mean' "Average x2" `x2_plusd' "+1 SD") ///
title("Aiken & West (1991)-style plot for interactions") ///
ytitle("Predicted y") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(longdash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
plot3opts(lpattern(shortdash)) ///
legend(subtitle(x1) ///
order(4 "- 1 SD" 5 "Average x1" 6 "+ 1 SD" 2 "95% CI")) ///
name(lines, replace)
// For some reason, marginsplot ignores the -label option-, therefore
// -order- is used here
// Calculations for contour plot
qui regress y c.x1##c.x2
predict y_hat
// Plot
qui twoway contour y_hat x1 x2, ///
title("Contour plot for interactions") ///
ztitle("Predicted y") ///
name(contours, replace)
// Combine figures
graph combine lines contours, col(1) ysize(8)
Mar 17, 2014
Random graphs (19): Simple interaction plots with -marginsplot- for continuous variables
sysuse auto, clear
regress mpg c.weight##c.gear_ratio
// Obtain means and standard deviation for plots
qui sum weight
local w_minsd = r(mean) - r(sd)
local w_mean = r(mean)
local w_plusd = r(mean) + r(sd)
qui sum gear_ratio
local g_minsd = r(mean) - r(sd)
local g_mean = r(mean)
local g_plusd = r(mean) + r(sd)
// Calculate predicted values for plot
margins, at(c.gear_ratio = (`g_minsd' `g_mean' `g_plusd') ///
c.weight = (`w_minsd' `w_mean' `w_plusd')) vsquish
// Plot
marginsplot, recastci(rarea) ciopts(color(gs10)) ///
title("Weight moderates the gear ratio{char 0150}mileage association") ///
xlabel(`g_minsd' "-1 SD" `g_mean' "Average gear ratio" `g_plusd' "+1 SD") ///
ytitle("Predicted mileage") ///
xtitle("Gear ratio") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(longdash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
plot3opts(lpattern(shortdash)) ///
legend(subtitle(Weight) ///
order(4 "- 1 SD" 5 "Mean" 6 "+ 1 SD" 2 "95% CI"))
// For some reason, marginsplot ignores the -label option-, therefore
// -order- is used here
gr export plot.png, as(png) replace
Feb 5, 2014
Random graphs (18): Simple interaction plots with -marginsplot- for categorical and continuous variables
sysuse auto, clear
reg mpg c.weight##i.foreign
qui sum weight
local minsd = r(mean) - r(sd)
local mean = r(mean)
local plusd = r(mean) + r(sd)
margins foreign, at(weight = (`minsd' `mean' `plusd'))
marginsplot, recast(line) noci scheme(lean1) ///
name(first, replace)
gr export first.png, as(png) replace
// Impossible to tell the difference between
// domestic and foreign lines because of the lean scheme
marginsplot, recast(line) noci ///
plot1opts(lpattern(dash)) ///
scheme(lean1) ///
name(second, replace)
gr export second.png, as(png) replace
// That fixes the problem
marginsplot, recast(line) recastci(rarea) ciopts(color(gs10)) ///
plot1opts(lpattern(dash)) ///
scheme(lean1) ///
name(third, replace)
gr export third.png, as(png) replace
// That's a nice way of plotting CI's
marginsplot, recast(line) recastci(rarea) ciopts(color(gs10)) ///
plot1opts(lpattern(dash)) ///
ytitle("Predicted mileage") ///
title("Car type moderates the weight{char 0150}mileage association") ///
xlabel(`minsd' "-1 SD" `mean' "Average weight (lbs.)" `plusd' "+1 SD") ///
xtitle("") ///
scheme(lean1) ///
name(fourth, replace)
gr export fourth.png, as(png) replace
// That labels the graph nicely
Mar 31, 2013
Random graphs (13): Interaction plots using -marginsplot-
estimates use m2 // Open estimates saved via -estimates save-
xtcloglog // Show results
// Set esample() as explained in the help-file to -estimates save-:
mark miss_mod1
markout miss_mod1 fail2 female age1519 age2529 age3034 ///
hatlev1d loggdp hatfield_main
estimates esample: if miss_mod1
// Calculate predicted probabilities:
margins female#hatlev1d, pred(pu0)
// Simple default marginsplot:
marginsplot, xtitle("Gender") ytitle("Predicted probability") ///
caption("{it:Source:} EU-LFS 2009 AHM", span) ///
title("") name(sexXeduc1, replace)
// As a bar graph:
marginsplot, xtitle("") ytitle("Predicted probability") title("") ///
recast(bar) /// // Draw bar graph instead of default
xdimension(female hatlev1d) ///
/* Get both interacting variables onto x-axis */ ///
caption("{it:Source:} EU-LFS 2009 AHM", span) ///
name(sexXeduc2, replace)
// As a bar graph that better separates the different variable combinations:
marginsplot, xtitle("") ytitle("Predicted probability") ///
recast(bar) ///
xdimension(female) ///
bydimension(hatlev1d, ///
elabels(1 "Low education" ///
2 "Medium education" ///
3 "High education")) ///
/* Use -bydimension- to better structure the x-axis */ ///
name(sexXeduc3, replace)
// As a proper-looking graph:
marginsplot, xtitle("") ytitle("Predicted probability") ///
recast(bar) ///
xdimension(female) ///
bydimension(hatlev1d, elabels(1 "Low education" ///
2 "Medium education" ///
3 "High education")) ///
byopts(row(1) /// All plots in one row
noiyaxes /// // Don't display individual y-axes
imargin(zero) /// // No margin between by() dimensions
title("") /// // Suppress title
caption("{it:Note:} Error bars are 95 % CI's" ///
"{it:Source:} EU-LFS 2009 AHM", span)) ///
/* Caption here to suppress displaying it three times */ ///
plotregion(lwidth(none)) ///
/* Removes lines around separate by dimension plots */
subtitle(, pos(6)) /// // Place label of by dimensions below plot
plotopts(fcolor(gs14) lcolor(b)) ///
/* Options of the recast() type can go here as well */ ///
name(sexXeduc4, replace)
Subscribe to:
Posts (Atom)