Standard deviation of binary variables

clear

twoway (function y = sqrt(x * (1 - x)), range(0 1)), ///
        ytitle("Sample standard deviation") xtitle("Proportion") ///
        xlabel(0 (.25) 1, format(%6.2f)) ///
        name(figure1, replace)

Random graphs (138): Functions

twoway (function y = 1 + 1*x, range(0 5)), ///
        ylabel("") xlabel("") xsize(4) ysize(4) ///
        title(Linear function) name(figure9a, replace)

twoway (function y = exp(x) / (1 + exp(x)), range(-5 5)), ///
        ylabel("") xlabel("") xsize(4) ysize(4) ///
        title(Logistic function) name(figure9b, replace)

graph combine figure9a figure9b, col(2) xsize(8) ysize(4) ///
              altshrink name(figure9, replace)

Random graphs (135): Scatterplot with OLS regression

clear

// Generate data
set seed 1
set obs 50
gen hours = rnormal(3, 1) // Number of hours studied
gen e = rnormal(1,1)
gen questions = 2 + 2*hours + 1*e // Questions answered correctly

// 1) Basic scatterplot
twoway (scatter questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Number of questions answered correctly") ///
        name(figure1, replace) ylabel(0 (5) 15) xsize(4) ysize(4)

// 2) Scatterplot with regression line
twoway (scatter questions hours) ///
       (lfit questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 (5) 15) legend(off) name(figure2, replace) xsize(4) ysize(4)

// 3) Scatterplot with regression line and equation
regress questions hours
local intercept = round(_b[_cons], .1)
local x = round(_b[hours], .1)

twoway (scatter questions hours) ///
       (lfit questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(14 2 "y = `intercept' + `x' + e", size(large)) ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 (5) 15) legend(off) name(figure3, replace) xsize(4) ysize(4)

// 4) Scatterplot with regression line and equation and labels for components
twoway (scatter questions hours) ///
       (lfit questions hours, lpattern(solid) range(0 5)) ///
       (function y = 1.6 + 2.5, range(1 2)) ///
       (function y = 2, range(4.1 6.6) horizontal lpattern(solid) lcolor(red)) ///
       (scatteri 1.6 0 (3) "Intercept", msymbol(o) mlabcolor(red)) ///
       (scatteri 5 2 "Slope", msymbol(i) mlabcolor(red)), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(14 2 "y = `intercept' + `x' x + e", size(large)) ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 1.6 5 10 15) legend(off) name(figure4, replace) xsize(4) ysize(4)
    
graph combine figure1 figure2 figure3 figure4, col(2) xsize(8) ysize(8) altshrink name(figures14, replace)

Random graphs (134): Chi-squared function

twoway (function y = (x^(1  -1)*exp(-1*x/2))/((2^(1.0))*exp(lngamma(1.0))), range(0 20)) ///
       (function y = (x^(1.5-1)*exp(-1*x/2))/((2^(1.5))*exp(lngamma(1.5))), range(0 20)) ///
       (function y = (x^(2.5-1)*exp(-1*x/2))/((2^(2.5))*exp(lngamma(2.5))), range(0 20)) ///
       (function y = (x^(5-1)*exp(-1*x/2))/((2^(5))*exp(lngamma(5))), range(0 20)) ///
       (scatteri .5   0.1 "{it:df} = 2", msymbol(i) mlabpos(3)) ///
       (scatteri .24  1.7 "{it:df} = 3", msymbol(i) mlabpos(3)) ///
       (scatteri .15  5.0 "{it:df} = 5", msymbol(i) mlabpos(3)) ///    
       (scatteri .10 10.0 "{it:df} = 10", msymbol(i) mlabpos(3)) /// 
      , xtitle(Chi-squared) ytitle("") ylabel(, format(%6.1f)) legend(off) name(figu, replace) //ysize(6) xsize(6)

// Also search chidemo

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)

Random graphs (105): Scatterplot

// Open ISSP data
use C_ALPHAN SEX V51 WEIGHT using ZA4350_v2-0-0.dta, clear

// Fix country variable
replace C_ALPHAN = "BE" if C_ALPHAN == "BE-FLA"
replace C_ALPHAN = "DE" if C_ALPHAN == "DE-E"
replace C_ALPHAN = "DE" if C_ALPHAN == "DE-W"
replace C_ALPHAN = "UK" if C_ALPHAN == "GB-GBN"

// Recode job satisfaction
recode V51 (7 = 0 "Completely dissatisfied") ///
           (6 = 1 "Very dissatisfied") ///
           (5 = 2 "Fairly dissatisfied") ///
           (4 = 3 "Neither satisfied nor dissatisfied") ///
           (3 = 4 "Fairly satisfied") ///
           (2 = 5 "Very satisfied") ///
           (1 = 6 "Completely satisfied"), gen(jobsat)
label var jobsat "Job satisfaction"

// Recode sex
gen female = (SEX == 2) if !missing(SEX)
label define female 0 "Male" 1 "Female"
label val female female
label var female "Sex"
drop if missing(female)

// Get rid of old variables
drop SEX V51

// Collapse and reshape data
collapse jobsat [pw = WEIGHT], by(C_ALPHAN female)
reshape wide jobsat, i(C_ALPHAN) j(female)

label var jobsat0 "Job satisfaction: men"
label var jobsat1 "Job satisfaction: women"

// Plot
twoway (scatter jobsat0 jobsat1 if jobsat0 > jobsat1, mlabel(C_ALPHAN) mlabpos(9)) ///
       (scatter jobsat0 jobsat1 if jobsat0 < jobsat1, mlabel(C_ALPHAN) mlabpos(6)) ///
       (function y = x, range(3.5 5.1))  ///
     , xlabel(3.5 (.5) 5.0, format(%6.2f)) ///
       ylabel(3.5 (.5) 5.0, format(%6.2f)) ///
       ytitle("Job satisfaction: men") ///
       xtitle("Job satisfaction: women") ///
       title("Cross-national variation in the gender gap in job satisfaction", span) ///
       note(" " "{it:Source:} International Social Survey Program (ISSP) 2005, doi:10.4232/1.11648", span) ///
       legend(off) name(figure1, replace)

Dropping one group at a time

// Interaction plot excluding one country at a time
foreach n of numlist 1/35 {  /// Countries are numbered from 1 to 35

  preserve
  qui drop if country == `n'      // Drop one country
  local naam: label country `n'   // Save name in a local
  *di "`naam'"

  qui mixed happiness c.workfamconf##c.gdppc || country: workfamconf, cov(uns)
  predict u1 u0, reffect

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

  gen wfceffect = u1 + _b[workfamconf] + _b[c.workfamconf#c.gdppc] * gdppc

  qui sum gdppc // 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) mlabsize(*.8)) ///
         (function y = _b[workfamconf] + _b[c.workfamconf#c.gdppc] * x, range(gdppc)) ///
        , ytitle("WFC coefficient", size(*.6)) ///
          xlabel(`gdpmin2sd' "-2 SD" `gdpminsd' "-1 SD" ///
                 `gdpmean' `"Avg."' ///
                 `gdpplusd' "+1 SD" `gdpplu2sd' "+2 SD", labsize(*.8)) ///
          title("Excluding `naam'") ///
          legend(off) ///
          name(figurewo`n', replace) xsize(3) ysize(2.5) nodraw
  restore
}

graph combine figurewo1  figurewo2 ///
              figurewo3  figurewo4 ///
              figurewo5  figurewo6 ///
              figurewo7  figurewo8 ///
              figurewo9  figurewo10 ///
              figurewo11 figurewo12 ///
              figurewo13 figurewo14 ///
              figurewo15 figurewo16 , col(4) row(4) xsize(12) ysize(10)
// And so on ...

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