Random graphs (148): Dynamite plot

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)

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))

Random graphs (121): Overlaid prediction plots

// Open cumulated Allbus 2004-2014
use if year >= 2004 using "ZA4584_v1-0-0.dta", clear

// Select data
keep if inrange(v729, 25, 65)  // Only those between 25 and 65
drop if v495 == 0 // Health split

// Poor self-rated health
generate poorhealth = inlist(v495, 4, 5) if !inlist(v495, 0, 9)

// Top-bottom
recode v152 (0 99 = .), gen(sss)
label define sss 1 `""1" "Bottom""' 2 "2" 3 "3" 4 "4" 5 "5" 6 "6" 7 "7" 8 "8" 9 "9" 10 `""10" "Top""'
label val sss sss
label var sss "Subjective SES"

// Education
recode v767 (94 99 = .), gen(educ)
label define educ 1 "Basic" 2 "Lower secondary" 3 "Upper secondary" 4 "Post-secondary" 5 "Tertiary"
label val educ educ
label var educ "Education"

// ISEI 1988
generate isei = v789 if !inlist(v789, 0, 99)
replace  isei = v825 if missing(isei) & !inlist(v825, 0, 99)
label var isei "Occupational status (ISEI)"

// EGP
recode v784 (0 = .a) (10004 = .b) (10009 = .c) (1 2 = .d) // ISCO-88
recode v820 (0 = .a) (10004 = .b) (10009 = .c) (1 2 = .d) // Last job ISCO-88
clonevar isco88 = v784
replace  isco88 = v820 if missing(isco88) & !missing(v820)
recode v804 (0 = .a "Not available") (9999 = .b "No answer"), gen(supervisor) // Number of people supervised
recode v773 (0 = .a "Not available") (99 = .b "No answer") ///
            (10/24 = 1 "Self-employed") (40/65 = 0 "Employed") ///
            (30 70/74 = .c "Not employed"), gen(selfemployed)  // Self-employed
iskoegp egp, isko(isco88) supvis(supervisor) sempl(selfemployed)

// Household income
generate income = v924+1 if !inlist(v924, 99996, 99997, 99999)
generate lninc  = log(income)
label var lninc "Household income (logged, equivalized)"

// Controls
generate state = v1374
generate female = (v731 == 2)
generate age = v729 if v729 != 999

// Macro for controls
local control i.female c.age##c.age i.state##i.year

// Fit some models
eststo clear
qui eststo: logit poorhealt sss                     `control', robust cluster(year) 
qui margins, over(sss) post
est store model1
eststo: logit poorhealt sss i.educ                  `control', robust cluster(year) 
qui margins, over(sss) post
est store model2
eststo: logit poorhealt sss i.educ isei             `control', robust cluster(year) 
qui margins, over(sss) post
est store model3
eststo: logit poorhealt sss i.educ isei i.egp       `control', robust cluster(year) 
qui margins, over(sss) post
est store model4
eststo: logit poorhealt sss i.educ isei i.egp lninc `control', robust cluster(year) 
qui margins, over(sss) post
est store model5

// Plot findings
coefplot model1 model2 model3 model4 model5, vertical ///
         legend(order(2 "Age, federal state, survey year" 4 "plus Education" ///
                      6 "plus occupational status (ISEI)" ///
                      8 "plus social class (EGP)" ///
                      10 "plus household income") pos(2) ring(0) col(1) ///
                title(Model accounting for)) ///
         scheme(plotplainblind) ///
         xtitle(Subjective socioeconomic status) ///
         ytitle(Predicted probability of poor health)

Random graphs (119): Line plots and dot plots

// Prepare ESS round 8
use cntry essround wrkctra dweight using "ESS8e01.dta", clear
recode wrkctra (6 = .a) (7 = .b) (8 = .c) (9 = .d)

// Add ESS rounds 1-7
append using "ESS1-7e01.dta", keep(cntry essround wrkctr wrkctra dweight)

  // Prepare variables
generate nocontract = (wrkctra == 3) if !missing(wrkctra)
kountry cntry, from(iso2c)
rename NAMES_STD country

  // Prepare files for post commands
tempname foo1
tempname foo2
postfile `foo1' str20 country essround ll nocontract ul using `foo2', replace

qui levelsof country, local(country)

 // Loop
foreach x of local country {
 foreach i of numlist 1/8 {
  capture logit nocontract [pw = dweight] if country == "`x'" & essround == `i'
  if _rc == 0 {
   qui margins [pw = dweight]
   matrix prevs = r(table)
   local prev = prevs[1,1] * 100
   local ll = prevs[5,1] * 100
   local ul = prevs[6,1] * 100
   *di "`x'" _skip(5) `i' _skip(5)  `ll' _skip(5) `prev' _skip(5) `ul'
   post `foo1' ("`x'") (`i') (`ll') (`prev') (`ul')
  }
 }
}
postclose `foo1'

// Plot estimates
use `foo2', clear

  // Fix value label
label define essround 1 "2002" 2 "2004" 3 "2006" 4 "2008" ///
                      5 "2010" 6 "2012" 7 "2014" 8 "2016", modify
label val essround essround

// First plot
sort country essround
twoway (rarea ll ul essround, lcolor(white)) ///
       (connected nocontract essround), ///
        by(country, ///
           note("") ///
           legend(off)) ///
        xlabel(1/8, val ang(90)) ///
        xtitle("") ytitle("% without job contract") ///
        name(figure2a, replace) ysize(9)

// Second plot
reshape wide ll ul nocontract, i(country) j(essround)    

scores average =  mean(nocontract*)
egen order_ = rank(average), unique
labmask order_, value(country)
    
twoway (dot nocontract1 order_, horizontal) ///
       (dot nocontract2 order_, horizontal) ///
       (dot nocontract3 order_, horizontal) ///
       (dot nocontract4 order_, horizontal) ///
       (dot nocontract5 order_, horizontal) ///
       (dot nocontract6 order_, horizontal) ///
       (dot nocontract7 order_, horizontal) ///
       (dot nocontract8 order_, horizontal) ///
       (rspike ul1 ll1 order_, horizontal) ///
       (rspike ul2 ll2 order_, horizontal) ///
       (rspike ul3 ll3 order_, horizontal) ///
       (rspike ul4 ll4 order_, horizontal) ///
       (rspike ul5 ll5 order_, horizontal) ///
       (rspike ul6 ll6 order_, horizontal) ///
       (rspike ul7 ll7 order_, horizontal) ///
       (rspike ul8 ll8 order_, horizontal), ///
        ylabel(1/32, val) ytitle("") xscale(alt) ///
        xtitle("% without job contract") ///
        legend(order(1 "2002" 2 "2004" 3 "2006" 4 "2008" ///
                     5 "2010" 6 "2012" 7 "2014" 8 "2016") ///
               pos(5) ring(0)) ///
        name(figure2b, replace) ysize(9)
   
graph combine figure2a figure2b, ///
      col(2) name(figure2, replace) ///
      note(" " "{it:Source:} European Social Survey 2002-16, weighted data. {it:Note:} Error bands/spikes denote 95% confidence intervals", span  size(*.7)) 

Random graphs (115): Bar graphs

// Data
use ESS7e02_1.dta, clear

// Drinking variable
recode alcfreq (7 = 0 "Never") (6 = 1 "Less than once a month") ///
               (5 = 2 "Once a month") (4 = 3 "2{c 150}3 times a month") ///
               (3 = 4 "Once a week") (2 = 5 "Several times a week") ///
               (1 = 6 "Every day") (77 88 99 = .), gen(alcohol)
label var alcohol "Alcohol consumption"

// Dichotomize alcohol variable
generate drinking = inlist(alcohol, 5, 6) if !missing(alcohol)      

// Country name variable
kountry cntry, from(iso2c) marker
rename NAMES_STD country
  
// Education
recode eisced (1 2   = 0 "Lower") ///
              (3 4 5 = 1 "Medium") ///
     (6 7   = 2 "High") ///
     (55/99 = .), gen(education)
label var education "Education"

// Alcohol consumption by country
histogram alcohol, percent disc horizontal ///
                   by(country, title("{bf:A} Frequency of drinking alcohohl", ///
                         justification(left) bexpand span) ///
                               note("")) ///
                   ylabel(0(1)6, val) ytitle("") ///
                   ysize(8) name(figure1, replace)

// Alcohol consumption by education by country
  // Collect predicted probabilities
capture program drop my_logit
program define my_logit, eclass
    syntax[if]
    marksample touse
    logit drinking if `touse'
    margins if `touse', post
    exit
end

statsby point = _b[_cons] se = _se[_cons], by(country education) clear: my_logit

 // Calculate stuff
replace point = point * 100
replace se = se * 100

gen lb = point - 1.96 * se
gen ub = point + 1.96 * se

  // Plot
twoway (bar point point education) ///
       (rspike lb ub education), ///
    by(country, legend(off) ///
                   title("{bf:B} Drinking by educational attainment", ///
                         justification(left) bexpand span) ///
                   note("")) ///
       xlabel(0 1 2, val ang(v)) ylabel(0 (10) 50) ///
       ytitle("% drinking more than once a week") ///
    ysize(8) ///
    name(figure2, replace)

graph combine figure1 figure2, col(2) ///
                               note(" " "{it:Source:} European Social Survey 2014", ///
                                    span justification(right) bexpand size(*.8))

Comparing AME's and the LPM

use doi hs01 age sex educ mstat using ZA5250_v2-0-0.dta, clear

// Poor health
recode hs01 (1 2 3 = 0 "Non-poor health") ///
            (  4 5 = 1 "Poor health") ///
            (-9 = .), gen(poorhealth)
label var poorhealth "Poor health"

// Female sex
recode sex (2 = 1 "Female") ///
           (1 = 0 "Male"), gen(female)
label var female "Female sex"

// Age
recode age (18/24 = 0 "18-24 y.") ///
           (25/34 = 1 "25-34 y.") ///
           (35/44 = 2 "35-44 y.") ///
           (45/54 = 3 "45-54 y.") ///
           (55/64 = 4 "55-64 y.") ///
           (65/74 = 5 "65-74 y.") ///
           (75/84 = 6 "75-84 y.") ///
           (85/97 = 7 "85-97 y.") ///
           (-32   = . ), gen(agecat)
label var agecat "Age"

// Education
recode educ (1 2 = 0 "Hauptschule or less") ///
            (  3 = 1 "Mittlere Reife") ///
            (4 5 = 2 "(Fach)Hochschulreife or more") ///
            (-41 -9 6 7 = .), gen(education)
label var education "Education"

// Marital status
recode mstat (1 6     = 2 "Married/cohabiting") ///
             (2 3 4 9 = 1 "Divorced, widowed etc.") ///
             (5       = 0 "Never married") ///
             (-9      = .), gen(married)
label var married "Marital status"

// Model 
eststo clear
  // AME
logit poorhealth i.female i.agecat i.education i.married 
margins, dydx(*) post
eststo
  // LPM
eststo: regress poorhealth i.female i.agecat i.education i.married, robust
esttab using amelpm.tex, drop(_cons) mtitle("AME" "LPM") label b(2) se(2) nonumbers booktabs ///
                         title("Predictors of poor self-rated health, Germany 2016 \label{tab1}") ///
                         addnote("\emph{Source}: Allbus 2016, doi: 10.4232/1.12754"  ///
                                 "AME: Average marginal effects, LPM: Linear probability model") ///
                         refcat(0.female "\emph{Sex}" ///
                                0.agecat "\emph{Age}" ///
                                0.education "\emph{Education}" ///
                                0.married "\emph{Marital status}", nolabel) ///
                         alignment(D{.}{.}{-1}) width(0.9\hsize) replace

Random graphs (94): Predicted probabilities and their differences

use 2010_ah.dta

// Prepare variables
recode ahm2010_varhours (1 = 0 "Inflexible") (2 3 4 5 = 1 "Flexible"), gen(flexible)
label variable flexible "Flexible working hours"
decode country, gen(cntry)
gen female = (sex == 2) if !missing(sex) 

// Set up loop for posting results
preserve
levelsof cntry, local(country)

tempname foo
tempname foo2
postfile `foo' str20 cntry sexdiffer lb ub using `foo2', replace

foreach x of local country {
      // Estimate model
  qui logit flexible i.female if cntry == "`x'"             
      // Predict probabilities, the r operator gives
      // differences from the reference (base) level
  qui capture margins r.female if cntry == "`x'", post
  local differ      = 100 * _b[r1vs0.female]
  local differ_loci = 100 * (_b[r1vs0.female] + (1.96 * _se[r1vs0.female]))
  local differ_hici = 100 * (_b[r1vs0.female] - (1.96 * _se[r1vs0.female]))
  post `foo' ("`x'") (`differ') (`differ_loci') (`differ_hici')
}
postclose `foo'

// Plot results
use `foo2', clear
   // Sort by size
egen order_ = rank(-sexdiffer), unique
labmask order_, value(cntry)

twoway (rcap sexdiffer sexdiffer order_, horizontal dsymbol(x)) ///
       (rspike ub lb order_, horizontal ) , ///
        xline(0) ylabel(1/30, val ang(h)) ///
        ytitle("") ///
        xtitle("Gender gap in flexible hours" "(Women minus men)") ///
        xscale(alt) ///
        legend(off) ///
        name(by_sex, replace) 
restore 

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).")

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)

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)

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)

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

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.

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)

Random Graphs (50): Dot graphs with confidence intervals

// Create temporary object

tempname bqlrer

// Define postfile

postfile `bqlrer' count str100 commandline str100 descri str10 entity ///
                  prev prev_loci prev_hici differ differ_loci differ_hici ///
                  loed mide hied ratio_ ratio_loci ratio_hici n ///
                  using "results\results", replace

   // Count variable
local count = 0

  // Levels for loop
levelsof entity, local(entity_levels)
*di `entity_levels'

  // Loop: 1 round per entity
foreach Y of local entity_levels  {
     local entity_level = "`Y'"
     di "`entity_level'"
  *qui sum poorhealth if entity == "`Y'"
  *scalar proportion = r(mean)
  
  // Model 1: OV: PSRH lowest2, adjusted -- Prevalence rate only, education not in model
  qui logit poorhealth c_age c_female if entity == "`Y'"
  di _rc
  if _rc == 0 {
  *prev prev_loci prev_hici
  qui margins if entity == "`Y'", post
     scalar prev = 100 * _b[_cons]
     scalar prev_loci = 100 * (_b[_cons] - 1.96 * _se[_cons])
     scalar prev_hici = 100 * (_b[_cons] + 1.96 * _se[_cons])
  di `prev' 
  di `prev_loci' 
  di `prev_hici'
  }
  
  // Model 2: OV: PSRH lowest2, adjusted
     local count = `count' + 1
     local desc "DV: Poor health (bottom 2), IV: Education, Age, Sex"
     qui logit poorhealth i.education c_age c_female if entity == "`Y'"
  local commandline = e(cmdline)
  di _rc
  if _rc == 0 {
  
  // Differences Model 2
  qui margins r.education if entity == "`Y'"
  scalar differ      = -10 * _b[r3vs1.education]
  scalar differ_loci = -10 * (_b[r3vs1.education] + (1.96 * _se[r3vs1.education]))
  scalar differ_hici = -10 * (_b[r3vs1.education] - (1.96 * _se[r3vs1.education]))

  // Predicted probabilities Model 2
  qui margins i.education if entity == "`Y'", post
  matrix preds = r(table)
  scalar loed = preds[1,1]
     scalar mide = preds[1,2]
     scalar hied = preds[1,3]
     matrix drop preds
  
  // Ratios Model 2
  qui nlcom _b[1.education]/_b[3.education], post
  scalar ratio_     = _b[_nl_1]
  scalar ratio_loci = _b[_nl_1] - 1.96 * _se[_nl_1]
  scalar ratio_hici = _b[_nl_1] + 1.96 * _se[_nl_1]
  
  }
  else {
  scalar differ = -99
  scalar differ_loci = -99
  scalar differ_hici = -99
  scalar ratio_ = -99
   scalar ratio_loci = -99
     scalar ratio_hici = -99
  scalar loed = -99
  scalar mide = -99
  scalar hied = -99  
  }
     post `bqlrer' (`count') (`"`commandline'"') ("`desc'") ("`Y'") ///
                   (prev) (prev_loci) (prev_hici) ///
                   (differ) (differ_loci) (differ_hici) ///
                   (loed) (mide) (hied) ///
                   (ratio_) (ratio_loci) (ratio_hici) (e(N))
}

// Fit one model for entire EVS
// Model 1: OV: PSRH lowest2, adjusted -- Prevalence rate only, education not in model
qui logit poorhealth c_age c_female [pw = eu_weight] if survey == "EVS"
qui margins if survey == "EVS", post
scalar prev = 100 * _b[_cons]
scalar prev_loci = 100 * (_b[_cons] - 1.96 * _se[_cons])
scalar prev_hici = 100 * (_b[_cons] + 1.96 * _se[_cons])
di `prev' 
di `prev_loci' 
di `prev_hici'

// Model 2: OV: PSRH lowest2, adjusted
local count = `count' + 1
local desc "DV: Poor health (bottom 2), IV: Education, Age, Sex"
qui logit poorhealth i.education c_age c_female /*[pw = eu_weight]*/ if survey == "EVS"
local commandline = e(cmdline)

// Differences Model 2
qui margins r.education if survey == "EVS"
scalar differ      = -10 * _b[r3vs1.education]
scalar differ_loci = -10 * (_b[r3vs1.education] + (1.96 * _se[r3vs1.education]))
scalar differ_hici = -10 * (_b[r3vs1.education] - (1.96 * _se[r3vs1.education]))

  // Predicted probabilities Model 2
qui margins i.education if survey == "EVS", post
matrix preds = r(table)
scalar loed = preds[1,1]
scalar mide = preds[1,2]
scalar hied = preds[1,3]
matrix drop preds
  
  // Ratios Model 2
qui nlcom _b[1.education]/_b[3.education], post
scalar ratio_     = _b[_nl_1]
scalar ratio_loci = _b[_nl_1] - 1.96 * _se[_nl_1]
scalar ratio_hici = _b[_nl_1] + 1.96 * _se[_nl_1]

local Y = "EVS"

post `bqlrer' (`count') (`"`commandline'"') ("`desc'") ("`Y'") ///
              (prev) (prev_loci) (prev_hici) ///
              (differ) (differ_loci) (differ_hici) ///
              (loed) (mide) (hied) ///
              (ratio_) (ratio_loci) (ratio_hici) (e(N))
  
postclose `bqlrer'

// Create data set with FIPS codes and state names
preserve
tempfile codescheme
egen pickone = tag(entity)
keep if pickone
keep entity entity_num
save `codescheme'
restore

// Merge FIPS codes to results
use results\results, clear
merge 1:1 entity using `codescheme', keepusing(entity entity_num)
drop _merge

label var count "Generic counter"
label var commandline "Command line"
label var descri "Description of model fitted"
label var entity "Entity (string)"
label var prev   "Proportion poor health"
label var prev_loci "Proportion 95 % CI (lower)"
label var prev_hici "Proportion 95 % CI (higher)"
label var loed   "Proportion poor health lower educ."
label var mide   "Proportion poor health mid educ."
label var hied   "Proportion poor health higher educ."
label var differ "Difference lower educated - higher"
label var differ_loci  "Difference 95 % CI (lower)"
label var differ_hici  "Difference 95 % CI (higher)"
label var ratio_      "Ratio lower educated / higher"
label var ratio_loci  "Ratio 95 % CI (lower)"
label var ratio_hici  "Ratio 95 % CI (higher)"
label var n           "Size of entity sample"

// Save results, also as csv-file
save results\results, replace
outsheet using "results\results.csv", comma replace nolabel


use results\results, replace

// Create US identifier
gen usa = 0 
replace usa = 1 if regexm(entity, "US-")

replace entity = "{bf:Europe}" if entity == "EVS"

egen order_ratio = rank(-ratio_), unique 
labmask order_ratio, value(entity)

egen order_differ = rank(-differ), unique 
labmask order_differ, value(entity)

egen order_prev = rank(-prev), unique
labmask order_prev, value(entity)

twoway (dot prev order_prev if usa == 0, horizontal ndots(20)) ///
       (dot prev order_prev if usa == 1, horizontal ndots(20)) ///
    (rspike prev_loci prev_hici order_prev, horizontal) ///
    , ///
    legend(label(1 "Europe") label(2 "US") label(3 "95% CI") pos(1) ring(0)) ///
    ylabel(1/97, valuelabels ang(h) labsize(*.65)) ///
    ytitle("") ///
    xtitle("Prevalence poor health") ///
    name(prevalence, replace) ///
    xsize(2.7) ysize(10)

generate ratio_loci_p = ratio_loci               // Shorten CI's to make graph fit better
replace  ratio_loci_p =  0 if ratio_loci_p <  0  
generate ratio_hici_p = ratio_hici
replace  ratio_hici_p = 6 if ratio_hici_p > 6
 
twoway (dot ratio_ order_ratio if usa == 0, horizontal ndots(20) msymbol(smplus)) ///
       (dot ratio_ order_ratio if usa == 1, horizontal ndots(20) msymbol(smx)) ///
    (rspike ratio_loci_p ratio_hici_p order_ratio, horizontal) ///
    , ///
    legend(label(1 "Europe") label(2 "US") label(3 "95 % CI") pos(1) ring(0)) ///
    ylabel(1/97, valuelabels ang(h) labsize(*.65)) ///
    xscale(range(0 6)) xlabel(0 1 3 5)  /// 
    ytitle("") ///
    xtitle("Relative" "inequalities") ///
    xline(1) ///
    xline(1.81019, lpattern(dash)) ///
    name(relative, replace) ///
    xsize(2.7) ysize(10)

generate differ_loci_p = differ_loci              // Shorten CI's to make graph fit better
replace  differ_loci_p = -1 if differ_loci_p < -1
    
twoway (dot differ order_differ if usa == 0, horizontal ndots(20) msymbol(smplus)) ///
       (dot differ order_differ if usa == 1, horizontal ndots(20) msymbol(smx)) ///
       (rspike differ_loci_p differ_hici order_differ, horizontal) ///
    , ///
    legend(label(1 "Europe") label(2 "US") label(3 "95 % CI") pos(1) ring(0)) ///
    xscale(range(0 30)) xlabel(0(10)30)  /// 
    ylabel(1/97, valuelabels ang(h) labsize(*.65)) ///
    ytitle("") ///
    xtitle("Absolute" "inequalities") ///
    xline(6.8754, lpattern(dash)) ///
    name(absolute, replace) ///
    xsize(2.7) ysize(10)

twoway (dot prev order_prev if usa == 0, horizontal ndots(20) msymbol(smplus)) ///
       (dot prev order_prev if usa == 1, horizontal ndots(20) msymbol(smx)) ///
    (rspike prev_loci prev_hici order_prev, horizontal) ///
    , ///
    legend(label(1 "Europe") label(2 "US") label(3 "95% CI") row(1) pos(12) size(vsmall) region(lwidth(vthin) lcolor(black)) bmargin(tiny) colgap(*.3)) ///
    ylabel(1/97, valuelabels ang(h) labsize(*.65)) ///
    ytitle("") ///
    xtitle("Prevalence" "poor health") ///
    xline(11.26938, lpattern(dash)) ///
    name(prevalence_leg, replace) ///
    xsize(2.7) ysize(10)    
    
// 8.27 × 11.69
*graph combine prevalence absolute relative, xsize(8.27) ysize(10) col(3) name(comb, replace)

grc1leg prevalence_leg absolute relative, xsize(8.27) ysize(10) col(3) imargin(small) legendfrom(prevalence_leg) name(oneleg, replace) pos(6) span // ignores size-commands

graph display oneleg, xsize(6.27) ysize(9.69)  // Thus redraw so that size commands take effect

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)

Random graphs (46): Plotting the contours of interactions

Next to the classical layout for interaction plots, it is also possible to make use of contour plots for visualizing 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)

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

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

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)