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.1
twoway (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.2
generate 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.4
eststo 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.2
estimates 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.5
use "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.7
eststo 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.3
estimates 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.5
estimates 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