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)

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

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

Dropping one group at a time

// Interaction plot excluding one country at a time
foreach n of numlist 1/35 {  /// Countries are numbered from 1 to 35

  preserve
  qui drop if country == `n'      // Drop one country
  local naam: label country `n'   // Save name in a local
  *di "`naam'"

  qui mixed happiness c.workfamconf##c.gdppc || country: workfamconf, cov(uns)
  predict u1 u0, reffect

  capture drop pickone
  egen pickone = tag(country)
  qui keep if pickone

  gen wfceffect = u1 + _b[workfamconf] + _b[c.workfamconf#c.gdppc] * gdppc

  qui sum gdppc // 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) mlabsize(*.8)) ///
         (function y = _b[workfamconf] + _b[c.workfamconf#c.gdppc] * x, range(gdppc)) ///
        , ytitle("WFC coefficient", size(*.6)) ///
          xlabel(`gdpmin2sd' "-2 SD" `gdpminsd' "-1 SD" ///
                 `gdpmean' `"Avg."' ///
                 `gdpplusd' "+1 SD" `gdpplu2sd' "+2 SD", labsize(*.8)) ///
          title("Excluding `naam'") ///
          legend(off) ///
          name(figurewo`n', replace) xsize(3) ysize(2.5) nodraw
  restore
}

graph combine figurewo1  figurewo2 ///
              figurewo3  figurewo4 ///
              figurewo5  figurewo6 ///
              figurewo7  figurewo8 ///
              figurewo9  figurewo10 ///
              figurewo11 figurewo12 ///
              figurewo13 figurewo14 ///
              figurewo15 figurewo16 , col(4) row(4) xsize(12) ysize(10)
// And so on ...

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