Specification curve analysis

use "ZA4612_v1-0-1.dta", clear
do "ZA4612_patch_v1-0-1.do" // Some patch from data provider

keep v399-v404 v933 v827 v301 v298 

// Outcomes
mvdecode v399-v404, mv(9 = .a )

generate stress  = 5 - v399
generate depress = 5 - v400
generate calm    = v401 - 1
generate energy  = v402 - 1
generate pain    = 5 - v403
generate lonely  = 5 - v404

// Predictors
  // Topbot
mvdecode v933 v827, mv(96 = .a \ 99 = .b) 
generate topbot = v933
replace  topbot = v827 if missing(topbot)
  // Age
mvdecode v301, mv(999 = .a)
rename v301 age
  // Female
generate female = (v298 == 2)

drop v399-v404 v933 v827 v298 // Drop unnecessary variables

// Post definition
tempname foo
postfile `foo' str50 spec k1 k2 k3 b se using deleteme.dta, replace

// Loop
forvalues k1 = 1/6 {
  if `k1' == 1 local y "stress"
  if `k1' == 2 local y "depress"
  if `k1' == 3 local y "pain"
  if `k1' == 4 local y "calm"
  if `k1' == 5 local y "energy"
  if `k1' == 6 local y "lonely"
    forvalues k2 = 1/4 {
      if `k2' == 1 local agecontrol "age"
      if `k2' == 2 local agecontrol "c.age##c.age"
      if `k2' == 3 local agecontrol "c.age##c.age##c.age"    
      if `k2' == 4 local agecontrol "c.age##c.age##c.age##c.age"
        forvalues k3 = 1/3 {
          if `k3' == 1 local ifs " "
          if `k3' == 2 local ifs "if female == 1"
          if `k3' == 3 local ifs "if female == 0"
 
        local spec regress `y' topbot `agecontrol' female `ifs'

        qui `spec'
        local  b =  _b[topbot]
        local se = _se[topbot]

        post `foo' ("`spec'") (`k1') (`k2') (`k3') (`b') (`se')
        }
    }
}
postclose `foo'

// Plot specification curve
use deleteme, clear

// Generate ranked analytical choice variable
sort b
generate sk = _n
label var sk "Specification (sorted by coefficient size)"

// Calculate CI's
generate ub = b + 1.96 * se
generate lb = b - 1.96 * se

// Remind yourself what the variables mean
label var k1 "Outcome"
label var k2 "Age control"
label var k3 "Subsample"

// Stack indicators
generate k3c = k3
generate k2c = k2 + 1 + 3 // 3 because K3 has 3 categories, 1 for title
generate k1c = k1 + 1 + 3 + 4 + 1  // K2 has 4 categories

// Calculate some things for size of second axis
qui summarize b
global brange = r(max) - r(min)
global bmin = r(min)
global bmax = r(max)
global from_y = $bmin - (4.5 * $brange)

// Plot
twoway (scatter k1c k2c k3c sk, msymbol(o o o) msize(vsmall vsmall vsmall) ///
        yscale(range(1 22)) ///        
        ylabel( 1 "Full sample" 2 "Females only" 3 "Males only" 4 "{bf:Subsample}" ///
                5 "Age" 6 "Age squared" 7 "Age cubed" 8 "Age quartic" 9 "{bf:Age control}" ///
               10 "Stress"   11 "Depression" 12 "Pain" ///
               13 "Calmness" 14 "Energy"     15 "Loneliness" 16 "{bf:Outcome}", tstyle(notick) axis(1)))  ///
       (rcap b b sk, yaxis(2) yscale(range($from_y $bmax) axis(2)) ///
                     ylab(, format(%2.1g) axis(2)) ///
                     ytitle("{bf:Coefficient size}", axis(2) placement(north)) ///
                     yline(0, axis(2))) ///
       (rspike ub lb sk, yaxis(2)), ///
        xlabel(none)  ///
        legend(off) name(curve, replace) 

Reference

Simonsohn, Uri, Joseph P. Simmons, and Leif D. Nelson. 2015. Specification Curve. Descriptive and Inferential Statistics on All Reasonable Specifications. University of Pennsylvania. doi: 10.2139/ssrn.2694998

Random graphs (142): Boxplot using twoway commands with overlaid data

use v358 v392 v585 v829 if inrange(v829, 1, 5) using "data allbus\ZA4612_v1-0-1.dta", clear

// Subjective social mobility
recode v829 (5 = 1 "Much lower status") ///
            (4 = 2 "Lower status") ///
            (3 = 3 "About equal") ///
            (2 = 4 "Higher status") ///
            (1 = 5 "Much higher status") ///
            (95 96 98 99 = .), gen(ssm)
label var ssm "Status of own job compared to father's"

// Respondent's social status ISEI
clonevar isei = v358
replace  isei =    . if inlist(isei, 0, 99)
replace  isei = v392 if missing(isei)
replace  isei =    . if inlist(isei, 0, 99)
label var isei "Respondent's ISEI"

// Father's social status ISEI
clonevar fisei = v585 
replace  fisei =    . if inlist(fisei, 0, 99)
label var fisei "Father's ISEI"

// Objective status mobility
gen osm = isei - fisei
label var osm "Objective status mobility"



// Boxplot using -graph- command
graph box osm, over(ssm, descending) horizontal ///
               yline(0) l1title("Status of own job compared to father's" "(self-assessment)") ///
               name(figure1, replace) 

// Boxplot using -twoway- command

// Calculate stuff for boxplot
qui bys ssm: sum osm, detail
qui bys ssm: egen min = pctile(osm), p(25) // Interquartile range
qui bys ssm: egen max = pctile(osm), p(75) // Interquartile range
qui bys ssm: egen med = median(osm)

// Plot
twoway (scatter ssm osm, msymb(o) jitter(2) mcolor(red%10)) ///
       (scatter ssm med, msymb(oh) mcolor(black)) ///
       (rspike  min max ssm, horizontal lcolor(black)), ///
        ylabel(1/5, valuelabel) xline(0) ///
        ytitle("Status of own job compared to father's" "(self-assessment)") ///
        legend(order(2 "Median" 3 "Interquartile range") ring(0) pos(11)) ///
        xtitle(Objective social mobility) name(figure2, replace)
drop min max med // Drop stuff for boxplot

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 (96): Categorical variables

use S003A S002 X001 C041 C039 C038 C037 if inlist(S002, 4, 5) ///
    using WVS_Longitudinal_1981-2014_stata_v_2014_11_25.dta, clear

// Declare missing values
mvdecode C041 C039 C038 C037, mv(-5 -4 -3 -2 -1)

// Factor analysis
factor C041 C039 C038 C037, pcf
alpha C041 C039 C038 C037, item
local alph = round(`r(alpha)', .01)

// Prepare variables
scores nonworkethic = mean(C041 C039 C038 C037), nv(3) 
generate workethic = 5 - nonworkethic 
gen male = (X001 == 1) if X001 != -2
drop nonworkethic X001


// Create descriptive plot    
foreach var of varlist C041 C039 C038 C037 {
 recode `var' (1 = 4 "Strongly agree") ///
              (2 = 3 "Agree") ///
              (3 = 2 "Neither agree nor disagree") ///
              (4 = 1 "Disagree") ///
              (5 = 0 "Strongly disagree"), gen(`var'_rec)
 numlabel `var'_rec, add mask("# ")
 twoway histogram `var'_rec, discrete horizontal yla(0/4, valuelabel) percent ///
                             title("`: variable label `var''") ///
                             ytitle("") ///
                             name(`var', replace) nodraw
 drop `var'_rec
}
graph combine C041 C039 C038 C037, col(2) row(2) altshrink title("Work ethic is measured as the average of four items:") ///
      caption("In a PCA, all items load on one dimension; explained variance = 49%, Cronbach's alpha = `alph'", span) note("{it:Source:} WVS 1999-2009, pooled", span) name(figure1, replace)
drop C041 C039 C038 C037


// Calculate country differences
preserve   
  statsby mean_ = _b[_cons] ///
          loci  = (_b[_cons] - 1.96 * _se[_cons]) ///
          hici  = (_b[_cons] + 1.96 * _se[_cons]) ///
        , by(S003A) clear: ///
          regress workethic

  // Sort coefficients by size
  egen order_ = rank(mean_), unique
  labmask order_, value(S003A) decode

  // Plot
  twoway (rcap mean_ mean_ order_, horizontal) ///
         (rspike loci hici order_, horizontal) ///
        , legend(off) ylabel(1/63, valuelabels ang(h) labsize(*.8)) ///
          xlabel(0 (1) 4, grid format(%6.1f)) name(all, replace)  ///
          xmtick(0 (.5) 4) ///
          ytitle("") xtitle("Work ethic across countries") ///
          title("{bf:A}", justification(left) bexpand span) ///
          xscale(alt) ysize(10) nodraw
restore  

// Country differences--men only
preserve   
  statsby mean_ = _b[_cons] ///
          loci  = (_b[_cons] - 1.96 * _se[_cons]) ///
          hici  = (_b[_cons] + 1.96 * _se[_cons]) ///
        , by(S003A) clear: ///
          regress workethic if male == 1

  // Sort coefficients by size
  egen order_ = rank(mean_), unique
  labmask order_, value(S003A) decode

  // Plot
  twoway (rcap mean_ mean_ order_, horizontal) ///
         (rspike loci hici order_, horizontal) ///
        , legend(off) ylabel(1/63, valuelabels ang(h) labsize(*.8)) ///
          xlabel(0 (1) 4, grid format(%6.1f)) name(men, replace)  ///
          xmtick(0 (.5) 4) ///
          ytitle("") xtitle("Men's work ethic across countries") ///
          title("{bf:B}", justification(left) bexpand span) ///
          xscale(alt) ysize(10) nodraw
restore  

// Calculate gender gap
preserve   
  statsby mean_ = _b[male] ///
          loci  = (_b[male] - 1.96 * _se[male]) ///
          hici  = (_b[male] + 1.96 * _se[male]) ///
        , by(S003A) clear: ///
          regress workethic male

  // Sort coefficients by size
  egen order_ = rank(mean_), unique
  labmask order_, value(S003A) decode

  // Plot
  twoway (rcap mean_ mean_ order_, horizontal) ///
         (rspike loci hici order_, horizontal) ///
        , legend(off) ylabel(1/63, valuelabels ang(h) labsize(*.8)) ///
          xlabel(-.25 (.25) .5, grid format(%6.2f)) name(gendergap, replace)  ///
          xmtick(-.25 (.1) .5) ///
          text(63 .5  "Men" "higher", place(sw)) ///
          text( 1 -.25 "Women" "higher", place(ne)) ///    
          ytitle("") xtitle("Gender gap in work ethic") ///
          title("{bf:C}", justification(left) bexpand span) ///
          xscale(alt) ysize(10) nodraw
restore  

graph combine all men gendergap, note(" " "{it:Source:} WVS 1999-2009, pooled", span) col(3) ysize(12) xsize(18) altshrink ///
                                 title("Work ethic in cross-national comparison", span) name(figure2, 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 

Generating a country–year variable

use "C:\ess 1-7\ESS1-7e01.dta", clear

// Generate country-year variable
egen cyear = group(cntry essround)

// Fit three-level model
mixed happy || cntry: || cyear:, variance
estat icc

capture drop u1* u0*
predict u1 u0, reffects 
predict u1se u0se, reses 


// Plot country-level variation
preserve
egen pickone = tag(cntry)
keep if pickone
   
egen order_ = rank(-u1), unique
labmask order_, value(cntry)

gen high = u1 + (1.96 * u1se)
gen low  = u1 - (1.96 * u1se)

twoway (rcap u1 u1 order_, dsymbol(x)) ///
       (rspike high low order_) , ///
        yline(0) xlabel(1/32, val ang(v)) ///
        xtitle("") ///
        ytitle("Country-level residuals") ///
        legend(off) ///
        name(cntry, replace) 
restore 


// Plot country-year level variation
preserve
egen pickone = tag(cyear)
keep if pickone
   
egen order_ = rank(-u0), unique
labmask order_, value(cyear)

gen high = u0 + (1.96 * u0se)
gen low  = u0 - (1.96 * u0se)

twoway (rcap u0 u0 order_, dsymbol(x) horizontal) ///
       (rspike high low order_, horizontal) , ///
        xline(0) ylabel(none) ///
        xscale(alt) ///
        ytitle("") ///
        xtitle("Country{c 150}year-level residuals") ///
        legend(off) ///
        name(cyear, replace) ysize(8)
restore 

Random graphs (91): Regression coefficients

// ESS round 1
use ESS1e06_4.dta, clear

// Prepare variables
  // Country variable
  encode cntry, gen(country)

  // Trade union membership
  recode trummb (0 = 0 "Non-member") (1 = 1 "Trade union member") (. = .), generate(unionmember)
  label var unionmember "Trade union member"

  // Age
  generate age = agea if agea != 999
  label var age "Age"

  // Sex
  gen female = (gndr == 2) if gndr != 9
  label var female "Female sex"
  label define female 1 "Female" 0 "Male"
  label val female female

  // Happiness
  recode happy (77 88 99 = .)

  // Center age and sex
  center age female
  
// First set of models
  preserve   
  statsby mean_ = _b[unionmember] ///
          loci  = (_b[unionmember] - 1.96 * _se[unionmember]) ///
          hici  = (_b[unionmember] + 1.96 * _se[unionmember]) ///
        , by(country) clear: ///
          regress happy unionmember

  // Sort coefficients by size
  egen order_ = rank(mean_), unique
  labmask order_, value(country) decode

  // Plot
  twoway (rcap mean_ mean_ order_, horizontal) ///
         (rspike loci hici order_, horizontal) ///
        , legend(off) ylabel(1/20, valuelabels ang(h) labsize(*.8)) ///
          xlabel(-.5 (.5) 1.5, grid format(%6.1f)) name(unadjusted, replace)  ///
          xmtick(-.5 (.1) 1.5) ///
          ytitle("") xtitle("Happiness advantage of trade union membership" " ") ///
          xscale(alt) ysize(4)
restore  

// Second set of models
  preserve   
  statsby mean_ = _b[unionmember] ///
          loci  = (_b[unionmember] - 1.96 * _se[unionmember]) ///
          hici  = (_b[unionmember] + 1.96 * _se[unionmember]) ///
        , by(country) clear: ///
          regress happy unionmember c_female c_age

  // Sort coefficients by size
  egen order_ = rank(mean_), unique
  labmask order_, value(country) decode

  // Plot
  twoway (rcap mean_ mean_ order_, horizontal) ///
         (rspike loci hici order_, horizontal) ///
        , legend(off) ylabel(1/20, valuelabels ang(h) labsize(*.8)) ///
          xlabel(-.5 (.5) 1.5, grid format(%6.1f)) name(adjusted, replace)  ///
          xmtick(-.5 (.1) 1.5) ///
          ytitle("") xtitle("Happiness advantage of trade union membership" "(adjusted for age and sex)") ///
          xscale(alt) ysize(4)
  restore  
  
// Plot both underneath one another
graph combine unadjusted adjusted, row(2) ysize(8) xcommon name(combined, replace) ///
          note(" " ///
             "{it:Source:} ESS round 1, own calculations. {it:Note:} Error bars denote 95% confidence intervals." , span size(*.8))  

Merging the Demographic and Health Surveys in Stata

// Unzip files
clear
cd "C:\dhs"

local filelist : dir . files "*.zip"  // Create local with all filenames ending with ".zip"
di `filelist'

local first : word 1 of `filelist'      // Identify first file
di "`first'"

local total_ : word count `filelist' // Identify total number of files
di `total_'        

forvalues x = 1/`total_' {
    di `x'
    local y : word `x' of `filelist'
    unzipfile "`y'", replace
}

// Append all files
clear
cd "C:\dhs"
local directorylist : dir . dirs  "*ir*"  // Create local with all directory names that contain women's data

di `directorylist'

local firstdir : word 1 of `directorylist'      // Identify first directory
di "`firstdir'"

local total_ : word count `directorylist' // Identify total number of directories
di `total_'

local firstfile = strupper("`firstdir'")
local firstfile = subinstr("`firstfile'","DT","FL",.)
di "`firstfile'"

use caseid v000 v005 v007 v012 v106 v107 v155 v191 v201 v212 v525 v531 using "`firstdir'/`firstfile'", clear                    
capture decode v106, gen(v106s)
drop v106

save testfile, replace

forvalues x = 2/`total_' {
    local y : word `x' of `directorylist'
 local filename = strupper("`y'")
    local filename = subinstr("`filename'", "DT", "FL", .)
 use "`y'/`filename'", clear
 
 // Source: https://stackoverflow.com/questions/17056016/stata-how-to-keep-a-list-of-variables-given-some-of-them-may-not-exist
 local masterlist "caseid v000 v005 v007 v012 v106 v107 v155 v191 v201 v212 v525 v531"
    local keeplist = ""

    foreach i of local masterlist  {
    capture confirm variable `i'
        if !_rc {
            local keeplist "`keeplist' `i'"
        }
     }
    keep `keeplist'
 capture decode v106, gen(v106s)
    capture drop v106
 
 tempfile new
 save `new', replace
 use testfile, clear
 append using `new', force
 save testfile, replace
}

// Prepare variables
use testfile, clear

//Country and wave identifiers
replace v000 = "VN3" if v000 == "VNT"
generate cntry = substr(v000,1,2)
generate wavex  = substr(v000,3,1)
replace wavex = "1" if wavex == ""
encode wavex, gen(wave)
drop wavex

  // Fix unusual country abbreviations
replace cntry = "BI" if cntry == "BU"
replace cntry = "IN" if cntry == "IA"
replace cntry = "KZ" if cntry == "KK"
replace cntry = "BI" if cntry == "BU"
replace cntry = "MD" if cntry == "MB"
replace cntry = "NA" if cntry == "NM"
replace cntry = "DO" if cntry == "DR"

kountry cntry, from(iso2c)
encode NAMES_STD, gen(country)
drop NAMES_STD

// Select last wave
keep if wave == 6

// Prepare variables
  // Age at first intercourse
generate age1stintercourse = .
replace  age1stintercourse = v531 if inrange(v531, 1, 63)
replace  age1stintercourse = .a   if v525 == 0
replace  age1stintercourse = .b   if v525 == 95
replace  age1stintercourse = .c   if v525 == 97
replace  age1stintercourse = .d   if v525 == 98
replace  age1stintercourse = .e   if v525 == 99
label define age1stintercourse .a "Not had intercourse" ///
                               .b "95?" ///
                               .c "inconsistent" ///
                               .d "don't know" ///
                               .e "99?" 
label val age1stintercourse age1stintercourse
label var age1stintercourse "Age at first intercourse"

  // Age at first birth
rename v212 afb
label var afb "Age at first birth"

// Calculate correlation
preserve
statsby mean_ = _b[age1stintercourse] ///
        loci  = (_b[age1stintercourse] - 1.96 * _se[age1stintercourse]) ///
        hici  = (_b[age1stintercourse] + 1.96 * _se[age1stintercourse]) ///
      , by(country) total clear: ///
        regress afb age1stintercourse

// Label total value
replace country = 1000 if country == .
label define country 1000 "{bf: Total}", modify

egen order_ = rank(mean_), unique
labmask order_, value(country) decode

// Plot
twoway (rcap mean_ mean_ order_, horizontal) ///
       (rspike loci hici order_, horizontal) ///
      , legend(off) ylabel(1/38, valuelabels ang(h) labsize(*.8)) ///
        xlabel(.7 (.1) 1.0, grid format(%6.1f)) name(ols, replace)  ///
        xmtick(.7 (.05) 1.0) ///
        ytitle("") xtitle("Association between" ///
                          "age at first intercourse" ///
                          "and age at first birth") xscale(alt) ysize(8) ///
        note(" " ///
             "{it:Source:} DHS VI, own calculations" , span size(*.8))
restore

Random graphs (79): Means with confidence intervals

use ZA5900_v3-0-0.dta, replace

renvars, lower // Switch variable names to lower case

// Fix country variable
generate cntry = c_alphan
replace  cntry = "GB" if cntry == "GB-GBN"
replace  cntry = "DE" if cntry == "DE-E"   | cntry == "DE-W"
replace  cntry = "BE" if cntry == "BE-BRU" | cntry == "BE-WAL" | cntry == "BE-FLA" 

// Generate country name variable
kountry cntry, from(iso2c)
encode NAMES_STD, gen(country)
drop NAMES_STD

// Happiness variable
recode v55 (1 = 6) (2 = 5) (3 = 4) (4 = 3) (5 = 2) (6 = 1) (7 = 0) (0 8 9 = .), gen(lsat)

label define lsat 0 "Completely unhappy" ///
                  1 "Very unhappy" ///
                  2 "Fairly unhappy" ///
                  3 "Neither happy nor unhappy" ///
                  4 "Fairly happy" ///
                  5 "Very happy" ///
                  6 "Completely happy"
label val lsat lsat
label var lsat "Happiness"

// Plot means by country
preserve
statsby mean_ = _b[_cons] ///
        loci  = (_b[_cons] - 1.96 * _se[_cons]) ///
        hici  = (_b[_cons] + 1.96 * _se[_cons]) ///
      , by(country) total clear: ///
        regress lsat

replace country = 1000 if country == .
label define country 1000 "{bf: Total}", modify

egen order_ = rank(mean_), unique
labmask order_, value(country) decode

twoway (rcap mean_ mean_ order_, horizontal) ///
       (rspike loci hici order_, horizontal) ///
      , legend(off) ylabel(1/41, valuelabels ang(h) labsize(*.8)) ///
        xlabel(3.5 (.5) 5.0, grid format(%6.1f)) name(lsat, replace)  ///
        xmtick(3.5 (.25) 5.0, grid) ///
        ytitle("") xtitle("Average happiness") xscale(alt) ysize(8) ///
        note(" " ///
             "{it:Note:} Happiness ranges from 0 ('Completely unhappy') to 6 ('Completely happy')" ///
             "{it:Source:} ISSP 2012, doi:10.4232/1.12339" , span size(*.8))
restore

Random graphs (72): Dot plots with the by() option

clear
input str20 cntry bindirect blbindirect ulbindirect bdirect blbdirect ulbdirect
Austria 0.002 -0.002 0.005 0.008 -0.016 0.032
Germany 0.006 0.000 0.012 0.040 0.017 0.063
Sweden 0.001 -0.005 0.007 0.054 0.031 0.076
Netherlands 0.000 -0.006 0.006 0.039 0.017 0.061
Spain 0.005 0.002 0.009 0.016 -0.001 0.033
Italy 0.000 -0.005 0.005 0.029 0.011 0.047
France 0.001 -0.003 0.006 0.043 0.025 0.061
Denmark 0.001 0.000 0.003 0.015 0.001 0.029
Greece 0.000 -0.003 0.003 0.041 0.026 0.057
Switzerland 0.002 -0.002 0.006 0.003 -0.020 0.026
Belgium 0.000 -0.003 0.004 0.033 0.016 0.050
"Czech Republic" 0.004 -0.001 0.010 0.073 0.031 0.115
Poland -0.004 -0.010 0.002 0.034 -0.003 0.071
end

// Reshape to long format
reshape long b blb ulb, j(fx "indirect" "direct") i(cntry) string

// Encode variables
encode cntry, gen(country)
encode fx, gen(effect)
label define effect 1 "Direct effect" 2 "Indirect effect", modify

// Sort by biggest direct effect size
egen order_ = rank(b) if effect == 1, unique
labmask order_, val(country) decode
bysort country (order_): replace order_ = order_[1] // Copy value to all cases 

// xrescale makes sure that x-axis fits both by() plots
twoway (dot b order_, horizontal by(effect, legend(off) note(" ") xrescale)) ///
       (rspike blb ulb order_, horizontal by(effect)) ///
      , ylabel(1 (1) 13, val) xscale(alt) ///
        ytitle(" ") xline(0) 
    

Random graphs (71): Proportions with confidence intervals

use sharew5_rel1-0-0_it.dta, clear

// Prepare variables
recode it004_ (1 = 1 "Yes") (5 = 0 "No") (-1 -2 . = .), generate(internetuser)
label var internetuser "Internet use in last 7 days"
decode country, gen(cntry)

// Define stuff
tempname foo
postfile `foo' str25 cntry perc perc_ll perc_ul using "C:\Windows\Temp\test.dta", replace
levelsof cntry, local(levels)

// Calculate proportions by country
foreach country of local levels {

    capture proportion internetuser if cntry == "`country'"
    matrix fcoefs  = r(table)
    local fperc    = fcoefs[1,2] * 100
    local fperc_ll = fcoefs[5,2] * 100
    local fperc_ul = fcoefs[6,2] * 100

    di "`country'"  _skip(2)  `fperc_ll' _skip(2) `fperc' _skip(2) `fperc_ul'

    if "`fperc'" != "" {    // Make sure that the loop doesn't break if empty
    post `foo' ("`country'") (`fperc') (`fperc_ll') (`fperc_ul')
    }
}
postclose `foo'

// Use posted data set
use "C:\Windows\Temp\test.dta", clear

// Sort countries by size
egen order_ = rank(perc), unique
labmask order_, value(cntry) 

// Plot
twoway (dot perc order_, horizontal) ///
       (rspike perc_ll perc_ul order_, horizontal), ///
       ylabel(1/15, val) legend(off) xscale(range(20 80) alt) ///
       ytitle("") xtitle("% Internet users (in last 7 days)") ///
       xlabel(20(10)80, grid) ///
       note(" " "{it:Source:} SHARE wave 5, doi:10.6103/SHARE.w5.100", span)

Random graphs (69): ICC's with confidence intervals

use sharew1_rel2-6-0_gv_isced.dta, clear

// Reshape data into long format
reshape long iscedy_c, i(mergeid) j(child)
// Fix variable of interest
recode iscedy_c (-7 = .a "not yet coded (temporary)") ///
                (-2 = .b "refusal") ///
                (-1 = .c "don't know") ///
                (95 = .d "still in school") ///
                (97 = .e "other") ///
                ( . = .f "missing") ///
               , gen(years)
label var years "Years of education"
    
// Israel doesn't provide a ISCED-to-years conversion, thus it's dropped here
drop if country == 25

// Fix country variable
decode country, gen(cntry)

// Calculate ICC's per country
tempname foo
postfile `foo' str100 commandline_str str20 cntry icc icclb iccub N N_groups using "C:\Windows\Temp\test.dta", replace    
    
levelsof cntry, local(country)

foreach x of local country {
  *di "`x'"
  qui mixed years || mergeid: if cntry == "`x'", reml
  qui estat icc
  
  matrix groups = e(N_g)
  scalar n_g = groups[1,1]
  matrix cis =  r(ci2)
  scalar icclb = cis[1,1]
  scalar iccub = cis[1,2]
  
  post `foo' (e(cmdline)) ("`x'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)
  
  matrix stuff1 = (r(icc2), icclb, iccub, e(N), n_g)
  if "`x'" == "Austria" matrix table1 = stuff1
  else matrix table1 = (table1\stuff1)
}

postclose `foo'

matrix rownames table1 = `country'
matrix colnames table1 = "ICC" "CI lower" "CI upper" "Children" "Families"
esttab matrix(table1, fmt(2 2 2 0 0))

use "C:\Windows\Temp\test.dta", clear

egen order_ = rank(-icc), unique
labmask order_, value(cntry)

twoway (rcap   icc icc order_, horizontal) ///
       (rspike icclb iccub order_, horizontal) ///
      , ylabel(1/11, val) legend(off) ytitle("") xscale(alt) ///
        xtitle("Sibling correlations in educational attainment") ///
        note(" " "{it:Source:} SHARE wave 1, doi:10.6103/SHARE.w1.260", span)

erase "C:\Windows\Temp\test.dta"

Random graphs (67): Dot plots with confidence intervals

clear
input str20 cntry str30 outcome lb icc ub
"Finland"   "Attainment" .35 .36 .38
"Germany"   "Attainment" .43 .51 .60
"Norway"   "Attainment" .41 .41 .42
"Sweden"   "Attainment" .41 .41 .42
"United Kingdom" "Attainment" .34 .44 .54
"United States"  "Attainment" .49 .51 .53
"Finland"   "GPA"        .   .   .
"Germany"   "GPA"        .19 .23 .28
"Norway"   "GPA"        .48 .48 .48
"Sweden"   "GPA"        .52 .52 .52
"United Kingdom" "GPA"        .   .   .
"United States"  "GPA"        .37 .41 .46
"Finland"   "Cognitive skills"  .   .   .
"Germany"   "Cognitive skills"  .39 .46 .52
"Norway"   "Cognitive skills"  .   .   .
"Sweden"   "Cognitive skills"  .50 .50 .51
"United Kingdom" "Cognitive skills"  .   .   .
"United States"  "Cognitive skills"  .55 .57 .60
end

// Defining the label first allows to determine order of -encode- categories
label define country 1 "Finland" 2 "Germany" 3 "Norway" 4 "Sweden" ///
                     5 "United Kingdom" 6 "United States"
encode cntry, gen(country) label(country)

twoway (rcap icc icc country, by(outcome, note("") row(1) legend(off)) horizontal) ///
       (rspike ub lb country, by(outcome) horizontal) ///
      , ylabel(1/6, val)  ytitle("") xscale(alt) yscale(reverse) 

Random graphs (66): Dot plots with confidence intervals

clear
input str50 variable str4 data icc lb ub
"Overall" DE .46 .40 .53
"Overall" SE .49 .49 .50
"Overall" PSID .57 .55 .6
"Low father education" DE .42 .34 .5
"Low father education" SE .44 .43 .45
"Low father education" PSID .42 .31 .53
"High father education" DE .4 .29 .53
"High father education" SE .47 .46 .47
"High father education" PSID .58 .55 .62
"Low mother education" DE .45 .38 .53
"Low mother education" SE .45 .44 .46
"Low mother education" PSID .56 .53 .59
"High mother education" DE .33 .21  .48
"High mother education" SE .45 .44 .46
"High mother education" PSID .56 .53 .59 
"Low parental occupation" DE .51 .43 .60
"Low parental occupation" SE .43 .42 .44
"Low parental occupation" PSID .59 .51 .67
"High parental occupation" DE .30 .21 .41
"High parental occupation" SE .39 .37 .40
"High parental occupation" PSID .54 .51 .57
"No migrant background" DE .40 .33 .49 
"No migrant background" SE .48 .48 .49
"No migrant background" PSID .57 .53 .61
"Migrant background" DE .48 .35 .61
"Migrant background" SE .54 .52 .55
"Migrant background" PSID .56 .44 .68
"Family size 2-3 children" DE .48 .41 .54
"Family size 2-3 children" SE .48 .47 .48
"Family size 2-3 children" PSID . . .
"Family size 4 and more children" DE .31 .12 .61
"Family size 4 and more children" SE .55 .54 .56
"Family size 4 and more children" PSID . . .
end

generate indicator =  1 if variable == "Overall"                  & data == "DE"
replace  indicator =  2 if variable == "Overall"                  & data == "SE"
replace  indicator =  3 if variable == "Overall"                  & data == "PSID"
replace  indicator =  4 if variable == "Low father education"     & data == "DE"
replace  indicator =  5 if variable == "Low father education"     & data == "SE"
replace  indicator =  6 if variable == "Low father education"     & data == "PSID"
replace  indicator =  7 if variable == "High father education"    & data == "DE"
replace  indicator =  8 if variable == "High father education"    & data == "SE"
replace  indicator =  9 if variable == "High father education"    & data == "PSID"
replace  indicator = 10 if variable == "Low mother education"     & data == "DE"
replace  indicator = 11 if variable == "Low mother education"     & data == "SE"
replace  indicator = 12 if variable == "Low mother education"     & data == "PSID"
replace  indicator = 13 if variable == "High mother education"    & data == "DE"
replace  indicator = 14 if variable == "High mother education"    & data == "SE"
replace  indicator = 15 if variable == "High mother education"    & data == "PSID"
replace  indicator = 16 if variable == "Low parental occupation"  & data == "DE"
replace  indicator = 17 if variable == "Low parental occupation"  & data == "SE"
replace  indicator = 18 if variable == "Low parental occupation"  & data == "PSID"
replace  indicator = 19 if variable == "High parental occupation" & data == "DE"
replace  indicator = 20 if variable == "High parental occupation" & data == "SE"
replace  indicator = 21 if variable == "High parental occupation" & data == "PSID"
replace  indicator = 22 if variable == "No migrant background" & data == "DE"
replace  indicator = 23 if variable == "No migrant background" & data == "SE"
replace  indicator = 24 if variable == "No migrant background" & data == "PSID"
replace  indicator = 25 if variable == "Migrant background" & data == "DE"
replace  indicator = 26 if variable == "Migrant background" & data == "SE"
replace  indicator = 27 if variable == "Migrant background" & data == "PSID"
replace  indicator = 28 if variable == "Family size 2-3 children" & data == "DE"
replace  indicator = 29 if variable == "Family size 2-3 children" & data == "SE"
replace  indicator = 30 if variable == "Family size 2-3 children" & data == "PSID"
replace  indicator = 31 if variable == "Family size 4 and more children" & data == "DE"
replace  indicator = 32 if variable == "Family size 4 and more children" & data == "SE"
replace  indicator = 33 if variable == "Family size 4 and more children" & data == "PSID"

label define indicator 1 "Overall" ///
                       4 "Low father education" ///
                       7 "High father education" ///
                      10 "Low mother education" ///
                      13 "High mother education" ///
                      16 "Low parental occupation" ///
                      19 "High parental occupation" ///
                      22 "No migrant background" ///
                      25 "Migrant background" ///
                      28 "Family size 2-3 children" ///
                      31 "Family size 4 and more children"
label val indicator indicator
  
twoway (dot icc indicator if data == "DE", horizontal msymbol(o)) ///
       (rspike ub lb indicator if data == "DE", horizontal) ///
       (dot icc indicator if data == "SE", horizontal msymbol(t)) ///
       (rspike ub lb indicator if data == "SE", horizontal) ///
       (dot icc indicator if data == "PSID", horizontal msymbol(s)) ///
       (rspike ub lb indicator if data == "PSID", horizontal) ///
      , legend(order(1 "Germany" 3 "Sweden" 5 "United States (PSID)") pos(2)) ///
        yscale(reverse) ylabel(1 (3) 31, val) ///
        xscale(alt) ytitle("") xtitle("Sibling correlations in cognitive skills") ///
        note("{it:Note:} Germany: {it:N} = 1,989 individuals in 1,431 families" ///
             "Sweden: {it:N} = 652,940 individuals in 536,224 families" ///
             "United States: {it:N} = 2,868 individuals in 1,968 families" , span)

Shrinkage (partial pooling) in multilevel models

unzipfile "output7717034280080927434.zip", replace
use cntry essround happy using ESS1-6e01_0_F1, clear

// Drastically reduce sample size
bysort cntry essround: keep if _n < 5

// Encode country variable
encode cntry, gen(country)

// Fit random intercept model
mixed happy || country: 
predict u0, reffects
predict u0se, reses

// Calculate posterior intercept
generate eb_happy = u0 + _b[_cons]

// Plot posterior intercept
preserve
egen pickone = tag(country)
keep if pickone
   
egen order_ = rank(-u0), unique
labmask order_, value(country) decode

gen high = u0 + _b[_cons] + (1.96 * u0se)
gen low  = u0 + _b[_cons] - (1.96 * u0se)
local avg =  _b[_cons]

twoway (rcap eb_happy eb_happy order_, dsymbol(x)) ///
       (rspike high low order_) , ///
        yline(`avg') xlabel(1/32, val ang(v)) ///
        xtitle("") ///
        ylabel(4/10) ///
        ytitle("Country-level residuals + intercept") ///
        legend(off) ///
        title("Random intercept model") ///
        name(emp_bayes, replace) nodraw
restore

// Fit OLS models
regress happy
local avg =  _b[_cons]
regress happy i.country
predict avg_happy
predict se_happy, stdp

// Plot OLS estimates
preserve
egen pickone = tag(country)
keep if pickone
   
egen order_ = rank(-avg_happy), unique
labmask order_, value(country) decode

gen high = avg_happy + (1.96 * se_happy)
gen low  = avg_happy - (1.96 * se_happy)

twoway (rcap avg_happy avg_happy order_, dsymbol(x)) ///
       (rspike high low order_) , ///
        yline(`avg') xlabel(1/32, val ang(v)) ///
        xtitle("") ///
        ylabel(4/10) ///
        ytitle("Country-level means") ///
        legend(off) ///
        title("OLS regression model") ///
        name(plain, replace) nodraw
restore

graph combine emp_bayes plain, col(1) ysize(8) name(combo, replace)

// Plot EB and OLS means against one another to illustrate 
// shrinkage in multilevel modeling

preserve
egen pickone = tag(country)
keep if pickone
twoway (scatter eb_happy avg_happy) ///
       (function y = x, range(0 10)), ///
    legend(off) xtitle("OLS estimates") ytitle("Posterior intercept") ///
    title("Shrinkage")
restore

Generating a country–year variable

unzipfile "output7717034280080927434.zip", replace
use cntry essround happy using ESS1-6e01_0_F1, clear


// Generate country-year variable
levelsof cntry, local(levels1)
levelsof essround, local(levels2)
gen cyear = ""

foreach country of local levels1 {
  foreach round of local levels2 {
     replace cyear = "`country'" + "`round'" if cntry == "`country'" & essround == `round'
  di "`country'" "`round'"
  }
}
label var cyear "Country-Year"

// Fit three-level model
mixed happy || cntry: || cyear:
estat icc

capture drop u1* u0*
predict u1 u0, reffects 
predict u1se u0se, reses   

  // Plot country-level variation
preserve
egen pickone = tag(cntry)
keep if pickone
   
egen order_ = rank(-u1), unique
labmask order_, value(cntry)

gen high = u1 + (1.96 * u1se)
gen low  = u1 - (1.96 * u1se)

twoway (rcap u1 u1 order_, dsymbol(x)) ///
       (rspike high low order_) , ///
        yline(0) xlabel(1/32, val ang(v)) ///
        xtitle("") ///
        ytitle("Country-level residuals") ///
        legend(off) ///
        name(cntry, replace) 
restore 

  // Plot country-year level variation
preserve
egen pickone = tag(cyear)
keep if pickone
   
egen order_ = rank(-u0), unique
labmask order_, value(cyear)

gen high = u0 + (1.96 * u0se)
gen low  = u0 - (1.96 * u0se)

twoway (rcap u0 u0 order_, dsymbol(x) horizontal) ///
       (rspike high low order_, horizontal) , ///
        xline(0) ylabel(none) ///
  xscale(alt) ///
        ytitle("") ///
        xtitle("Country{c 150}year-level residuals") ///
        legend(off) ///
        name(cyear, replace) ysize(8)
restore 

Random graphs (53): Visualizing multilevel models

use cntry happy mnactic agea gndr using ESS1e06_4.dta, clear

// Prepare variables
recode agea (999 = .), gen(age)
recode gndr (1 = 0 "Male") (2 = 1 "Female") (9 = .), gen(female)
recode mnactic (77 88 99 = .), gen(active)
label var female "Sex (ref. Male)"
label var age    "Age in decades"
label var happy  "Happiness"
label var active "Labor market status"
label val active mnactic

center happy age // Center outcome to make intercept smaller for plot
replace c_age = c_age / 10

// Listwise deletion
drop if missing(female, c_age, c_happy, active, cntry)
drop happy gndr agea mnactic // Unecessary variables can go

// Fit model
mixed c_happy c_age i.female i.active || cntry: c_age, ml cov(uns)

// Save random part
local var_age = round(exp(_b[lns1_1_1:_cons])^2, .0001)
local var_int = round(exp(_b[lns1_1_2:_cons])^2, .0001)
local covaria = round(tanh(_b[atr1_1_1_2:_cons]) * ///
                       exp(_b[lns1_1_1:_cons]) * ///
        exp(_b[lns1_1_2:_cons]), .0001)

// Plot fixed part
coefplot, xline(0) ///
xtitle(" " "Estimates and 95% CI's") ///
scheme(s1mono) ///
ciopts(recast(rcap)) ///
coeflabels(_cons = "{bf:Intercept}" ///
           c_age = "{bf:Age} in decades" ///
     8.active = "Housework") ///  // Shorten label
headings(c_age = " " ///
         1.female  = "{bf:Sex} ({it:ref.} Male)" ///
         2.active = `""{bf:Labor market status}" "({it:ref.} Paid work)""' ///
   _cons = " ") ///
mlabel format(%9.2f) mlabposition(12) mlabgap(*1.5) ///
xscale(alt) msymbol(x) ///
ysize(8) xsize(4) ///
title("Fixed part", span) ///
name(fixed_part, replace) nodraw

/* The default sizes of the available area are -ysize(4)- and -xsize(5.5)-,
   by the way. Letter size is -ysize(11)- and -xsize(8.5)-*/

// Plot random part   
  // Estimate residuals
capture drop u1* u0*
predict u1 u0, reffects 
predict u1se u0se, reses   

  // Plot intercept variation
preserve
egen pickone = tag(cntry)
keep if pickone
   
egen order_ = rank(-u0), unique
labmask order_, value(cntry)

gen high = u0 + (1.96 * u0se)
gen low  = u0 - (1.96 * u0se)

twoway (rcap u0 u0 order_, dsymbol(x)) ///
       (rspike high low order_) , ///
        yline(0) xlabel(1/22, val ang(v)) ///
        title("Intercept variance = `var_int'") ///
        xtitle("") ///
        ytitle("Random intercept residuals") ///
        legend(off) ///
        name(rand_int, replace) nodraw
restore 

  // Plot slope variation
preserve
egen pickone = tag(cntry)
keep if pickone

egen order_ = rank(-u1), unique
labmask order_, value(cntry)

gen high = u1 + (1.96 * u1se)
gen low  = u1 - (1.96 * u1se)

twoway (rcap u1 u1 order_, dsymbol(x)) ///
       (rspike high low order_) , ///
        yline(0) xlabel(1/22, val ang(v)) ///
        xtitle("") ///
        title("Slope variance = `var_age'") ///
        ytitle("Random slope residuals") ///
        ylabel(-1 (.5) 1) ///  // Make y-axis identical to other plot,
        legend(off) ///        // otherwise x-axis looks different, too
        name(rand_slope, replace) nodraw
restore 

  // Plot intercept-slope covariance
preserve

gen predRandomSlope= (_b[_cons] + u0) + ((_b[c_age] + u1) * c_age)
qui sum c_age
local hi1 =  1*r(sd)
local hi2 =  2*r(sd)
local hi3 =  3*r(sd)
local lo1 = -1*r(sd)
local lo2 = -2*r(sd)
      
sort cntry c_age
twoway (line predRandomSlope c_age, connect(ascending)), ///
        ytitle("Predicted Happiness") /// 
        xtitle("") ///
        title("Intercept{c 150}slope covariance = `covaria'") ///
        xlabel(`lo1' "-1 SD" 0 "Average age" `hi1' "+1 SD" `hi2' "+2 SD" `hi3' "+3 SD") ///
        ylabel(,format(%6.1f)) ///
        name(covariance, replace) nodraw
restore

// Combine Figures
graph combine rand_int rand_slope covariance, col(1) ysize(8) title("Random part") name(random_part, replace) nodraw
graph combine fixed_part random_part, col(2) ysize(8) xsize(8) altshrink title("Random coefficient model")

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 (11): Visualizing multiple group differences

// Prepare data

// Formal care use by income quintiles, households with a child
input str5 geo inc1 inc3 inc5
AT 10 7 9
BE 17 38 57
BG 0 10 15
CH 9 27 53
CY 16 19 30
CZ 3 3 4
DE 21 22 23
DK 87 72 83
EE 16 22 14
EL 6 11 12
ES 29 30 45
FI 18 27 41
FR 15 60 64
HR 6 8 13
HU 7 14 15
IE 8 12 34
IS 37 44 34
IT 17 26 28
LT 2 16 10
LU 23 34 56
LV 7 24 11
MT 0 16 15
NL 27 55 70
NO 34 57 53
PL 0 2 4
PT 14 44 36
RO 5 5 13
SE 44 56 32
SI 41 39 38
SK 2 5 0
UK 20 50 53
EU-27 17 34 36
end

// Encode country variable
encode geo, gen(country)

// Create country variable sortet by the size of the income gap
gen diff = inc5 - inc1 // Income gap
egen order = rank(-diff), unique  // Create rank variable
labmask order, value(country) decode // Assign value
  // labels of country to variable order based on its values

// Run-of-the-mill dot plot
twoway dot inc1 inc3 inc5 order, vertical ///
      legend(label(1 "1{sup:st} income quintile (poorest)") ///
   label(2 "3{sup:rd} income quintile") ///
   label(3 "5{sup:th} income quintile (richest)") ///
   order(3 2 1) ring(0) pos(12)) ///
   msymbol(th oh t) ///
      xlabel(1/32, valuelabels ang(v)) /// // turn on labels
      xtitle("") ///
   ytitle("% formal care use of households with" ///
          "a child younger than 3 years of age") ///
   caption("Note: Countries sorted by the size of the difference */
            /* between the 5{sup:th} and 1{sup:st} quintile" ///
           "Source: EU-SILC 2010", span) ///
   name(incomegap1, replace)

// Without dots
scatter inc1 inc3 inc5 order,   ///
   legend(label(1 "1{sup:st} income quintile (poorest)") ///
          label(2 "3{sup:rd} income quintile") ///
          label(3 "5{sup:th} income quintile (richest)") ///
         order(3 2 1) ring(0) pos(12)) ///
   msymbol(th oh t) ///
   xlabel(1/32, valuelabels ang(v)) /// // turn on labels
   xtitle("") ///
   ytitle("% formal care use of households with" ///
          "a child younger than 3 years of age") ///
   caption("Note: Countries sorted by the size of the difference /*
           */between the 5{sup:th} and 1{sup:st} quintile" ///
           "Source: EU-SILC 2010", span) ///
   name(incomegap2, replace)

// With numbers instead of symbols
   // Not sure whether this is the most elegant move
reshape long inc, i(geo) j(j)

twoway (dot inc order if j == 1, mlabel(j) mlabpos(0) msymbol(i))  ///
           (dot inc order if j == 3, mlabel(j) mlabpos(0) msymbol(i))  ///
           (dot inc order if j == 5, mlabel(j) mlabpos(0) msymbol(i)),  ///
      legend(label(1 "1 1{sup:st} income quintile (poorest)") ///
             label(2 "3 3{sup:rd} income quintile") ///
             label(3 "5 5{sup:th} income quintile (richest)") ///
             order(3 2 1) ring(0) pos(12)) ///
   xlabel(1/32, valuelabels ang(v)) /// // turn on labels
      xtitle("") ///
   ytitle("% formal care use of households with" ///
          "a child younger than 3 years of age") ///
   caption("Note: Countries sorted by the size of the difference /*
           */ between the 5{sup:th} and 1{sup:st} quintile" ///
           "Source: EU-SILC 2010", span) ///
   name(incomegap3, replace)

// With ordinal numbers and without dots
label define j 1 "1{sup:st}" 3 "3{sup:rd}" 5 "5{sup:th}" 
label value j j 

twoway (dot inc order if j == 1, mlabel(j) mlabpos(0) msymbol(i) ndots(0))  ///
       (dot inc order if j == 3, mlabel(j) mlabpos(0) msymbol(i) ndots(0))  ///
       (dot inc order if j == 5, mlabel(j) mlabpos(0) msymbol(i) ndots(0)),  ///
      legend(label(1 "1{sup:st} income quintile (poorest)") ///
             label(2 "3{sup:rd} income quintile") ///
             label(3 "5{sup:th} income quintile (richest)") ///
             order(3 2 1) ring(0) pos(12)) ///
   xlabel(1/32, valuelabels ang(v)) /// // turn on labels
      xtitle("") ///
   ytitle("% formal care use of households with" ///
          "a child younger than 3 years of age") ///
   caption("Note: Countries sorted by the size of the difference /*
           */ between the 5{sup:th} and 1{sup:st} quintile" ///
           "Source: EU-SILC 2010", span) ///
   name(incomegap4, replace)

// With arrows
reshape wide // Get data back into old format

twoway (pcarrow inc5 order inc1 order) ///
       (scatter inc5 order) ///
 ,   legend(label(2 "5{sup:th} income quintile (richest)") ///
           label(1 "Difference between the 5{sup:th} /*
                           */ and 1{sup:st} income quintile") ///
     order(1 2) ring(0) pos(12)) ///
     xlabel(1/32, valuelabels ang(v)) /// // turn on labels
      xtitle("") ///
   ytitle("% formal care use of households with" ///
          "a child younger than 3 years of age") ///
   caption("Note: Countries sorted by the size of the /*
           */ difference between the 5{sup:th} and 1{sup:st} quintile" ///
           "Source: EU-SILC 2010", span) ///
   name(incomegap5, replace)

// With a line
twoway (rspike inc5 inc1 order) ///
       (scatter inc5 order) ///
       (scatter inc1 order, msymbol(oh)) ///    
       (scatter inc3 order, msymbol(th)) ///    
 ,   legend(label(3 "1{sup:st} income quintile (poorest)") ///
            label(2 "5{sup:th} income quintile (richest)") ///
            label(1 "Difference between the 5{sup:th} /*
                           */ and 1{sup:st} income quintile") ///
      label(4 "3{sup:rd} income quintile") ///
      order(2 4 3) ring(0) pos(12)) ///
     xlabel(1/32, valuelabels ang(v)) /// // turn on labels
      xtitle("") ///
   ytitle("% formal care use of households with" ///
          "a child younger than 3 years of age") ///
   caption("Note: Countries sorted by the size of the difference /*
           */ between the 5{sup:th} and 1{sup:st} quintile" ///
           "Source: EU-SILC 2010", span) ///
   name(incomegap6, replace)