Random graphs (123): Descriptive statistics

// Open Allbus 2010
use v11 v13 v15 v17 v19 using ZA4612_v1-0-1.dta, clear

// Prepare variables
recode v11 v13 v15 v17 v19 (99 = .)
label define x 1 `""1" "Unimportant""' 7 `""7 Very" "important"'

label var v11 "v11: Job security"
label var v13 "v13: High income"
label var v15 "v15: Opportunities for promotion"
label var v17 "v17: Respectable occupation"
label var v19 "v19: Lots of leisure time"

// Plot histograms
foreach x of varlist v11 v13 v15 v17 v19 {
  label val `x' x
  histogram `x', frequency ///
                 xlabel(1 (1) 7, val) ///
                 ylabel(0 (500) 1500) ///
                 ytick(0 (250) 1750) ///
                 name(`x', replace) nodraw
}

// Calculate correlation matrix
cor v11-v19
matrix correlations = r(C)

// Plot correlation matrix
plotmatrix, mat(correlations) split(-1(0.10)1) color(sienna) ///
            freq formatcells(%6.2f) nodiag legend(title(Pearson's {it:r}) col(1)) ///
            name(cormatrix, replace) nodraw 

// Combine plots
graph combine v11 v13 v15 v17 v19 cormatrix, col(2) ysize(8) xsize(6) altshrink

Random graphs (116): Histograms

use "ESS8e01.dta", clear

kountry cntry, from(iso2c) marker
rename NAMES_STD country

recode lrscale (77 88 99 = .)

histogram lrscale, by(country, note("")) ///
                   percent discrete ///
                   xlabel(0 (2) 10, val) ///
                   ylabel(0 (10) 50) ///
                   xtitle("Placement on left-right scale") ///
                   name(figure1, replace)

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

OLS regression: overfitting and dichotomization

Babyak (2004) demonstrates a couple of aspects of OLS regression, making use of the following simulations:

// Figure 1
set seed 1

// Create file to store simulation results
tempname foo
postfile `foo' b using clt, replace

// Simulate analyses
forvalues i = 1/10000 {
    drop _all
 qui set obs 100
    generate x = rnormal()
 generate e = rnormal()
    generate y = .4 * x + e
    qui regress y x
 local b = _b[x]
 
 post `foo' (`b') 
}
postclose `foo'

// Open results from simulations and plot
use clt, clear
histogram b, freq ytitle("Frequency of b value") xtitle("Values of b") name(figure1, replace)



// Figure 2
set seed 1

// Create file to store simulation results
tempname foo
postfile `foo' r2_50 r2_100 r2_150 r2_200 using overfitting, replace

// Simulate analyses
forvalues i = 1/10000 {
    drop _all
 set obs 10000
    generate y = rnormal()
    foreach i of numlist 1/15 {
      generate x`i' = rnormal()
    }
 
 foreach j of numlist 50 100 150 200 {
  preserve
  sample `j', count
  qui reg y x*
  local r2_`j' = e(r2)
     restore
 }
    post `foo' (`r2_50') (`r2_100') (`r2_150') (`r2_200')
}
postclose `foo'

// Open results from simulations
use overfitting, clear

// Plot
twoway (kdensity r2_200)  ///
       (kdensity r2_150) ///
    (kdensity r2_100) ///
    (kdensity r2_50) ///
    , ytitle("Percent of samples") ///
         xtitle("R-square value from regression model") ///
   xlabel(0 (.1) .6) ///
         ylabel(0 (2) 20) ///
         legend(order(1 "ca. 13 cases/predictor ({it:N} = 200)" ///
                      2 "10 cases/predictor ({it:N} = 150)" ///
                      3 "ca. 7 cases/predictor ({it:N} = 100)" ///
                      4 "ca. 3 cases/predictor ({it:N} = 50)") ///
       pos(2) ring(0)) name(figure2, replace)
// Figure 2
set seed 1

// Create file to store simulation results
tempname foo
postfile `foo' r2_50 r2_100 r2_150 r2_200 using overfitting, replace

// Simulate analyses
forvalues i = 1/10000 {
    drop _all
 set obs 10000
    generate y = rnormal()
    foreach i of numlist 1/15 {
      generate x`i' = rnormal()
    }
 
 foreach j of numlist 50 100 150 200 {
  preserve
  sample `j', count
  qui reg y x*
  local r2_`j' = e(r2)
     restore
 }
    post `foo' (`r2_50') (`r2_100') (`r2_150') (`r2_200')
}
postclose `foo'

// Open results from simulations
use overfitting, clear

// Plot
twoway (kdensity r2_200)  ///
       (kdensity r2_150) ///
       (kdensity r2_100) ///
       (kdensity r2_50) ///
      , ytitle("Percent of samples") ///
        xtitle("R-square value from regression model") ///
        xlabel(0 (.1) .6) ///
        ylabel(0 (2) 20) ///
        legend(order(1 "ca. 13 cases/predictor ({it:N} = 200)" ///
                     2 "10 cases/predictor ({it:N} = 150)" ///
                     3 "ca. 7 cases/predictor ({it:N} = 100)" ///
                     4 "ca. 3 cases/predictor ({it:N} = 50)") ///
        pos(2) ring(0)) name(figure2, replace)



// Figure 4 (actually Table 1)
clear 
set seed 1
// Create file to store simulation results
tempname foo
postfile `foo' n correlation typei using dichotomization, replace

// Simulate analyses
forvalues i = 1/10000 {
    drop _all
    foreach j of numlist 50 100 200 {
   foreach k of numlist 0 .3 .5 .7 {
     qui drawnorm x1 x2, ///
              n(`j') ///
              corr(1, `k', 1) cstorage(lower) ///
              clear
        generate e = rnormal()
        generate y = .5*x1 + 0*x2 + e
        qui sum x1, detail
        generate x1s = (x1 > r(p50))
        qui sum x2, detail
        generate x2s = (x2 > r(p50))
        qui reg y x1s x2s
  local typei = _b[x2s]/_se[x2s]
  local sig   = (abs(`typei') > 1.96) 
  *di _b[x2s] _skip(5) _se[x2s] _skip(5) `typei' _skip(5) `sig'
  post `foo' (`j') (`k') (`sig') 
   }
 }
    
}
postclose `foo'

// Open results from simulations
use dichotomization, clear

collapse typei, by(n correlation)

graph hbar typei, over(n, relabel(1 "{it:N} = 50" 2 "{it:N} = 100" 3 "{it:N} = 200")) ///
                 over(correlation, relabel(1 "{it:Corr(x{sub:1}, x{sub:2})} = 0" ///
                                           2 "{it:Corr(x{sub:1}, x{sub:2})} = .3" ///
                                           3 "{it:Corr(x{sub:1}, x{sub:2})} = .5" ///
                                           4 "{it:Corr(x{sub:1}, x{sub:2})} = .7")) ///
                 ytitle("Type I error rate") yscale(alt) ylabel(, format(%6.2f)) ///
                 name(figure4, replace)

Reference

Babyak, Michael A. 2004. "What You See May Not Be What You Get. A Brief, Nontechnical Introduction to Overfitting in Regression-Type Models." Psychosomatic Medicine 66(3):411-421. doi: 10.1097/01.psy.0000127692.23278.a9

Chapters 9 and 10 of Singer and Willett's (2003) book on longitudinal data analysis

// Figure 9.1


use "C:\singer willett (2003)\teachers.dta", clear

label define censor 0 "Not censored" 1 "Censored"
label val censor censor

histogram t, by(censor, legend(off) note("")) freq xlabel(1 (1) 12) ylabel(0 (100) 500) ///
             addlabels discrete xtitle("Years of teaching") name(figure91, replace)

// Table 10.1


generate event = !censor
qui ltable t event, noadjust saving(lifetable, replace)

preserve
use lifetable, clear
list t0 t1 start deaths lost hazard survival, sep(0) noobs
erase lifetable.dta
restore

// Figure 10.1


stset t, failure(event)
sts generate h = h
twoway scatter h t, msymbol(i) connect(l) ylabel(0 (.05) .15) xlabel(0 (1) 13) ///
                    sort xtitle("Years in teaching") ///
                    ytitle("Estimated hazard probability") ///
                    name(figure101a, replace) nodraw
ltable t event, noadjust notab graph noconf xlabel(0 (1) 13) ylabel(0 (.5) 1) ///
                    yline(.5) xline(7.6) ///
                    xtitle("Years in teaching") ///
                    ytitle("Proportion surviving") ///
                    name(figure101b, replace) nodraw
graph combine figure101a figure101b, col(1) name(figure101, replace) ysize(8)

// Figure 10.2


use "C:\singer willett (2003)\relapse_days.dta", clear
generate weeks = int(days / 7) + 1 
generate event = !censor 
stset weeks, failure(event)
sts generate h = h
twoway (scatter h weeks, connect(l) msymbol(i)), xlabel(0 (1) 12) ylabel(, format(%6.2f)) ///
                    ytitle("Estimated hazard probability") ///
                    xtitle("Weeks after release") ///
                    name(figure102a1, replace) nodraw
ltable weeks event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) ///
                    xlabel(0(1)12) yline(0.5) ///
                    ytitle("Estimated survival probability") ///
                    xtitle("Weeks after release") ///
                    name(figure102a2, replace) nodraw
graph combine figure102a1 figure102a2, col(2) ///
      title("{bf:A} Time to cocaine relapse", pos(11)) name(figure102a, replace) nodraw

use "C:\singer willett (2003)\firstsex.dta", clear
generate event = !censor 
stset time, failure(event)
sts generate h = h
twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(6 (1) 12) ylabel(, format(%6.2f)) ///
                    ytitle("Estimated hazard probability") ///
                    xtitle("Grade") ///
                    name(figure102b1, replace) nodraw
ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) ///
                    xlabel(6 (1) 12) yline(0.5) ///
                    ytitle("Estimated survival probability") ///
                    xtitle("Grade") ///
                    name(figure102b2, replace) nodraw
graph combine figure102b1 figure102b2, col(2) ///
      title("{bf:B} Age at first intercourse for males", pos(11)) name(figure102b, replace) nodraw

use "C:\singer willett (2003)\suicide_orig.dta", clear
generate event = !censor 
stset time, failure(event)
sts generate h = h
twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(5 (2) 21) ylabel(, format(%6.2f)) ///
                    ytitle("Estimated hazard probability") ///
                    xtitle("Age") ///
                    name(figure102c1, replace) nodraw
ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) ///
                    xlabel(5 (2) 21) yline(0.5) ///
                    ytitle("Estimated survival probability") ///
                    xtitle("Age") ///
                    name(figure102c2, replace) nodraw
graph combine figure102c1 figure102c2, col(2) ///
      title("{bf:C} Age at first suicide ideation", pos(11)) name(figure102c, replace) nodraw

use "C:\singer willett (2003)\congress_orig.dta", clear
generate event = !censor 
stset time, failure(event)
sts generate h = h
twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(0 (1) 8) ylabel(0 (.1) .3, format(%6.2f)) ///
                    ytitle("Estimated hazard probability") ///
                    xtitle("Terms in office") ///
                    name(figure102d1, replace) nodraw
ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) ///
                    xlabel(0 (1) 8) yline(0.5) ///
                    ytitle("Estimated survival probability") ///
                    xtitle("Terms in office") ///
                    name(figure102d2, replace) nodraw
graph combine figure102d1 figure102d2, col(2) ///
      title("{bf:D} Duration of congressional careers for females", pos(11)) name(figure102d, replace) nodraw

graph combine figure102a figure102b figure102c figure102d, row(4) ysize(11) name(figure102, replace)

// Table 10.2
use "C:\singer willett (2003)\teachers.dta", clear

generate event = !censor
qui ltable t event, noadjust saving(lifetable, replace)

preserve
use lifetable, clear
list t0 hazard sehazard survival sesurvival, noobs sep(0)
erase lifetable.dta
restore

// Figure 10.4
use "C:\singer willett (2003)\teachers.dta", clear

list id t censor if inlist(id, 20, 126, 129), noobs sep(0)

expand t
bysort id: generate period = _n  // Period identifier
bysort id: generate event = !censor & _n == _N // Calculate outcome for 
                                               // discrete-time analysis

list id period event if inlist(id, 20, 126, 129), noobs sepby(id)

// Table 10.3
table period, c(sum event n event mean event)

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

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