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 (118): Dot plot

clear
input str5 cntry eqls2003 eqls2007 ewcs2005 ewcs2010 eusilc2014
SE 0.021949 0.030043 0.00213  0.000918 .
FI 0.022075 0.013452 0.02246  0.011247 0.00012
HR .   0.018228 0.017975 0.013302 .
SK 0.011168 0.035384 0.013425 0.019773 0.003893
CZ 0.005486 0.064484 0.008917 0.005968 0.000863
NO .   0.0226  0.021383 0.018745 0.012686
CH .   .   0.028731 .   .
HU 0.018815 0.060927 0.044574 0.013885 0.019517
SI 0.025373 0.030885 0.095199 0.019257 0.001327
NL 0.035021 0.048701 0.021501 0.034217 .
BE 0.039831 0.088688 0.025436 0.006923 0.01386
EE 0.027082 0.048187 0.052969 0.052337 0.006435
FR 0.030219 0.074314 0.040906 0.01435  .
LT 0.049985 0.077502 0.05434  0.032847 0.006778
BG 0.037858 0.063697 0.065784 0.032313 0.039127
PL 0.023636 0.075742 0.06099  0.037228 .
LU 0.106816 0.102864 0.004992 0.006159 0.032854
DE 0.090262 0.040957 0.038157 0.035452 .
RO 0.063563 0.072842 0.070762 0.022769 .
LV 0.09474  0.077167 0.054927 0.034465 0.038258
IT 0.040173 0.128064 0.06421  0.036467 .
AT 0.057465 0.132496 0.095462 0.05458  .
ES 0.112064 0.127856 0.086716 0.076224 0.038867
DK 0.129036 0.103272 0.107504 0.040896 .
UK 0.120762 0.124342 0.150114 0.083238 0.024998
ME .   .   .   0.101533 .
PT 0.173582 0.179532 0.073722 0.113514 0.091105
RS-KM .   .   .   0.160115 .
MK .   0.194579 .   0.128292 .
IE 0.270422 0.359196 0.293468 0.240577 0.03082
AL .   .   .   0.275311 .
GR 0.404795 0.53023  0.286437 0.292037 0.144072
MT 0.302118 0.465496 0.387637 0.277371 .
CY 0.322423 0.54425  0.436497 0.395157 .
TR 0.564689 0.496668 0.689807 0.643683 .
end

scores average = mean(eqls2003 eqls2007 ewcs2005 ewcs2010)

kountry cntry, from(iso2c) marker
rename NAMES_STD country
replace country = "Kosovo" if country == "rskm"

egen order_ = rank(average), unique
labmask order_, value(country)

twoway (dot eqls2003 order_, horizontal) ///
       (dot eqls2007 order_, horizontal) ///
       (dot ewcs2005 order_, horizontal) ///
       (dot ewcs2010 order_, horizontal) ///
       (dot eusilc2014 order_, horizontal) ///
      , ylabel(1/35, val) ysize(9) ///
        ytitle("") xscale(alt) ///
        xtitle("Proportion without job contract") ///
        legend(order(1 "EQLS 2003" ///
                     2 "EQLS 2003" ///
                     3 "EWCS 2005" ///
                     4 "EWCS 2010" ///
                     5 "EU-SILC 2014") title("Data source") pos(5) ring(0))

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

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 (74): Arrow plot

use 7363_F1.dta, clear

// Prepare wave variable
label define wavex 1 "1991" 2 "1995" 3 "2000" 4 "2005" 5 "2010"
label val wave wavex

// Using the internet
fre y10_q_16
tab  y10_q_16 wave

keep if wave > 3  // Drop waves where variable is missing

// Create variable of interest
generate nointernetuse = .
replace  nointernetuse = (y10_q_16 == 7) if !missing(y10_q_16)

// Collapse data
collapse (mean) nointernetuse, by(countid wave)
replace nointernetuse = nointernetuse * 100

reshape wide nointernetuse, i(countid) j(wave)

// Sort variable labels
egen order_ = rank(-nointernetuse5), unique
labmask order_, val(countid) decode
label define order_ 8 "Macedonia", modify

twoway (pcarrow nointernetuse4 order_ nointernetuse5 order_, horizontal) ///
       (scatter order_ nointernetuse5 if missing(nointernetuse4)) ///
      , ylabel(1/34, val) legend(order(1 "Decrease 2005 to 2010" ///
                                       2 "No data for 2005") pos(7) ring(0) size(*.8)) ///
        ytitle("") ysize(8) xscale(alt) xtitle("% of labor force" ///
                                               "never using internet/email" ///
                                               "for professional reasons") ///
        xlabel(10 (10) 80, grid) ///
        note(" " "{it:Source:} European Working Conditions Survey 2005 and 2010", span)

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 (61): Dot plot for two groups

tempname foo
tempname foox
postfile `foo' str2 cntry perc_f perc_f_ll perc_f_ul ///
                          perc_m perc_m_ll perc_m_ul using `foox' , replace

levelsof cntry, local(levels1)

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

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

    capture proportion nonesep_m if cntry == "`country'"
    *matrix list r(table)
    matrix mcoefs  = r(table)
    local mperc    = mcoefs[1,2] * 100
    local mperc_ll = mcoefs[5,2] * 100
    local mperc_ul = mcoefs[6,2] * 100 
 
    di "`country'"  _skip(2) `fperc_ll' _skip(2) `fperc' _skip(2) `fperc_ul'
    di "`country'"  _skip(2) `mperc_ll' _skip(2) `mperc' _skip(2) `mperc_ul'

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

postclose `foo'

// Use posted data set
use `foox', clear

// Reverse sign of one variable
replace perc_m    = -perc_m 
replace perc_m_ll = -perc_m_ll
replace perc_m_ul = -perc_m_ul

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

// Plot
twoway (dot perc_m order_, horizontal) ///
    (dot perc_f order_, horizontal) ///
      , ylabel(1/31, val) legend(order(1 "Men" 2 "Women") pos(6) row(1)) ///
        ytitle("") xtitle("% of men/women pooling all income") ///
  xline(0) ///
  xlabel(-100 "100" -50 "50" 0 "0" 50 "50" 100 "100") ///
        name(somefigure, replace) ysize(8)

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 (45): Re-creating -graph dot- using -twoway dot-

This plot recreates this one, but adds confidence intervals.

preserve

// Calculate country-specific means to be plotted
tempname foo
tempname idealage
postfile `foo' quintile incocomp incocomp_lb incocomp_ub using `foo2', replace

levelsof(quintile), local(quintile)        
        
foreach x of local quintile {
 qui reg incocomp if quintile == `x', cluster(cntry)
 local incocomp   = _b[_cons]
 local incocomplb = _b[_cons] - (1.96 * _se[_cons])
 local incocompub = _b[_cons] + (1.96 * _se[_cons])
 post `foo' (`x') (`incocomp') (`incocomplb') (`incocompub')
}

postclose `foo'

// Plot country-specific means
use `foo2', clear
    
twoway (dot incocomp quintile) ///  -twoway scatter- doesn't have the lines
       (rcap incocomp_lb incocomp_ub quintile) ///
      , yscale(range(1.5 3.0)) ///
        ylabel(1.5 (.25) 3.0, format(%6.2f)) ///
 legend(off)
        title("Income quintiles")
        xlabel(1 `""1" "(poorest)""' /// Double-compound quotes allow for the line break
               2 3 4 5 `""5" "(richest)"') ///
        fxsize(40) /// Forces x-size to be 40 percent of its original size
        xtitle("") ///
        xscale(range(0.2 5.8)) /// Making this a bit bigger makes room for the labels
 name(byquintile_hor, replace) 

restore

// Calculate country-specific means to be plotted 2
preserve
tempname foo
tempname idealage
postfile `foo' str2 cntry incocomp incocomp_lb incocomp_ub using `foo2', replace


levelsof(cntry), local(cntry)        
        
foreach x of local cntry {
 qui reg incocomp if cntry == "`x'"
 local incocomp   = _b[_cons]
 local incocomplb = _b[_cons] - (1.96 * _se[_cons])
 local incocompub = _b[_cons] + (1.96 * _se[_cons])
 post `foo' ("`x'") (`incocomp') (`incocomplb') (`incocompub')
}

postclose `foo'

use `foo2', clear

// Create sorted-by-size variable
egen order_ = rank(-incocomp), unique
labmask order_, value(cntry)

// Plot country-specific means 2      
twoway (dot incocomp order_) ///
       (rcap incocomp_lb incocomp_ub order_) ///
      , yscale(off) /// Not necessary, as we're using the one from the other graph
        xlabel(1/23, val alt) ///
        title("Countries") ///
        xtitle("") ///
        legend(off) ///
        /* graphregion(margin(b=16)) */ /// Could help as -ycommon- isn't accounting for the different label sizes
        name(bycntry_hor, replace)

// Combine plots
graph combine byquintile_hor bycntry_hor, ///
      l1title("Mean income comparison orientation") ///
      ycommon row(1) ///
      imargin(zero) // Sets margin between both plots to 0

restore

Random graphs (43): Means with confidence intervals

use "ESS3e03_5.dta", clear

// Restrict to respondents 25-42 y
keep if agea >= 25 & agea <= 42 

// Generate variables of interest

// Split ballot identifier
*fre icsbfm

// Sex
generate female = (gndr == 2) if gndr != .a
drop if female == .

// Country
replace cntry = "UK" if cntry == "GB"

// iagpnt "In your opinion, what is the ideal age for a XXX 
//         to become a mother/father?"
recode iagpnt (  0 = .a "No ideal age") ///
              (777 = .b "Refusal") ///
              (888 = .c "Don't know") ///
              (999 = .d "No answer") ///
              (998 = .e "Split ballot") ///
              , gen(idealageparent)

clonevar idealagefather = idealageparent
replace  idealagefather = .e if icsbfm == 1

clonevar idealagemother = idealageparent
replace  idealagemother = .e if icsbfm == 2

// tochld "After what age would you say a woman/man is generally too old to
//        "consider having any more children?"
recode tochld (  0 = .a "Never too old") ///
              (777 = .b "Refusal") ///
              (888 = .c "Don't know") ///
              (999 = .d "No answer") ///
              (998 = .e "Split ballot") ///
              (997 = .f "Wrong age group") ///
              , gen(toooldforchild)

clonevar  toooldforchildf = toooldforchild
replace   toooldforchildf = .e if icsbfm == 1
label var toooldforchildf "Man too old for a(nother) child"

clonevar  toooldforchildm = toooldforchild
replace   toooldforchildm = .e if icsbfm == 2
label var toooldforchildm "Woman too old for a(nother) child"

// Set outliers and "never too old" to country-specific 99th percentile
levelsof(cntry), local(country)
foreach x of varlist toooldforchildf toooldforchildm {
  foreach y of local country {
   qui sum `x' if cntry == "`y'", detail
  *di `x' _skip(2) "`y'" _skip(2) r(p95) _skip(2) r(p99) 
   replace `x' = r(p99) if `x'   == .a ///
                         & cntry == "`y'"
   replace `x' = r(p99) if `x'    > r(p99) ///
                         & !missing(`x') ///
                         & cntry == "`y'"
 }
}

// Create temporary files and postfile
tempname foo
tempname idealage
postfile `foo' str2 cntry idealagem idealagemlb idealagemub ///
                          idealagef idealageflb idealagefub ///
                          toooldforchildfm toooldforchildflb toooldforchildfub  ///
                          toooldforchildmm toooldforchildmlb toooldforchildmub ///
                          using `idealage', replace

levelsof(cntry), local(country)
foreach x of local country {
 qui reg idealagemother if cntry == "`x'"
 local idealagem   = _b[_cons]
 local idealagemlb = _b[_cons] - (1.96 * _se[_cons])
 local idealagemub = _b[_cons] + (1.96 * _se[_cons])

 qui reg idealagefather if cntry == "`x'"
 local idealagef = _b[_cons]
 local idealageflb = _b[_cons] - (1.96 * _se[_cons])
 local idealagefub = _b[_cons] + (1.96 * _se[_cons])

 qui reg toooldforchildf if cntry == "`x'"
 local toooldforchildfm   = _b[_cons]
 local toooldforchildflb = _b[_cons] - (1.96 * _se[_cons])
 local toooldforchildfub = _b[_cons] + (1.96 * _se[_cons])

 qui reg toooldforchildm if cntry == "`x'"
 local toooldforchildmm = _b[_cons]
 local toooldforchildmlb = _b[_cons] - (1.96 * _se[_cons])
 local toooldforchildmub = _b[_cons] + (1.96 * _se[_cons])

post `foo' ("`x'") (`idealagem') (`idealagemlb') (`idealagemub') ///
                   (`idealagef') (`idealageflb') (`idealagefub') ///
                   (`toooldforchildfm') (`toooldforchildflb') (`toooldforchildfub') ///
                   (`toooldforchildmm') (`toooldforchildmlb') (`toooldforchildmub') 
}
postclose `foo'

use `idealage', clear

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

twoway (scatter toooldforchildmm order_) ///
       (rcap toooldforchildmub toooldforchildmlb order_) ///
       (scatter toooldforchildfm order_) ///
       (rcap toooldforchildfub toooldforchildflb order_) ///       
       , legend(label(1 "... women") ///
                label(3 "... men") ///
                order(3 1) pos(1) ring(0)) ///
         xlabel(1/23, val alt) ///
         ylabel(40(5)60) ///
         xtitle(" ") ytitle("Age in years") ///
         title("Age when one is too old to have a(nother) child for ...") ///
         name(tooold, replace)
   
drop order_

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

twoway (scatter idealagem order_) ///
       (rcap idealagemub idealagemlb order_) ///
       (scatter idealagef order_) ///       
       (rcap idealagefub idealageflb order_) ///
       , legend(label(1 "... mother") ///
                label(3 "... father") ///
                order(3 1) pos(1) ring(0)) ///
         xlabel(1/23, val alt) ///
         ylabel(20(5)40) ///
         xtitle(" ") ytitle("Age in years") ///
         title("Ideal age to become a ...") ///
         name(idealage, replace)

graph combine idealage tooold, ///
          col(1) ysize(8) /// 
          note("{it:Source:} European Social Survey Round 3, own calculations" ///
               "{it:Notes:} Respondents age 25{c 150}42 y only. Error bars denote 95 % CI's", span size(small))

Random graphs (37): Proportions with confidence intervals

use acountry asex aage ayear a602 a612 aweight using GGS_Wave1_V.4.2.dta, clear

// Create outcome variable
recode a612 (1 2 1601 = 1 "(Likely) infertile") ///
            (3 4 1501 1602 = 0 "(Likely) fertile") ///
            (2001 . .a .b .c = .x "Don't know"), ///
            gen(infertility)

replace infertility = 0 if a602 == 1 // Also fertile when currently pregnant
   
// Create country identifier
gen cntry = ""
replace cntry = "BG" if acountry == 11
replace cntry = "RU" if acountry == 12
replace cntry = "GE" if acountry == 13
replace cntry = "DE" if acountry == 14
replace cntry = "FR" if acountry == 15
replace cntry = "HU" if acountry == 16
replace cntry = "IT" if acountry == 17
replace cntry = "NL" if acountry == 18
replace cntry = "RO" if acountry == 19
replace cntry = "NO" if acountry == 20
replace cntry = "AT" if acountry == 21
replace cntry = "EE" if acountry == 22
replace cntry = "BE" if acountry == 23
replace cntry = "AU" if acountry == 24
replace cntry = "LT" if acountry == 25
replace cntry = "PL" if acountry == 26
replace cntry = "CZ" if acountry == 28


// Recode gender variable
gen sex = ""
replace sex = "Men" if asex == 1
replace sex = "Women" if asex == 2

// Drop unnecessary countries
drop if inlist(cntry, "IT", "AU", "NL")

// Prevalence plot by country

preserve

// Drop unnecesary age groups
keep if aage <= 39
keep if aage >= 34

svyset [pweight = aweight] // Declare weight

// Create temporary files and postfile
tempname foo
tempname infert
postfile `foo' str2 cntry str5 sex perc_infert perc_infert_ll perc_infert_ul using `infert' , replace

levelsof cntry, local(levels1)
levelsof sex, local(sex1)


// Calculate proportions by country and sex
foreach sex of local sex1 {
foreach country of local levels1 {
      // Calculate proportions
    capture proportion infertility if cntry == "`country'" & sex == "`sex'"
    *matrix list r(table)
    matrix coefs  = r(table)
    local perc    = coefs[1,2] * 100
    local perc_ll = coefs[5,2] * 100
    local perc_ul = coefs[6,2] * 100

    di "`country'"  _skip(2) `perc_ll' _skip(2) `perc' _skip(2) `perc_ul'

    if "`perc'" != "" {    // Make sure that the loop doesn't break if empty
    post `foo' ("`country'") ("`sex'") (`perc') (`perc_ll') (`perc_ul') 
  }
}
}
postclose `foo'

// Use posted data set
use `infert', clear

// Sort countries by size for women
egen order_ = rank(-perc_infert) if sex == "Women", unique
labmask order_, value(cntry)

// Plot for women
twoway (dot perc_infert order_ if sex == "Women", horizontal) ///
       (rcap perc_infert_ll perc_infert_ul order_ if sex == "Women", horizontal) ///
      , ylabel(1/14, val) legend(off) ///
        ytitle("") xtitle("% of {bf:women} aged 35{c 0150}39 reporting (suspected) infertility") ///
        /*caption("{it:Note:} Error bars denote 95 % confidence intervals", size(small) span)*/ ///
        name(infert3539women, replace)

// Drop women's order
drop order_

// Sort countries by size for men
egen order_ = rank(-perc_infert) if sex == "Men", unique
labmask order_, value(cntry)

// Plot for men
twoway (dot perc_infert order_ if sex == "Men", horizontal) ///
       (rcap perc_infert_ll perc_infert_ul order_ if sex == "Men", horizontal) ///
      , ylabel(1/13, val) legend(off) ///
        ytitle("") xtitle("% of {bf:men} aged 35{c 0150}39 reporting (suspected) infertility") ///
        /*caption("{it:Note:} Error bars denote 95 % confidence intervals", size(small) span)*/ ///
        name(infert3539men, replace)

// Combine graphs
graph combine infert3539women infert3539men ///
      , ysize(8) col(1) xcommon ///
        caption("{it:Source:} Generations and Gender Survey Wave 1 (2005{c 0150}2011}" ///
                "{it:Note:} Error bars denote 95 % confidence intervals", size(vsmall) span)
restore




// Plot infertility by age

// Recode age into groups
recode aage (20/24 = 1 "20–24") ///
            (25/29 = 2 "25–29") ///
            (30/34 = 3 "30–34") ///
            (35/39 = 4 "35–39") ///
            (40/45 = 5 "40–45") ///
            (1/19  = .a "Below 20 y") ///
            (46/85 = .b "46+ y") ///
          , gen(age)

// Save valuelabels in temporary do-file
tempfile valuelabels
label save age using `valuelabels'

svyset [pweight = aweight] // Declare weight

// Create temporary files and postfile
tempname foo
tempname infertage
postfile `foo' str2 cntry str5 sex age perc_infert perc_infert_ll perc_infert_ul using `infertage' , replace

levelsof cntry, local(levels1)
levelsof sex, local(sex1)

// Calculate proportion by age group, country, and sex
foreach sex of local sex1 {
foreach country of local levels1 {
  forvalues x = 1/5 {
      // Calculate proportions
    capture proportion infertility if age == `x' & cntry == "`country'" & sex == "`sex'"
    *matrix list r(table)
    matrix coefs  = r(table)
    local perc    = coefs[1,2] * 100
    local perc_ll = coefs[5,2] * 100
    local perc_ul = coefs[6,2] * 100

    di "`country'" _skip(2) `x' _skip(2) `perc_ll' _skip(2) `perc' _skip(2) `perc_ul'

    if "`perc'" != "" {    // Make sure that the loop doesn't break if 
    post `foo' ("`country'") ("`sex'") (`x') (`perc') (`perc_ll') (`perc_ul') 
 }
  }
}

}
postclose `foo'

// Use posted values
use `infertage', clear         

drop if perc_infert == .   // Drop empty rows

// Define and set saved value labels again
do `valuelabels'
label val age age

// Plot for women by country and age
twoway (rarea perc_infert_ll perc_infert_ul age if sex == "Women", fcolor(gs15) cmissing(no)) ///
       (line perc_infert age if sex == "Women", cmissing(no)) ///
      , by(cntry, cols(3) legend(off) note("")  ///
        caption("{it:Source:} Generations and Gender Survey Wave 1 (2005{c 0150}2011)" ///
                "{it:Note:} Gray areas denote 95 % confidence intervals", size(vsmall) span)) ///
        ytitle("% of {bf:women} reporting (suspected) infertility") xtitle("Age group") ///
        xlabel(1/5, val alternate) ///
        ysize(8) name(infertage, replace)

Random graphs (26): Plot means by country

use cntry stflife using "ESS1e06.3_F1.dta", clear

encode cntry, gen(country)  // De-string country identifier

preserve

statsby mean_ = _b[_cons] ///
        loci  = (_b[_cons] - 1.96 * _se[_cons]) ///
        hici  = (_b[_cons] + 1.96 * _se[_cons]) ///
      , by(country) total clear: /// // -total- fits model for total sample
        regress stflife

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

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

twoway (dot mean_ order_, msize(vsmall)) ///
       (rcap loci hici order_) ///
     , legend(off) xlabel(1/23, valuelabels ang(v)) ///
       ylabel(,format(%6.1f)) name(lifesat, replace)  ///
       xtitle("") ytitle("Life satisfaction (avg.)")

restore

Random graphs (21): Confidence interval plots

clear

input str33 fiel prop_women prop_women_se ci_low ci_hi str3 test_
"Education" 0.7365 0.0029 0.7308 0.7422 No
"Social and behavioral sciences" 0.5251 0.0163 0.4931 0.5569 Yes 
"Psychology" 0.8282 0.0183 0.7921 0.8639 No
"Nursing and caring" 0.9298 0.0046 0.9208 0.9388 No
"Therapy and rehabilitation" 0.8598 0.0336 0.7939 0.9257 No
"Child care and youth service" 0.9517 0.0178 0.9168 0.9866 No
"Social work and counseling" 0.8419 0.0238 0.7953 0.8885 No
end

// Transform from proportion to percentage
replace prop_women = prop_women * 100
replace ci_low     = ci_low * 100
replace ci_hi      = ci_hi * 100

// Sort educational fields by % females
egen order_ = rank(-prop_women), unique
labmask order_, val(fiel)

twoway (dot prop_women order_, horizontal) ///
       (rcap ci_low ci_hi order_, horizontal) ///
      , legend(off) ylabel(1/7, valuelabels) ///
        ytitle("Educational fields") ///
        xtitle("% Females") ///
        xline(50) xlabel(40 (10) 100, format(%6.0f)) ///
        note("{it: Source:} European Labor Force Survey 2009" ///
             "{it: Note:} Error bars denote 95% CI's", span)

clear
input str44 fiel prop_women prop_women_se ci_low ci_hi hard
"Science, mathematics, and computing" 0.445   0.0114 0.4226 0.4673 1
"Physical science" 0.383    0.0069 0.3694 0.3964 1
"Physics"  0.390    0.0288 0.3336 0.4464 1
"Mathematics"  0.571    0.0193 0.5331 0.6089 1
"Computer science"  0.252    0.0051 0.1460 0.3580 1
"Engineering, manufacturing, and construction"    0.161 0.0012 0.1586 0.1634 1
"Engineering and engineering trades"              0.199 0.0182 0.1633 0.2347 1
"Humanities"  0.717   0.0096   0.698 0.7358 0
"Social and behavioral science" 0.525    0.0163 0.4931 0.5569 0
"Psychology"  0.828   0.0183 0.7921 0.8639 0
"Sociology and cultural studies" 0.663   0.0354 0.5936 0.7323 0
end

// Transform from proportion to percentage
replace prop_women = prop_women * 100
replace ci_low     = ci_low * 100
replace ci_hi      = ci_hi * 100

// Sort educational fields by % females and by hard/soft
egen orderhard = rank(-prop_women) if hard, unique
labmask orderhard, val(fiel)

egen ordersoft = rank(-prop_women) if hard == 0, unique
labmask ordersoft, val(fiel)

twoway (dot prop_women orderhard, horizontal) ///
       (rcap ci_low ci_hi orderhard, horizontal) ///
      , legend(off) ylabel(1/7, valuelabels) ///
        ytitle("Hard educational fields") ///
        xtitle("% Females") ///
        xline(50) xlabel(20 (10) 100, format(%6.0f)) ///
 xscale(off) /// // remove x-axis
 name(hard, replace)

twoway (dot prop_women ordersoft, horizontal) ///
       (rcap ci_low ci_hi ordersoft, horizontal) ///
      , legend(off) ylabel(1/4, valuelabels) ///
        ytitle("Soft educational fields") ///
        xtitle("% Females") ///
        graphregion(margin(l=28)) /// Account for different y-axis label length
        xline(50) xlabel(20 (10) 100, format(%6.0f)) ///
 name(soft, replace)

graph combine hard soft,  col(1) imargin(b = 2 t = 1) ///
        note("{it: Source:} European Labor Force Survey 2009" ///
             "{it: Note:} Error bars denote 95% CI's", span)

graph export Graph.png 

Random graphs (16): Plots with confidence intervals

preserve

reg srh i.country [pweight = weight]

predict means, xb
predict semeans, stdp

gen loci = means - 1.96 * semeans
gen hici = means + 1.96 * semeans

keep if pickone
 
egen order_ = rank(-means), unique
labmask order_, value(country) decode

twoway (dot means order_) (rcap loci hici order_) ///
        , legend(off) xlab(1/29, valuelabels ang(v)) ///
         ylabel(, format(%6.1f)) ///
         ytitle("Self-rated health") ///
         xtitle("") ///
         note(`"Note: SRH ranges from 0 ("poor") to 4 ("excellent")"', span)

restore

preserve

logit srh_bin i.country [pweight = weight]

predict lr_index, xb
predict se_index, stdp
generate p_hat = exp(lr_index)/(1+exp(lr_index)) // Predicted probabilities
 // (following http://www.stata.com/support/faqs/...)
gen lb = lr_index - invnormal(0.975)*se_index
gen ub = lr_index + invnormal(0.975)*se_index
gen plb = exp(lb)/(1+exp(lb))
gen pub = exp(ub)/(1+exp(ub))

replace p_hat = p_hat * 100
replace plb = plb * 100
replace pub = pub * 100
keep if pickone
 
egen order_ = rank(p_hat), unique
labmask order_, value(country) decode

twoway (dot p_hat order_) (rcap plb pub order_) ///
       , legend(off) xlab(1/29, valuelabels ang(v)) ///
         ylabel(, format(%6.0f)) ///
         ytitle("% Poor health") xtitle("")

Random graphs (14): Combining dot plots

collapse (mean) incocomp, by(quintile)

graph dot incocomp, ///
      over(quintile, relabel(1 `""1" "(poorest)""'  5 `""5" "(richest)""')) ///
      vertical /// // Labels with line break
      ytitle("Income comparison orientation") ///
      yscale(range(1.5 3.0)) ///
      ylabel(1.5 (.25) 3.0, format(%6.2f)) /// // Format axis numbering
      exclude0 ///
      b1title("Income quintile averages") /// // Label over() axis
      name(byquintile, replace) ///
      fxsize(50) // Reduce size for combining

restore 
collapse (mean) incocomp, by(cntry)

// Create neatly ordered variable for x-axis
egen order = rank(-incocomp), unique
encode cntry, gen(geo)
labmask order, value(geo) decode

graph dot incocomp, over(order, label(alternate)) vertical ///
          ytitle("Income comparison orientation") ///
          name(bycountry, replace) ///
          ylabel(1.5 (.25) 3.0, format(%6.2f)) exclude0 yscale(off) ///
   b1title("Country averages") fxsize(100)

graph combine byquintile bycountry, xsize(6) ycommon imargin(zero)

restore
collapse (mean) incocomp, by(quintile cntry)

// Get data into proper shape
sort cntry quintile
reshape wide incocomp, i(cntry) j(quintile)

// Create ordered variable for x-axis labels
gen diff = incocomp5 - incocomp1 // Income comparison gap
egen order = rank(incocomp5), unique
encode cntry, gen(geo)
labmask order, value(geo) decode
sort order

// Create variables as marker labels
gen incocomp1l = 1
gen incocomp3l = 3
gen incocomp5l = 5

twoway (pcspike incocomp1 order incocomp5 order, lcolor(gs14)) ///
       (scatter incocomp1 order, mlabel(incocomp1l) mlabpos(0) msymbol(none)) ///
       (scatter incocomp5 order, mlabel(incocomp5l) mlabpos(0) msymbol(none)) ///
        , xlab(1/23, valuelabels ang(v)) ///
   ytitle("Average income comparison orientation") ///
   ylabel(, format(%6.1f)) ///
   xtitle("") ///
   name(bycountry, replace) ///
   legend(label(1 "Test") ///
          label(2 "1 First quintile (poorest)") ///
   label(3 "5 Fifth quintile (richest)") ///
   order(2 3) ring(0) pos(5))

Random graphs (12): Stacked bar graphs

use ahm2009_wide, clear
 
replace country_str = "EL" if country_str == "GR"
 
quietly tab contracttype, gen(contr) 
ren contr1 perma 
ren contr2 tempo 
ren contr3 selfe 
ren contr4 famil
 
collapse perma tempo selfe famil [iw=eu29weight], by(country_str)
 
replace perma = perma * 100 
replace tempo = tempo * 100 
replace selfe = selfe * 100 
replace famil = famil * 100
 
egen order = rank(-perma), unique 
labmask order, value(country_str)
 
graph hbar perma tempo selfe famil ///
   , stack /// Stack bars on top of one another
   over(order, label(labsize(*0.9))) /// // Reduce size of country abbreviations
   ytitle("Percentage first contract") /// 
   legend(label(1 "Permanent") /// 
          label(2 "Temporary") /// 
          label(3 "Self-employed") /// 
          label(4 "Family worker") /// 
          col(2) pos(6)) ///  
   xsize(8.25) ysize(10) /// 
   note("{it:Source:} EU-LFS AHM 2009 (weighted)", span) /// 
   name(contract, replace)