// 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 mixed. Show all posts
Showing posts with label mixed. Show all posts
Aug 2, 2018
Meta-regression
Jul 27, 2017
Random graphs (106): Interaction plot with overlaid density
// Open Allbus 2008
use V154 V151 V156 V760 V754 V755 V5 V767 using ZA4602_v1-0-0.dta, clear
// Age
generate age = V154 if V154 != 999
// Sex
generate female = (V151 == 2)
// Migrant
generate migrant = (V156 == 2)
// Interviewer ID
rename V760 id
// Interviewer sex
generate femaleinterviewer = (V754 == 2)
// Interviewer age
generate interviewerage = V755
// Attractiveness rating by interviewer before and after interview
qui factor V5 V767, pcf
predict attractiveness
// Calculate Spearman-Brown for two-item scales according
// to Eisinga et al. (https://doi.org/10.1007/s00038-012-0416-3):
spearman V5 V767
// Fit model with respondents clustered in interviewers
mixed attractiveness c.age##i.female i.migrant i.femaleinterviewer##c.interviewerage || id:
// Calculate margins and plot
qui margins, at(age = (18 (10) 98) female = (0 1))
marginsplot, recast(line) recastci(rarea) ciopt(color(gs14)) ///
plot1opts(lpattern(dash)) ///
legend(ring(0) pos(2)) ///
title("Female beauty premium disappears with age", span) ///
ytitle("Predicted attractiveness", axis(1)) ///
xtitle("Age") ///
xlabel(20 (10) 100) ///
addplot(histogram age, discrete yaxis(2) ///
lcolor(white) ///
ylabel(0 0.01 0.02, format(%6.2f) axis(2)) ///
yscale(alt range(0 0.1) axis(2)) ///
ytitle("Age density", axis(2)) ///
legend(order(3 "Men" 4 "Women"))) ///
note(" " "{it:Source:} German General Social Survey Allbus 2008, doi:10.4232/1.12345", span) ///
name(figure1, 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
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
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 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 30, 2016
Random graphs (69): ICC's with confidence intervals
use sharew1_rel2-6-0_gv_isced.dta, clear
// Reshape data into long format
reshape long iscedy_c, i(mergeid) j(child)
// Fix variable of interest
recode iscedy_c (-7 = .a "not yet coded (temporary)") ///
(-2 = .b "refusal") ///
(-1 = .c "don't know") ///
(95 = .d "still in school") ///
(97 = .e "other") ///
( . = .f "missing") ///
, gen(years)
label var years "Years of education"
// Israel doesn't provide a ISCED-to-years conversion, thus it's dropped here
drop if country == 25
// Fix country variable
decode country, gen(cntry)
// Calculate ICC's per country
tempname foo
postfile `foo' str100 commandline_str str20 cntry icc icclb iccub N N_groups using "C:\Windows\Temp\test.dta", replace
levelsof cntry, local(country)
foreach x of local country {
*di "`x'"
qui mixed years || mergeid: if cntry == "`x'", reml
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]
post `foo' (e(cmdline)) ("`x'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)
matrix stuff1 = (r(icc2), icclb, iccub, e(N), n_g)
if "`x'" == "Austria" matrix table1 = stuff1
else matrix table1 = (table1\stuff1)
}
postclose `foo'
matrix rownames table1 = `country'
matrix colnames table1 = "ICC" "CI lower" "CI upper" "Children" "Families"
esttab matrix(table1, fmt(2 2 2 0 0))
use "C:\Windows\Temp\test.dta", clear
egen order_ = rank(-icc), unique
labmask order_, value(cntry)
twoway (rcap icc icc order_, horizontal) ///
(rspike icclb iccub order_, horizontal) ///
, ylabel(1/11, val) legend(off) ytitle("") xscale(alt) ///
xtitle("Sibling correlations in educational attainment") ///
note(" " "{it:Source:} SHARE wave 1, doi:10.6103/SHARE.w1.260", span)
erase "C:\Windows\Temp\test.dta"
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 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 27, 2015
Preacher et al. (2008): Latent Growth Curve Modeling in Stata
The website accompanying Preacher et al.'s (2008) book on latent growth curve modeling provide syntax files for Lisrel, Mplus, and Mx, but not for Stata. These are the Stata commands for replicating the first four models presented in the book.
version 13 // Chapter 2 clear // Read in data from clsn_cov.dat ssd init clsn1 clsn3 clsn4 clsn5 clsn6 sex ssd set observations 851 ssd set means 37.9542 37.2785 37.0463 36.5696 36.1363 0.49 #delimit ; ssd set cov 6.3944 \ 3.2716 7.5282 \ 4.1435 6.0804 10.7290 \ 3.7058 5.1597 6.5672 10.2920 \ 4.1286 5.7608 7.2365 7.6463 12.9085 \ -0.0940 -0.0390 -0.1521 -0.1104 -0.1469 0.2502 ; #delimit cr // Table 2.2 ssd list // Model 0: The null model // Only mean intercept and residual variance are estimated
sem (Intercept@1 -> clsn1, ) /// // Set loadings to 1
(Intercept@1 -> clsn3, ) /// // for estimating constant
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain means to be
(clsn3 <- _cons@a, ) /// // the same
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
, latent(Intercept ) ///
cov(Intercept@0 /// Constrain intercept variance to 0
e.clsn1@b e.clsn3@b /// // Constrain residual variances
e.clsn4@b e.clsn5@b /// // to be the same
e.clsn6@b) ///
nocapslatent
estimates store m0
estat gof, stats(chi2 rmsea indices residuals)
// Non-normed fit index (NNFI) is called
// Tucker-Lewis index (TLI) in Stata
// Model 1: Random intercept model (Table 2.3)
// Only mean intercept, intercept variance, and residual variance are
// estimated
sem (Intercept@1 -> clsn1, ) /// // Set loadings to 1
(Intercept@1 -> clsn3, ) /// // for estimating constant
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain means to be
(clsn3 <- _cons@a, ) /// // the same
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
, latent(Intercept ) ///
cov(e.clsn1@b e.clsn3@b /// // Constrain residual variances
e.clsn4@b e.clsn5@b /// // to be the same
e.clsn6@b) ///
nocapslatent
estimates store m1
di 5.27 / (5.27 + 4.68) // Random intercept model allows calculating an ICC
estat gof, stats(chi2 rmsea indices residuals)
// Poor model fit according to all tests
// Likelihood ratio test:
lrtest m1 m0 // Massive improvement in fit for Model 1, though
// Model 2: Fixed intercept, fixed slope model (Table 2.4)
// Only mean intercept, intercept variance, and residual variance are
// estimated
sem (Intercept@1 -> clsn1, ) /// // Constrain paths to be 1
(Intercept@1 -> clsn3, ) ///
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain intercepts to be the same
(clsn3 <- _cons@a, ) ///
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
(Slope@0 -> clsn1, ) /// // Determine temporal structure
(Slope@2 -> clsn3, ) ///
(Slope@3 -> clsn4, ) ///
(Slope@4 -> clsn5, ) ///
(Slope@5 -> clsn6, ) ///
, covstruct(_lexogenous, diagonal) ///
latent(Intercept Slope ) ///
cov(Intercept@0 /// // Set intercept variance to 0
Slope@0 /// // Set slope variance to 0
e.clsn1@b ///
e.clsn3@b ///
e.clsn4@b ///
e.clsn5@b ///
e.clsn6@b) ///
means(Slope) /// // Estimate slope
nocapslatent
estimates store m2
estat gof, stats(chi2 rmsea indices residuals)
// Poor fit, even worse than Model 1
// Model 3: Random intercept, fixed slope (Table 2.5)
sem (Intercept@1 -> clsn1, ) /// // Constrain paths to be 1
(Intercept@1 -> clsn3, ) ///
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain intercepts to be the same
(clsn3 <- _cons@a, ) ///
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
(Slope@0 -> clsn1, ) /// // Determine temporal structure
(Slope@2 -> clsn3, ) ///
(Slope@3 -> clsn4, ) ///
(Slope@4 -> clsn5, ) ///
(Slope@5 -> clsn6, ) ///
, covstruct(_lexogenous, diagonal) ///
latent(Intercept Slope ) ///
cov(Slope@0 /// // Set slope variance to 0
e.clsn1@b ///
e.clsn3@b ///
e.clsn4@b ///
e.clsn5@b ///
e.clsn6@b) ///
means(Slope) /// // Estimate slope
nocapslatent
estimates store m3
estat gof, stats(chi2 rmsea indices residuals)
// Likelihood ratio test:
lrtest m3 m2 // Improvement in fit for Model 3
// Model 4: Random intercept, random slope (Table 2.6)
sem (Intercept@1 -> clsn1, ) /// // Constrain paths to be 1
(Intercept@1 -> clsn3, ) ///
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain intercepts to be the same
(clsn3 <- _cons@a, ) ///
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
(Slope@0 -> clsn1, ) /// // Determine temporal structure
(Slope@2 -> clsn3, ) ///
(Slope@3 -> clsn4, ) ///
(Slope@4 -> clsn5, ) ///
(Slope@5 -> clsn6, ) ///
, covstruct(_lexogenous, diagonal) ///
latent(Intercept Slope ) ///
cov(Intercept*Slope /// // Include intercept-slope covariance
e.clsn1@b ///
e.clsn3@b ///
e.clsn4@b ///
e.clsn5@b ///
e.clsn6@b) ///
means(Slope) /// // Estimate slope
nocapslatent
estimates store m4
estat gof, stats(chi2 rmsea indices residuals)
// Likelihood ratio test:
lrtest m4 m3 // Improvement in fit for Model 4
// Model 4 estimated as a multilevel model (Table 4.1):
// Read in data from clsn_cov.dat, this time using -corr2data-
#delimit ;
matrix input C = (
6.3944,
3.2716, 7.5282,
4.1435, 6.0804, 10.7290,
3.7058, 5.1597, 6.5672, 10.2920,
4.1286, 5.7608, 7.2365, 7.6463, 12.9085,
-0.0940, -0.0390, -0.1521, -0.1104, -0.1469, 0.2502
)
;
#delimit cr
corr2data clsn1 clsn3 clsn4 clsn5 clsn6 sex, n(851) ///
means(37.9542 37.2785 37.0463 36.5696 36.1363 0.49) ///
cov(C) cstorage(lower) clear
correlate, cov // Seems to work
gen id = _n // Create individual identifier
reshape long clsn, i(id) j(grade) // Convert to long format
mixed clsn grade || id: grade, var cov(uns)
Reference
Preacher, Kristopher J., Aaron L. Wichman, Robert C. MacCallum, and Nancy E. Briggs. 2008. Latent Growth Curve Modeling. Sage. doi: 10.4135/9781412984737
Aug 24, 2015
Obtaining estimates from non-converging models
// Create model that doesn't converge use http://www.stata-press.com/data/r10/pig.dta, clear gen week2 = week * week gen week3 = week * week * week gen week4 = week * week * week * week mixed weight week week2 week3 week4 || id: week week2 week3 week4, var cov(uns) // Gives up after 389 iterations // Stop one interation before to get last estimates mixed weight week week2 week3 week4 || id: week week2 week3 week4, var cov(uns) iterate(388)
Oct 5, 2012
Random graphs (4): Plotting intercept and slope variation
use V3 V51 ISCO88 using "C:\Users\User\work\data sets\issp - work orientations\work orientations iii (2005)\ZA4350_F1.dta", clear // Generate variables gen jobsatisf = 7 - V51 iskoisei isei, isko(ISCO88) drop V51 ISCO88 preserve// Calculating intercept and slope variation using OLS statsby inter = _b[_cons] /// slope = _b[isei] /// , by(V3) /// saving(ols, replace): /// regress jobsatisf isei merge m:1 V3 using ols drop _merge // Visualizing intercept and slope variation gen yhat_ols = inter + slope*isei separate jobsatisf, by(V3) separate yhat_ols, by(V3) twoway (line yhat_ols1-yhat_ols43 isei, sort(V3 isei)) /// (lfit jobsatisf isei, clwidth(vvthick) clcolor(black)) /// , legend(off) ytitle("Job satisfaction") xtitle("ISEI") /// xlabel(16 25 50 75 90) ylabel(,format(%6.1f)) /// caption("{it:Source:} ISSP 2005 (Work Orientations III), own calculations", span) /// name(one, replace) restore// Calculating intercept and slope variation using a random effects model center isei, inplace mixed jobsatisf isei || V3: isei, var predict u1 u0, reffects // Visualizing intercept and slope variation gen predRandomSlope = (_b[_cons] + u0) + ((_b[isei] + u1) * isei) twoway (line predRandomSlope isei, connect(ascending) sort(V3 isei)), /// ytitle("Job satisfaction") /// xtitle("ISEI (centered)") /// xlabel(-25 0 25 50) /// ylabel(,format(%6.1f)) /// caption("{it:Source:} ISSP 2005 (Work Orientations III), own calculations", span) /// name(two, replace)
Subscribe to:
Posts (Atom)