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 

Random graphs (63): Area plots

generate home = .
replace  home = 1 if homewk == 1 | homewk == 2
replace  home = 0 if homewk == 3

generate usually = .
replace  usually = 1 if homewk == 1
replace  usually = 0 if homewk == 2 | homewk == 3

generate sometimes = .
replace  sometimes = 1 if homewk == 2
replace  sometimes = 0 if homewk == 1 | homewk == 3

collapse (mean) home usually sometimes, by(country refyear)

replace home = home * 100
replace usually = usually * 100
replace sometimes = sometimes * 100

drop if country == 14 & refyear < 2006  

replace sometimes = sometimes + usually

twoway (area sometimes refyear, by(country, note("{it:Source:} EU Labor Force Survey 1998{c 150}2013, HR 2004/5 removed", span))) ///
       (area usually refyear, by(country)), ///
        ytitle("% of labor force working from home") ///
        xlabel(1998 2000 2005 2010 2013, ang(h) alternate) xtitle("") ///
        legend(order(2 "% usually" 1 "% sometimes") row(1) pos(9)) ///
        name(graph, replace)

Appending all files in a directory

cd "E:\eu-lfs\converted files"

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

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

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

use "`first'", clear                    

forvalues x = 2/`total_' {
    di `x'
    local y : word `x' of `filelist'
    append using "`y'", force
}

save eulfs, replace

Random graphs (47): Complex scatterplot

import excel "C:\table+graph field match countries.xlsx", sheet("Sheet2") cellrange(A1:I27) clear firstrow

renvars, lower // Switch variable names to lower case

list // Check data

// Bring variable names into systematic shape
foreach x of varlist healthwelfare education socialsciencesbusinesslaw ///
                     artshumanities services sciencemathcomputing ///
                     agricultureveterinary engineeringmanufacturing {
  ren `x' var_`x'
}

list // Check data

// Reshape data
reshape long var_, i(country) j(field_str) string

// Convert proportion into percentage

replace var_ = var_ * 100

// Get field_str variable into numerical format
encode field_str, gen(field1)

label define field1 1 "Agriculture and Veterinary" ///
                    2 "Arts and Humanities" ///
                    3 "Education" ///
                    4 "Engineering and Manufacturing" ///
                    5 "Health and Welfare" ///
                    6 "Science, Math, and Computing" ///
                    7 "Services" ///
                    8 "Social Sciences, Business, and Law", modify

// Sort by average segregation per field
list field1 var if country == "EU average"  
recode field1 (5 = 1 "Health and Welfare") ///
              (3 = 2 "Education") ///
              (8 = 3 "Social Sciences, Business, and Law") ///
              (2 = 4 "Arts and Humanities") ///     
              (7 = 5 "Services") ///          
              (6 = 6 "Science, Math, and Computing") ///
              (1 = 7 "Agriculture and Veterinary") ///
              (4 = 8 "Engineering and Manufacturing") ///
              , gen(field)
       
// Plot
sort field var_
twoway (scatter field var_ if country != "EU average") ///
       (scatter field var_ if country == "EU average", connect(l) lwidth(thick)) ///
      , ylabel(1/8, val) xscale(alt) ///
        xtitle("Percentage women by field of study" "in European countries") ///
        ytitle("") note(" " "{it:Source:} European Labor Force Survey 2011, own calculations.", span) ///
        legend(order(2) label(2 "Average" "across" "countries") ring(0) pos(1))

Random graphs (29): Line plots with confidence bands

use mainstat age country refyear coeff if refyear == 2011 using eulfs.dta, replace

drop refyear // Drop unnecessary variable
generate retired = (mainstat == 4) if mainstat != . // Create dummy variable for retired

svyset [pweight = coeff] // Declare weight

// Create temporary files and postfile
tempname foo
tempname retage
postfile `foo' cntry age perc_retired perc_retired_ll perc_retired_ul using `retage' , replace

levelsof country, local(levels1)

foreach country of local levels1 {
  forvalues x = 57 (5) 77 {
      // Calculate proportions
    capture proportion retired if age == `x' & country == `country'
    *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

    qui 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') (`x') (`perc') (`perc_ll') (`perc_ul') 
 }
  }
}

postclose `foo'

use `retage', clear         

drop if perc_retired == .   // Drop empty rows

  // Label variables
label define age 57 "55{c 150}59 y" 62 "60{c 150}64 y" 67 "65{c 150}69 y" ///
                 72 "70{c 150}74 y" 77 "75{c 150}79 y", modify
label value age age
label define cntry  1 "AT"  2 "BE"  3 "BG"  4 "CH"  5 "CY" ///
                    6 "CZ"  7 "DE"  8 "DK"  9 "EE" 10 "ES" ///                    
                   11 "FI" 12 "FR" 13 "GR" 14 "HU" 15 "IE" ///
                   16 "IS" 17 "IT" 18 "LT" 19 "LU" 20 "LV" ///
                   21 "MT" 22 "NL" 23 "NO" 24 "PL" 25 "PT" ///
                   26 "RO" 27 "SE" 28 "SI" 29 "SK" 30 "UK" ///
                 , modify
label value cntry cntry

  // Plot plots     
twoway (rarea perc_retired_ll perc_retired_ul age, fcolor(gs14)) ///
       (line perc_retired age) ///   
      , by(cntry, cols(3) legend(off) note("") ///
                  caption("{it:Source:} EU-LFS 2011" ///
                          "{it:Notes:} Gray areas denote 95 % confidence intervals"))) ///
        ytitle("% retired") xtitle("Age group") ///
        xlabel(57 (5) 77, val ang(v)) ///
        ysize(10) 

/* ,by(, note("")) suppresses the default "Graphs by" note) */
/* legend(off) needs to be in by(, legend(off)) to work */

Random graphs (25): Line plot

use country refyear coeff age wstator hwactual exist2j hwactua2 if refyear == 2011 using eulfs.dta, clear

// Handle missing values:
replace hwactual = 0 if age < 15
replace hwactual = 0 if inlist(wstator, 3, 4, 5)
replace hwactua2 = 0 if inlist(exist2j, 1, .)
replace hwactua2 = 0 if age < 15

// Sum up working hours from first and second job
gen hours = hwactual + hwactua2

// Drop unnecessary variables
drop refyear wstator exist2j

// Collapse to calculate means
collapse (mean) hours [pweight = coeff], by(country age)

// Label variables
label define age2 02  "0{c 150}4"  07  "5{c 150}9" ///
                 12 "10{c 150}14" 17 "15{c 150}19" ///
                 22 "20{c 150}24" 27 "25{c 150}29" ///
                 32 "30{c 150}34" 37 "35{c 150}39" ///
     42 "40{c 150}44" 47 "45{c 150}49" ///
                 52 "50{c 150}54" 57 "55{c 150}59" ///
     62 "60{c 150}64" 67 "65{c 150}69" ///
                 72 "70{c 150}74" 77 "75{c 150}79" ///
                 82 "80{c 150}84" 87 "85{c 150}89" ///
                 92 "90{c 150}94" 97 "95+", modify
label value age age2
label define country  1 "AT"  2 "BE"  3 "BG"  4 "CH"  5 "CY" ///
                      6 "CZ"  7 "DE"  8 "DK"  9 "EE" 10 "ES" ///                    
                     11 "FI" 12 "FR" 13 "GR" 14 "HU" 15 "IE" ///
                     16 "IS" 17 "IT" 18 "LT" 19 "LU" 20 "LV" ///
                     21 "MT" 22 "NL" 23 "NO" 24 "PL" 25 "PT" ///
                     26 "RO" 27 "SE" 28 "SI" 29 "SK" 30 "UK" ///
                 , modify
label value country country

// Plot

twoway (line hours age if country == 30) ///
       (line hours age if country ==  7) ///
       (line hours age if country == 22) ///
       (line hours age if country == 10) ///
      , caption("{it:Source:} EU Labor Force Survey 2011, own calculations", span) ///
        ytitle("Hours worked per week per capita") xtitle("Age group") ///
        xlabel( 2 (5) 97, val ang(45)) ylabel(0 (5) 30) ///
        legend(label(1 "United Kingdom") ///
               label(2 "Germany") ///
               label(3 "The Netherlands") ///
               label(4 "Spain") ring(0) pos(1))

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 (13): Interaction plots using -marginsplot-

 
estimates use m2   // Open estimates saved via -estimates save-

xtcloglog  // Show results
 
// Set esample() as explained in the help-file to -estimates save-:
mark miss_mod1
markout miss_mod1 fail2 female age1519 age2529 age3034 /// 
                  hatlev1d loggdp hatfield_main
estimates esample:  if miss_mod1
 
// Calculate predicted probabilities:
margins female#hatlev1d, pred(pu0) 

  
// Simple default marginsplot:
marginsplot, xtitle("Gender") ytitle("Predicted probability") /// 
             caption("{it:Source:} EU-LFS 2009 AHM", span) /// 
             title("") name(sexXeduc1, replace) 

  
// As a bar graph:
marginsplot, xtitle("") ytitle("Predicted probability") title("") /// 
             recast(bar) /// // Draw bar graph instead of default 
             xdimension(female hatlev1d) ///  
             /* Get both interacting variables onto x-axis */ ///
             caption("{it:Source:} EU-LFS 2009 AHM", span) /// 
             name(sexXeduc2, replace) 

 
// As a bar graph that better separates the different variable combinations:
marginsplot, xtitle("") ytitle("Predicted probability") ///  
             recast(bar)  ///
             xdimension(female) ///  
             bydimension(hatlev1d, ///
                         elabels(1 "Low education" ///
                                 2 "Medium education" /// 
                                 3 "High education")) /// 
      /* Use -bydimension- to better structure the x-axis */ ///
             name(sexXeduc3, replace) 

 
// As a proper-looking graph:
marginsplot, xtitle("") ytitle("Predicted probability") /// 
             recast(bar)  /// 
       xdimension(female) /// 
             bydimension(hatlev1d, elabels(1 "Low education" ///
                                           2 "Medium education" /// 
                                           3 "High education")) /// 
             byopts(row(1) /// All plots in one row
                    noiyaxes /// // Don't display individual y-axes
      imargin(zero) /// // No margin between by() dimensions  
                    title("") /// // Suppress title 
      caption("{it:Note:} Error bars are 95 % CI's" ///
      "{it:Source:} EU-LFS 2009 AHM", span)) /// 
                    /* Caption here to suppress displaying it three times */ ///
             plotregion(lwidth(none)) /// 
             /* Removes lines around separate by dimension plots */
      subtitle(, pos(6)) /// // Place label of by dimensions below plot 
      plotopts(fcolor(gs14) lcolor(b)) /// 
             /* Options of the recast() type can go here as well */ ///
      name(sexXeduc4, replace)  

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) 

Random graphs (8): Coefficient plot by hand

While the result certainly looks cool, such a graph can be extremely tedious to make.

 
// Requires parmest, ingap, and eclplot
 
use eulfsahm2010, clear
 
// Create and label dummy variables for regression model
fre age
quietly tab age, gen(age)
label var age1 "15-19 years"
label var age2 "20-24 years"
label var age3 "25-29 years"
label var age4 "30-34 years"
label var age5 "35-39 years"
label var age6 "40-44 years"
label var age7 "45-49 years"
label var age8 "50-54 years"
label var age9 "55-59 years"
label var age10 "60-64 years"
 
fre marstat
quietly tab marstat, gen(marstat)
label var marstat1 "Widowed, divorced, or separated" 
label var marstat2 "Single" 
label var marstat3 "Married"
 
fre education
quietly tab education, gen(education)
label var education1 "Low education" 
label var education2 "Medium education" 
label var education3 "Tertiary education" 
 
fre parttime
quietly tab parttime, gen(parttime)
label var parttime1 "Marginal part-time (< 20 hrs.)" 
label var parttime2 "Substantial part-time (20-34 hrs.)" 
label var parttime3 "Full-time (35-80 hrs.)" 

fre occgroups 
quietly tab occgroups, gen(occgroups) 
label var occgroups1 "Unskilled job" 
label var occgroups2 "Skilled manual job" 
label var occgroups3 "Skilled routine services job" 
label var occgroups4 "High-skilled services" 
 
fre industry 
quietly tab industry , gen(industry) 
label var industry1  "Agriculture" 
label var industry2  "Industry and construction" 
label var industry3  "Market services" 
label var industry4  "Non-market services" 
 
// Re-label some more variables here. 
// Variable labels HAVE to be identical to those specified below. 
// Formatting is possible here.
label var female     "{bf:Female sex}"
label var smallfirm  "{bf:Small firm} (< 10 employees)" 
label var fixedterm  "{bf:Fixed-term contract}" 
label var jobtenure  "{bf:Job tenure} (in years)"
label var supervisor "{bf:Supervisor status}" 

// Estimate regression model
ologit posstend female age1-age3 age5-age10 marstat1 marstat3 /// 
       education1 education3 parttime1 parttime2 smallfirm fixedterm /// 
       jobtenure supervisor occgroups1 occgroups3 occgroups4 /// 
       industry1 industry3 industry4, cluster(cntry) 
 
// Save regression parameters in data set
parmest, norestore eform label  
 
// Add gaps to add some structure to the regression plots
ingap 2 11 13 15 21 24, rowlabel(label) /// 
                        growlabel("Age ({it:Ref.} 30-34 years)" /// 
                        "Marital status ({it:Ref.} Single)" /// 
   "Education ({it:Ref.} Medium)" ///  
   "Working hours ({it:Ref.} Full-time (35-80 hrs)" ///
   "Occup. group ({it:Ref.} Skilled manual job)" /// 
   "Industry ({it:Ref.} Industry and construction)")
           
// Define label by hand to be able to determine order of predictors on axis 
  // Labels must be exactly the same as above 
  // 
  // Formatting axis label text seems 
  // to be working here as well 
  // 
  // An easier solution might be the one here: 
  // http://www.survey-design.com.au/Stata%20Graphs.html 
  // ("Odds ratio graph - selecting order of bars")
#delimit ;
label define label1    1 "{bf:Female sex}"                            
                       2 "{bf:Age} ({it:Ref.} 30-34 years)"           
                       3 "15-19 years" 
                       4 "20-24 years" 
                       5 "25-29 years"
                       6 "35-39 years" 
                       7 "40-44 years" 
                       8 "45-49 years"
                       9 "50-54 years"  
                      10 "55-59 years" 
                      11 "60-64 years" 
                      12 "{bf:Marital status} ({it:Ref.} Single)"
                      13 "Widowed, divorced, or separated" 
                      14 "Married" 
                      15 "{bf:Education} ({it:Ref.} Medium)" 
                      16 "Low education" 
                      17 "Tertiary education" 
                      18 "{bf:Working hours} ({it:Ref.} Full-time (35-80 hrs)" 
                      19 "Marginal part-time (< 20 hrs.)"  
                      20 "Substantial part-time (20-34 hrs.)"
                      21 "{bf:Small firm} (< 10 employees)"
                      22 "{bf:Fixed-term contract}" 
                      23 "{bf:Job tenure} (in years)" 
                      24 "{bf:Supervisor status}" 
                      25 "{bf:Occupational group} ({it:Ref.} Skilled manual job)" 
                      26 "Unskilled job" 
                      27 "Skilled routine services job" 
                      28 "High-skilled services" 
                      29 "{bf:Industry} ({it:Ref.} Industry and construction)"
                      30 "Agriculture" 
                      31 "Market services" 
                      32 "Non-market services" 
                      33 "Constant"  
                      34 "Constant" ; 
#delimit cr 

encode label, gen(label1) label(label1)
  // Specifying a label in the -encode- command makes sure that -encode- 
  // does not sort labels in alphabetical order

eclplot estimate min95 max95 label1 /// 
        if label1 <= 32                /// Don't plot the cut values
      , horizontal ylabel(1(1)32, labsize(vsmall)) /// Reduce font size on axis 
        xline(1)                       /// Add vertical line
        xtitle("Odds ratio") ///
 ytitle("Predictors") /// 
 title("Possibility to vary start and/or end of working day", span) /// 
 caption("Source: EU-LFS AHM 2010, authors' calculations." /* 
        */ "Ordered logit regression, /* 
        */ 95% CI's adjusted for clustering in countries", span) 

Random graphs (3): Kaplan–Meier plots by subsample

This approach of graphing a Kaplan–Meier survival curve is an alternative to -sts graph-, which does not seem to allow for combining overlaying and separating Kaplan–Meier curves of subgroups.

stset ahm2009_jobsearch, failure(fail) /* Declares survival data */

sts generate surv_sex = s, by(country sex)

/* s refers to Kaplan-Meier product-limit estimate of the survivor function,
   -sts generate- is also able to produce other estimates and their standard
   errors. */

twoway  (line surv_sex _t if sex == 1, sort) ///
        (line surv_sex _t if sex == 2, sort) ///
       , by(country, title("Kaplan–Meier survival estimates") col(3)) ///
         legend(label(1 "Males") label(2 "Females") col(2)) ///
         xtitle("Time in months") ytitle("Survival probability") ///
         xsize(5.5) ysize(10)

/* _t has been created by -stset- and refers to the analysis time when
   record ends. */

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