// 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-regressionesttab, 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)
Showing posts with label Multilevel modeling. Show all posts
Showing posts with label Multilevel modeling. Show all posts
Aug 2, 2018
Meta-regression
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)
Jul 4, 2017
Chapter 5 of Singer and Willett's (2003) book on longitudinal data analysis
There is another take on the chapter here, but I like mine better.
// Table 5.1 use "C:\singer willett (2003)\reading_pp.dta", clear format age %6.2f list id wave agegrp age piat if inlist(id, 4, 27, 31, 33, 41, 49, 69, 77, 87), sep(0) noobs // Figure 5.1twoway (scatter piat age) /// (lfit piat age) /// (scatter piat agegrp) /// (lfit piat agegrp) /// if inlist(id, 4, 27, 31, 33, 41, 49, 69, 77, 87) /// , by(id, note("")) xtitle("{it:AGE} or {it:AGEGRP}") /// ylabel(0 (20) 80) xlabel(6 (1) 12, format(%6.0f)) ytitle("{it:PIAT}") /// legend(order(1 "Age" 2 "Linear fit age" /// 3 "Target age" 4 "Linear fit target age") col(2)) /// name(figure51, replace) ysize(8) // Table 5.2generate agegrp_65 = agegrp - 6.5 generate age_65 = age - 6.5 capture program drop randomslopetable program define randomslopetable eststo `1' qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance exp([lns1_1_1]_b[_cons]) * /// exp([lns1_1_2]_b[_cons]) qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance end eststo clear mixed piat agegrp_65 || id: agegrp_65, var cov(uns) mle eststo agegrp randomslopetable agegrp mixed piat age_65 || id: age_65, var cov(uns) mle eststo age randomslopetable age esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(agegrp_65 age_65) /// // Align coefficients coeflabels(_cons "Intercept" age_65 "Change") /// // Label coefficients order(_cons age_65) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(2 2 2 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("AGEGRP - 6.5" "AGE - 6.5") /// keep(piat:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.3 use "C:\singer willett (2003)\wages_pp.dta", clear list id exper lnw black hgc uerate if inlist(id, 206, 332, 1028), sep(0) noobs // Table 5.4eststo clear mixed lnw exper || id: exper, var mle cov(uns) eststo modela randomslopetable modela mixed lnw c.exper##c.hgc_9 c.exper##black || id: exper, var mle cov(uns) eststo modelb randomslopetable modelb mixed lnw exper hgc_9 c.exper#black || id: exper, var mle cov(uns) eststo modelc randomslopetable modelc esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// 1.black "BLACK" exper "Change" /// c.exper#c.hgc_9 "Change x (HGC - 9)" /// 1.black#c.exper "Change x BLACK") /// // Label coefficients order(_cons hgc_9 1.black /// exper c.exper#c.hgc_9 1.black#c.exper) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.2estimates restore modelc margins, at(exper = (0 (1) 10) black = (0 1) hgc_9 = (0 3)) marginsplot, ytitle("Predicted log hourly wage") title("") /// recastci(rarea) ciopt(color(gs14)) /// xtitle("{it:EXPER}") /// addplot(scatteri 2.35 10 "White", msymbol(none) mlabpos(0) /// || scatteri 2.23 10 "White", msymbol(none) mlabpos(0) /// || scatteri 2.185 10 "Black", msymbol(none) mlabpos(0) /// || scatteri 2.07 10 "Black", msymbol(none) mlabpos(0) /// || scatteri 1.75 0 "9th grade", msymbol(none) mlabpos(0) /// || scatteri 1.88 0 "12th grade", msymbol(none) mlabpos(0)) /// legend(off) ylabel(, format(%6.1f)) /// name(figure52, replace) // Table 5.5use "C:\singer willett (2003)\wages_small_pp.dta", clear eststo clear mixed lnw hgc_9 exper c.exper#black || id: exper, var mle cov(uns) eststo modela randomslopetable modela mixed lnw hgc_9 exper c.exper#black || id: , var mle cov(uns) eststo modelc qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(agegrp_65 age_65) /// // Align coefficients coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// 1.black "BLACK" exper "Change" /// 1.black#c.exper "Change x BLACK") /// // Label coefficients order(_cons hgc_9 1.black /// exper 1.black#c.exper) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.6 use "C:\singer willett (2003)\unemployment_pp.dta", clear list id months cesd unemp if inlist(id, 7589, 55697, 67641, 65441, 53782), sepby(id) noobs // Table 5.7eststo clear mixed cesd months || id: months, var mle cov(uns) eststo modela randomslopetable modela mixed cesd months i.unemp || id: months, var mle cov(uns) eststo modelb randomslopetable modelb mixed cesd c.months##i.unemp || id: months, var mle cov(uns) eststo modelc randomslopetable modelc version 10 // This one is a challenge to fit generate unempXmonths = unemp * months xtmixed cesd unemp unempXmonths || id: unemp unempXmonths, var mle cov(uns) eststo modeld qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v4 = exp(2*[lns1_1_2]_b[_cons]) // Unemp variance qui estadd scalar v3 = exp(2*[lns1_1_3]_b[_cons]) // Unemp x months variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance version 14 esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(unemp 1.unemp unempXmonths 1.unemp#c.months) /// coeflabels(_cons "Intercept" months "Change" /// 1.unemp "UNEMP" /// 1.unemp#c.months "Change x UNEMP") /// // Label coefficients order(_cons months /// 1.unemp 1.unemp#c.months) /// // Order coefficients stats(v_e v1 v2 v3 v4 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Var(UNEMP)" /// "Var(UNEMP x TIME)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C" "Model D") /// keep(cesd:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.3estimates restore modelb margins, at(months = (0 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("{it:Predicted CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("") /// name(figure53, replace) ** All other plots of Figure 5.3 are only variants of this one ** // Figure 5.4 estimates restore modelb margins, at(months = (0 (1) 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("Model B") subtitle("Main effects of" /// "{it:UNEMP} and {it:TIME}") /// name(figure54A, replace) estimates restore modelc margins, at(months = (0 (1) 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("Model C") subtitle("Interaction between" /// "{it:UNEMP} and {it:TIME}") /// name(figure54B, replace) estimates restore modeld // not sure how to do this with marginsplot predict pd twoway (line pd months if unemp == 0, c(L)) /// (line pd months if unemp == 1, c(L)) /// (scatteri 11.5 15 "Employed", msymbol(none) mlabpos(9)) /// (scatteri 15 15 "Unemployed", msymbol(none) mlabpos(9)) /// , ylabel(5 (5) 20) legend(off) /// xlabel(0 (2) 14) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// title("Model D") subtitle("Constraining the effect of {it:TIME}" /// "among the re-employed") /// name(figure54C, replace) graph combine figure54A figure54B figure54C, row(1) xsize(11) // Table 5.8 use "C:\singer willett (2003)\wages_pp.dta", clear eststo clear mixed lnw hgc_9 ue_7 exper c.exper#black || id: exper, mle cov(un) var eststo modela randomslopetable modela mixed lnw hgc_9 ue_mean ue_person_centered exper c.exper#black || id: exper, mle cov(un) var eststo modelb randomslopetable modelb mixed lnw hgc_9 ue1 ue_centert1 exper c.exper#black || id: exper, mle cov(un) var eststo modelc randomslopetable modelc esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(ue_mean ue_7 ue1 ue_7 ue_centert1 ue_person_centered) /// coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// ue_7 "UERATE" ue_person_centered "Deviation UERATE" /// exper "Change" 1.exper#black "Change x BLACK") /// // Label coefficients order(_cons hgc_9 ue_7 ue_person_centered /// exper c.exper#black) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.9 use "C:\singer willett (2003)\medication_pp.dta", clear list wave day timeofday time time333 time667 in 1/11, noobs sep(0) // Table 5.10 eststo clear mixed pos i.treat##c.time || id: time, mle cov(uns) var eststo modela randomslopetable modela mixed pos i.treat##c.time333 || id: time333, mle cov(uns) var eststo modelb randomslopetable modelb mixed pos i.treat##c.time667 || id: time667, mle cov(uns) var eststo modelc randomslopetable modelc esttab /// , b(2) se(2) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(time333 time time667 time /// 1.treat#c.time333 1.treat#c.time /// 1.treat#c.time667 1.treat#c.time) /// coeflabels(_cons "Intercept" time "Change" /// 1.treat "TREAT" 1.treat#c.time "Change x TREAT") /// // Label coefficients order(_cons 1.treat time 1.treat#c.time) /// // Order coefficients stats(v_e v1 v2 cov dev aic bic nc N, /// Add variance components to table fmt(2 2 2 2 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Cov(Init., Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(pos:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.5estimates restore modela margins, at(time = (0 (1) 7) treat = (1 0)) marginsplot, recastci(rarea) ciopts(color(gs14)) legend(off) /// title("") xtitle("Days") ytitle("Predicted {it:POS}") /// addplot(scatteri 187 7 "Treatment", msymbol(none) mlabpos(11) /// || scatteri 153 7 "Control", msymbol(none) mlabpos(11) /// xlabel(0 (1) 7)) /// name(figure55, replace)
Reference
Singer, Judith D., and John B. Willett. 2003. Applied Longitudinal Data Analysis. Modeling Change and Event Occurrence. Oxford University Press. doi: 10.1093/acprof:oso/9780195152968.001.0001
Jun 21, 2017
Chapter 4 of Singer and Willett's (2003) book on longitudinal data analysis
There is another take on the chapter here, but I like mine better.
use alcohol1_pp.dta, clear // Figure 4.1twoway (scatter alcuse age) (lfit alcuse age) /// if inlist(id, 4, 14, 23, 32, 41, 56, 65, 82), /// by(id, row(2) legend(off) note("")) /// xlabel(13 (1) 17) ytitle({it:ALCUSE}) xtitle({it:AGE}) /// name(figure41, replace) // Figure 4.2preserve bysort id: sample 32, count qui sum peer generate hipeer = (peer >= r(mean)) qui regress alcuse i.id##c.age predict predicted_alcuse label var predicted_alcuse "Predicted {it:ALCUSE}" xtline predicted_alcuse if coa == 0, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("{it:COA} = 0") ylabel(-1 (1) 4) /// name(coa0, replace) nodraw xtline predicted_alcuse if coa == 1, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("{it:COA} = 1") ylabel(-1 (1) 4) /// name(coa1, replace) nodraw xtline predicted_alcuse if hipeer == 0, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("Low {it:PEER}") ylabel(-1 (1) 4) /// name(hipeer0, replace) nodraw xtline predicted_alcuse if hipeer == 1, overlay t(age) i(id) legend(off) /// xtitle({it:AGE}) xlabel(13 (1) 17) /// title("High {it:PEER}") ylabel(-1 (1) 4) /// name(hipeer1, replace) nodraw graph combine coa0 coa1 hipeer0 hipeer1, col(2) ysize(8) name(figure42, replace) restore // Table 4.1generate time = age - 14 // Create program that gets all the random slope parameters for esttab capture program drop randomslopetable program define randomslopetable eststo `1' qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance exp([lns1_1_1]_b[_cons]) * /// exp([lns1_1_2]_b[_cons]) qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance capture drop p1 // For R-squared based on observed and predicted outcomes qui predict p1 // For R-squared based on observed and predicted outcomes qui regress alcuse p1 // For R-squared based on observed and predicted outcomes local orsq = e(r2) // For R-squared based on observed and predicted outcomes estimates drop . // For R-squared based on observed and predicted outcomes estimates restore `1' // For R-squared based on observed and predicted outcomes qui estadd scalar orsq = `orsq' // For R-squared based on observed and predicted outcomes scalar v_e = exp(2*[lnsig_e]_b[_cons]) // For R-squared based on residual variance reduction scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar rsqe = (v_eumm - v_e) / v_eumm // For R-squared based on residual variance reduction of unconditional means model if !missing(`v2_ugm') { qui estadd scalar rsqv2 = (v2_ugm - v2) / v2_ugm // For R-squared based on residual variance reduction of unconditional means model qui estadd scalar rsqv1 = (v1_ugm - v1) / v1_ugm // For R-squared based on residual variance reduction of unconditional means model } end // This was based on: https://www.statalist.org/forums/forum/general-stata-discussion/general/1309801-saving-commands-as-macro-variables-local-and-global-to-make-do-files-shorter eststo clear mixed alcuse || id: , variance mle // Model A, unconditional means model eststo modela qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance scalar v_eumm = exp(2*[lnsig_e]_b[_cons]) // Residual variance for R-squared mixed alcuse time || id: time, cov(un) variance mle // Model B, unconditional growth model randomslopetable modelb scalar v2_ugm = exp(2*[lns1_1_1]_b[_cons]) // Slope variance scalar v1_ugm = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance mixed alcuse i.coa##c.time || id: time, cov(un) variance mle // Model C randomslopetable modelc mixed alcuse i.coa##c.time c.peer##c.time || id: time, cov(un) variance mle // Model D randomslopetable modeld mixed alcuse i.coa c.peer##c.time || id: time, cov(un) variance mle // Model E randomslopetable modele mixed alcuse i.coa c.cpeer##c.time || id: time, cov(un) variance mle // Model F randomslopetable modelf mixed alcuse ccoa c.cpeer##c.time || id: time, cov(un) variance mle // Model G randomslopetable modelg esttab /// , b(3) se(3) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define starts rename(ccoa 1.coa cpeer peer c.cpeer#c.time c.peer#c.time) /// // Align coefficients coeflabels(1.coa "COA" 1.coa#c.time "COA x time" peer "PEER" /// c.peer#c.time "PEER x time" _cons "Intercept" time "Time") /// // Label coefficients order(_cons 1.coa peer ) /// // Order coefficients stats(v_e v1 v2 cov orsq rsqe rsqv1 rsqv2 dev aic bic nc N, /// Add variance components to table fmt(3 3 3 3 3 3 3 3 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Cov(Init., Change)" /// "R-squared obs./pred." /// "R-squared var(Residual)" /// "R-squared var(Initial)" /// "R-squared var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C" "Model D" "Model E" "Model F" "Model G") /// keep(alcuse:) // Drop variance components in weird shapes // ICC of p. 96: estimates restore modela estat icc // Figure 4.3estimates restore modelb qui margins, at(time = (0 1 2)) marginsplot, xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") ylabel(0 .5 1 1.5 2) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Unconditional growth model", span) name(modelb, replace) nodraw estimates restore modelc qui margins, at(time = (0 1 2) coa = (0 1)) marginsplot, xlabel(-1 "13" 0 "14" .5 " " 1 "15" 1.5 " " 2 "16" 3 "17", ) ylabel(0 (1) 2, ) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Uncontrolled effects of {it:COA}", span) legend(off) /// addplot(scatteri .9 2 "{it:COA} = 0" 1.5 2 "{it:COA} = 1", msymbol(none)) /// name(modelc, replace) nodraw estimates restore modele margins, at(time = (0 1 2) coa = (0 1) peer = (.655 1.381)) marginsplot, xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17", format(%6.0f)) ylabel(0 1 2, ) ytitle("Predicted {it:ALCUSE}") /// recastci(rarea) ciopts(color(gs14)) /// title("Controlled effects of {it:COA}", span bexpand) legend(off) /// addplot(scatteri .95 2 "{it:COA} = 0" /// 1.5 2 "{it:COA} = 1" /// .14 -0.85 "Low {it:PEER}" /// .64 -0.85 "High {it:PEER}" /// .72 -0.85 "Low {it:PEER}" /// 1.21 -0.85 "High {it:PEER}", msymbol(none)) /// name(modele, replace) nodraw graph combine modelb modelc modele, ycommon xcommon col(3) xsize(12) ysize(6) altshrink name(figure43, replace) // Figure 4.4preserve statsby, by(id) saving(temp, replace): regress alcuse time merge m:1 id using temp qui cor _b_cons coa local r = round(r(rho), .01) twoway (scatter _b_cons coa) /// (scatteri 4 .5 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(-.5 " " 0 1 1.5 " ") ytitle(Intercept) /// xtitle("{it:COA}") legend(off) name(g1, replace) nodraw qui cor _b_cons peer local r = round(r(rho), .01) twoway (scatter _b_cons peer) /// (scatteri 4 1 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(0 1 2 3) ytitle(Intercept) /// xtitle("{it:PEER}") legend(off) name(g2, replace) nodraw qui cor _b_time coa local r = round(r(rho), .01) twoway (scatter _b_time coa) /// (scatteri 4 .5 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(-.5 " " 0 1 1.5 " ") ytitle(Change) /// xtitle("{it:COA}") legend(off) name(g3, replace) nodraw qui cor _b_time peer local r = round(r(rho), .01) twoway (scatter _b_time peer) /// (scatteri 4 1 "{it:r} = `r'", msymbol(none) mlabpos(0)) /// , xlabel(0 1 2 3) ytitle(Change) /// xtitle("{it:PEER}") legend(off) name(g4, replace) nodraw graph combine g1 g2 g3 g4, name(figure44, replace) restore // Figure 4.5estimates restore modelf predict e, resid predict e_time e_cons, reffects relevel(id) qnorm e, yline(0) ytitle(Var(residual)) name(g1, replace) nodraw qnorm e_cons, yline(0) ytitle(Var(intercept)) name(g2, replace) nodraw qnorm e_time, yline(0) ytitle(Var(change)) name(g3, replace) nodraw egen ze = std(e) egen ze_time = std(e_time) egen ze_cons = std(e_cons) scatter ze id, yline(0) ytitle("Standardized" "var(residual)") name(g4, replace) nodraw scatter ze_cons id, yline(0) ytitle("Standardized" "var(intercept)") name(g5, replace) nodraw scatter ze_time id, yline(0) ytitle("Standardized" "var(change)") name(g6, replace) nodraw graph combine g1 g2 g3 g4 g5 g6, cols(2) colfirst name(figure45, replace) ysize(8) // Figure 4.6twoway (scatter e time), ylabel(-2(1)2) ytitle("Var(residual)") /// xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") /// xtitle("{it:TIME}") /// yline(0) name(g1, replace) nodraw twoway (scatter e_cons coa), ylabel(-2(1)2) ytitle("Var(intercept)") /// xlabel(-.5 " " 0 1 1.5 " ") xtitle("{it:COA}") /// yline(0) name(g2, replace) nodraw twoway (scatter e_cons peer), ylabel(-2(1)2) ytitle("Var(intercept)") /// xlabel(0 1 2 3) xtitle("{it:PEER}") /// yline(0) name(g3, replace) nodraw twoway (scatter e_time coa), ylabel(-2(1)2) ytitle("Var(change)") /// xlabel(-.5 " " 0 1 1.5 " ") xtitle("{it:COA}") /// yline(0) name(g4, replace) nodraw twoway (scatter e_time peer), ylabel(-2(1)2) ytitle("Var(change)") /// xlabel(0 1 2 3) xtitle("{it:PEER}") /// yline(0) name(g5, replace) nodraw graph combine g1 g2 g3 g4 g5, cols(2) hole(2) name(figure46, replace) ysize(8) // Figure 4.7qui reg alcuse time coa cpeer c.cpeer#c.time predict pa estimates restore modelf gen bayes0 = _b[_cons] + _b[1.coa]*coa + _b[cpeer]*cpeer + e_cons gen bayes1 = _b[time] + _b[c.cpeer#c.time]*cpeer + e_time gen bayes = bayes0 + bayes1 * time twoway (scatter alcuse time) /// (lfit alcuse time) /// (line pa time, sort) /// (line bayes time, sort) /// if inlist(id, 4, 14, 23, 32, 41, 56, 65, 82), by(id, col(4) note("")) /// ytitle({it:ALCUSE}) xlabel(-1 "13" 0 "14" 1 "15" 2 "16" 3 "17") /// legend(order(1 "Data points" 2 "OLS" /// 3 "Population average" 4 "Bayes estimate") col(4)) /// xtitle("{it:TIME}") name(figure47, replace)
Reference
Singer, Judith D., and John B. Willett. 2003. Applied Longitudinal Data Analysis. Modeling Change and Event Occurrence. Oxford University Press. doi: 10.1093/acprof:oso/9780195152968.001.0001
May 25, 2017
Simulating multilevel data
clear set obs 200 // Number of level-2 units gen j = _n // ID for level-2 units gen c_j = rnormal(0,1) // Level-2 covariate gen u_j = rnormal(0,1) // Level-2 error term expand 100 // Number of level-1 units per level-2 unit bysort j: gen i = _n // ID for level-1 units gen x_ij = rnormal(0,1) // Level-1 covariate gen e_ij = rnormal(0,1) // Level-1 error term gen y_ij = 1 + 1 * x_ij + 1 * c_j + u_j + e_ij // Regression equation mixed y_ij x_ij c_j
Dec 27, 2016
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
Aug 2, 2016
Random graphs (90): Regression coefficients
version 14
use 7348_F1.dta, clear
// Prepare variables:
// Country variable
rename Y11_Country ctry
decode ctry, gen(country)
kountryadd "Macedonia (FYROM)" to "Macedonia" add
kountry country, from(other) stuck
rename _ISO3N_ geo
kountry geo, from(iso3n) to(iso2c)
rename _ISO2C_ cntry
replace cntry = "XK" if country == "Kosovo"
drop geo
label define Y11_Country 30 "Macedonia", modify
// Wave identifier
rename Wave wave
// Gender
generate female = (Y11_HH2a == 2)
// Work-family conflict
*factor Y11_Q12a Y11_Q12b Y11_Q12c, pcf
*tab Y11_Q12a wave, mis
*tab Y11_Q12b wave, mis
*tab Y11_Q12c wave, mis
*alpha Y11_Q12a Y11_Q12b Y11_Q12c, item
scores wfb = mean(Y11_Q12a Y11_Q12b Y11_Q12c)
generate wfc = 5 - wfb
drop wfb
eststo clear
levelsof wave, local(wave)
foreach x of local wave {
eststo: mixed wfc female || ctry: female if wave == `x', cov(uns)
preserve
// Get no. of countries
matrix groups = e(N_g)
local n_g = groups[1,1]
// Predict residuals
predict u1 u0, reffects
predict u1se u0se, reses
// Calculate posterior slope
generate eb_female = u1 + _b[female]
// Plot posterior slope
egen pickone = tag(country) // Keep one case per countr
keep if pickone
egen order_ = rank(-u1), unique
labmask order_, value(ctry) decode
gen high = u1 + _b[female] + (1.96 * u0se)
gen low = u1 + _b[female] - (1.96 * u0se)
local avg = _b[female]
twoway (rcap eb_female eb_female order_, horizontal) ///
(rspike high low order_, horizontal) , ///
xline(`avg') ylabel(1/`n_g', val ang(h)) ///
ytitle("") ///
xscale(alt) ///
xlabel(-.2 (.1) .4) ///
xtitle("Female WFC disadvantage") ///
legend(off) ///
title(`: label (wave) `x'') ///
name(emp_bayes_`x', replace) xsize(3) nodraw
restore
}
graph combine emp_bayes_1 emp_bayes_2 emp_bayes_3, col(3) xsize(6) name(variation, replace) altshrink ///
note("{it:Source:} European Quality of Life Surveys, 2003-11. {it:Notes:} Horizontal line shows average coefficient, country-specific" ///
"estimates indicate deviation from this average. Error bars are 95% CI's based on random-effects models.")
coefplot est1, bylabel(EQLS 2003) || ///
est2, bylabel(EQLS 2007) || ///
est3, bylabel(EQLS 2011) || ///
, xlabel(0 (.05) .2) drop(_cons) xscale(alt) baselevel coeflabel(female = "Female") ///
xtitle("Female WFC disadvantage") xline(0) byopts(row(1)) ciopts(recast(rcap))
coefplot est1 est2 est3, xlabel(0 (.05) .2) drop(_cons) xscale(alt) baselevel coeflabel(female = "Female") ///
xtitle("Female WFC disadvantage") xline(0) byopts(row(1)) ciopts(recast(rcap)) ///
legend(order(2 "EQLS 2003" 4 "EQLS 2007" 6 "EQLS 2011"))
Apr 12, 2016
Snijders and Bosker's (2012) chapter on multivariate multilevel analysis using -runmlwin-
version 13.1
clear
// Data set for Example 16.1
// mlbook2, data set with pupils having missings on ses or IQ_verb or (langPOST and aritPOST) excluded.
input schoolnr pupilNR_new langPOST aritPOST ses IQ_verb IQ_perf Minority denomina sch_ses sch_iqv sch_min
1.00 1.00 . 12.00 -17.73 -1.37 -3.75 1.00 1.00 -14.035 -1.4039 0.630
1.00 3.00 46.00 24.00 -4.73 3.13 1.25 0.00 1.00 -14.035 -1.4039 0.630
1.00 4.00 45.00 19.00 -17.73 2.63 -1.08 1.00 1.00 -14.035 -1.4039 0.630
1.00 5.00 33.00 24.00 -12.73 -2.37 -0.08 0.00 1.00 -14.035 -1.4039 0.630
1.00 6.00 46.00 26.00 -4.73 -0.87 -1.08 0.00 1.00 -14.035 -1.4039 0.630
1.00 7.00 20.00 9.00 -17.73 -3.87 -4.41 0.00 1.00 -14.035 -1.4039 0.630
1.00 8.00 30.00 13.00 -17.73 -2.37 -2.08 1.00 1.00 -14.035 -1.4039 0.630
// .. Skip a few cases ..
end
// 1. State Mlwin path:
global MLwiN_path C:\Program Files (x86)\MLwiN v2.31\i386\mlwin.exe
// 2. Create constant for models:
generate cons = 1
// 3. Fit models:
// Results of Table 16.1:
runmlwin (langPOST cons, eq(1)) (aritPOST cons, eq(2)), ///
level2(schoolnr: (cons, eq(1)) (cons, eq(2))) ///
level1(pupilNR_new: (cons, eq(1)) (cons, eq(2))) ///
nopause
// Population correlation coefficients:
di [RP2]cov(cons_1\cons_2) / sqrt([RP2]var(cons_1)*[RP2]var(cons_2)) // School-level
di [RP1]cov(cons_1\cons_2) / sqrt([RP1]var(cons_1)*[RP1]var(cons_2)) // Student-level
// This one seems a bit off in the book:
di ([RP2]cov(cons_1\cons_2) + [RP1]cov(cons_1\cons_2)) / sqrt( ([RP2]var(cons_1) + [RP1]var(cons_1)) ///
* ([RP2]var(cons_2) + [RP1]var(cons_2)))
// Correlation between group means for groups of a hypothetical size n = 30
di ([RP2]cov(cons_1\cons_2) + [RP1]cov(cons_1\cons_2) /30) / sqrt( ([RP2]var(cons_1) + [RP1]var(cons_1) /30) ///
* ([RP2]var(cons_2) + [RP1]var(cons_2) /30))
// Results of Table 16.2:
// - Create interaction terms (-runmlwin- doesn't accept factor variables):
generate iqXses = IQ_verb * ses
generate schiqXschses = sch_iqv * sch_ses
// - Fit model:
runmlwin (langPOST cons IQ_verb ses sch_iqv sch_ses iqXses schiqXschses, eq(1)) ///
(aritPOST cons IQ_verb ses sch_iqv sch_ses iqXses schiqXschses, eq(2)), ///
level2(schoolnr: (cons, eq(1)) (cons, eq(2))) ///
level1(pupilNR_new: (cons, eq(1)) (cons, eq(2))) ///
nopause
References
Leckie, George, and Chris Charlton. 2013. "runmlwin. A Program to Run the MLwiN Multilevel Modeling Software from within Stata." Journal of Statistical Software 52(11):1-40.Snijders, Tom, and Roel Boskers. 2012. Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling, 2nd ed. Sage.
Apr 11, 2016
Snijders and Bosker's (2012) chapter on multivariate multilevel analysis using -mixed-
Download the complete do-file including data here.
version 13.1
clear
// Data set for Example 16.1
// mlbook2, data set with pupils having missings on ses or IQ_verb or (langPOST and aritPOST) excluded.
input schoolnr pupilNR_new langPOST aritPOST ses IQ_verb IQ_perf Minority denomina sch_ses sch_iqv sch_min
1.00 1.00 . 12.00 -17.73 -1.37 -3.75 1.00 1.00 -14.035 -1.4039 0.630
// .. Skip a couple of thousand observations
end
// Rename outcome variables for -reshape-:
rename langPOST score1
rename aritPOST score2
// Convert data into long format:
reshape long score, i(pupilNR_new) j(outcome)
// Create dummy variables that identify outcome variable:
quietly tab outcome, gen(outcome)
// Results in Table 16.1:
mixed score outcome1 outcome2, nocons ///
|| schoolnr: outcome1 outcome2, nocons cov(un) ///
|| pupilNR_new: , nocons cov(un) residuals(un, t(outcome))
// Results in Table 16.2:
mixed score outcome1 outcome2 ///
outcome#c.IQ_verb outcome#c.ses ///
outcome#c.sch_iqv outcome#c.sch_ses ///
outcome#c.IQ_verb#c.ses outcome#c.sch_iqv#c.sch_ses, nocons ///
|| schoolnr: outcome1 outcome2, nocons cov(un) ///
|| pupilNR_new:, nocons cov(un) residuals(un, t(outcome))
Reference
Snijders, Tom, and Roel Boskers. 2012. Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling, 2nd ed. Sage.
Mar 21, 2016
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 ...
Mar 19, 2016
Avoding the "initial values not feasible" (r(1400)) error with -melogit-
melogit home time || country: time, cov(uns) oryields an error message "initial values not feasible." The same error occurs with the much simpler:
melogit home || country:, orIn order to further simplify estimation, the number of integration points can be reduced (the default is 7):
melogit home || country:, intpoints(2) orDrawing on the estimates of this model as starting values, the desired model can then be estimated:
melogit home || country:, intpoints(2) or mat a = e(b) melogit home time || country:, intpoints(2) or from(a, skip) mat b = e(b) melogit home time || country: time, cov(uns) intpoints(2) or from(b, skip)
Mar 6, 2016
Shrinkage (partial pooling) in multilevel models
unzipfile "output7717034280080927434.zip", replace
use cntry essround happy using ESS1-6e01_0_F1, clear
// Drastically reduce sample size
bysort cntry essround: keep if _n < 5
// Encode country variable
encode cntry, gen(country)
// Fit random intercept model
mixed happy || country:
predict u0, reffects
predict u0se, reses
// Calculate posterior intercept
generate eb_happy = u0 + _b[_cons]
// Plot posterior intercept
preserve
egen pickone = tag(country)
keep if pickone
egen order_ = rank(-u0), unique
labmask order_, value(country) decode
gen high = u0 + _b[_cons] + (1.96 * u0se)
gen low = u0 + _b[_cons] - (1.96 * u0se)
local avg = _b[_cons]
twoway (rcap eb_happy eb_happy order_, dsymbol(x)) ///
(rspike high low order_) , ///
yline(`avg') xlabel(1/32, val ang(v)) ///
xtitle("") ///
ylabel(4/10) ///
ytitle("Country-level residuals + intercept") ///
legend(off) ///
title("Random intercept model") ///
name(emp_bayes, replace) nodraw
restore
// Fit OLS models
regress happy
local avg = _b[_cons]
regress happy i.country
predict avg_happy
predict se_happy, stdp
// Plot OLS estimates
preserve
egen pickone = tag(country)
keep if pickone
egen order_ = rank(-avg_happy), unique
labmask order_, value(country) decode
gen high = avg_happy + (1.96 * se_happy)
gen low = avg_happy - (1.96 * se_happy)
twoway (rcap avg_happy avg_happy order_, dsymbol(x)) ///
(rspike high low order_) , ///
yline(`avg') xlabel(1/32, val ang(v)) ///
xtitle("") ///
ylabel(4/10) ///
ytitle("Country-level means") ///
legend(off) ///
title("OLS regression model") ///
name(plain, replace) nodraw
restore
graph combine emp_bayes plain, col(1) ysize(8) name(combo, replace)
// Plot EB and OLS means against one another to illustrate // shrinkage in multilevel modeling
preserve
egen pickone = tag(country)
keep if pickone
twoway (scatter eb_happy avg_happy) ///
(function y = x, range(0 10)), ///
legend(off) xtitle("OLS estimates") ytitle("Posterior intercept") ///
title("Shrinkage")
restore
Feb 26, 2016
Generating a country–year variable
unzipfile "output7717034280080927434.zip", replace
use cntry essround happy using ESS1-6e01_0_F1, clear
// Generate country-year variable
levelsof cntry, local(levels1)
levelsof essround, local(levels2)
gen cyear = ""
foreach country of local levels1 {
foreach round of local levels2 {
replace cyear = "`country'" + "`round'" if cntry == "`country'" & essround == `round'
di "`country'" "`round'"
}
}
label var cyear "Country-Year"
// Fit three-level model
mixed happy || cntry: || cyear:
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
Jan 30, 2016
Random graphs (56): Intraclass Correlation Coefficients with confidence intervals
webuse productivity, clear
tempname foo
postfile `foo' str100 commandline_str str100 description_str icc icclb iccub N N_groups using "C:\Windows\Temp\test.dta", replace
qui mixed gsp || region: if year <= 1973
qui estat icc
matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis = r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1970 to 1973"
post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)
qui mixed gsp || region: if year >= 1974 | year <= 1977
qui estat icc
matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis = r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1974 to 1977"
post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)
qui mixed gsp || region: if year >= 1978 | year <= 1981
qui estat icc
matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis = r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1978 to 1981"
post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)
qui mixed gsp || region: if year >= 1981 | year <= 1984
qui estat icc
matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis = r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1981 to 1984"
post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)
qui mixed gsp || region: if year >= 1985 | year <= 1986
qui estat icc
matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis = r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1985 to 1986"
post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)
postclose `foo'
use "C:\Windows\Temp\test.dta", clear
list
encode description_str, gen(description)
twoway (dot icc description, horizontal) ///
(rcap icclb iccub description, horizontal) ///
, ylabel(1/5, val) ytitle("") legend(off) ///
xlabel(0 (.1) 1) xtitle("ICC") name(figure, replace)
erase "C:\Windows\Temp\test.dta"
Oct 6, 2015
Random graphs (54): Cross-level interactions
use deleteme.dta
// Cross-level interaction model
qui mixed lsat c.wfc##c.gdp || country: wfc , cov(uns)
estimates store cross_lvl
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Create table esttab cross_lvl /// , se /// stats(v1 v2 v_e cov dev nc N, /// Add variance components to table labels("Var(Intercept)" /// "Var(Slope)" /// "Var(Residual)" /// "Cov(Int., Slope)" /// "Deviance" /// "No. clusters" /// "No. individuals")) /// coeflabel(wfc "Work-family conflict" /// gdp "GDP" /// c.wfc#c.gdp "Work-family conflict X GDP" /// _cons "Intercept") /// mtitles("Cross-level interaction") /// nonumbers varwidth(28) modelwidth(24) /// keep(lsat:) // Drop variance components in weird shapes // Interaction plot (A) preserve estimates restore cross_lvl qui sum wfc // get standard deviation and mean local wfcminsd = r(mean) - r(sd) local wfcmean = r(mean) local wfcplusd = r(mean) + r(sd) qui sum gdp // 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) margins, at(gdp = (`gdpmin2sd' `gdpminsd' `gdpmean' `gdpplusd' `gdpplu2sd') /// wfc = (`wfcminsd' `wfcmean' `wfcplusd')) vsquish marginsplot, x(gdp) noci /// xlabel(`gdpmin2sd' "-2 SD" /// `gdpminsd' "-1 SD" /// `gdpmean' `""Average" "GDP""' /// `gdpplusd' "+1 SD" /// `gdpplu2sd' "+2 SD") /// xtitle("") /// ylabel(, format(%6.1f)) /// ytitle("Predicted happiness") /// plotopts(msymbol(none)) /// // Turn off markers plot1opts(lpattern(longdash)) /// // Define line types here plot2opts(lpattern(solid)) /// plot3opts(lpattern(shortdash)) /// legend(subtitle("Work{c 150}family conflict" , size(small)) /// order(1 "- 1 SD" 2 "Mean" 3 "+ 1 SD") size(small) /// pos(11) ring(0)) /// title("(A) Interaction plot") /// name(xlvl_plot1, replace) restore // Interaction plot (B) preserve estimates restore cross_lvl predict u1 u0, reffect capture drop pickone egen pickone = tag(country) keep if pickone gen wfceffect = u1 + _b[wfc] + _b[c.wfc#c.gdp] * gdp qui sum gdp // 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)) /// (function y = _b[wfc] + _b[c.wfc#c.gdp] * x, range(gdp)) /// , ytitle("Work{c 150}family conflict coefficient") /// xlabel(`gdpmin2sd' "-2 SD" `gdpminsd' "-1 SD" /// `gdpmean' `""Average" "GDP""' /// `gdpplusd' "+1 SD" `gdpplu2sd' "+2 SD") /// title("(B) Interaction plot") /// legend(order(2 "Work{c 150}family conflict coefficient by GDP") /// pos(7) ring(0) size(small)) /// name(xlvl_plot2, replace) restore // Combine graphs graph combine xlvl_plot1 xlvl_plot2, row(1) ysize(3) xsize(5.5) altshrink
Sep 30, 2015
Random graphs (53): Visualizing multilevel models
use cntry happy mnactic agea gndr using ESS1e06_4.dta, clear
// Prepare variables
recode agea (999 = .), gen(age)
recode gndr (1 = 0 "Male") (2 = 1 "Female") (9 = .), gen(female)
recode mnactic (77 88 99 = .), gen(active)
label var female "Sex (ref. Male)"
label var age "Age in decades"
label var happy "Happiness"
label var active "Labor market status"
label val active mnactic
center happy age // Center outcome to make intercept smaller for plot
replace c_age = c_age / 10
// Listwise deletion
drop if missing(female, c_age, c_happy, active, cntry)
drop happy gndr agea mnactic // Unecessary variables can go
// Fit model
mixed c_happy c_age i.female i.active || cntry: c_age, ml cov(uns)
// Save random part
local var_age = round(exp(_b[lns1_1_1:_cons])^2, .0001)
local var_int = round(exp(_b[lns1_1_2:_cons])^2, .0001)
local covaria = round(tanh(_b[atr1_1_1_2:_cons]) * ///
exp(_b[lns1_1_1:_cons]) * ///
exp(_b[lns1_1_2:_cons]), .0001)
// Plot fixed part
coefplot, xline(0) ///
xtitle(" " "Estimates and 95% CI's") ///
scheme(s1mono) ///
ciopts(recast(rcap)) ///
coeflabels(_cons = "{bf:Intercept}" ///
c_age = "{bf:Age} in decades" ///
8.active = "Housework") /// // Shorten label
headings(c_age = " " ///
1.female = "{bf:Sex} ({it:ref.} Male)" ///
2.active = `""{bf:Labor market status}" "({it:ref.} Paid work)""' ///
_cons = " ") ///
mlabel format(%9.2f) mlabposition(12) mlabgap(*1.5) ///
xscale(alt) msymbol(x) ///
ysize(8) xsize(4) ///
title("Fixed part", span) ///
name(fixed_part, replace) nodraw
/* The default sizes of the available area are -ysize(4)- and -xsize(5.5)-,
by the way. Letter size is -ysize(11)- and -xsize(8.5)-*/
// Plot random part
// Estimate residuals
capture drop u1* u0*
predict u1 u0, reffects
predict u1se u0se, reses
// Plot intercept variation
preserve
egen pickone = tag(cntry)
keep if pickone
egen order_ = rank(-u0), unique
labmask order_, value(cntry)
gen high = u0 + (1.96 * u0se)
gen low = u0 - (1.96 * u0se)
twoway (rcap u0 u0 order_, dsymbol(x)) ///
(rspike high low order_) , ///
yline(0) xlabel(1/22, val ang(v)) ///
title("Intercept variance = `var_int'") ///
xtitle("") ///
ytitle("Random intercept residuals") ///
legend(off) ///
name(rand_int, replace) nodraw
restore
// Plot slope 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/22, val ang(v)) ///
xtitle("") ///
title("Slope variance = `var_age'") ///
ytitle("Random slope residuals") ///
ylabel(-1 (.5) 1) /// // Make y-axis identical to other plot,
legend(off) /// // otherwise x-axis looks different, too
name(rand_slope, replace) nodraw
restore
// Plot intercept-slope covariance
preserve
gen predRandomSlope= (_b[_cons] + u0) + ((_b[c_age] + u1) * c_age)
qui sum c_age
local hi1 = 1*r(sd)
local hi2 = 2*r(sd)
local hi3 = 3*r(sd)
local lo1 = -1*r(sd)
local lo2 = -2*r(sd)
sort cntry c_age
twoway (line predRandomSlope c_age, connect(ascending)), ///
ytitle("Predicted Happiness") ///
xtitle("") ///
title("Intercept{c 150}slope covariance = `covaria'") ///
xlabel(`lo1' "-1 SD" 0 "Average age" `hi1' "+1 SD" `hi2' "+2 SD" `hi3' "+3 SD") ///
ylabel(,format(%6.1f)) ///
name(covariance, replace) nodraw
restore
// Combine Figures
graph combine rand_int rand_slope covariance, col(1) ysize(8) title("Random part") name(random_part, replace) nodraw
graph combine fixed_part random_part, col(2) ysize(8) xsize(8) altshrink title("Random coefficient model")
Sep 1, 2015
Making -mixed- tables via -esttab-
// Fit a couple of models
// Model 1
eststo clear
eststo: xtmixed wellbeing || cluster: , ml var
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Model 2
eststo: xtmixed wellbeing hourscgm || cluster: hourscgm, ml var cov(uns)
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Model 3
eststo: xtmixed wellbeing hourscwc || cluster: hourscwc, ml var cov(uns)
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Model 4
eststo: xtmixed wellbeing sizecgm hourscgm || cluster: hourscgm, ml var cov(uns)
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Model 5
eststo: xtmixed wellbeing sizecgm hourscwc || cluster: hourscwc, ml var cov(uns)
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Create table via -esttab-
esttab est1 est2 est3 est4 est5 ///
, se ///
stats(v1 v2 v_e cov dev nc N, /// Add variance components to table
labels("Var(Intercept)" ///
"Var(Slope)" ///
"Var(Residual)" ///
"Cov(Int., Slope)" ///
"Deviance" ///
"No. clusters" ///
"No. individuals")) ///
label /// Use variable labels
keep(wellbeing:) // Drop variance components in weird shapes
// Give it another go
esttab est1 est2 est4 est3 est5 ///
, se ///
stats(v1 v2 v_e cov dev nc N, ///
labels("Var(Intercept)" ///
"Var(Slope)" ///
"Var(Residual)" ///
"Cov(Int., Slope)" ///
"Deviance" ///
"No. clusters" ///
"No. individuals")) ///
rename(hourscwc hours hourscgm hours) /// Match up both hour variables
varlabels(hours "Work hours (CGM/CWC)" /// Relabel variables
sizecgm "Workgroup size (CGM)" ///
_cons "Constant") ///
wrap /// Wrap variable labels
mgroups("Null model" ///
"Work hours CGM" /// Group models
"Work hours CWC", pattern(1 1 0 1 0) span) ///
eqlabels("") /// Suppress equation label
nomtitle /// Suppress model titles
keep(wellbeing:) // Drop weird variance components
Aug 31, 2015
Modeling variance functions
version 12.1
set more off
// Open National Longitudinal Survey (Young women 14-26 years of age in 1968)
webuse nlswork, clear
// Level 1: annual measurement
// Level 2: young women
// Null model
xtmixed ln_w || id: , var
xtmrho // VPC in random intercept case w/out covariates equals ICC
// Prepare level-1 covariate
center tenure, replace
label var tenure "Job tenure (years)"
// Add level-1 covariate
xtmixed ln_w tenure || id: , var
xtmrho // VPC after adjusting for a covariate
// 1) Constant variance
// Create variables: constant variance function at level 2
gen lev2var = e(var_u1)
label var lev2var "Variance at level 2"
// Create variables: constant variance function at level 1
gen lev1var = e(var_e)
label var lev1var "Variance at level 1"
// Create variable: predicted values
predict fit, xb
// Create variables: constant bounds for between-neighborhood variation
gen lev2cbhi = fit + 1.96*(lev2var^.5)
gen lev2cblo = fit - 1.96*(lev2var^.5)
label var fit "Fixed relationship"
label var lev2cbhi "Upper bound between-person variation"
label var lev2cblo "Lower bound between-person variation"
// Graphs: Heterogeneity as a constant function
qui twoway (line lev2var tenure, sort) ///
(line lev1var tenure, sort) ///
, legend(ring(0) pos(5)) ///
ylabel(.085 (.01) .12,format(%6.3f)) ///
ytitle("Variance") ///
name(variancefunction1, replace)
qui twoway (line fit lev2cbhi lev2cblo tenure, sort) ///
, legend(ring(0)) ///
ylabel(1 (.5) 3,format(%6.1f)) ///
ytitle("Predicted log wage") ///
name(betweenperson1, replace)
graph combine variancefunction1 between1
, col(1) ysize(8) xcommon ///
title("Random intercept model") ///
name(constantvariance, replace)
// Clean up
drop lev2var lev1var fit lev2cbhi lev2cblo
// Add random slope to model
xtmixed ln_w tenure || id: tenure, var cov(uns)
xtmrho // VPC after adjusting for a covariate and random slope
// 2) Quadratic variance
// Getting at the covariance parameter
// and slope variation as
// -xtmrho- hasn't stored it in a better format
local cov = tanh([atr1_1_1_2]_b[_cons])* exp([lns1_1_1]_b[_cons])* exp([lns1_1_2]_b[_cons])
di `cov'
local var_tenure = exp([lns1_1_1]_b[_cons])^2
// In order to get the other variance components
// without -xtmrho-, one needs:
*local var_cons = exp([lns1_1_2]_b[_cons])^2 // Intercept variance
*local var_res = exp([lnsig_e]_b[_cons])^2 // Residual variance
// Create variables: quadratic variance function at level 2
gen lev2var = e(var_u1) ///
+ 2 * `cov' * tenure ///
+ `var_tenure' * tenure * tenure
label var lev2var "Variance function at level 2"
// Create variable: predicted values
predict fit, xb
// Create variables: quadratic bound bounds for between-neighborhood variation
gen lev2cbhi = fit + 1.96*(lev2var^.5)
gen lev2cblo = fit - 1.96*(lev2var^.5)
label var fit "Fixed relationship"
label var lev2cbhi "Upper bound between-person variation"
label var lev2cblo "Lower bound between-person variation"
// Graphs: Heterogeneity as a quadratic function
qui twoway (line lev2var tenure, sort) ///
, legend(ring(0)) ///
ytitle("Variance") ///
ylabel(, format(%6.3f)) ///
legend(ring(0) pos(11)) ///
name(variancefunction2, replace)
qui twoway (line fit lev2cbhi lev2cblo tenure, sort) ///
, legend(ring(0) pos(11)) ///
name(betweenperson2, replace)
qui graph combine variancefunction2 betweenperson2, ///
, col(1) xcommon name(quadraticvariance, replace) ///
title("Random coefficient model")
graph combine constantvariance quadraticvariance, row(1) xsize(12) ysize(8)
// Clean up
drop lev2var fit lev2cbhi lev2cblo
Aug 11, 2014
Random graphs (27): Between and within regression
// Case A
clear
set seed 1
set obs 300
generate e = 0 + (200 - 0) * runiform() // To generate random variates over the
generate x = 0 + (600 - 0) * runiform() // interval [a,b), a+(b-a)*runiform()
generate y = e // Relationship between X and Y random
generate id = 1
replace y = (e + 200) if x > 200 & x <= 400
replace id = 2 if x > 200 & x <= 400
replace y = (e + 400) if x > 400
replace id = 3 if x > 400
bys id: egen ym = mean(y) // Generate cluster means for Y and X
bys id: egen xm = mean(x)
twoway (scatter y x if id == 1, msymbol(oh)) ///
(scatter y x if id == 2, msymbol(dh)) ///
(scatter y x if id == 3, msymbol(sh)) ///
(scatter ym xm, msymbol(S)) ///
(lfit ym xm, lpattern(solid)) ///
(lfit y x if id == 1, lpattern(dash)) ///
(lfit y x if id == 2, lpattern(dash)) ///
(lfit y x if id == 3, lpattern(dash)) ///
, /*legend(order(1 "A" 2 "B" 3 "C") title(Region))*/ legend(off) ///
xtitle("Component 1") ytitle("Component 2") ///
title("Case A:" "Regional correlation," "individual orthogonality") ///
xlabel(none) ylabel(none) name(casea, replace)
// Case B
clear
set seed 1
set obs 4 // Generate four clusters
generate id = _n
generate u_i = 300 if id == 1 // Generate cluster-specific intercepts
replace u_i = 0 if id == 2
replace u_i = 0 if id == 3
replace u_i = -300 if id == 4
expand 100 // Generate units in clusters
bysort id: generate x = -100 + (600 + 100) * runiform()
generate e_ij = rnormal(0,70) // Generate level-1 error term
generate y = x + u_i + e_ij
drop if y < 0 | x < 0
drop if id == 2 & x >= 300
drop if id == 3 & x < 300
bys id: egen ym = mean(y)
bys id: egen xm = mean(x)
twoway (scatter ym xm, msymbol(S)) ///
(scatter y x if id == 1, msymbol(oh)) ///
(scatter y x if id == 2, msymbol(dh)) ///
(scatter y x if id == 3, msymbol(sh)) ///
(scatter y x if id == 4, msymbol(th)) ///
(lfit y x if id == 1, lpattern(dash)) ///
(lfit y x if id == 2, lpattern(dash)) ///
(lfit y x if id == 3, lpattern(dash)) ///
(lfit y x if id == 4, lpattern(dash)) ///
(lfit ym xm, lpattern(solid)) ///
, legend(order(2 "A" 3 "B" 4 "C" 5 "D" 1 "Means") size(small) title(Region:, size(small)) row(1)) ///
xtitle("Component 1") ytitle("Component 2") ///
title("Case B:" "Regional orthogonality," "individual correlation") ///
xlabel(none) ylabel(none) name(caseb, replace)
// Case C
clear
set seed 1
set obs 3
generate id = _n
generate u_i = 300 if id == 1
replace u_i = 0 if id == 2
replace u_i = -300 if id == 3
expand 150
bysort id: generate x = -100 + (600 + 100) * runiform()
generate e_ij = rnormal(0,70) // Generate level-1 error term
generate y = x + u_i + e_ij
drop if y < 0 | x < 0
bys id: egen ym = mean(y)
bys id: egen xm = mean(x)
twoway (scatter y x if id == 1, msymbol(oh)) ///
(scatter y x if id == 2, msymbol(dh)) ///
(scatter y x if id == 3, msymbol(sh)) ///
(scatter ym xm, msymbol(S)) ///
(lfit y x if id == 1, lpattern(dash)) ///
(lfit y x if id == 2, lpattern(dash)) ///
(lfit y x if id == 3, lpattern(dash)) ///
(lfit ym xm, lpattern(solid)) ///
, /*legend(order(1 "A" 2 "B" 3 "C") title(Region))*/ legend(off) ///
xtitle("Component 1") ytitle("Component 2") ///
title("Case C:" "Negative regional correlation," "positive individual correlation") ///
xlabel(none) ylabel(none) ///
name(casec, replace)
grc1leg casea caseb casec, legendfrom(caseb) span pos(6) row(1) name(combined, replace)
Jul 1, 2014
Random graphs (23): Between and within regression
clear
input j i xij x_j yij y_j
// Data from Snijders and Bosker (1999, Table 3.2):
1 1 1 2 5 6
1 2 3 2 7 6
2 1 2 3 4 5
2 2 4 3 6 5
3 1 3 4 3 4
3 2 5 4 5 4
4 1 4 5 2 3
4 2 6 5 4 3
5 1 5 6 1 2
5 2 7 6 3 2
end
// Plot as Figure 3.4 of Snijders ans Bosker (1999):
twoway (scatter yij xij) ///
(lfit yij xij) ///
(lfit y_j x_j) ///
(lfit yij xij if j == 1, lpattern(dash_dot)) ///
(lfit yij xij if j == 2, lpattern(dash_dot)) ///
(lfit yij xij if j == 3, lpattern(dash_dot)) ///
(lfit yij xij if j == 4, lpattern(dash_dot)) ///
(lfit yij xij if j == 5, lpattern(dash_dot)) ///
, legend(label(2 "Total regression") ///
label(3 "Between regression") ///
label(4 "Within regression") ///
order(2 3 4) ///
pos(1) ring(0)) ///
xtitle(X) ytitle(Y) ///
title("Within, between, and total relations") ///
xlabel(none) ylabel(none) ///
name(wbt_legend, replace)
// Lines directly labeled
twoway (scatter yij xij) ///
(lfit yij xij) ///
(lfit y_j x_j) ///
(lfit yij xij if j == 1, lpattern(dash_dot)) ///
(lfit yij xij if j == 2, lpattern(dash_dot)) ///
(lfit yij xij if j == 3, lpattern(dash_dot)) ///
(lfit yij xij if j == 4, lpattern(dash_dot)) ///
(lfit yij xij if j == 5, lpattern(dash_dot)) ///
, legend(off) ///
text(1.8 6.5 "Between") ///
text(3.0 7.5 "Total") ///
text(7.0 3.5 "Within") ///
text(5.0 5.5 "Within") ///
xtitle(X) ytitle(Y) ///
title("Within, between, and total relations") ///
xlabel(none) ylabel(none) ///
name(wbt_nolegend, replace)
graph combine wbt_legend wbt_nolegend, col(1) ysize(8) xcommon ycommon
// ysize(4) and xsize(5.5) are defaults
Reference
Snijders, Tom, and Roel Boskers. 1999. Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling. Sage.
Subscribe to:
Posts (Atom)