Mummolo & Peterson (2018): Fixed-effects regression

set scheme lean1

webuse nlswork, clear
xtset idcode year

// Model of interest:
xtreg ln_w ttl_exp age c.age#c.age tenure c.tenure#c.tenure not_smsa south i.year, fe
local b : display %4.2f _b[ttl_exp]

// Isolate relevant variation in the treatment:
reghdfe ttl_exp, absorb(idcode year) residuals(relevantvariation)

// Identify a plausible counterfactual shift in X given the data:
center ttl_exp

sum relevantvariation
local sd1 : display %4.2f r(sd)

sum ttl_exp
local sd2 : display %4.2f r(sd)

twoway (kdensity c_ttl_exp) ///
       (kdensity relevantvariation), ///
       ytitle("Density") xtitle(Treatment (centered around 0)) ///
       title("{bf:A} Treatment distribution", span bexpand justification(left)) ///
       legend(order(1 "Original distribution, {it:SD} = `sd2'" ///
                    2 "Distribution after FE, {it:SD} = `sd1'" ) ///
                    pos(2) ring(0)) ///
       name(figure1A, replace)

      // Calculate range within individuals
bysort idcode: egen max = max(ttl_exp)
bysort idcode: egen min = min(ttl_exp)
generate range = max - min 
drop max min

qui sum range, detail
local mean : display %4.2f r(mean)
local p95 :  display %4.2f r(p95)

twoway (histogram range, freq) ///
       (scatteri 2000 `mean' (3) "Mean = `mean'", msymbol(i)) ///
       (scatteri 2000 `p95' (3) "95th percentile = `p95'" , msymbol(i)), ///
       title("{bf:B} Within-individual range of treatment", span bexpand justification(left)) ///
       ylabel(, format(%7.0gc)) ytitle(Frequency) ///
       xline(`mean') xline(`p95') ///
       xtitle(Within-individual range) legend(off) ///
       name(figure1B, replace)

graph combine figure1A figure1B, col(2)

// Calculate % of cases where treatment does not vary
egen pickone = tag(idcode)

fre range if pickone
  // In 12% of cases the treatment does not vary

// Multiply the estimated coefficient of interest by the revised standard deviation

display "Estimated coefficient of interest = `b'"
display "Within-standard deviation = `sd1'"
display "Product = " _skip(1) `sd1' * `b'

Reference Mummolo, Jonathan, and Erik Peterson. 2018. "Improving the Interpretation of Fixed Effects Regression Results." Political Science Research and Methods 6(4):829-835. doi: 10.1017/psrm.2017.44

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)    

Random graphs (106): Interaction plot with overlaid density

// Open Allbus 2008
use V154 V151 V156 V760 V754 V755 V5 V767 using ZA4602_v1-0-0.dta, clear

// Age
generate age = V154 if V154 != 999

// Sex
generate female = (V151 == 2)

// Migrant
generate migrant = (V156 == 2)

// Interviewer ID
rename V760 id

// Interviewer sex
generate femaleinterviewer = (V754 == 2)

// Interviewer age
generate interviewerage = V755

// Attractiveness rating by interviewer before and after interview
qui factor V5 V767, pcf
predict attractiveness

// Calculate Spearman-Brown for two-item scales according
// to Eisinga et al. (https://doi.org/10.1007/s00038-012-0416-3):
spearman V5 V767

// Fit model with respondents clustered in interviewers
mixed attractiveness c.age##i.female i.migrant i.femaleinterviewer##c.interviewerage || id: 

// Calculate margins and plot
qui margins, at(age = (18  (10) 98) female = (0 1))
marginsplot, recast(line) recastci(rarea) ciopt(color(gs14)) ///
             plot1opts(lpattern(dash)) ///
             legend(ring(0) pos(2)) ///
             title("Female beauty premium disappears with age", span) ///
             ytitle("Predicted attractiveness", axis(1)) ///
             xtitle("Age") ///
             xlabel(20 (10) 100) ///
             addplot(histogram age, discrete yaxis(2) ///
                                    lcolor(white) ///
                                    ylabel(0 0.01 0.02, format(%6.2f) axis(2)) ///
                                    yscale(alt range(0 0.1) axis(2)) ///
                                    ytitle("Age density", axis(2)) ///
                                    legend(order(3 "Men" 4 "Women"))) ///
             note(" " "{it:Source:} German General Social Survey Allbus 2008, doi:10.4232/1.12345", span) ///
             name(figure1, replace)

Random graphs (96): Categorical variables

use S003A S002 X001 C041 C039 C038 C037 if inlist(S002, 4, 5) ///
    using WVS_Longitudinal_1981-2014_stata_v_2014_11_25.dta, clear

// Declare missing values
mvdecode C041 C039 C038 C037, mv(-5 -4 -3 -2 -1)

// Factor analysis
factor C041 C039 C038 C037, pcf
alpha C041 C039 C038 C037, item
local alph = round(`r(alpha)', .01)

// Prepare variables
scores nonworkethic = mean(C041 C039 C038 C037), nv(3) 
generate workethic = 5 - nonworkethic 
gen male = (X001 == 1) if X001 != -2
drop nonworkethic X001


// Create descriptive plot    
foreach var of varlist C041 C039 C038 C037 {
 recode `var' (1 = 4 "Strongly agree") ///
              (2 = 3 "Agree") ///
              (3 = 2 "Neither agree nor disagree") ///
              (4 = 1 "Disagree") ///
              (5 = 0 "Strongly disagree"), gen(`var'_rec)
 numlabel `var'_rec, add mask("# ")
 twoway histogram `var'_rec, discrete horizontal yla(0/4, valuelabel) percent ///
                             title("`: variable label `var''") ///
                             ytitle("") ///
                             name(`var', replace) nodraw
 drop `var'_rec
}
graph combine C041 C039 C038 C037, col(2) row(2) altshrink title("Work ethic is measured as the average of four items:") ///
      caption("In a PCA, all items load on one dimension; explained variance = 49%, Cronbach's alpha = `alph'", span) note("{it:Source:} WVS 1999-2009, pooled", span) name(figure1, replace)
drop C041 C039 C038 C037


// Calculate country differences
preserve   
  statsby mean_ = _b[_cons] ///
          loci  = (_b[_cons] - 1.96 * _se[_cons]) ///
          hici  = (_b[_cons] + 1.96 * _se[_cons]) ///
        , by(S003A) clear: ///
          regress workethic

  // Sort coefficients by size
  egen order_ = rank(mean_), unique
  labmask order_, value(S003A) decode

  // Plot
  twoway (rcap mean_ mean_ order_, horizontal) ///
         (rspike loci hici order_, horizontal) ///
        , legend(off) ylabel(1/63, valuelabels ang(h) labsize(*.8)) ///
          xlabel(0 (1) 4, grid format(%6.1f)) name(all, replace)  ///
          xmtick(0 (.5) 4) ///
          ytitle("") xtitle("Work ethic across countries") ///
          title("{bf:A}", justification(left) bexpand span) ///
          xscale(alt) ysize(10) nodraw
restore  

// Country differences--men only
preserve   
  statsby mean_ = _b[_cons] ///
          loci  = (_b[_cons] - 1.96 * _se[_cons]) ///
          hici  = (_b[_cons] + 1.96 * _se[_cons]) ///
        , by(S003A) clear: ///
          regress workethic if male == 1

  // Sort coefficients by size
  egen order_ = rank(mean_), unique
  labmask order_, value(S003A) decode

  // Plot
  twoway (rcap mean_ mean_ order_, horizontal) ///
         (rspike loci hici order_, horizontal) ///
        , legend(off) ylabel(1/63, valuelabels ang(h) labsize(*.8)) ///
          xlabel(0 (1) 4, grid format(%6.1f)) name(men, replace)  ///
          xmtick(0 (.5) 4) ///
          ytitle("") xtitle("Men's work ethic across countries") ///
          title("{bf:B}", justification(left) bexpand span) ///
          xscale(alt) ysize(10) nodraw
restore  

// Calculate gender gap
preserve   
  statsby mean_ = _b[male] ///
          loci  = (_b[male] - 1.96 * _se[male]) ///
          hici  = (_b[male] + 1.96 * _se[male]) ///
        , by(S003A) clear: ///
          regress workethic male

  // Sort coefficients by size
  egen order_ = rank(mean_), unique
  labmask order_, value(S003A) decode

  // Plot
  twoway (rcap mean_ mean_ order_, horizontal) ///
         (rspike loci hici order_, horizontal) ///
        , legend(off) ylabel(1/63, valuelabels ang(h) labsize(*.8)) ///
          xlabel(-.25 (.25) .5, grid format(%6.2f)) name(gendergap, replace)  ///
          xmtick(-.25 (.1) .5) ///
          text(63 .5  "Men" "higher", place(sw)) ///
          text( 1 -.25 "Women" "higher", place(ne)) ///    
          ytitle("") xtitle("Gender gap in work ethic") ///
          title("{bf:C}", justification(left) bexpand span) ///
          xscale(alt) ysize(10) nodraw
restore  

graph combine all men gendergap, note(" " "{it:Source:} WVS 1999-2009, pooled", span) col(3) ysize(12) xsize(18) altshrink ///
                                 title("Work ethic in cross-national comparison", span) name(figure2, replace)

Random graphs (62): Histograms

qui summarize challenging
local m = round(r(mean), .1)
local sd = round(r(mean), .1)
local n = r(N)
numlabel challenging, add mask("# ")
twoway histogram challenging, discrete yla(, valuelabel) horizontal gap(20)  ///
                              freq ytitle("") start(1) xlabel(0/4, grid) ///
                              title("How" "{bf:challenging}" "was this paper?", span) ///
                              xsize(3) note("Mean = `m', {it:SD} = `sd', {it:N} = `n'", span) ///
                              name(challenging, replace)

qui summarize interesting
local m = round(r(mean), .1)
local sd = round(r(mean), .1)
local n = r(N)
numlabel interesting, add mask("# ")
twoway histogram interesting, discrete yla(, valuelabel) gap(20) horizontal ///
                              freq ytitle("") start(1) xlabel(0/4, grid) ///
                              title("How" "{bf:interesting}" "was this paper?", span) ///
                              xsize(3) note("Mean = `m', {it:SD} = `sd', {it:N} = `n'", span) ///
                              name(interesting, replace)

qui summarize recommendable
local m = round(r(mean), .1)
local sd = round(r(mean), .1)
local n = r(N)
numlabel recommendable, add mask("# ")
twoway histogram recommendable, discrete yla(, valuelabel) gap(20) horizontal ///
                                freq ytitle("") start(1) xlabel(0/4, grid) ///
                                title("How likely are you to" "{bf:recommend}" "this paper to a friend?", span) ///
                                xsize(3) note("Mean = `m', {it:SD} = `sd', {it:N} = `n'", span) ///
                                name(recommendable, replace)
         
graph combine challenging interesting recommendable, row(1)

Random graphs (58): Descriptives of a categorical variable

numlabel health, add mask("[#] ")
twoway histogram health, discrete yla(, valuelabel) gap(20) horizontal percent 
numlabel health, remove mask("[#] ")

Random graphs (30): Log-log histograms

Healy and Moody (2014) recommend log-log histograms as an alternative to regular histograms. Given the fact that important variables in Sociology are often highly skewed, log-log histograms reveal greater details at the upper end of the distribution

// Membership variables of the first round (2002) of the ESS
use sptcmmb cltommb trummb prfommb ///
    cnsommb hmnommb epaommb rlgommb ///
    prtymmb setommb sclcmmb othvmmb using "ESS1e06.3_F1.dta", clear

// Generate count variable of memberships
scores memberships = total(sptcmmb cltommb trummb prfommb ///
                           cnsommb hmnommb epaommb rlgommb ///
                                 prtymmb setommb sclcmmb othvmmb)
                           
// Standard histogram
twoway (histogram memberships, percent     discrete), ///
        xtitle("Total no. of memberships in" "voluntary associations in last 12 mo.") ///
          ytitle("Respondents in ESS round 1 (%)") xlabel(0/12) name(a, replace)
         
// Log-log histogram
twoway__histogram_gen memberships, gen(y x) width(1) fraction
replace y = y * 100 // Convert fractions to %

scatter y x, yscale(log) xscale(log) ///
        xtitle("Total no. of memberships in" "voluntary associations in last 12 mo.") ///
        ytitle("Respondents in ESS round 1 (%)") xlabel(0/12) ///
  ylab(.1 1 10 20 40) ///
  name(b, replace)

graph combine a b    

Reference 

Healy, Kieran, and James Moody. 2014. "Data Visualization in Sociology." Annual Review of Sociology 40:105-128. doi: 10.1146/annurev-soc-071312-145551