Random graphs (137): Logistic regression

clear

// Generate data
set seed 1
set obs 50
gen hours = rnormal(3, 1) // Number of hours studied
gen e = rnormal(1,1)
gen questions = 2 + 2*hours + 1*e // Questions answered correctly
qui sum questions, detail
generate pass = (questions >= r(p75)) // Passing the exam

// 1) Histogram of outcome variable
twoway (histogram pass, discrete percent), ///
        xlabel(0 "[0] Failed" 1 "[1] Passed") xtitle("") ///
        ytitle("Percent of students") xsize(4) ysize(4) name(figure5, replace)

// 2) Scatterplot
twoway (scatter pass hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Exam success") ///
        ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
        xsize(4) ysize(4) name(figure6, replace)    

// 3) Scatterplot with regression line
regress pass hours
local intercept = round(_b[_cons], .01)
local x = round(_b[hours], .01)      
       
twoway (scatter pass hours) ///
       (lfit pass hours, lpattern(solid) range(1 5)), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(.8 2 "y = `intercept' + `x' x + e", size(large)) ///
        ytitle("Exam success") ///
        ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
        xsize(4) ysize(4) name(figure7, replace)

// 4) Logit curve
logit pass hours
predict yhat   

twoway (scatter pass hours) ///
       (line yhat hours, lpattern(solid) sort), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Exam success") ///
        ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
        xsize(4) ysize(4) name(figure8, replace)

graph combine figure5 figure6 figure7 figure8, ///
              col(2) xsize(8) ysize(8) altshrink name(figures58, replace)

Random graphs (136): Linear probability model

clear

// Generate data
set seed 1
set obs 50
gen hours = rnormal(3, 1) // Number of hours studied
gen e = rnormal(1,1)
gen questions = 2 + 2*hours + 1*e // Questions answered correctly
qui sum questions, detail
generate pass = (questions >= r(p75)) // Passing the exam

// 1) Histogram of outcome variable
twoway (histogram pass, discrete percent), ///
        xlabel(0 "[0] Failed" 1 "[1] Passed") xtitle("") ///
        ytitle("Percent of students") xsize(4) ysize(4) name(figure5, replace)

// 2) Scatterplot
twoway (scatter pass hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Exam success") ///
        ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
        xsize(4) ysize(4) name(figure6, replace)    

// 3) Scatterplot with regression line
regress pass hours
local intercept = round(_b[_cons], .01)
local x = round(_b[hours], .01)      
       
twoway (scatter pass hours) ///
       (lfit pass hours, lpattern(solid) range(1 5)), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(.8 2 "y = `intercept' + `x' x + e", size(large)) ///
        ytitle("Exam success") ///
        ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
        xsize(4) ysize(4) name(figure7, replace)

// 4) Heteroskedastic residuals
regress pass hours
predict resid, resid

twoway (scatter resid hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        yline(0) name(figure7a, replace)

graph combine figure5 figure6 figure7 figure7a, ///
              col(2) xsize(8) ysize(8) altshrink name(figures57, replace)    

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

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

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

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

Calculate inflection point of curvilinear relationships

// Generate data
clear
set seed 1
set obs 500
generate  e = rnormal(0,10)     
generate  x = rnormal(0,2)     
generate  y  = (-2*x + 2*(x*x) + e)
regress y c.x##c.x

matrix list e(b) // Find out names of coefficients

// Plot regression coefficients
coefplot, xline(0) ///
xtitle(" " "Regression coefficients and 95% CI's") ///
drop(_cons) scheme(s1mono) ///
coeflabels(x = "x" c.x#c.x = "x{char 178}") ///
ciopts(recast(rcap)) ///
headings(x = "U-shaped effect of x") ///
name(coeffs, replace)

// Calculate bend
local bend = round(-_b[x]/(2*_b[c.x#c.x]), .01)
display "b1 = " round(_b[x], .01) _newline ///
        "b2 = " round(_b[c.x#c.x], .01) _newline ///
        "-b1/(2*b2) = " `bend'

predict yhat  // Generate predicted values

// Plot results
twoway (line yhat x, sort(x)) ///
       (scatter y x) ///
   , legend(off) ///
    xline(`bend') ///
    ytitle("y") ///
    xtitle("x") ///
    name(scatter, replace) ///
    title("Inflection point at `bend'")
    

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

Random graphs (16): Plots with confidence intervals

preserve

reg srh i.country [pweight = weight]

predict means, xb
predict semeans, stdp

gen loci = means - 1.96 * semeans
gen hici = means + 1.96 * semeans

keep if pickone
 
egen order_ = rank(-means), unique
labmask order_, value(country) decode

twoway (dot means order_) (rcap loci hici order_) ///
        , legend(off) xlab(1/29, valuelabels ang(v)) ///
         ylabel(, format(%6.1f)) ///
         ytitle("Self-rated health") ///
         xtitle("") ///
         note(`"Note: SRH ranges from 0 ("poor") to 4 ("excellent")"', span)

restore

preserve

logit srh_bin i.country [pweight = weight]

predict lr_index, xb
predict se_index, stdp
generate p_hat = exp(lr_index)/(1+exp(lr_index)) // Predicted probabilities
 // (following http://www.stata.com/support/faqs/...)
gen lb = lr_index - invnormal(0.975)*se_index
gen ub = lr_index + invnormal(0.975)*se_index
gen plb = exp(lb)/(1+exp(lb))
gen pub = exp(ub)/(1+exp(ub))

replace p_hat = p_hat * 100
replace plb = plb * 100
replace pub = pub * 100
keep if pickone
 
egen order_ = rank(p_hat), unique
labmask order_, value(country) decode

twoway (dot p_hat order_) (rcap plb pub order_) ///
       , legend(off) xlab(1/29, valuelabels ang(v)) ///
         ylabel(, format(%6.0f)) ///
         ytitle("% Poor health") xtitle("")

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.

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)