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 (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 (27): Between and within regression

// Case A
clear
set seed 1
set obs 300

generate  e = 0 + (200 - 0) * runiform()  // To generate random variates over the
generate  x = 0 + (600 - 0) * runiform()  // interval [a,b), a+(b-a)*runiform()

generate  y = e // Relationship between X and Y random
generate id = 1
replace   y = (e + 200) if x > 200 & x <= 400
replace  id = 2         if x > 200 & x <= 400
replace   y = (e + 400) if x > 400
replace  id = 3         if x > 400

bys id: egen ym = mean(y)  // Generate cluster means for Y and X
bys id: egen xm = mean(x)

twoway (scatter y x if id == 1, msymbol(oh)) ///
       (scatter y x if id == 2, msymbol(dh)) ///
       (scatter y x if id == 3, msymbol(sh)) ///
       (scatter ym xm,          msymbol(S)) ///
       (lfit ym xm,             lpattern(solid)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
      , /*legend(order(1 "A" 2 "B" 3 "C") title(Region))*/  legend(off) ///
        xtitle("Component 1") ytitle("Component 2") ///
        title("Case A:" "Regional correlation," "individual orthogonality") ///
        xlabel(none) ylabel(none) name(casea, replace)

// Case B
clear
set seed 1

set obs 4 // Generate four clusters

generate id = _n
generate u_i = 300 if id == 1  // Generate cluster-specific intercepts
replace u_i = 0 if id == 2
replace u_i = 0 if id == 3
replace u_i = -300 if id == 4

expand 100  // Generate units in clusters

bysort id: generate x = -100 + (600 + 100) * runiform()
generate e_ij = rnormal(0,70)                // Generate level-1 error term

generate y = x + u_i + e_ij
drop if y < 0 | x < 0
drop if id == 2 & x >= 300
drop if id == 3 & x  < 300

bys id: egen ym = mean(y)
bys id: egen xm = mean(x)

twoway (scatter ym xm, msymbol(S)) ///
       (scatter y  x  if id == 1, msymbol(oh)) ///
       (scatter y  x  if id == 2, msymbol(dh)) ///
       (scatter y  x  if id == 3, msymbol(sh)) ///
       (scatter y  x  if id == 4, msymbol(th)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
       (lfit y x if id == 4,    lpattern(dash)) ///
       (lfit ym xm, lpattern(solid)) ///
      , legend(order(2 "A" 3 "B" 4 "C" 5 "D" 1 "Means") size(small) title(Region:, size(small)) row(1)) ///
        xtitle("Component 1") ytitle("Component 2") ///
        title("Case B:" "Regional orthogonality," "individual correlation") ///
        xlabel(none) ylabel(none) name(caseb, replace)


// Case C
clear
set seed 1
set obs 3

generate  id = _n
generate u_i = 300 if id == 1
replace u_i = 0 if id == 2
replace u_i = -300 if id == 3

expand 150

bysort id: generate x = -100 + (600 + 100) * runiform()
generate e_ij = rnormal(0,70)                // Generate level-1 error term

generate y = x + u_i + e_ij
drop if y < 0 | x < 0

bys id: egen ym = mean(y)
bys id: egen xm = mean(x)

twoway (scatter y  x  if id == 1, msymbol(oh)) ///
       (scatter y  x  if id == 2, msymbol(dh)) ///
       (scatter y  x  if id == 3, msymbol(sh)) ///
       (scatter ym xm, msymbol(S)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
       (lfit ym xm, lpattern(solid)) ///
       , /*legend(order(1 "A" 2 "B" 3 "C") title(Region))*/ legend(off) ///
       xtitle("Component 1") ytitle("Component 2") ///
       title("Case C:" "Negative regional correlation," "positive individual correlation") ///
       xlabel(none) ylabel(none) ///
       name(casec, replace)

grc1leg casea caseb casec, legendfrom(caseb) span pos(6) row(1) name(combined, replace)

Random graphs (28): Plotting categorical variables

foreach x of varlist V39 V40 V41 {

    // Recode each variable:
  recode `x' (1/5 = 2 "Valid answer") ///
             ( .c = 1 "Can't choose") ///
             ( .n = 0 "No answer") ///
    , gen(`x'_mis)
  label var `x'_mis "Recode of `x'"
  
    // Plot each variable:
  catplot `x'_mis, over(C_ALPHAN, label(ang(v)))
                   stack asyvars perc(C_ALPHAN) ///
                   legend(pos(1) row(1)) recast(bar) ///
     ytitle("Percentage") b1title("") ///
     title("Recode of `x'") ///
     name(`x'_mis, replace)
}

grc1leg V39_mis V40_mis V41_mis, col(1) name(combined, replace)
graph display combined, ysize(11) xsize(7)

Random graphs (20): Bar graph

clear

// Table from Vaupel and Loichinger 2006, 
// "Redistributing Work in Aging Europe,"
// Science 312(5782):1911-1913.
// doi: 10.1126/science.1127487
// Table 1

input str2 country r2005 r2025 rchange h2005 h2025 hchange
DE  1.27  1.47  16  16.28  14.95  -8
DK  0.97  1.12  15  17.46  16.11  -8
FR  1.43  1.69  18  15.09  13.63  -10
IT  1.59  1.86  17  15.19  13.48  -11
NL  1.01  1.20  19  15.31  13.88  -9
UK  1.09  1.19  9  17.32  16.34  -6
US  1.09  0.99  -9  18.71  18.29  -2
end

drop rchange hchange // drop unnecessary columns

graph bar r2005 r2025, over(country) bargap(-30) ///
      ytitle("Ratio nonworkers per worker") ///
      legend(label(1 "2005") label(2 "2025") row(1)) ///
      ylabel(,format(%6.1f)) ///  Format axis labels
      name(r, replace)

graph bar h2005 h2025, over(country) bargap(-30) ///
      ytitle("Hours worked per week per capita") ///
      legend(label(1 "2005") label(2 "2025")) ///
      ylabel(,format(%6.1f)) ///
      name(h, replace)
   
grc1leg r h, xcommon row(1) ///
             note("{it: Source:} Vaupel and Loichinger 2006, p. 1912" ///
                  "{it: Note:} The values for 2025 assume change in the" ///
                  "population pyramid but no change in labor force " ///
                  "participation or effort" "by age and sex.", span)

Random graphs (22): Graphing functions

clear

twoway (function y = 2-.01*x, range(0 100))  ///
       (function y = 1-.01*x, range(0 100)), ///
    text(2 20 "High SES") ///
    text(1 20 "Low SES") ///
    xtitle(Age) ytitle(Health) ylabel(0(1)2.5, nolabel noticks) /// 
    /* ylabel(none) doesn't allow to control axis, thus
       this is a helpful workaround */ ///
    xlabel(0 100, nolabel noticks) ///
    legend(off) ///
    /*legend(label(1 "High SES") label(2 "Low SES") pos(6) col(2))*/ ///
    title("Status maintenance", box bexpand) ///
    name(maint, replace) 
twoway (function y = 2-.004*x, range(0 100))  ///
       (function y = 1-.01*x, range(0 100)), ///
    text(2 20 "High SES") ///
    text(1 20 "Low SES") ///
    xtitle(Age) ytitle(Health) ylabel(0(1)2.5, nolabel  noticks) ///
    xlabel(0 100, nolabel noticks) ///
    legend(off) ///
    /*legend(label(1 "High SES") label(2 "Low SES")) */ ///
    title("Cumulative (dis-)advantage", box bexpand) ///
    name(cumulat, replace) 
twoway (function y = 2-.015*x, range(0 100))  ///
       (function y = 1-.01*x, range(0 100)), ///
    text(2 20 "High SES") ///
    text(1 20 "Low SES") ///
    xtitle(Age) ytitle(Health) ylabel(0(1)2.5, nolabel noticks) ///
    xlabel(0 100, nolabel noticks) ///
    legend(off) ///
    /*legend(label(1 "High SES") label(2 "Low SES")) */ ///
    title("Age as leveler", box bexpand) ///
    name(leveler, replace) 
graph combine cumulat maint leveler, xcommon col(1) xsize(3) ysize(8)
graph export Graph.png, replace

// Same graph with one legend at the bottom:
twoway (function y = 2-.01*x, range(0 100))  ///
       (function y = 1-.01*x, range(0 100)), ///
    xtitle(Age) ytitle(Health) ylabel(0(1)2.5, nolabel noticks) /// 
    /* ylabel(none) doesn't allow to control axis, thus
       this is a helpful workaround */ ///
    xlabel(0 100, nolabel noticks) ///
    legend(label(1 "High SES") label(2 "Low SES") pos(6) col(2)) ///
    title("Status maintenance", box bexpand) ///
    name(maint, replace) 
twoway (function y = 2-.004*x, range(0 100))  ///
       (function y = 1-.01*x, range(0 100)), ///
    xtitle(Age) ytitle(Health) ylabel(0(1)2.5, nolabel  noticks) ///
    xlabel(0 100, nolabel noticks) ///
    legend(off) ///
    title("Cumulative (dis-)advantage", box bexpand) ///
    name(cumulat, replace) 
twoway (function y = 2-.015*x, range(0 100))  ///
       (function y = 1-.01*x, range(0 100)), ///
    xtitle(Age) ytitle(Health) ylabel(0(1)2.5, nolabel noticks) ///
    xlabel(0 100, nolabel noticks) ///
    legend(off) ///
    title("Age as leveler", box bexpand) ///
    name(leveler, replace) 
grc1leg cumulat maint leveler, xcommon col(1) ///
xsize(3) ysize(8) name(combined, replace) legendfrom(maint)
// Problem: -grc1leg- appears to ignore -xsize()- and -ysize()-
graph display combined, xsize(3) ysize(8)
// Solution: Redraw so that size commands take effect
graph export Graphwlegend.png, replace