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