Specification curve analysis

use "ZA4612_v1-0-1.dta", clear
do "ZA4612_patch_v1-0-1.do" // Some patch from data provider

keep v399-v404 v933 v827 v301 v298 

// Outcomes
mvdecode v399-v404, mv(9 = .a )

generate stress  = 5 - v399
generate depress = 5 - v400
generate calm    = v401 - 1
generate energy  = v402 - 1
generate pain    = 5 - v403
generate lonely  = 5 - v404

// Predictors
  // Topbot
mvdecode v933 v827, mv(96 = .a \ 99 = .b) 
generate topbot = v933
replace  topbot = v827 if missing(topbot)
  // Age
mvdecode v301, mv(999 = .a)
rename v301 age
  // Female
generate female = (v298 == 2)

drop v399-v404 v933 v827 v298 // Drop unnecessary variables

// Post definition
tempname foo
postfile `foo' str50 spec k1 k2 k3 b se using deleteme.dta, replace

// Loop
forvalues k1 = 1/6 {
  if `k1' == 1 local y "stress"
  if `k1' == 2 local y "depress"
  if `k1' == 3 local y "pain"
  if `k1' == 4 local y "calm"
  if `k1' == 5 local y "energy"
  if `k1' == 6 local y "lonely"
    forvalues k2 = 1/4 {
      if `k2' == 1 local agecontrol "age"
      if `k2' == 2 local agecontrol "c.age##c.age"
      if `k2' == 3 local agecontrol "c.age##c.age##c.age"    
      if `k2' == 4 local agecontrol "c.age##c.age##c.age##c.age"
        forvalues k3 = 1/3 {
          if `k3' == 1 local ifs " "
          if `k3' == 2 local ifs "if female == 1"
          if `k3' == 3 local ifs "if female == 0"
 
        local spec regress `y' topbot `agecontrol' female `ifs'

        qui `spec'
        local  b =  _b[topbot]
        local se = _se[topbot]

        post `foo' ("`spec'") (`k1') (`k2') (`k3') (`b') (`se')
        }
    }
}
postclose `foo'

// Plot specification curve
use deleteme, clear

// Generate ranked analytical choice variable
sort b
generate sk = _n
label var sk "Specification (sorted by coefficient size)"

// Calculate CI's
generate ub = b + 1.96 * se
generate lb = b - 1.96 * se

// Remind yourself what the variables mean
label var k1 "Outcome"
label var k2 "Age control"
label var k3 "Subsample"

// Stack indicators
generate k3c = k3
generate k2c = k2 + 1 + 3 // 3 because K3 has 3 categories, 1 for title
generate k1c = k1 + 1 + 3 + 4 + 1  // K2 has 4 categories

// Calculate some things for size of second axis
qui summarize b
global brange = r(max) - r(min)
global bmin = r(min)
global bmax = r(max)
global from_y = $bmin - (4.5 * $brange)

// Plot
twoway (scatter k1c k2c k3c sk, msymbol(o o o) msize(vsmall vsmall vsmall) ///
        yscale(range(1 22)) ///        
        ylabel( 1 "Full sample" 2 "Females only" 3 "Males only" 4 "{bf:Subsample}" ///
                5 "Age" 6 "Age squared" 7 "Age cubed" 8 "Age quartic" 9 "{bf:Age control}" ///
               10 "Stress"   11 "Depression" 12 "Pain" ///
               13 "Calmness" 14 "Energy"     15 "Loneliness" 16 "{bf:Outcome}", tstyle(notick) axis(1)))  ///
       (rcap b b sk, yaxis(2) yscale(range($from_y $bmax) axis(2)) ///
                     ylab(, format(%2.1g) axis(2)) ///
                     ytitle("{bf:Coefficient size}", axis(2) placement(north)) ///
                     yline(0, axis(2))) ///
       (rspike ub lb sk, yaxis(2)), ///
        xlabel(none)  ///
        legend(off) name(curve, replace) 

Reference

Simonsohn, Uri, Joseph P. Simmons, and Leif D. Nelson. 2015. Specification Curve. Descriptive and Inferential Statistics on All Reasonable Specifications. University of Pennsylvania. doi: 10.2139/ssrn.2694998

Merging the Demographic and Health Surveys in Stata

// Unzip files
clear
cd "C:\dhs"

local filelist : dir . files "*.zip"  // Create local with all filenames ending with ".zip"
di `filelist'

local first : word 1 of `filelist'      // Identify first file
di "`first'"

local total_ : word count `filelist' // Identify total number of files
di `total_'        

forvalues x = 1/`total_' {
    di `x'
    local y : word `x' of `filelist'
    unzipfile "`y'", replace
}

// Append all files
clear
cd "C:\dhs"
local directorylist : dir . dirs  "*ir*"  // Create local with all directory names that contain women's data

di `directorylist'

local firstdir : word 1 of `directorylist'      // Identify first directory
di "`firstdir'"

local total_ : word count `directorylist' // Identify total number of directories
di `total_'

local firstfile = strupper("`firstdir'")
local firstfile = subinstr("`firstfile'","DT","FL",.)
di "`firstfile'"

use caseid v000 v005 v007 v012 v106 v107 v155 v191 v201 v212 v525 v531 using "`firstdir'/`firstfile'", clear                    
capture decode v106, gen(v106s)
drop v106

save testfile, replace

forvalues x = 2/`total_' {
    local y : word `x' of `directorylist'
 local filename = strupper("`y'")
    local filename = subinstr("`filename'", "DT", "FL", .)
 use "`y'/`filename'", clear
 
 // Source: https://stackoverflow.com/questions/17056016/stata-how-to-keep-a-list-of-variables-given-some-of-them-may-not-exist
 local masterlist "caseid v000 v005 v007 v012 v106 v107 v155 v191 v201 v212 v525 v531"
    local keeplist = ""

    foreach i of local masterlist  {
    capture confirm variable `i'
        if !_rc {
            local keeplist "`keeplist' `i'"
        }
     }
    keep `keeplist'
 capture decode v106, gen(v106s)
    capture drop v106
 
 tempfile new
 save `new', replace
 use testfile, clear
 append using `new', force
 save testfile, replace
}

// Prepare variables
use testfile, clear

//Country and wave identifiers
replace v000 = "VN3" if v000 == "VNT"
generate cntry = substr(v000,1,2)
generate wavex  = substr(v000,3,1)
replace wavex = "1" if wavex == ""
encode wavex, gen(wave)
drop wavex

  // Fix unusual country abbreviations
replace cntry = "BI" if cntry == "BU"
replace cntry = "IN" if cntry == "IA"
replace cntry = "KZ" if cntry == "KK"
replace cntry = "BI" if cntry == "BU"
replace cntry = "MD" if cntry == "MB"
replace cntry = "NA" if cntry == "NM"
replace cntry = "DO" if cntry == "DR"

kountry cntry, from(iso2c)
encode NAMES_STD, gen(country)
drop NAMES_STD

// Select last wave
keep if wave == 6

// Prepare variables
  // Age at first intercourse
generate age1stintercourse = .
replace  age1stintercourse = v531 if inrange(v531, 1, 63)
replace  age1stintercourse = .a   if v525 == 0
replace  age1stintercourse = .b   if v525 == 95
replace  age1stintercourse = .c   if v525 == 97
replace  age1stintercourse = .d   if v525 == 98
replace  age1stintercourse = .e   if v525 == 99
label define age1stintercourse .a "Not had intercourse" ///
                               .b "95?" ///
                               .c "inconsistent" ///
                               .d "don't know" ///
                               .e "99?" 
label val age1stintercourse age1stintercourse
label var age1stintercourse "Age at first intercourse"

  // Age at first birth
rename v212 afb
label var afb "Age at first birth"

// Calculate correlation
preserve
statsby mean_ = _b[age1stintercourse] ///
        loci  = (_b[age1stintercourse] - 1.96 * _se[age1stintercourse]) ///
        hici  = (_b[age1stintercourse] + 1.96 * _se[age1stintercourse]) ///
      , by(country) total clear: ///
        regress afb age1stintercourse

// Label total value
replace country = 1000 if country == .
label define country 1000 "{bf: Total}", modify

egen order_ = rank(mean_), unique
labmask order_, value(country) decode

// Plot
twoway (rcap mean_ mean_ order_, horizontal) ///
       (rspike loci hici order_, horizontal) ///
      , legend(off) ylabel(1/38, valuelabels ang(h) labsize(*.8)) ///
        xlabel(.7 (.1) 1.0, grid format(%6.1f)) name(ols, replace)  ///
        xmtick(.7 (.05) 1.0) ///
        ytitle("") xtitle("Association between" ///
                          "age at first intercourse" ///
                          "and age at first birth") xscale(alt) ysize(8) ///
        note(" " ///
             "{it:Source:} DHS VI, own calculations" , span size(*.8))
restore

Selection of regression predictors

clear
set seed 1
set obs 1000

// Generate random variable y
generate y = rnormal()

// Generate 50 random variables x
forvalues i = 01/50 {
 generate x`i' = rnormal()
}

// Model 1: Regress y on the x's
quietly regress y x1-x50
estimates store model1
coefplot model1, xline(0) xscale(alt) ylabel(, labsize(*.7)) ///
                 xtitle("Regression weights and 90% CI's", size(*.8)) levels(90) ///
                 xlabel(, format(%6.2f) labsize(*.8)) ///
                 title("Model 1") ysize(8) xsize(3) ///
                 drop(_cons) msymbol(o) name(model1, replace)

// Identify variables significant at 10% level
forvalues i = 1/50 {
 local t = _b[x`i'] / _se[x`i']
 local p = 2 * ttail(e(df_r), abs(`t'))
 if `p' <= .10 {
   local significant10 `significant10' x`i'
 }
 di "x`i'" _skip(5) `t' _skip(5) `p' _skip(5) "`significant10'"
}

// Model 2: Regress y on the x's significant at the 10 per cent level
quietly regress y `significant10'
estimates store model2
coefplot model2, xline(0) xscale(alt) ylabel(, labsize(*.7)) ///
                 xtitle("Regression weights and 90% CI's", size(*.8)) levels(90) ///
                 xlabel(, format(%6.2f) labsize(*.8)) ///
                 title("Model 2") ysize(8) xsize(3) ///
                 drop(_cons) msymbol(o) name(model2, replace)
    

// Identify variables significant at 25% level
estimates restore model1
forvalues i = 1/50 {
 local t = _b[x`i'] / _se[x`i']
 local p = 2 * ttail(e(df_r), abs(`t'))
 if `p' <= .25 {
   local significant25 `significant25' x`i'
 }
 di "x`i'" _skip(5) `t' _skip(5) `p' _skip(5) "`significant25'"
}

// Model 3: Regress y on the x's significant at the 25 per cent level
quietly regress y `significant25'
estimates store model3
coefplot model3, xline(0) xscale(alt) ylabel(, labsize(*.7)) ///
                 xtitle("Regression weights and 90% CI's", size(*.8)) levels(90) ///
                 xlabel(, format(%6.2f) labsize(*.8)) ///
                 title("Model 3") ysize(8) xsize(3)  ///
                 drop(_cons) msymbol(o) name(model3, replace)

graph combine model1 model2 model3, row(1) xcommon

Appending all files in a directory

cd "E:\eu-lfs\converted files"

local filelist : dir . files "*_y.dta"  // Create local wit all filenames ending with "_y.dta
di `filelist'

local first : word 1 of `filelist'      // Identify first file
di "`first'"

local total_ : word count `filelist' // Identify total number of files
di `total_'

use "`first'", clear                    

forvalues x = 2/`total_' {
    di `x'
    local y : word `x' of `filelist'
    append using "`y'", force
}

save eulfs, replace

Random graphs (29): Line plots with confidence bands

use mainstat age country refyear coeff if refyear == 2011 using eulfs.dta, replace

drop refyear // Drop unnecessary variable
generate retired = (mainstat == 4) if mainstat != . // Create dummy variable for retired

svyset [pweight = coeff] // Declare weight

// Create temporary files and postfile
tempname foo
tempname retage
postfile `foo' cntry age perc_retired perc_retired_ll perc_retired_ul using `retage' , replace

levelsof country, local(levels1)

foreach country of local levels1 {
  forvalues x = 57 (5) 77 {
      // Calculate proportions
    capture proportion retired if age == `x' & country == `country'
    *matrix list r(table)
    matrix coefs  = r(table)
    local perc    = coefs[1,2] * 100
    local perc_ll = coefs[5,2] * 100
    local perc_ul = coefs[6,2] * 100

    qui di `country' _skip(2) `x' _skip(2) `perc_ll' _skip(2) `perc_' _skip(2) `perc_ul'

    if "`perc'" != "" {    // Make sure that the loop doesn't break if 
    post `foo' (`country') (`x') (`perc') (`perc_ll') (`perc_ul') 
 }
  }
}

postclose `foo'

use `retage', clear         

drop if perc_retired == .   // Drop empty rows

  // Label variables
label define age 57 "55{c 150}59 y" 62 "60{c 150}64 y" 67 "65{c 150}69 y" ///
                 72 "70{c 150}74 y" 77 "75{c 150}79 y", modify
label value age age
label define cntry  1 "AT"  2 "BE"  3 "BG"  4 "CH"  5 "CY" ///
                    6 "CZ"  7 "DE"  8 "DK"  9 "EE" 10 "ES" ///                    
                   11 "FI" 12 "FR" 13 "GR" 14 "HU" 15 "IE" ///
                   16 "IS" 17 "IT" 18 "LT" 19 "LU" 20 "LV" ///
                   21 "MT" 22 "NL" 23 "NO" 24 "PL" 25 "PT" ///
                   26 "RO" 27 "SE" 28 "SI" 29 "SK" 30 "UK" ///
                 , modify
label value cntry cntry

  // Plot plots     
twoway (rarea perc_retired_ll perc_retired_ul age, fcolor(gs14)) ///
       (line perc_retired age) ///   
      , by(cntry, cols(3) legend(off) note("") ///
                  caption("{it:Source:} EU-LFS 2011" ///
                          "{it:Notes:} Gray areas denote 95 % confidence intervals"))) ///
        ytitle("% retired") xtitle("Age group") ///
        xlabel(57 (5) 77, val ang(v)) ///
        ysize(10) 

/* ,by(, note("")) suppresses the default "Graphs by" note) */
/* legend(off) needs to be in by(, legend(off)) to work */