Showing posts with label twoway rarea. Show all posts
Showing posts with label twoway rarea. Show all posts

Nov 17, 2017

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


Oct 11, 2017

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)

Sep 7, 2017

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)


Aug 12, 2017

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)

Jun 20, 2015

Random graphs (49): Means with confidence intervals


tempname foo
tempname cogscore
postfile `foo' time ART cog coglb cogub using `cogscore', replace

levelsof(time), local(timepoint)
levelsof(ART), local(art)
foreach y of local art {
  foreach x of local timepoint {
    qui reg COG if time == `x' & ART == `y'
    local cog   = _b[_cons]
    local coglb  = _b[_cons] - (1.96 * _se[_cons])
    local cogub = _b[_cons] + (1.96 * _se[_cons])
 di `timepoint' _skip(5) `art' _skip(5) `cog' _skip(5) `coglb' _skip(5) `cogub'
    post `foo' (`x') (`y') (`cog') (`coglb') (`cogub')
  }
}
postclose `foo'

use `cogscore', clear

label define time 1 "Age 3" 2 "Age 5" 3 "Age 7" 4 "Age 11"
label val time time

twoway (rarea coglb cogub time if ART == 1) ///
       (rarea coglb cogub time if ART == 0) ///
       (line cog time if ART == 0) ///
       (line cog time if ART == 1) ///
      , legend(label(3 "Planned natural conception") ///
               label(4 "ART conception") col(1) ring(0) pos(11) order(4 3)) ///
        ytitle("Average cognitive development") ///
        xtitle(" ") xlabel(1(1)4, valuelabels) ///
        note("{it:Note:} Shaded areas denote 95 % confidence intervals." ///
             "{it:Source:} Millennium Cohort Study.", span) name(bands, replace)

Nov 30, 2014

Random graphs (37): Proportions with confidence intervals

use acountry asex aage ayear a602 a612 aweight using GGS_Wave1_V.4.2.dta, clear // Create outcome variable recode a612 (1 2 1601 = 1 "(Likely) infertile") /// (3 4 1501 1602 = 0 "(Likely) fertile") /// (2001 . .a .b .c = .x "Don't know"), /// gen(infertility) replace infertility = 0 if a602 == 1 // Also fertile when currently pregnant // Create country identifier gen cntry = "" replace cntry = "BG" if acountry == 11 replace cntry = "RU" if acountry == 12 replace cntry = "GE" if acountry == 13 replace cntry = "DE" if acountry == 14 replace cntry = "FR" if acountry == 15 replace cntry = "HU" if acountry == 16 replace cntry = "IT" if acountry == 17 replace cntry = "NL" if acountry == 18 replace cntry = "RO" if acountry == 19 replace cntry = "NO" if acountry == 20 replace cntry = "AT" if acountry == 21 replace cntry = "EE" if acountry == 22 replace cntry = "BE" if acountry == 23 replace cntry = "AU" if acountry == 24 replace cntry = "LT" if acountry == 25 replace cntry = "PL" if acountry == 26 replace cntry = "CZ" if acountry == 28 // Recode gender variable gen sex = "" replace sex = "Men" if asex == 1 replace sex = "Women" if asex == 2 // Drop unnecessary countries drop if inlist(cntry, "IT", "AU", "NL") // Prevalence plot by country preserve // Drop unnecesary age groups keep if aage <= 39 keep if aage >= 34 svyset [pweight = aweight] // Declare weight // Create temporary files and postfile tempname foo tempname infert postfile `foo' str2 cntry str5 sex perc_infert perc_infert_ll perc_infert_ul using `infert' , replace levelsof cntry, local(levels1) levelsof sex, local(sex1) // Calculate proportions by country and sex foreach sex of local sex1 { foreach country of local levels1 { // Calculate proportions capture proportion infertility if cntry == "`country'" & sex == "`sex'" *matrix list r(table) matrix coefs = r(table) local perc = coefs[1,2] * 100 local perc_ll = coefs[5,2] * 100 local perc_ul = coefs[6,2] * 100 di "`country'" _skip(2) `perc_ll' _skip(2) `perc' _skip(2) `perc_ul' if "`perc'" != "" { // Make sure that the loop doesn't break if empty post `foo' ("`country'") ("`sex'") (`perc') (`perc_ll') (`perc_ul') } } } postclose `foo' // Use posted data set use `infert', clear // Sort countries by size for women egen order_ = rank(-perc_infert) if sex == "Women", unique labmask order_, value(cntry) // Plot for women twoway (dot perc_infert order_ if sex == "Women", horizontal) /// (rcap perc_infert_ll perc_infert_ul order_ if sex == "Women", horizontal) /// , ylabel(1/14, val) legend(off) /// ytitle("") xtitle("% of {bf:women} aged 35{c 0150}39 reporting (suspected) infertility") /// /*caption("{it:Note:} Error bars denote 95 % confidence intervals", size(small) span)*/ /// name(infert3539women, replace) // Drop women's order drop order_ // Sort countries by size for men egen order_ = rank(-perc_infert) if sex == "Men", unique labmask order_, value(cntry) // Plot for men twoway (dot perc_infert order_ if sex == "Men", horizontal) /// (rcap perc_infert_ll perc_infert_ul order_ if sex == "Men", horizontal) /// , ylabel(1/13, val) legend(off) /// ytitle("") xtitle("% of {bf:men} aged 35{c 0150}39 reporting (suspected) infertility") /// /*caption("{it:Note:} Error bars denote 95 % confidence intervals", size(small) span)*/ /// name(infert3539men, replace) // Combine graphs graph combine infert3539women infert3539men /// , ysize(8) col(1) xcommon /// caption("{it:Source:} Generations and Gender Survey Wave 1 (2005{c 0150}2011}" /// "{it:Note:} Error bars denote 95 % confidence intervals", size(vsmall) span) restore
// Plot infertility by age // Recode age into groups recode aage (20/24 = 1 "20–24") /// (25/29 = 2 "25–29") /// (30/34 = 3 "30–34") /// (35/39 = 4 "35–39") /// (40/45 = 5 "40–45") /// (1/19 = .a "Below 20 y") /// (46/85 = .b "46+ y") /// , gen(age) // Save valuelabels in temporary do-file tempfile valuelabels label save age using `valuelabels' svyset [pweight = aweight] // Declare weight // Create temporary files and postfile tempname foo tempname infertage postfile `foo' str2 cntry str5 sex age perc_infert perc_infert_ll perc_infert_ul using `infertage' , replace levelsof cntry, local(levels1) levelsof sex, local(sex1) // Calculate proportion by age group, country, and sex foreach sex of local sex1 { foreach country of local levels1 { forvalues x = 1/5 { // Calculate proportions capture proportion infertility if age == `x' & cntry == "`country'" & sex == "`sex'" *matrix list r(table) matrix coefs = r(table) local perc = coefs[1,2] * 100 local perc_ll = coefs[5,2] * 100 local perc_ul = coefs[6,2] * 100 di "`country'" _skip(2) `x' _skip(2) `perc_ll' _skip(2) `perc' _skip(2) `perc_ul' if "`perc'" != "" { // Make sure that the loop doesn't break if post `foo' ("`country'") ("`sex'") (`x') (`perc') (`perc_ll') (`perc_ul') } } } } postclose `foo' // Use posted values use `infertage', clear drop if perc_infert == . // Drop empty rows // Define and set saved value labels again do `valuelabels' label val age age // Plot for women by country and age twoway (rarea perc_infert_ll perc_infert_ul age if sex == "Women", fcolor(gs15) cmissing(no)) /// (line perc_infert age if sex == "Women", cmissing(no)) /// , by(cntry, cols(3) legend(off) note("") /// caption("{it:Source:} Generations and Gender Survey Wave 1 (2005{c 0150}2011)" /// "{it:Note:} Gray areas denote 95 % confidence intervals", size(vsmall) span)) /// ytitle("% of {bf:women} reporting (suspected) infertility") xtitle("Age group") /// xlabel(1/5, val alternate) /// ysize(8) name(infertage, replace)

Sep 14, 2014

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 */