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)
Showing posts with label xtmixed postestimation. Show all posts
Showing posts with label xtmixed postestimation. Show all posts
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.
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"
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
Jul 22, 2013
Random graphs (15): Caterpillar plots
use http://www.ats.ucla.edu/stat/examples/imm/imm23, clear // Data from Kreft and De Leeuw (1998: 58) xtmixed math || schid:, var predict means1, reffects predict means1se, reses egen pickone = tag(schid) sort means1 egen means1rank = rank(means1) if pickone == 1 serrbar means1 means1se means1rank if pickone == 1, /// scale(1.96) yline(0) xlabel(none) /// xtitle("Schools") ytitle("School effects on math scores") /// title("Null model") /// name(means1, replace) xtmixed math meanses public ratio || schid:, var predict means2, reffects predict means2se, reses serrbar means2 means2se means1rank if pickone == 1, /// scale(1.96) yline(0) xlabel(none) /// xtitle("Schools") ytitle("School effects on math scores") /// title("After controlling for average SES, school type, /// and student–teacher ratio") /// name(means2, replace) graph combine means1 means2, ycommon xcommon col(1)
Reference
Kreft, Ita, and Jan De Leeuw. 1998. Introducing Multilevel Modeling. Sage.
Subscribe to:
Posts (Atom)