Random graphs (48): Interaction plot with overlaid data points

// Simulate data
clear
set seed 1
set obs 500

generate  e  = 0 + (500 - 0) * runiform()     // To generate random variates over the
generate  x1 = 0 + (800 - 0) * runiform()     // interval [a,b), a+(b-a)*runiform()

generate  x2 = round(0 + (1 - 0) * runiform()) // Binary variable
generate  y  = (x1 + x2 + (x1*x2) + e) / 1000

// Calculations for Aiken & West (1991) style plot
regress y c.x1##x2

qui sum x1
local x1_minsd = r(mean) - (2*r(sd))
local x1_mean  = r(mean)
local x1_plusd = r(mean) + (2*r(sd))


// Calculate predicted values for plot
margins, at(c.x1 = (`x1_minsd' `x1_mean' `x1_plusd') ///
            c.x2 = (0 1)) vsquish
// Plot
qui marginsplot, recastci(rarea) ciopts(color(gs10)) ///
   title("Interaction plot with overlaid data points") ///
   ytitle("y") ///  
   ylabel(, format(%6.1f)) ///
   xtitle("") ///
   plotopts(msymbol(none)) ///        // Turn off markers
   plot1opts(lpattern(longdash)) ///  // Define line types here
   plot2opts(lpattern(solid)) ///
   addplot(scatter y x1 ///
         , symbol(o) ///
           xlabel(`x1_minsd' "-2 SD" ///       // The addplot seems to override
                  `x1_mean' "Average x1" ///   // the regular axis label 
                  `x1_plusd' "+2 SD") ///      // command
           legend(subtitle(x2) ///             // And the legend command
                  order(3 "x2 = 0" 4 "x2 = 1" 2 "95 % CI"))) ///    
   legend(pos(5) ring(0)) ///                  // But not completely
   name(plot, replace)

Random graphs (34): Country scatterplots

use ART.dta, replace

// Fix missing values
recode avgemtrans (-1 -3 = .)
recode emhum      (-1 = .)

// Calculate r-squared
qui regress avgemtrans emhum
local r2 = round(e(r2), .01) // Round coeff.

// Create first panel
twoway (scatter avgemtrans emhum, mlab(countrycode) mlabpos(0) msymbol(i)) ///
       (lfit avgemtrans emhum) ///
    , xlabel(, format(%6.1f)) ylabel(, format(%6.1f)) ///
    ytitle("Avg. no. of fresh non-donor embryos" "per fresh embryo transfer cycle") ///
    xtitle("Avg. level disagreement that" "embryo is a human being") ///
       legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(8)) ///
    name(figure2a, replace)

// Fix missing values
recode afford clonehelpinf (-1 = .)

// Listwise deletion to control y-axis range
preserve
keep if !missing(afford, clonehelpinf)

// Calculate r-squared    
regress afford clonehelpinf
local r2 = round(e(r2), .01) // Round coeff.

// Create second panel
twoway (scatter afford clonehelpinf, mlab(countrycode) mlabpos(0) msymbol(i)) ///
       (lfit afford clonehelpinf) ///
    , xlabel(, format(%6.1f)) ylabel(0 (5) 25, format(%6.1f)) ///
    ytitle("Affordability (net cost of a fresh ART cycle" "as % of annual disposable income)") ///
    xtitle("Avg. level disagreement with" "cloning to help infertile couples") ///
       legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(8)) ///
    name(figure2b, replace)
restore

// Create final Figure
graph combine figure2a figure2b, row(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 

Calculate coefficient alpha by country

use "ZA5900_v1-0-0.dta", replace

renvars, lower // Switch variable names to lower case

kountry v4, from(iso3n) to(iso2c)  // Generate country abbreviation variable
encode _ISO2C_, gen(country)

fre v51-v54 // Work-family conflict items

// Factor analysis yields one-dimensional solution:
factor v51-v54, pcf


// Calculate Cronbach's alpha for each country of the sample:

preserve

statsby alpha_ = r(alpha) ///
      , by(country) clear total: /// -total- adds row for total sample
    alpha v51-v54

replace country = 25 if country == .  // Label row for total sample
label define country 25 "All", modify // Label row for total sample    
    
sort alpha_         // Make list 
list, sep(0)

  // Create graph
egen order_ = rank(-alpha_), unique
labmask order_, value(country) decode

twoway (dot alpha_ order_) ///
       , yline(.70) ///
      xlab(1/25, valuelabels ang(v)) xtitle("Country") ///
   ylabel(, format(%6.2f)) ytitle("Cronbach's alpha of work{c 150}family conflict scale") ///
   note("{it: Note:} Horizontal line denotes conventional cut-off value for Cronbach's alpha", span)

restore

Random graphs (20): Bar graph

clear

// Table from Vaupel and Loichinger 2006, 
// "Redistributing Work in Aging Europe,"
// Science 312(5782):1911-1913.
// doi: 10.1126/science.1127487
// Table 1

input str2 country r2005 r2025 rchange h2005 h2025 hchange
DE  1.27  1.47  16  16.28  14.95  -8
DK  0.97  1.12  15  17.46  16.11  -8
FR  1.43  1.69  18  15.09  13.63  -10
IT  1.59  1.86  17  15.19  13.48  -11
NL  1.01  1.20  19  15.31  13.88  -9
UK  1.09  1.19  9  17.32  16.34  -6
US  1.09  0.99  -9  18.71  18.29  -2
end

drop rchange hchange // drop unnecessary columns

graph bar r2005 r2025, over(country) bargap(-30) ///
      ytitle("Ratio nonworkers per worker") ///
      legend(label(1 "2005") label(2 "2025") row(1)) ///
      ylabel(,format(%6.1f)) ///  Format axis labels
      name(r, replace)

graph bar h2005 h2025, over(country) bargap(-30) ///
      ytitle("Hours worked per week per capita") ///
      legend(label(1 "2005") label(2 "2025")) ///
      ylabel(,format(%6.1f)) ///
      name(h, replace)
   
grc1leg r h, xcommon row(1) ///
             note("{it: Source:} Vaupel and Loichinger 2006, p. 1912" ///
                  "{it: Note:} The values for 2025 assume change in the" ///
                  "population pyramid but no change in labor force " ///
                  "participation or effort" "by age and sex.", span)

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