Random graphs (146): Bar graph

clear

// Install rcall command
*github install haghish/rcall, stable

// Clear memory of R session
rcall clear

//#install.packages("essurvey");
rcall: library(essurvey)

// Download data
rcall: download_rounds(c(9), ///
                       ess_email = "REGISTERED E-MAIL ADDRESS HERE", ///
                       output_dir = "data", ///
                       format = 'stata')
					   
// Open data set
local vars cntry gndr polintr dweight
use `vars' using "data/ESS9/ESS9e01_2.dta", clear

// Variables
generate interest = (inlist(polintr, 1, 2))
replace  interest =  . if missing(polintr)
generate female = (gndr == 2)

// Calculate gender gap per country
tempname foo
postfile `foo' str50 cntry men women gap gap_se using deleteme.dta, replace
levelsof cntry, local(cntry)
foreach i of local cntry {
 qui regress interest i.female [pw = dweight] if cntry == `"`i'"'
 qui margins i.female, post
 local   men    = _b[0.female] * 100
 local women    = _b[1.female] * 100
 qui regress interest i.female [pw = dweight] if cntry == `"`i'"'
 qui margins r.female, post
 local   gap    = _b[r1vs0.female] * -100
 local   gap_se = _se[r1vs0.female] * -100
 post `foo' ("`i'") (`men') (`women') (`gap') (`gap_se')
}
postclose `foo'

// Open data set with gender gap per country
use deleteme.dta, clear

kountry cntry, from(iso2c)
egen country = rank(gap), unique
labmask country, value(NAMES_STD)

generate ub = gap + 1.96 * gap_se
generate lb = gap - 1.96 * gap_se

generate y1 = -women
generate y2 = 0 
generate y3 = -women + men

generate x1 = country - 0.2
generate x2 = country + 0.2

mylabels 0(20)60, myscale(-@) local(show)

local spaces = 30 * " "

twoway (scatter y1 country, msymbol(p) yaxis(1 2)) ///
       (rbar y1 y2 x1, blcolor(gs10%40) bfcolor(gs10%40) barwidth(0.4)) ///
       (rbar y1 y3 x2, blcolor(gs5%40) bfcolor(gs5%40) barwidth(0.4)) ///
       (rbar gap y2 x2, blcolor(gs3%40) bfcolor(gs3%40) barwidth(0.4)) ///
       (rspike ub lb x2) ///
     , ylabel(`show',  grid ang(h)) /// 
       ylabel(0(10)40, grid ang(h) axis(2)) /// 
       ytitle("% interested in politics") ///
       ytitle("`spaces' Gender gap", axis(2)) ///
       xlabel(1/19, noticks ang(45) valuelabel labsize(small))  ///
       legend(order(2 "Men" 3 "Women" 4 "Gender gap" 5 "Gender gap 95% CI") ///
              col(4) pos(12) ring(0)) name(gendergap, replace)

Random graphs (144): Automatically labeling figures with letters

// Open ESS round 5
use agea stflife cntry if inlist(cntry, "DE", "GB", "NL", "SE") using "ESS5e03_4.dta", clear

// Country variable
kountry cntry, from(iso2c)
rename NAMES_STD country

// Restrict to complete cases
keep if !missing(agea, stflife, cntry)

// Center age
qui sum agea, detail
replace agea = `r(p99)' if agea > `r(p99)' & !missing(agea)
bys country: center agea, gen(age)
drop agea
label var age "Age (centered)"

// Letters
tokenize "`c(ALPHA)'"
local i 1

levelsof country, local(country)

foreach z of local country {
  qui sum age if country == `"`z'"', detail 
  foreach x of numlist 10 25 50 75 90 {
    local x`x': di %9.1f r(p`x')
    di `x`x''
  }
  qui regress stflife c.age##c.age if country == `"`z'"'
  
  qui margins, at(age = (`x10' `x25' `x50' `x75' `x90'))
  marginsplot, title(`"{bf:``i''} `z'"') ///
               ytitle("Life satisfaction" "(predicted)") ///
               xtitle("Age (percentiles)") ///
               xlabel(`x10' "10th" ///
                      `x25' "25th" ///
                      `x50' "50th" ///
                      `x75' "75th" ///
                      `x90' "90th") ///
               ylabel(, format(%6.1f)) ///
               recastci(rarea) ciopts(color(gs10)) ///
               name("``i''", replace) nodraw
  local ++i
}

graph combine A B C D, col(2) ycommon

Downloading the European Social Survey from within Stata

clear

// Install rcall command
*github install haghish/rcall, stable

// Clear memory of R session
rcall clear

//#install.packages("essurvey");
rcall: library(essurvey)

// Download data
rcall: download_rounds(c(1, 2, 3, 4, 5, 6, 7, 8, 9), ///
                       ess_email = "REGISTERED E-MAIL ADDRESS HERE", ///
                       output_dir = "C:/Users/User/Desktop/ess/data", ///
                       format = 'stata')

Downloading the European Social Survey from within Stata

clear

rsource, terminator(END_OF_R) rpath("Path of Rterm.exe") roptions(`"--vanilla"')

//#install.packages("essurvey");
library(essurvey)
//# Set the working directory
setwd("Where you want to download the data to");

download_rounds(c(1, 2, 3, 4, 5, 6, 7, 8),
                ess_email = "Registered email address",
                output_dir = "data",
                format = 'stata');
    
END_OF_R

Random graphs (132): Distributions

// Open ESS 2006 data
use ESS3e03_6.dta, clear

// Prepare variables
kountry cntry, from(iso2c) 
rename NAMES_STD country

recode agea (999 = .)
recode ygcdbyr (6666/9999 = .)
recode inwyye (9999 = .)
recode gndr (9 = .)

// Generate relevant variables
generate yearborn = inwyye - agea
gen grandparent =  ygcdbyr - yearborn

//Plot figure
twoway kdensity grandparent if gndr == 1, by(country, title(Men  , justification(left) span) note("")) xline(50) ylabel(0 .05 .10, format(%6.2f)) ytitle("Density") xtitle("") /*xtitle("Age at grandparenthood")*/ name(men, replace) nodraw
twoway kdensity grandparent if gndr == 2, by(country, title(Women, justification(left) span) note("")) xline(50) ylabel(0 .05 .10, format(%6.2f)) ytitle("Density") xtitle("Age at grandparenthood") name(women, replace) nodraw

graph combine men women, col(1) ysize(8) note("{it:Note:} Vertical line indicates age 50", size(vsmall))

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 (116): Histograms

use "ESS8e01.dta", clear

kountry cntry, from(iso2c) marker
rename NAMES_STD country

recode lrscale (77 88 99 = .)

histogram lrscale, by(country, note("")) ///
                   percent discrete ///
                   xlabel(0 (2) 10, val) ///
                   ylabel(0 (10) 50) ///
                   xtitle("Placement on left-right scale") ///
                   name(figure1, replace)

Random graphs (115): Bar graphs

// Data
use ESS7e02_1.dta, clear

// Drinking variable
recode alcfreq (7 = 0 "Never") (6 = 1 "Less than once a month") ///
               (5 = 2 "Once a month") (4 = 3 "2{c 150}3 times a month") ///
               (3 = 4 "Once a week") (2 = 5 "Several times a week") ///
               (1 = 6 "Every day") (77 88 99 = .), gen(alcohol)
label var alcohol "Alcohol consumption"

// Dichotomize alcohol variable
generate drinking = inlist(alcohol, 5, 6) if !missing(alcohol)      

// Country name variable
kountry cntry, from(iso2c) marker
rename NAMES_STD country
  
// Education
recode eisced (1 2   = 0 "Lower") ///
              (3 4 5 = 1 "Medium") ///
     (6 7   = 2 "High") ///
     (55/99 = .), gen(education)
label var education "Education"

// Alcohol consumption by country
histogram alcohol, percent disc horizontal ///
                   by(country, title("{bf:A} Frequency of drinking alcohohl", ///
                         justification(left) bexpand span) ///
                               note("")) ///
                   ylabel(0(1)6, val) ytitle("") ///
                   ysize(8) name(figure1, replace)

// Alcohol consumption by education by country
  // Collect predicted probabilities
capture program drop my_logit
program define my_logit, eclass
    syntax[if]
    marksample touse
    logit drinking if `touse'
    margins if `touse', post
    exit
end

statsby point = _b[_cons] se = _se[_cons], by(country education) clear: my_logit

 // Calculate stuff
replace point = point * 100
replace se = se * 100

gen lb = point - 1.96 * se
gen ub = point + 1.96 * se

  // Plot
twoway (bar point point education) ///
       (rspike lb ub education), ///
    by(country, legend(off) ///
                   title("{bf:B} Drinking by educational attainment", ///
                         justification(left) bexpand span) ///
                   note("")) ///
       xlabel(0 1 2, val ang(v)) ylabel(0 (10) 50) ///
       ytitle("% drinking more than once a week") ///
    ysize(8) ///
    name(figure2, replace)

graph combine figure1 figure2, col(2) ///
                               note(" " "{it:Source:} European Social Survey 2014", ///
                                    span justification(right) bexpand size(*.8))

Random graphs (114): Line plot and choropleth

use ESS1-7e01, clear

// Country name variable
kountry cntry, from(iso2c) marker
rename NAMES_STD country

// Create No contract variable
gen nocontract = (wrkctra == 3) if !missing(wrkctra)

// Calculate per year and country
preserve
statsby cont = _b[_cons] contse = _se[_cons], by(country essround) clear: regress nocontract

replace cont = cont * 100             // Convert proportion into precentage
generate  lb = cont - (contse * 100)
generate  ub = cont + (contse * 100)

label define essround 1 "2002" 2 "2004" 3 "2006" 4 "2008" 5 "2010" 6 "2012" 7 "2014", modify

sort essround country
twoway (rarea lb ub essround) ///
       (connected cont essround), ///+
        by(country, legend(off) note(" " "{it:Source:} ESS 2002-2014", span) ///
        title("{bf:A}", justification(left) bexpand span)) ///
        xtitle("") xlabel(1/7, val ang(45)) ytitle("% no contract") ///
        xsize(6) ysize(6) ///
        name(figurea, replace) 
restore
  
// Calculate average % no contracts per country
collapse nocontract, by(cntry)
replace nocontract = nocontract * 100
format nocontract %6.1f
saveold cont, replace
  
// Read in map data
  // Source: shapefile from http://www.naturalearthdata.com/downloads/10m-cultural-vectors/
shp2dta using "maps\ne_10m_admin_0_countries", database(database)  coordinates(coordinates) genid(id) replace

  // Restrict coordinates to Europe
  // Source: https://en.wikipedia.org/wiki/Extreme_points_of_Europe
use coordinates, clear
replace _Y = . if _Y < 36
replace _Y = . if _Y > 71 & !missing(_Y)
replace _X = . if _X < -28 
replace _X = . if _X > 33 & !missing(_X)
saveold europecoordinates, replace

  // Create data set with Europe internal borders
clonevar id = _ID
merge m:1 id using database, nogenerate
keep if CONTINENT == "Europe"
saveold europe, replace

  // Merge map data with set with data to be plotted
use ISO_A2 NAME id using database, clear
rename ISO_A2 cntry
replace cntry = "FR" if NAME == "France" // Somehow not correct in map
replace cntry = "NO" if NAME == "Norway" // Somehow not correct in map
drop if cntry == "-99"
merge 1:1 cntry using cont, keep(matched) nogenerate

spmap nocontract using europecoordinates, id(id) xsize(6) ysize(6) ///
                                     polygon(data(europe)) legstyle(2) clnumber(9) ///
                                     legend(position (9) ring(0)) fcolor(Oranges) ///
                                     legorder(hilo) ocolor(none ..) ///
                                     title("{bf:B}", justification(left) bexpand span) ///
                                     legtitle("{bf:% no contract, 2002-2014 averages}") legjunction({c 150}) ///
                                     note(" " "{it:Source:} ESS 2002-2014", span) name(figureb, replace)
          
graph combine figurea figureb, row(1) xsize(12)

Random graphs (112): Line plot

clear 
// Open Eurostat data
unzipfile "macrodata\lfsi_pt_a.zip"
insheet using lfsi_pt_a_1_Data.csv

// Prepare variables
replace value = "." if value == ":"
destring value, gen(fixedterm)
rename time year
replace geo = "France" if geo == "France (metropolitan)"
kountry geo, from(other) stuck
rename _ISO3N_ country
kountry country, from(iso3n) to(iso2c)
rename _ISO2C_ cntry

// Select data and save
keep if year >= 2004
keep fixedterm geo year cntry
save lfs_fixedterm, replace

// Open ESS data
use essround agea mnact wrkctra pspwght cntry using ESS1-7e01, clear

// Select data
keep if essround >= 2
keep if inrange(agea, 20, 64)
keep if mnact == 1

// Prepare variables
generate ess_fixedterm = (wrkctra == 2)
replace  ess_fixedterm = (ess_fixedterm * 100)
generate year = 2004 if essround == 2
replace  year = 2006 if essround == 3
replace  year = 2008 if essround == 4
replace  year = 2010 if essround == 5
replace  year = 2012 if essround == 6
replace  year = 2014 if essround == 7

// Point estimates and standard errors
statsby ess_fixedterm = _b[_cons] se = _se[_cons], clear by(cntry year): regress ess_fixedterm

// Merge Eurostat data with ESS data
merge 1:1 cntry year using lfs_fixedterm

// Select countries 
drop if inlist(cntry, "RU", "IL", "UA", "LV", "MK", "MT", "RO")

// Generate country name variable
kountry cntry, from(iso2c) 
ren NAMES_STD country

// Calculate confidence intervals
generate lb = ess_fixedterm - 1.96 * se
generate ub = ess_fixedterm + 1.96 * se

// Plot figure
sort year
twoway (rarea ub lb year, lcolor(white)) ///
       (connected ess_fixedterm year) ///
       (line fixedterm year), by(country, note("")) ///
        xtitle("") ytitle("Percentage of total employed (20{c 150}64 y)" "on temporary contract") ///
        legend(order(2 "ESS" 1 "95% CI" 3 "Eurostat") row(1)) xlabel(2004 (2) 2014)

Random graphs (109): Dropping one country at a time

// Open ESS 6 data
use stflife eduyrs agea gndr cntry using "ESS6e02_2.dta", clear

// Life satisfaction
recode stflife (77 88 99 = .), gen(lifesat)

// Education
sum eduyrs if eduyrs < 77, detail
generate education = eduyrs if eduyrs < 77
replace  education = r(p99) if education >= r(p99) & !missing(education)

// Age
recode agea (999 = .), gen(age)

// Gender
generate female = (gndr == 2) if gndr != 9

// Fit model across all countries and save estimates
regress lifesat education c.age##c.age i.female, cluster(cntry)
local opointest = _b[education]
local olb       = _b[education] - 1.96 * _se[education]
local oub       = _b[education] + 1.96 * _se[education]

// Define temporary objects
tempname foo
tempname foox
postfile `foo' str3 cntry pointest lb ub using `foox', replace

// Drop one country at a time
levelsof cntry, local(country)
foreach i of local country {
  qui regress lifesat education c.age##c.age i.female if cntry != "`i'"
  local pointest = _b[education]
  local lb       = _b[education] - 1.96 * _se[education]
  local ub       = _b[education] + 1.96 * _se[education]
  post `foo' ("`i'") (`pointest') (`lb') (`ub')
}
postclose `foo'

// Open estimates
use `foox', clear

// Creats country variable
kountry cntry, from(iso2c) 
rename NAMES_STD geo
encode geo, generate(country)
label define country 29 "Kosovo", modify

// Plot estimates
sort country
twoway (rarea lb ub country, horizontal color(gs14)) ///
       (function y = `opointest', horizontal range(country) lpattern(solid)) ///
       (function y = `oub', horizontal range(country) lpattern(dash)) ///
       (function y = `olb', horizontal range(country) lpattern(dash)) ///    
       (dot pointest country, horizontal) ///
      , ylabel(1/29, val) ytitle("") xscale(alt) ///
        legend(order(5 "Point estimate after excluding country" ///
                     2 "Point estimate from complete sample" ///
                     1 "95% CI after excluding country" ///
                     3 "95% CI (cluster robust) from complete sample") pos(6) span) ///
        ysize(7) name(robustness, replace)

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 (88): Model predictions

use gndr fcldbrn yrbrn health agea ///
    fltdpr flteeff slprl wrhpp fltlnl enjlf fltsd cldgng cntry ///
    using ESS3e03_5.dta, clear

// Gender variable
gen female = (gndr == 2) if gndr !=9

// Calculate age at first birth
recode fcldbrn (6666 = .a "NA") ///
               (7777 = .b "Refusal") ///
               (8888 = .c "Don't know") ///
               (9999 = .d "No answer") ///
               , gen(yrfirstbirth)
recode yrbrn   (7777 = .b "Refusal") ///
               (8888 = .c "Don't know") ///
               (9999 = .d "No answer") ///
               , gen(yrbirth)
generate afb = yrfirstbirth - yrbirth

qui centile afb, centile(.5 99) // Truncate
replace afb = r(c_1) if afb < r(c_1)
replace afb = r(c_2) if afb >= r(c_2) & !missing(afb)

// Self-rated health variable
recode health ( 1 =  4 "Very good") ///
              ( 2 =  3 "Good") ///
              ( 3 =  2 "Fair") ///
              ( 4 =  1 "Bad") ///
              ( 5 =  0 "Very bad") ///
              ( 7 = .a "Refusal") ///
              ( 8 = .b "Don't know") ///
              ( 9 = .c "No answer") ///
              , generate(srh)

// Depression variable
recode fltdpr flteeff slprl wrhpp fltlnl ///
       enjlf fltsd cldgng ///
       (7 = .a) (8 = .b) (9 = .c)
recode wrhpp enjlf (1 = 4) (2 = 3) (3 = 2) (4 = 1) (. = .)
scores depression = mean(fltdpr flteeff slprl wrhpp fltlnl enjlf fltsd cldgng)

// Age variable
recode agea (999 = .a)

// Models and plots
regress srh i.afb i.agea if female == 1, cluster(cntry)
qui margins, over(afb)
marginsplot, ytitle("Predicted self-rated health") xtitle("Age at first birth") ///
             name(womensrh, replace) recastci(rarea) ciopts(color(gs12)) xsize(3) ///
             title("") ylabel(2 (.2) 3.2, format(%6.1f)) nodraw

regress depression i.afb i.agea if female == 1, cluster(cntry)
qui margins, over(afb)
marginsplot, ytitle("Predicted depression") xtitle("Age at first birth") ///
             name(womendep, replace) recastci(rarea) ciopts(color(gs12)) xsize(3) ///
             title("") ylabel(1.2 (.2) 2.4, format(%6.1f)) nodraw

regress srh i.afb i.agea if female == 0, cluster(cntry)
qui margins, over(afb)
marginsplot, ytitle("Predicted self-rated health") xtitle("Age at first birth") ///
             name(mensrh, replace) recastci(rarea) ciopts(color(gs12)) xsize(3) ///
             title("") ylabel(2 (.2) 3.2, format(%6.1f)) nodraw

regress depression i.afb i.agea if female == 0, cluster(cntry)
qui margins, over(afb)
marginsplot, ytitle("Predicted depression") xtitle("Age at first birth") ///
             name(mendep, replace) recastci(rarea) ciopts(color(gs12)) xsize(3) ///
             title("") ylabel(1.2 (.2) 2.4, format(%6.1f)) nodraw

graph combine womensrh womendep, col(1) ysize(8) name(a, replace) ///
      title("Women", span) nodraw
graph combine   mensrh   mendep, col(1) ysize(8) name(b, replace) ///
      title("Men", span) nodraw
graph combine a b, col(2) ysize(8) xsize(8) ycommon xcommon ///
      note("{it:Source:} European Social Survey, round 3." ///
           "{it:Note:} All models control for age as dummy variables. 95% CI's based on robust SE's", span size(*.9))

Random graphs (86): Interaction plots

version 14
use ESS1e06_4.dta, clear

// Age variable
replace agea = . if agea == 999
qui centile(agea), centile(99.9)
replace agea = r(c_1) if age > r(c_1) & !missing(agea)
rename agea age

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

// Education variable
recode eduyrs (77 88 99 = .)
qui centile eduyrs, centile(99.9)
replace eduyrs = r(c_1) if eduyrs > r(c_1) & !missing(eduyrs)

// Fit model
regress happy c.eduyrs##c.age##c.age, robust

// Caculate education slope as function of age
qui margins, dydx(eduyrs) over(age)
marginsplot, recastci(rarea) recast(line) ciopts(color(gs12)) ///
             ytitle("Effect of education on happiness") ylabel(, format(%6.2f)) ///
             xtitle("Age") xlabel(20(10)90) xmtick(14(1)93) ///
             title("") name(slope, replace) xsize(4) nodraw

// Predict happiness for different values of age and education 
qui margins, at(eduyrs=(9 12 15) age=(14 20 (5) 90 93))
marginsplot, xdimension(age) recastci(rarea) recast(connected) ciopts(color(gs12)) ///
             ytitle("Predicted happiness") ylabel(, format(%6.1f)) ///
             xtitle("Age") xlabel(20(10)90) xmtick(14(1)93) ///
             legend(order(6 "15 yrs" 5 "12 yrs." 4 "9 yrs.") ///
                    title(Education, size(*.8)) pos(11) ring(0)) ///
             title("") name(prediction, replace) xsize(4) nodraw

// Plot
graph combine slope prediction, xsize(8) row(1) note("{it:Source:} European Social Survey, round 1 (2002/2003).")

Random graphs (82): Quadratic by categorical interaction

unzipfile output7138637360988071942.zip

use agea gndr happy using ESS1-6e01_1_F1, clear

// Generate age variable, truncated at 99.9th percentile
qui centile(agea), centile(99.9)
generate age = agea
replace  age = r(c_1) if age > r(c_1) & !missing(age)

// Generate gender variable 
generate female =(gndr == 2) if !missing(gndr)

// Fit model and plot
regress happy c.age##c.age##i.female
qui margins, at(age=(15(5)90) female=(0 1))
marginsplot, recastci(rarea) ciopts(color(gs12)) ytitle("Predicted happiness") ///
             legend(order(3 "Men" 4 "Women") ring(0) pos(2)) xtitle("Age") ///
             ylab(, format(%6.1f)) title("") ///
             note(" " "{it:Source:} European Social Survey, round 1{c 150}6." ///
                  "{it:Note:} Gray-shaded areas denote 95% confidence intervals.", span)

Random graphs (81): Plot average and confidence interval for each value of other variable

use pspwght happy gndr agea if inrange(agea, 15, 86) using ESS1-6e01_1_F1.dta, clear

// Create gender variable
generate female = (gndr == 2) if !missing(gndr)
label define female 0 "Men" 1 "Women"
label val female female

// Fit model
regress happy i.agea##i.female [pweight = pspwght]

// Calculate means and confidence intervals
margins, over(agea female)

// Plot means and confidence intervals
marginsplot, recastci(rarea) ciopts(color(gs12)) recast(line) ylabel(, format(%6.1f)) ///
             xtitle("Age") ytitle("Average happiness") ///
             title("Happiness over the life course") ///
             legend(pos(2) ring(0)) xlabel(15 (10) 86) ///
             plot1opts(lpattern(dash)) plot2opts(lpattern(solid)) /// 
             note(" " "{it:Source:} European Social Survey, rounds 1{c 150}6, weighted data" ///
                  "{it:Note:} Gray areas denote 95% confidence intervals", span)

Random graphs (80): Line plots using the -by- option

unzipfile output961465243413697143.zip, replace

use ESS1-6e01_1_F1.dta, clear

// Label waves
label define essround 1 "2002" 2 "2004" 3 "2006" 4 "2008" 5 "2010" 6 "2012"
label val essround essround 

// Keep only countries that contribute six observations
egen pickone = tag(cntry essround)
tab cntry essround if pickone // Look at no. of waves per country
egen numrounds = total(pickone), by(cntry) // Calculate no. of waves per country
keep if numrounds == 6

// Identify migrants
generate migrant = 1 if (brncntr == 2) & !missing(brncntr)
replace  migrant = 2 if (brncntr == 1 & (facntr == 2 | mocntr == 2)) & !missing(brncntr, facntr, mocntr)
replace  migrant = 0 if (brncntr == 1 &  facntr == 1 & mocntr == 1) & !missing(brncntr, facntr, mocntr)

label define migrant 0 "Native" 1 "First-generation migrant" 2 "Second-generation migrant"
label value migrant migrant

// Generate country identifier
kountry cntry, from(iso2c)
rename NAMES_STD country
replace country = "Slovenia" if country == "si"

// By country, year, and migrant status
collapse (mean) imbgeco imueclt imwbcnt [pweight = dweight], by(country essround migrant)

drop if migrant == .

sort country migrant essround
twoway (connected imbgeco essround if migrant == 0, ///
         by(country, note(" " "{it:Source:} European Social Survey, rounds 1{c 150}6, weighted data", span size(*.8)))) ///
       (connected imbgeco essround if migrant == 1) ///
       (connected imbgeco essround if migrant == 2) ///
      , xlabel(1/6, val ang(v)) ytitle("Immigration good for economy") ysize(8) xtitle("") name(imbgecocym, replace) ///
        legend(order(1 "Natives" ///
                     2 "First-gen." ///
                     3 "Second-gen.") ring(0) row(1)) 

Random graphs (75): Line plots

use essround cntry netuse dweight using ESS1-6e01_0_F1.dta, clear

// Drop parts of data and label wave identifier
drop if essround == 6 // Doesn't have internet use variable
drop if inlist(cntry, "IS", "LT") // drop countries with only one measurement
label define essround 1 "2002" 2 "2004" 3 "2006" 4 "2008" 5 "2010"
label val essround essround 

// Generate variable of interest
generate dailyuser = (netuse == 7) if !missing(netuse)

// Generate values by country and round
collapse (mean) dailyuser [pweight = dweight], by(cntry essround)
replace dailyuser = dailyuser * 100

// Generate country identifier
kountry cntry, from(iso2c)
rename NAMES_STD country

// Plot
twoway (connected dailyuser essround, ///
        by(country, note(" " "{it:Source:} European Social Survey, rounds 1{c 150}5, weighted data", span size(*.8)) ///
     col(3))), xlabel(, val) ytitle("% daily internet users") ysize(8) xtitle("")

Random graphs (73): Density functions

use happy eduyrs using ESS1-6e01_0_F1.dta, clear

// Pull in outliers
qui sum eduyrs, detail
replace eduyrs = r(p99) if eduyrs >= r(p99) & !missing(eduyrs)

// Calculate density functions at different years of education
kdensity happy if eduyrs ==  4, n(5000) k(gauss) bw(.8) gen(h0 d0) nograph
kdensity happy if eduyrs ==  8, n(5000) k(gauss) bw(.8) gen(h1 d1) nograph
kdensity happy if eduyrs == 12, n(5000) k(gauss) bw(.8) gen(h2 d2) nograph
kdensity happy if eduyrs == 16, n(5000) k(gauss) bw(.8) gen(h3 d3) nograph

// Scale and place density functions
replace d0 = (d0*18) +  4
replace d1 = (d1*18) +  8
replace d2 = (d2*18) + 12
replace d3 = (d3*18) + 16

// Plot
twoway (line h0 d0 if inrange(h0,3,10)) ///
       (line h1 d1 if inrange(h1,3,10)) ///
       (line h2 d2 if inrange(h2,3,10)) ///
       (line h3 d3 if inrange(h3,3,10)) ///
       (lfit happy eduyrs) ///
      , legend(off) xline(4 8 12 16) xlabel(0 (1) 22) ///
        xtitle("Years of education") ytitle("Happiness") ///
        yscale(range(3 10)) ylabel(3(1)10)