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 (90): Regression coefficients

version 14
use 7348_F1.dta, clear

// Prepare variables:
// Country variable
rename Y11_Country ctry
decode ctry, gen(country)
kountryadd "Macedonia (FYROM)" to "Macedonia" add
kountry country, from(other) stuck
rename _ISO3N_ geo
kountry geo, from(iso3n) to(iso2c)
rename _ISO2C_ cntry
replace cntry = "XK" if country == "Kosovo"
drop geo 

label define Y11_Country 30 "Macedonia", modify

// Wave identifier
rename Wave wave

// Gender
generate female = (Y11_HH2a == 2)

// Work-family conflict
*factor Y11_Q12a Y11_Q12b Y11_Q12c, pcf
*tab Y11_Q12a wave, mis
*tab Y11_Q12b wave, mis
*tab Y11_Q12c wave, mis
*alpha Y11_Q12a Y11_Q12b Y11_Q12c, item
scores wfb = mean(Y11_Q12a Y11_Q12b Y11_Q12c)
generate wfc = 5 - wfb
drop wfb

eststo clear

levelsof wave, local(wave)

foreach x of local wave {
eststo: mixed wfc female || ctry: female if wave == `x', cov(uns)

preserve
// Get no. of countries
matrix groups = e(N_g)
local n_g = groups[1,1]

// Predict residuals
predict u1 u0, reffects
predict u1se u0se, reses

// Calculate posterior slope
generate eb_female = u1 + _b[female]

// Plot posterior slope
egen pickone = tag(country)  // Keep one case per countr
keep if pickone
   
egen order_ = rank(-u1), unique
labmask order_, value(ctry) decode

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

twoway (rcap eb_female eb_female order_, horizontal) ///
       (rspike high low order_, horizontal) , ///
        xline(`avg') ylabel(1/`n_g', val ang(h)) ///
        ytitle("") ///
  xscale(alt) ///
        xlabel(-.2 (.1) .4) ///
        xtitle("Female WFC disadvantage") ///
        legend(off) ///
        title(`: label (wave) `x'') ///
        name(emp_bayes_`x', replace) xsize(3) nodraw
restore
}

graph combine emp_bayes_1 emp_bayes_2 emp_bayes_3, col(3) xsize(6) name(variation, replace) altshrink ///
      note("{it:Source:} European Quality of Life Surveys, 2003-11. {it:Notes:} Horizontal line shows average coefficient, country-specific" ///
        "estimates indicate deviation from this average. Error bars are 95% CI's based on random-effects models.")



coefplot est1, bylabel(EQLS 2003) || ///
         est2, bylabel(EQLS 2007) || ///
         est3, bylabel(EQLS 2011) || ///
       , xlabel(0 (.05) .2) drop(_cons) xscale(alt) baselevel coeflabel(female = "Female") ///
         xtitle("Female WFC disadvantage") xline(0) byopts(row(1)) ciopts(recast(rcap))



coefplot est1 est2 est3, xlabel(0 (.05) .2) drop(_cons) xscale(alt) baselevel coeflabel(female = "Female") ///
         xtitle("Female WFC disadvantage") xline(0) byopts(row(1)) ciopts(recast(rcap)) ///
  legend(order(2 "EQLS 2003" 4 "EQLS 2007" 6 "EQLS 2011")) 

Random graphs (70): Bar graphs for interaction plots

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

// Select Round 5
keep if inlist(essround, 5)

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

// Satisfaction with WLB
rename stfjbot swlb

// Contracted working hours
recode wkhct (  0/20    = 1 "Marginal part-time") ///
             ( 21/34.75 = 2 "Substantial part-time") ///
             (35/220    = 0 "Full-time") ///
            , gen(parttime)
label var parttime "Working hours"

// Professional status
generate professional = .
replace  professional = (iscoco >= 1000 & iscoco <= 2999)
label var professional "Professional status"
label define professional 0 "Non-prof." 1 "Professional"
label val professional professional 

// Sex
recode gndr ( 2 =  1 "Female") ///
            ( 1 =  0 "Male") ///
            (.a = .a "No answer") ///
            , gen(female) label(female)
label var female "Gender"

// Select cases: Those in paid work and with a partner
keep if partner == 1 & mnactic == 1

// Are part-time workers more satisfied with their WLB than full-time employees?
qui regress swlb i.parttime i.country, cluster(country)
estimates store m1 

// Hypothesis 1b: The effect is stronger for marginal PT than substantial PT.
qui test 1.parttime == 2.parttime
local F = round(r(F), .001)
local p = round(r(p), .001)
local r = r(df)

qui margins, at(parttime=(0 1 2))
marginsplot, recast(bar) ///
             title("SWLB difference between full-time and part-time workers", size(large)) ///
             plotopts(fcolor(gs14) lcolor(black)) ytitle("Predicted SWLB") ylabel(6 (.5) 7.5, format(%6.1f) grid) ///
             name(noint, replace) xtitle("") ysize(3) ///
             note(" " ///
                  "Marginal and substantial part-time" ///
                  "differs significantly:" ///
                  "{it:F}(`r', `r(df_r)') = `F', {it:p} = `p'", ///
                  pos(11) ring(0) bmargin(small)) ///
             nodraw
    
// Are professional part-time workers less satisfied than non-professional part-time workers?
qui regress swlb i.parttime##i.professional i.country, cluster(country)
estimates store m2
qui margins, at(parttime=(0 1 2) professional=(0 1))
marginsplot, recast(bar) xdimension(professional) ///
             bydimension(parttime) byopts(row(1) noiyaxes imargin(zero) ///
    title("Interaction part-time status and professional status")) ///
    subtitle(, pos(6)) /// // Place label of by dimensions below plot 
             plotopts(fcolor(gs14) lcolor(black)) ytitle("Predicted SWLB") ylabel(6 (.5) 7.5, format(%6.1f)) ///
    name(int1, replace) xtitle("") ysize(3) nodraw
    
// Are women working part-time more satisfied with their SWLB than part-time working men?
qui regress swlb i.parttime##i.female i.country, cluster(country)
estimates store m3
qui margins, at(parttime=(0 1 2) female=(0 1))
marginsplot, recast(bar) xdimension(female) ///
             bydimension(parttime) byopts(row(1) noiyaxes imargin(zero) ///
    title("Interaction part-time status and gender")) ///
    subtitle(, pos(6)) /// // Place label of by dimensions below plot 
             plotopts(fcolor(gs14) lcolor(black)) ytitle("Predicted SWLB") ylabel(6 (.5) 7.5, format(%6.1f)) ///
    name(int2, replace) xtitle("") ysize(3) nodraw


// Output table and figure
esttab m1 m2 m3 using test.tex, compress replace se label nomtitles ///
                                indicate(Country dummies = *country) ///
                                varwidth(30) interaction(" X ") ///
                                title(Regression table\label{tab1}) ///
                                booktabs
graph combine noint int1 int2, col(1) ysize(9) ///

                               note("95% CI's based on cluster-robust standard errors")

Random graphs (62): Histograms

qui summarize challenging
local m = round(r(mean), .1)
local sd = round(r(mean), .1)
local n = r(N)
numlabel challenging, add mask("# ")
twoway histogram challenging, discrete yla(, valuelabel) horizontal gap(20)  ///
                              freq ytitle("") start(1) xlabel(0/4, grid) ///
                              title("How" "{bf:challenging}" "was this paper?", span) ///
                              xsize(3) note("Mean = `m', {it:SD} = `sd', {it:N} = `n'", span) ///
                              name(challenging, replace)

qui summarize interesting
local m = round(r(mean), .1)
local sd = round(r(mean), .1)
local n = r(N)
numlabel interesting, add mask("# ")
twoway histogram interesting, discrete yla(, valuelabel) gap(20) horizontal ///
                              freq ytitle("") start(1) xlabel(0/4, grid) ///
                              title("How" "{bf:interesting}" "was this paper?", span) ///
                              xsize(3) note("Mean = `m', {it:SD} = `sd', {it:N} = `n'", span) ///
                              name(interesting, replace)

qui summarize recommendable
local m = round(r(mean), .1)
local sd = round(r(mean), .1)
local n = r(N)
numlabel recommendable, add mask("# ")
twoway histogram recommendable, discrete yla(, valuelabel) gap(20) horizontal ///
                                freq ytitle("") start(1) xlabel(0/4, grid) ///
                                title("How likely are you to" "{bf:recommend}" "this paper to a friend?", span) ///
                                xsize(3) note("Mean = `m', {it:SD} = `sd', {it:N} = `n'", span) ///
                                name(recommendable, replace)
         
graph combine challenging interesting recommendable, row(1)

Random graphs (54): Cross-level interactions

use deleteme.dta

// Cross-level interaction model
qui mixed lsat c.wfc##c.gdp || country: wfc , cov(uns)
estimates store cross_lvl  
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g)          // Number of level 2 units
qui estadd scalar nc  = foo[1,1] // Number of level 2 units
qui estadd scalar v2  = exp(2*[lns1_1_1]_b[_cons])     // Slope variance
qui estadd scalar v1  = exp(2*[lns1_1_2]_b[_cons])     // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
                        exp([lns1_1_1]_b[_cons])  * ///
                        exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons])     // Residual variance

// Create table
esttab cross_lvl ///
     , se ///
       stats(v1 v2 v_e cov dev nc N, ///     Add variance components to table
             labels("Var(Intercept)" ///
                    "Var(Slope)" ///
                    "Var(Residual)" ///
                    "Cov(Int., Slope)" ///
                    "Deviance" ///
                    "No. clusters" ///
                    "No. individuals")) ///
       coeflabel(wfc "Work-family conflict" ///
                 gdp "GDP" ///
                 c.wfc#c.gdp "Work-family conflict X GDP" ///
                 _cons "Intercept") ///
       mtitles("Cross-level interaction") ///
       nonumbers varwidth(28)  modelwidth(24) ///
       keep(lsat:)                  // Drop variance components in weird shapes

// Interaction plot (A)
preserve
estimates restore cross_lvl

qui sum wfc // get standard deviation and mean
local wfcminsd = r(mean) - r(sd)
local wfcmean  = r(mean)
local wfcplusd = r(mean) + r(sd)

qui sum gdp // get standard deviation
local gdpmin2sd = r(mean) - 2*r(sd)
local gdpminsd  = r(mean) -   r(sd)
local gdpmean   = r(mean)
local gdpplusd  = r(mean) +   r(sd)
local gdpplu2sd = r(mean) + 2*r(sd)
  
margins, at(gdp = (`gdpmin2sd' `gdpminsd' `gdpmean' `gdpplusd' `gdpplu2sd') ///
            wfc = (`wfcminsd' `wfcmean' `wfcplusd')) vsquish
marginsplot, x(gdp) noci ///
             xlabel(`gdpmin2sd' "-2 SD" ///
                    `gdpminsd'  "-1 SD" ///
                    `gdpmean'  `""Average" "GDP""' ///
                    `gdpplusd'  "+1 SD" ///
                    `gdpplu2sd' "+2 SD") ///     
             xtitle("") ///
             ylabel(, format(%6.1f)) ///
             ytitle("Predicted happiness") ///
             plotopts(msymbol(none)) ///        // Turn off markers
             plot1opts(lpattern(longdash)) ///  // Define line types here
             plot2opts(lpattern(solid)) ///
             plot3opts(lpattern(shortdash)) ///
             legend(subtitle("Work{c 150}family conflict" , size(small)) ///
             order(1 "- 1 SD" 2 "Mean" 3 "+ 1 SD") size(small) ///
             pos(11) ring(0)) ///
             title("(A) Interaction plot") ///
             name(xlvl_plot1, replace)    
restore 
  
// Interaction plot (B)
preserve

estimates restore cross_lvl  

predict u1 u0, reffect

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

gen wfceffect = u1 + _b[wfc] + _b[c.wfc#c.gdp] * gdp

qui sum gdp // get standard deviation
local gdpmin2sd = r(mean) - 2*r(sd)
local gdpminsd  = r(mean) -   r(sd)
local gdpmean   = r(mean)
local gdpplusd  = r(mean) +   r(sd)
local gdpplu2sd = r(mean) + 2*r(sd)

twoway (scatter wfceffect gdp, mlabel(country) mlabpos(0) msymbol(none)) ///
       (function y = _b[wfc] + _b[c.wfc#c.gdp] * x, range(gdp)) ///
      , ytitle("Work{c 150}family conflict coefficient") ///
        xlabel(`gdpmin2sd' "-2 SD" `gdpminsd' "-1 SD" ///
               `gdpmean' `""Average" "GDP""' ///
               `gdpplusd' "+1 SD" `gdpplu2sd' "+2 SD") ///
        title("(B) Interaction plot") ///
        legend(order(2 "Work{c 150}family conflict coefficient by GDP") ///
               pos(7) ring(0)  size(small)) ///
        name(xlvl_plot2, replace)
restore    

// Combine graphs
graph combine xlvl_plot1 xlvl_plot2, row(1) ysize(3) xsize(5.5) altshrink

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 (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 (42): Stacked bar graph

use "ESS3e03_5.dta", clear

// Drop two remote countries

drop if cntry == "UA"
drop if cntry == "RU"

// Prepare variables of interest

  // Country variable
encode(cntry), gen(country)
label define country 12 "UK", modify

  // Keep only respondents asked about women (split ballot)
keep if icsbfm == 1

  /// Question about values
recode aftjbyc (1 2 = 1  "(Strongly) disapprove") ///
               (3   = 2  "Neutral") ///
               (4 5 = 3  "(Strongly) approve") ///
               (7/9 = .a "Missing/NA") ///
              , gen(value) label(value)
quietly tab value, gen(val)
ren val1 disapprove
ren val2 neutral
ren val3 approve

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

// Calculate average scores across all countries by gender
preserve
collapse (mean) disapprove neutral approve [pweight = dweight], by(female)
gen country = 50  // Assign some value to sample average
tempfile euaverage
save `euaverage', replace
restore

// Calculate average scores for each country by gender
collapse (mean) disapprove neutral approve [pweight = dweight], by(country female)

// Add average score
append using `euaverage'
label define country 50 "{bf:EU}", modify    // Add bold label for sample averages

// Plot
graph bar approve neutral disapprove if female == 1 ///
        , over(country, sort(approve) descending label(alternate)) stack percentages ///
          title("Popular (dis)approval of a full-time" "working woman with a child under 3 years of age") ///
          ytitle("% of women") yscale(range(0 100)) ylabel(0(20)100) legend(off) ///
          name(women, replace)

graph bar approve neutral disapprove if female == 0 ///
        , over(country, sort(approve) descending label(alternate)) stack percentages ///
          ytitle("% of men") yscale(range(0 100)) ylabel(0(20)100) ///
          legend(label(3 "(Strongly) disapprove") ///
                 label(2 "Neutral") ///
                 label(1 "(Strongly) approve") ///
                 order(1 2 3) pos(6) row(1)) ///
          caption(" " ///
                  "{it: Source:} European Social Survey 2006/07, own calculations." ///
                  `"{it: Note:} "EU" refers to average for the 20 EU member states in this Figure."', ///
                  span size(small)) ///
          name(men, replace)
  
graph combine women men, col(1) ysize(8)

Random graphs (39): Plotting on a log axis

Rather than looking at a variable in absolute terms, namely GDP per capita in the upper panel of the Figure, it is also sometimes helpful to look at it in terms of percentage increases, as in the lower panel. Equal distances on the x-axis refer to equal percentage increases in GDP per capita. Four equidistant points on the x-axis are labeled, each indicating a fourfold increase in GDP.

use wvs2005_v20090901a.dta, clear

// Country variable
// 1) Turn into string
decode v2, gen(ctry)

// 2) Abbreviate
kountry ctry, from(other) stuck marker
ren _ISO3N_ cntry
kountry cntry, from(iso3n) to(iso2c)
ren _ISO2C_ country

// Generate outcome: % in good health
generate goodhealth = 100 if v11 <= 2
replace  goodhealth =   0 if v11  > 2
replace  goodhealth =   . if v11 == .

// Collapse data set
collapse (mean) goodhealth [pw = v259], by(country)

// Generate year variable
gen year = 2006

// Preserve collapsed data set
preserve

// Get GDP from World Bank data base
wbopendata, language(en - English) country() topics() indicator(NY.GDP.PCAP.PP.CD) clear long

// Fix obtained data set
ren ny_gdp_pcap_pp_cd gdp
ren iso2code country
keep if year == 2006
keep gdp country
drop if country == ""

// Save obtained GDP data
tempfile gdp
save `gdp'

// Restore 
restore

// Merge GDP with collapsed data set
merge m:1 country using `gdp'
keep if _merge == 3
drop _merge

// Create labels with thousand separator
label define gdp 20000 "20,000" ///
                 40000 "40,000" ///
                 60000 "60,000"
label val gdp gdp

// Plot on unlogged axis
twoway (scatter goodhealth gdp, mlabel(country) mlabpos(0) msymbol(none)) ///
       (lfit goodhealth gdp) ///
      , legend(off) xtitle("GDP per capita, 2006, PPP in current international $") ///    
        ytitle("% in good health") ///
        xlabel(0(20000)60000, valuelabels) ///
        name(unlogged, replace)

// Generate logged variable
gen loggdp = log(gdp) * 1000  // Multiply by 1,000 because only integers can be labeled

// Generate numbers for labeling
//  Round them to three decimal digits, then multiply by 1,000 to get integers
local log1 = round(log(1000), .001) * 1000
local log2 = round(log(1000 * 4), .001) * 1000
local log3 = round(log(1000 * 4 * 4), .001) * 1000
local log4 = round(log(1000 * 4 * 4 * 4), .001) *1000

*di `log1' _skip(2) `log2' _skip(2) `log3' _skip(2) `log4' 

// Create labels for numbers to be labeled
label define loggdp `log1' "1,000" ///
                    `log2' "4,000" ///
                    `log3' "16,000" ///                    
                    `log4' "64,000"
label value loggdp loggdp

// Plot on log axis
twoway (scatter goodhealth loggdp, mlabel(country) mlabpos(0) msymbol(none)) ///
       (lfit goodhealth loggdp) ///
      , legend(off) xtitle("GDP per capita, 2006, PPP in current international $") ///   
        ytitle("% in good health") ///
        xlabel(`log1' `log2' `log3' `log4', valuelabels) ///
        name(logged, replace)

// Combine plots
graph combine unlogged logged, col(1) ysize(8)

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 (22): Graphing functions

clear

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

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

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 (5): Combining plots of means

set autotabgraphs on  // Recycle graph window

import excel "X:\Data\Desktop\means.xlsx", sheet("Sheet1") ///
             firstrow clear case(lower)
  // Read in aggregate data from Excel sheet:
  // - Variable name (variable) [string]
  // - Country name (country) [string]
  // - Average score (mean)
  // - Confidence bounds (lerr uerr)

// Create numerical country variable for x-axis:
// Fast alternative: Use -encode- if you do not
// want to manipulate the order of countries

generate cntry = 1 if country == "SE"
replace  cntry = 2 if country == "UK"
replace  cntry = 3 if country == "NL"
replace  cntry = 4 if country == "DE"
replace  cntry = 5 if country == "PT"

label define cntry 1 "SE" 2 "UK" 3 "NL" 4 "DE" 5 "PT"
label val cntry cntry
label var cntry "Country"

graph drop _all  // Drop graphs in memory, usefulfor re-running this file

twoway (scatter mean cntry if variable == "Job autonomy") ///
       (rcap lerr uerr cntry if variable == "Job autonomy") ///
       , legend(off) xlabel(, valuelabels) ///
         ytitle("Job autonomy") xscale(off) name(JA)

twoway (scatter mean cntry if variable == "WF culture") ///
       (rcap lerr uerr cntry if variable == "WF culture") ///
       , legend(off) xlabel(, valuelabels) ///
         ytitle("WF culture") xscale(off) name(WFC)

twoway (scatter mean cntry if variable == "WF supervisor support") ///
       (rcap lerr uerr cntry if variable == "WF supervisor support") ///
       , legend(off) xlabel(, valuelabels) ///
         ytitle("WF supervisor support") xscale(off) name(WFSS)

twoway (scatter mean cntry if variable == "WF co-worker support") ///
       (rcap lerr uerr cntry if variable == "WF co-worker support") ///
       , legend(off) xlabel(, valuelabels) ///
         ytitle("WF co-worker support") xscale(off) name(WFCS)

twoway (scatter mean cntry if variable == "FWA use") ///
       (rcap lerr uerr cntry if variable == "FWA use") ///
       , legend(off) xlabel(, valuelabels) ///
         ytitle("FWA use") name(FWA)

graph combine JA WFC WFSS WFCS FWA, col(1) xcommon imargin(b = 0 t = 0) ysize(10)