Making tables for logit models using -esttab-

// Open Allbus 2016
use "ZA5250_v2-0-0.dta", clear

// Prepare variables
recode hs01 (4 5 = 1 "Poor health") (3 2 1 = 0 "Good health") (-9 = .), gen(poorhealth)
label var poorhealth "Poor self-rated health"
recode sex  (2 = 1 "  Female") (1 = 0 "  Male"), gen(female)
label var female "Female sex"
recode age (-32 = .)
label var age "Age"
recode isced97 (1 2 = 0 "  Low") (3 4 = 1 "  Medium") (5 6 = 2 "  High") (-32 = .), gen(education)
label var education "Education"

// Fit model 1
eststo clear
eststo: qui logit poorhealth i.female i.education age
estadd expb
qui sum poorhealth
estadd scalar avg = r(mean) * 100

// Long table

esttab, cells(b(star fmt(2) label("B")) ///
              se(par fmt(2) label("(SE B)")) ///
              expb(par([ ]) label("[OR]"))) ///
        stats(N avg chi2 df_m, fmt(%8.0gc %8.1f %8.1gc 0) ///
                     label("Observations" ///
                           "% poor health" ///
                           "Chi-squared" ///
                           "df")) ///
        label varwidth(30) modelwidth(25) nonumber nomtitle varlabel(_cons "Intercept") ///
        eqlabel(" ") ///
        nobaselevel ///
        refcat(1.female "Sex (ref. male)" 1.education "Education (ref. low)", nol) ///
        addnote("* p<0.05, ** p<0.01, *** p<0.001") ///
        title("Poor self-rated health regressed on sex, education, and age. Logistic model")

// Fit models 2
eststo clear
eststo: qui logit poorhealth  i.education age if female == 0
qui estadd expb
qui sum poorhealth if female == 0
qui estadd scalar avg = r(mean) * 100
eststo: qui logit poorhealth age i.education if female == 1
qui estadd expb
qui sum poorhealth if female == 1
qui estadd scalar avg = r(mean) * 100
  
// Wide table

esttab, cell("b(fmt(2) label(B)) se(fmt(2) label(SE B)) expb(fmt(2) label(OR) star)") ///
        stats(N avg chi2 df_m, fmt(%8.0gc %8.1f %8.1gc 0) ///
                     label("Observations" ///
                           "% poor health" ///
                           "Chi-squared" ///
                           "df")) ///
        label varwidth(30) modelwidth(8) nonumber mtitle("Men" "Women") varlabel(_cons "Intercept") ///
        eqlabel(" ") nobaselevels ///
        refcat(1.education "Education (ref. low)", nol) ///
        addnote("* p<0.05, ** p<0.01, *** p<0.001") ///
        title("Poor self-rated health regressed on education and age, stratified by sex. Logistic model")

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 (119): Line plots and dot plots

// Prepare ESS round 8
use cntry essround wrkctra dweight using "ESS8e01.dta", clear
recode wrkctra (6 = .a) (7 = .b) (8 = .c) (9 = .d)

// Add ESS rounds 1-7
append using "ESS1-7e01.dta", keep(cntry essround wrkctr wrkctra dweight)

  // Prepare variables
generate nocontract = (wrkctra == 3) if !missing(wrkctra)
kountry cntry, from(iso2c)
rename NAMES_STD country

  // Prepare files for post commands
tempname foo1
tempname foo2
postfile `foo1' str20 country essround ll nocontract ul using `foo2', replace

qui levelsof country, local(country)

 // Loop
foreach x of local country {
 foreach i of numlist 1/8 {
  capture logit nocontract [pw = dweight] if country == "`x'" & essround == `i'
  if _rc == 0 {
   qui margins [pw = dweight]
   matrix prevs = r(table)
   local prev = prevs[1,1] * 100
   local ll = prevs[5,1] * 100
   local ul = prevs[6,1] * 100
   *di "`x'" _skip(5) `i' _skip(5)  `ll' _skip(5) `prev' _skip(5) `ul'
   post `foo1' ("`x'") (`i') (`ll') (`prev') (`ul')
  }
 }
}
postclose `foo1'

// Plot estimates
use `foo2', clear

  // Fix value label
label define essround 1 "2002" 2 "2004" 3 "2006" 4 "2008" ///
                      5 "2010" 6 "2012" 7 "2014" 8 "2016", modify
label val essround essround

// First plot
sort country essround
twoway (rarea ll ul essround, lcolor(white)) ///
       (connected nocontract essround), ///
        by(country, ///
           note("") ///
           legend(off)) ///
        xlabel(1/8, val ang(90)) ///
        xtitle("") ytitle("% without job contract") ///
        name(figure2a, replace) ysize(9)

// Second plot
reshape wide ll ul nocontract, i(country) j(essround)    

scores average =  mean(nocontract*)
egen order_ = rank(average), unique
labmask order_, value(country)
    
twoway (dot nocontract1 order_, horizontal) ///
       (dot nocontract2 order_, horizontal) ///
       (dot nocontract3 order_, horizontal) ///
       (dot nocontract4 order_, horizontal) ///
       (dot nocontract5 order_, horizontal) ///
       (dot nocontract6 order_, horizontal) ///
       (dot nocontract7 order_, horizontal) ///
       (dot nocontract8 order_, horizontal) ///
       (rspike ul1 ll1 order_, horizontal) ///
       (rspike ul2 ll2 order_, horizontal) ///
       (rspike ul3 ll3 order_, horizontal) ///
       (rspike ul4 ll4 order_, horizontal) ///
       (rspike ul5 ll5 order_, horizontal) ///
       (rspike ul6 ll6 order_, horizontal) ///
       (rspike ul7 ll7 order_, horizontal) ///
       (rspike ul8 ll8 order_, horizontal), ///
        ylabel(1/32, val) ytitle("") xscale(alt) ///
        xtitle("% without job contract") ///
        legend(order(1 "2002" 2 "2004" 3 "2006" 4 "2008" ///
                     5 "2010" 6 "2012" 7 "2014" 8 "2016") ///
               pos(5) ring(0)) ///
        name(figure2b, replace) ysize(9)
   
graph combine figure2a figure2b, ///
      col(2) name(figure2, replace) ///
      note(" " "{it:Source:} European Social Survey 2002-16, weighted data. {it:Note:} Error bands/spikes denote 95% confidence intervals", span  size(*.7)) 

Random graphs (115): Bar graphs

// Data
use ESS7e02_1.dta, clear

// Drinking variable
recode alcfreq (7 = 0 "Never") (6 = 1 "Less than once a month") ///
               (5 = 2 "Once a month") (4 = 3 "2{c 150}3 times a month") ///
               (3 = 4 "Once a week") (2 = 5 "Several times a week") ///
               (1 = 6 "Every day") (77 88 99 = .), gen(alcohol)
label var alcohol "Alcohol consumption"

// Dichotomize alcohol variable
generate drinking = inlist(alcohol, 5, 6) if !missing(alcohol)      

// Country name variable
kountry cntry, from(iso2c) marker
rename NAMES_STD country
  
// Education
recode eisced (1 2   = 0 "Lower") ///
              (3 4 5 = 1 "Medium") ///
     (6 7   = 2 "High") ///
     (55/99 = .), gen(education)
label var education "Education"

// Alcohol consumption by country
histogram alcohol, percent disc horizontal ///
                   by(country, title("{bf:A} Frequency of drinking alcohohl", ///
                         justification(left) bexpand span) ///
                               note("")) ///
                   ylabel(0(1)6, val) ytitle("") ///
                   ysize(8) name(figure1, replace)

// Alcohol consumption by education by country
  // Collect predicted probabilities
capture program drop my_logit
program define my_logit, eclass
    syntax[if]
    marksample touse
    logit drinking if `touse'
    margins if `touse', post
    exit
end

statsby point = _b[_cons] se = _se[_cons], by(country education) clear: my_logit

 // Calculate stuff
replace point = point * 100
replace se = se * 100

gen lb = point - 1.96 * se
gen ub = point + 1.96 * se

  // Plot
twoway (bar point point education) ///
       (rspike lb ub education), ///
    by(country, legend(off) ///
                   title("{bf:B} Drinking by educational attainment", ///
                         justification(left) bexpand span) ///
                   note("")) ///
       xlabel(0 1 2, val ang(v)) ylabel(0 (10) 50) ///
       ytitle("% drinking more than once a week") ///
    ysize(8) ///
    name(figure2, replace)

graph combine figure1 figure2, col(2) ///
                               note(" " "{it:Source:} European Social Survey 2014", ///
                                    span justification(right) bexpand size(*.8))

Analyzing natural policy experiments

Hu et al. 2017 simulate data of a natural policy experiment and show how to analyze the data with regression adjustment, propensity score matching, difference-in-differences, and fixed effects regression. (The paper also includes IV, regression discontinuity, and interrupted time series, but does not describe the data created for these analyses.)

clear
set seed 2

// Create data set with experimental conditions
input str4 educ str6 sex str9 treatmentstr str4 health1 number
      Low  Male   Exposed   Poor   333
      Low  Male   Exposed   Good   917
      Low  Male   Unexposed Poor  1000
      Low  Male   Unexposed Good  2750   
      Low  Female Exposed   Poor   500
      Low  Female Exposed   Good  3250  
      Low  Female Unexposed Poor   167 
      Low  Female Unexposed Good  1083
      High Male   Exposed   Poor    83
      High Male   Exposed   Good   542
      High Male   Unexposed Poor   584
      High Male   Unexposed Good  3791
      High Female Exposed   Poor   125 
      High Female Exposed   Good  1750
      High Female Unexposed Poor   208
      High Female Unexposed Good  2917
end

expand number  // Create full number of cases

// Transform strings to numerical variables
generate loeduc    = (educ         == "Low")
generate female    = (sex          == "Female")
generate treatment = (treatmentstr == "Exposed")
generate good1     = (health1      == "Good")
label define loeduc 0 "High" 1 "Low"
label val loeduc loeduc
label define female 0 "Male" 1 "Female"
label val female female
label define treatment 0 "Unexposed" 1 "Exposed"
label val treatment treatment
label define health 0 "Poor" 1 "Good"
label val good1 health

// Simulate outcome variable
generate good2 = good1
replace  good2 =       1 if good1 == 0 & loeduc == 1 & (runiform() <= .05)
replace  good2 =       1 if good1 == 0 & loeduc == 0 & (runiform() <= .20)
replace  good2 =       1 if good1 == 0 & treatment == 1 & (runiform() <= .30)
label val good2 health

// Transform some more and clean up
generate poor2 = (good2 == 0)
generate poor1 = (good1  == 0)
drop educ sex treatmentstr health1 number good1 good2

// Table 1-ish
table poor2, by(loeduc female treatment) contents(freq) 

// 1) Regression adjustment
logit poor2 treatment female if loeduc == 0, or
logit poor2 treatment female if loeduc == 1, or
logit poor2 treatment##loeduc female##loeduc, or

// 2) Propensity score matching
teffects nnmatch (poor2 female) (treatment) if loeduc == 0
teffects nnmatch (poor2 female) (treatment) if loeduc == 1

// 3) Difference in difference
  // Transform to long format
gen id = _n
reshape long poor, i(id) j(year)

logit poor treatment##c.year if loeduc == 0, or 
logit poor treatment##c.year if loeduc == 1, or 
logit poor loeduc##treatment##c.year, or 

// 4) Fixed effects model
replace treatment = 0 if year == 1

xtset id year
xtreg poor treatment year if loeduc == 0, fe 
xtreg poor treatment year if loeduc == 1, fe 
xtreg poor treatment##loeduc year##loeduc, fe 

Reference

Hu, Yannan, Frank J. van Lenthe, Rasmus Hoffmann, Karen van Hedel, and Johan P. Mackenbach. 2017. "Assessing the Impact of Natural Policy Experiments on Socioeconomic Inequalities in Health. How to Apply Commonly Used Quantitative Analytical Methods?" BMC Medical Research Methodology 17(1):68. doi: 10.1186/s12874-017-0317-5

Comparing AME's and the LPM

use doi hs01 age sex educ mstat using ZA5250_v2-0-0.dta, clear

// Poor health
recode hs01 (1 2 3 = 0 "Non-poor health") ///
            (  4 5 = 1 "Poor health") ///
            (-9 = .), gen(poorhealth)
label var poorhealth "Poor health"

// Female sex
recode sex (2 = 1 "Female") ///
           (1 = 0 "Male"), gen(female)
label var female "Female sex"

// Age
recode age (18/24 = 0 "18-24 y.") ///
           (25/34 = 1 "25-34 y.") ///
           (35/44 = 2 "35-44 y.") ///
           (45/54 = 3 "45-54 y.") ///
           (55/64 = 4 "55-64 y.") ///
           (65/74 = 5 "65-74 y.") ///
           (75/84 = 6 "75-84 y.") ///
           (85/97 = 7 "85-97 y.") ///
           (-32   = . ), gen(agecat)
label var agecat "Age"

// Education
recode educ (1 2 = 0 "Hauptschule or less") ///
            (  3 = 1 "Mittlere Reife") ///
            (4 5 = 2 "(Fach)Hochschulreife or more") ///
            (-41 -9 6 7 = .), gen(education)
label var education "Education"

// Marital status
recode mstat (1 6     = 2 "Married/cohabiting") ///
             (2 3 4 9 = 1 "Divorced, widowed etc.") ///
             (5       = 0 "Never married") ///
             (-9      = .), gen(married)
label var married "Marital status"

// Model 
eststo clear
  // AME
logit poorhealth i.female i.agecat i.education i.married 
margins, dydx(*) post
eststo
  // LPM
eststo: regress poorhealth i.female i.agecat i.education i.married, robust
esttab using amelpm.tex, drop(_cons) mtitle("AME" "LPM") label b(2) se(2) nonumbers booktabs ///
                         title("Predictors of poor self-rated health, Germany 2016 \label{tab1}") ///
                         addnote("\emph{Source}: Allbus 2016, doi: 10.4232/1.12754"  ///
                                 "AME: Average marginal effects, LPM: Linear probability model") ///
                         refcat(0.female "\emph{Sex}" ///
                                0.agecat "\emph{Age}" ///
                                0.education "\emph{Education}" ///
                                0.married "\emph{Marital status}", nolabel) ///
                         alignment(D{.}{.}{-1}) width(0.9\hsize) replace

Random graphs (94): Predicted probabilities and their differences

use 2010_ah.dta

// Prepare variables
recode ahm2010_varhours (1 = 0 "Inflexible") (2 3 4 5 = 1 "Flexible"), gen(flexible)
label variable flexible "Flexible working hours"
decode country, gen(cntry)
gen female = (sex == 2) if !missing(sex) 

// Set up loop for posting results
preserve
levelsof cntry, local(country)

tempname foo
tempname foo2
postfile `foo' str20 cntry sexdiffer lb ub using `foo2', replace

foreach x of local country {
      // Estimate model
  qui logit flexible i.female if cntry == "`x'"             
      // Predict probabilities, the r operator gives
      // differences from the reference (base) level
  qui capture margins r.female if cntry == "`x'", post
  local differ      = 100 * _b[r1vs0.female]
  local differ_loci = 100 * (_b[r1vs0.female] + (1.96 * _se[r1vs0.female]))
  local differ_hici = 100 * (_b[r1vs0.female] - (1.96 * _se[r1vs0.female]))
  post `foo' ("`x'") (`differ') (`differ_loci') (`differ_hici')
}
postclose `foo'

// Plot results
use `foo2', clear
   // Sort by size
egen order_ = rank(-sexdiffer), unique
labmask order_, value(cntry)

twoway (rcap sexdiffer sexdiffer order_, horizontal dsymbol(x)) ///
       (rspike ub lb order_, horizontal ) , ///
        xline(0) ylabel(1/30, val ang(h)) ///
        ytitle("") ///
        xtitle("Gender gap in flexible hours" "(Women minus men)") ///
        xscale(alt) ///
        legend(off) ///
        name(by_sex, replace) 
restore 

Firebaugh (1997): Analyzing Repeated Surveys in Stata (ch. 3)

set maxvar 6000
use GSS7214_R4.DTA, clear

fre sample
// Try to recreate sample
drop if sample == 4 // FP 1970 Black oversample
drop if sample == 5 // BFP 1980 Black oversample
drop if sample == 7 // FP 1980 Black oversample

// Retirement status
fre wrkstat

generate retired_others  = (wrkstat == 5)
generate retired_workers = (wrkstat == 5)
replace  retired_workers = . if inlist(wrkstat, 7, 6, 8) // Exclude homemakers, 
                                                         // students, and others
// Spending on education
gen ed_toolittle = 0
replace ed_toolittle = 1 if nateduc == 1 & !missing(nateduc)
replace ed_toolittle = 0 if nateduc  > 1 & !missing(nateduc)
replace ed_toolittle = . if missing(nateduc)

gen ed_toomuch = 0
replace ed_toomuch = 1 if nateduc == 3 & !missing(nateduc)
replace ed_toomuch = . if missing(nateduc)

// Spending on social security
gen ss_toolittle = 0
replace ss_toolittle = 1 if natsoc == 1 & !missing(natsoc)
replace ss_toolittle = 0 if natsoc  > 1 & !missing(natsoc)
replace ss_toolittle = . if missing(natsoc)

gen ss_toomuch = 0
replace ss_toomuch = 1 if natsoc == 3 & !missing(natsoc)
replace ss_toomuch = 0 if natsoc  < 3 & !missing(natsoc)
replace ss_toomuch = . if missing(natsoc)

eststo clear
// Table 3.1

preserve
keep if year >= 1973 & year <= 1993
gen trend = year - 1973

eststo: logit ed_toomuch     c.trend##retired_others
eststo: logit ed_toolittle   c.trend##retired_others
 
eststo: logit ed_toomuch     c.trend##retired_workers
eststo: logit ed_toolittle   c.trend##retired_workers

esttab est1 est2 est3 est4 using table3.1.tex, replace varwidth(50) ///
       booktabs label ///
       mgroups("Retirees vs. others" "Retirees vs. workers", pattern(1 0 1 0) ///
       prefix(\multicolumn{@span}{c}{) suffix(})   ///
       span erepeat(\cmidrule(lr){@span}))         ///
       alignment(D{.}{.}{-1}) page(dcolumn) nonumber ///
       b(%9.3f) not ///
       drop(0b.retired_others ///
            0b.retired_others#co.trend ///
            0b.retired_workers ///
            0b.retired_workers#co.trend ///
   _cons) ///
    rename(1.retired_others 1.retired_workers ///
              1.retired_others#c.trend 1.retired_workers#c.trend) ///
       coeflabels(1.retired_workers "Initial difference, retirees minus workers/others" ///
               trend "Trend for retirees" ///
      1.retired_workers#c.trend "Trend difference retirees minus workers/others" ///
                  1.retired_others "Initial difference, retirees minus others" ///
      1.retired_others#c.trend "Trend difference retirees minus workers/others") ///
       mlabels("\multicolumn{1}{c}{Too much}" ///
            "\multicolumn{1}{c}{Too little}" ///
      "\multicolumn{1}{c}{Too much}" ///
            "\multicolumn{1}{c}{Too little}") ///
       order(1.retired_workers trend 1.retired_workers#c.trend) ///
    title("Table 3.1: Trend analysis for spending on education, 1973-1993: Logit coefficients")

restore

// Table 3.2
preserve
keep if year >= 1984 & year <= 1993
gen trend = year - 1984

eststo: logit ss_toomuch   c.trend##retired_others
eststo: logit ss_toolittle c.trend##retired_others
eststo: logit ss_toomuch   c.trend##retired_workers
eststo: logit ss_toolittle c.trend##retired_workers

esttab est5 est6 est7 est8 using table32.tex, replace varwidth(50) ///
       booktabs label ///
       mgroups("Retirees vs. others" "Retirees vs. workers", pattern(1 0 1 0) ///
       prefix(\multicolumn{@span}{c}{) suffix(})   ///
       span erepeat(\cmidrule(lr){@span}))         ///
       alignment(D{.}{.}{-1}) page(dcolumn) nonumber ///
       b(%9.3f) not ///
       drop(0b.retired_others ///
            0b.retired_others#co.trend ///
            0b.retired_workers ///
            0b.retired_workers#co.trend ///
   _cons) ///
    rename(1.retired_others 1.retired_workers ///
              1.retired_others#c.trend 1.retired_workers#c.trend) ///
       coeflabels(1.retired_workers "Initial difference, retirees minus workers/others" ///
               trend "Trend for retirees" ///
      1.retired_workers#c.trend "Trend difference retirees minus workers/others" ///
                  1.retired_others "Initial difference, retirees minus others" ///
      1.retired_others#c.trend "Trend difference retirees minus workers/others") ///
       mlabels("\multicolumn{1}{c}{Too much}" ///
            "\multicolumn{1}{c}{Too little}" ///
      "\multicolumn{1}{c}{Too much}" ///
            "\multicolumn{1}{c}{Too little}") ///
       order(1.retired_workers trend 1.retired_workers#c.trend) ///
    title("Table 3.2: Trend analysis for spending on social security, 1984-1993: Logit coefficients")
restore

Reference

Firebaugh, Glenn. 1997. Analyzing Repeated Surveys. Sage. doi: 10.4135/9781412983396