Random graphs (139): Forest plot

clear

input str50 study year r n
"Duncan et al." 1972 .573 14347
"Chamberlain & Griliches" 1977 .508 584
"Hauser et al." 1982 
"Hauser"    1984 .404  518
"Hauser & Mossel"    1985 .404  518
"Hauser & Sewell"    1986 .459  532
"Hauser & Wong"    1989 .57  164
"Dronkers"  1993 .6 244
"Scarr & Weinberg" 1994 .32 105
"Hauser et al." 1999 .421 1100
"Ashenfelter & Zimmerman" 1997 .51 143
"Sieben & de Graaf" 2001 .515  2200
"Sieben & de Graaf" 2003 .342 152618
"Conley & Glauber"    2007 .576 1777
"Beenstock" 2008 .435 6926
"Conley"    2008 .578 1777
"Conley & Glauber"    2008 .627 1127
"Mazumder"  2008 .602 6097
"Olneck" 1977 .549 346
"Lee"       2009 .478 4561
"Ermisch & Pronzanto" 2010 .37 68957
"Bjorklund & Salvanes" 2011 .39 34476
"Lindahl" 2011 .411 9864
"Mazumder"  2011 .665 671
"Adermon" 2012 .424 338222
"Schnitzlein" 2014 .656 1480
"Andrade" 2016 .122 1632366
"Bredtman & Smith" 2016 .327 3087
"Pfeffer et al." 2016 .46 5215 
end

generate z = atanh(r) // r-to-z = inverse hyperbolic tangent
generate sez = sqrt(1/(n - 3))

admetan z sez

display _newline ///
" Pooled estimate of r = " tanh(r(eff)) _newline ///
"Lower Limit of 95% CI = " tanh(r(eff) - (1.96 * r(se_eff))) _newline ///
"Upper Limit of 95% CI = " tanh(r(eff) + (1.96 * r(se_eff)))

// Prepare data for -forestplot-
generate _USE = 1

  // Generate CI's for r
generate lb = tanh(_LCI)
generate ub = tanh(_UCI)

  // Generate study labels for -forestplot-
generate _LABELS = study + " (" + string(year, "%02.0f") + ")"

label var n "Sample size"

  // Add effect size to data set
local new = _N + 1
set obs `new'
replace _LABELS = "{bf:Overall}"                    if _n == _N
replace r       = tanh(r(eff))                      if _LABELS == "{bf:Overall}" 
replace lb      = tanh(r(eff) - (1.96 * r(se_eff))) if _LABELS == "{bf:Overall}" 
replace ub      = tanh(r(eff) + (1.96 * r(se_eff))) if _LABELS == "{bf:Overall}" 
replace _USE    = 5                                 if _LABELS == "{bf:Overall}" 

  // Forest plot
forestplot r lb ub, effect("Correlation") rcol(n) leftjustify nowt nonull xlabel(0 .1 .2 .3 .4 .5 .6 .7)

Meta-regression

// Open ISSP 2005
use COUNTRY WRKHRS V51 using "ZA4350_v2-0-0.dta", clear

// Prepare variables
rename COUNTRY country
  // Job satisfaction
generate jobsatis = 7 - V51
rename WRKHRS workhrs
  // Average working hours per country
bys country: egen avghrs = mean(workhrs)
  // Center working hours
center workhrs, inplace

// Fit models
eststo clear
eststo: regress jobsatis workhrs avghrs                  // OLS regression
eststo:   mixed jobsatis workhrs avghrs || country:      // RE regression
 qui matrix foo = e(N_g)          // Number of individuals
 qui estadd scalar nc  = foo[1,1] // Number of individuals

 // Prepare data for meta-regression
preserve  // Generate country-level data set of average working hours
 egen pickone = tag(country)
 keep if pickone
 keep country avghrs 
 tempfile temp
 save `temp', replace
restore

 // Average job satisfaction per country, controlling for individual working hours
statsby _b[_cons] _se[_cons] e(N), clear by(country): regress jobsatis workhrs
 // Merge with average working hours
merge 1:1 country using `temp'

eststo: regress _stat_1  avghrs                // OLS regression two-step
eststo: vwls _stat_1  avghrs, sd(_stat_2)      // FE meta-regression
eststo: metareg _stat_1  avghrs, wsse(_stat_2) // RE meta-regression



esttab, keep(main:) mtitles("OLS" "Multilevel" "OLS two-step" "FE meta-regression" "RE meta-regression") b(3) se(3) ///
        coeflabel(workhrs "Working hours" ///
                  avghrs "Average working hours" ///
                  _cons "Intercept") ///
         stats(nc N, labels("No. countries") fmt(%12.0gc %12.0gc)) varwidth(25) modelwidth(20) eqlabels("", none)