Showing posts with label local. Show all posts
Showing posts with label local. Show all posts

Jun 19, 2017

Random graphs (99): Scatterplots

regress treatments index
local r2 = round(e(r2), .01) // Round coeff.
twoway (scatter treatments index) ///
       (lfit treatments index) ///
       (scatter treatments index if inlist(cntry, "IT", "AT", "UA", "BE", "DK", "RU"), msymbol(o)) ///
      , xtitle(ART comprehensive availability) ytitle("ART treatments" "per million women 15-44 y.") ///
        title("{bf:B}", justification(left) bexpand span) ///
        xlabel(0 (1) 9) ///
        text(1500 6.6 "Denmark") ///   
        text(1500 8.6 "Belgium") ///   
        text(470  0.4 "Italy") ///
        text(400  1.6 "Austria") ///
        text(120  7.5 "Ukraine") ///
        text(210  8.0  "Kazakhstan", placement(west)) ///
        text(150  8.0 "Russia", placement(east)) ///
     legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(7)) name(policy, replace)

Mar 21, 2016

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

Mar 16, 2016

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

Mar 10, 2016

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)

Jan 30, 2016

Random graphs (56): Intraclass Correlation Coefficients with confidence intervals


webuse productivity, clear

tempname foo
postfile `foo' str100 commandline_str str100 description_str icc icclb iccub N N_groups using "C:\Windows\Temp\test.dta", replace

qui mixed gsp || region: if year <= 1973
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1970 to 1973"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

qui mixed gsp || region: if year >= 1974 | year <= 1977
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1974 to 1977"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

qui mixed gsp || region: if year >= 1978 | year <= 1981
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1978 to 1981"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

qui mixed gsp || region: if year >= 1981 | year <= 1984
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1981 to 1984"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

qui mixed gsp || region: if year >= 1985 | year <= 1986
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1985 to 1986"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

postclose `foo'

use "C:\Windows\Temp\test.dta", clear
list

encode description_str, gen(description)
twoway (dot icc description, horizontal) ///
       (rcap icclb iccub description, horizontal)  ///
  , ylabel(1/5, val) ytitle("") legend(off) ///
    xlabel(0 (.1) 1) xtitle("ICC") name(figure, replace)
  
  
erase "C:\Windows\Temp\test.dta"

Sep 27, 2015

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


Sep 3, 2015

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



Aug 18, 2015

Random Graphs (50): Dot graphs with confidence intervals


// Create temporary object

tempname bqlrer

// Define postfile

postfile `bqlrer' count str100 commandline str100 descri str10 entity ///
                  prev prev_loci prev_hici differ differ_loci differ_hici ///
                  loed mide hied ratio_ ratio_loci ratio_hici n ///
                  using "results\results", replace

   // Count variable
local count = 0

  // Levels for loop
levelsof entity, local(entity_levels)
*di `entity_levels'

  // Loop: 1 round per entity
foreach Y of local entity_levels  {
     local entity_level = "`Y'"
     di "`entity_level'"
  *qui sum poorhealth if entity == "`Y'"
  *scalar proportion = r(mean)
  
  // Model 1: OV: PSRH lowest2, adjusted -- Prevalence rate only, education not in model
  qui logit poorhealth c_age c_female if entity == "`Y'"
  di _rc
  if _rc == 0 {
  *prev prev_loci prev_hici
  qui margins if entity == "`Y'", post
     scalar prev = 100 * _b[_cons]
     scalar prev_loci = 100 * (_b[_cons] - 1.96 * _se[_cons])
     scalar prev_hici = 100 * (_b[_cons] + 1.96 * _se[_cons])
  di `prev' 
  di `prev_loci' 
  di `prev_hici'
  }
  
  // Model 2: OV: PSRH lowest2, adjusted
     local count = `count' + 1
     local desc "DV: Poor health (bottom 2), IV: Education, Age, Sex"
     qui logit poorhealth i.education c_age c_female if entity == "`Y'"
  local commandline = e(cmdline)
  di _rc
  if _rc == 0 {
  
  // Differences Model 2
  qui margins r.education if entity == "`Y'"
  scalar differ      = -10 * _b[r3vs1.education]
  scalar differ_loci = -10 * (_b[r3vs1.education] + (1.96 * _se[r3vs1.education]))
  scalar differ_hici = -10 * (_b[r3vs1.education] - (1.96 * _se[r3vs1.education]))

  // Predicted probabilities Model 2
  qui margins i.education if entity == "`Y'", post
  matrix preds = r(table)
  scalar loed = preds[1,1]
     scalar mide = preds[1,2]
     scalar hied = preds[1,3]
     matrix drop preds
  
  // Ratios Model 2
  qui nlcom _b[1.education]/_b[3.education], post
  scalar ratio_     = _b[_nl_1]
  scalar ratio_loci = _b[_nl_1] - 1.96 * _se[_nl_1]
  scalar ratio_hici = _b[_nl_1] + 1.96 * _se[_nl_1]
  
  }
  else {
  scalar differ = -99
  scalar differ_loci = -99
  scalar differ_hici = -99
  scalar ratio_ = -99
   scalar ratio_loci = -99
     scalar ratio_hici = -99
  scalar loed = -99
  scalar mide = -99
  scalar hied = -99  
  }
     post `bqlrer' (`count') (`"`commandline'"') ("`desc'") ("`Y'") ///
                   (prev) (prev_loci) (prev_hici) ///
                   (differ) (differ_loci) (differ_hici) ///
                   (loed) (mide) (hied) ///
                   (ratio_) (ratio_loci) (ratio_hici) (e(N))
}

// Fit one model for entire EVS
// Model 1: OV: PSRH lowest2, adjusted -- Prevalence rate only, education not in model
qui logit poorhealth c_age c_female [pw = eu_weight] if survey == "EVS"
qui margins if survey == "EVS", post
scalar prev = 100 * _b[_cons]
scalar prev_loci = 100 * (_b[_cons] - 1.96 * _se[_cons])
scalar prev_hici = 100 * (_b[_cons] + 1.96 * _se[_cons])
di `prev' 
di `prev_loci' 
di `prev_hici'

// Model 2: OV: PSRH lowest2, adjusted
local count = `count' + 1
local desc "DV: Poor health (bottom 2), IV: Education, Age, Sex"
qui logit poorhealth i.education c_age c_female /*[pw = eu_weight]*/ if survey == "EVS"
local commandline = e(cmdline)

// Differences Model 2
qui margins r.education if survey == "EVS"
scalar differ      = -10 * _b[r3vs1.education]
scalar differ_loci = -10 * (_b[r3vs1.education] + (1.96 * _se[r3vs1.education]))
scalar differ_hici = -10 * (_b[r3vs1.education] - (1.96 * _se[r3vs1.education]))

  // Predicted probabilities Model 2
qui margins i.education if survey == "EVS", post
matrix preds = r(table)
scalar loed = preds[1,1]
scalar mide = preds[1,2]
scalar hied = preds[1,3]
matrix drop preds
  
  // Ratios Model 2
qui nlcom _b[1.education]/_b[3.education], post
scalar ratio_     = _b[_nl_1]
scalar ratio_loci = _b[_nl_1] - 1.96 * _se[_nl_1]
scalar ratio_hici = _b[_nl_1] + 1.96 * _se[_nl_1]

local Y = "EVS"

post `bqlrer' (`count') (`"`commandline'"') ("`desc'") ("`Y'") ///
              (prev) (prev_loci) (prev_hici) ///
              (differ) (differ_loci) (differ_hici) ///
              (loed) (mide) (hied) ///
              (ratio_) (ratio_loci) (ratio_hici) (e(N))
  
postclose `bqlrer'

// Create data set with FIPS codes and state names
preserve
tempfile codescheme
egen pickone = tag(entity)
keep if pickone
keep entity entity_num
save `codescheme'
restore

// Merge FIPS codes to results
use results\results, clear
merge 1:1 entity using `codescheme', keepusing(entity entity_num)
drop _merge

label var count "Generic counter"
label var commandline "Command line"
label var descri "Description of model fitted"
label var entity "Entity (string)"
label var prev   "Proportion poor health"
label var prev_loci "Proportion 95 % CI (lower)"
label var prev_hici "Proportion 95 % CI (higher)"
label var loed   "Proportion poor health lower educ."
label var mide   "Proportion poor health mid educ."
label var hied   "Proportion poor health higher educ."
label var differ "Difference lower educated - higher"
label var differ_loci  "Difference 95 % CI (lower)"
label var differ_hici  "Difference 95 % CI (higher)"
label var ratio_      "Ratio lower educated / higher"
label var ratio_loci  "Ratio 95 % CI (lower)"
label var ratio_hici  "Ratio 95 % CI (higher)"
label var n           "Size of entity sample"

// Save results, also as csv-file
save results\results, replace
outsheet using "results\results.csv", comma replace nolabel


use results\results, replace

// Create US identifier
gen usa = 0 
replace usa = 1 if regexm(entity, "US-")

replace entity = "{bf:Europe}" if entity == "EVS"

egen order_ratio = rank(-ratio_), unique 
labmask order_ratio, value(entity)

egen order_differ = rank(-differ), unique 
labmask order_differ, value(entity)

egen order_prev = rank(-prev), unique
labmask order_prev, value(entity)

twoway (dot prev order_prev if usa == 0, horizontal ndots(20)) ///
       (dot prev order_prev if usa == 1, horizontal ndots(20)) ///
    (rspike prev_loci prev_hici order_prev, horizontal) ///
    , ///
    legend(label(1 "Europe") label(2 "US") label(3 "95% CI") pos(1) ring(0)) ///
    ylabel(1/97, valuelabels ang(h) labsize(*.65)) ///
    ytitle("") ///
    xtitle("Prevalence poor health") ///
    name(prevalence, replace) ///
    xsize(2.7) ysize(10)

generate ratio_loci_p = ratio_loci               // Shorten CI's to make graph fit better
replace  ratio_loci_p =  0 if ratio_loci_p <  0  
generate ratio_hici_p = ratio_hici
replace  ratio_hici_p = 6 if ratio_hici_p > 6
 
twoway (dot ratio_ order_ratio if usa == 0, horizontal ndots(20) msymbol(smplus)) ///
       (dot ratio_ order_ratio if usa == 1, horizontal ndots(20) msymbol(smx)) ///
    (rspike ratio_loci_p ratio_hici_p order_ratio, horizontal) ///
    , ///
    legend(label(1 "Europe") label(2 "US") label(3 "95 % CI") pos(1) ring(0)) ///
    ylabel(1/97, valuelabels ang(h) labsize(*.65)) ///
    xscale(range(0 6)) xlabel(0 1 3 5)  /// 
    ytitle("") ///
    xtitle("Relative" "inequalities") ///
    xline(1) ///
    xline(1.81019, lpattern(dash)) ///
    name(relative, replace) ///
    xsize(2.7) ysize(10)

generate differ_loci_p = differ_loci              // Shorten CI's to make graph fit better
replace  differ_loci_p = -1 if differ_loci_p < -1
    
twoway (dot differ order_differ if usa == 0, horizontal ndots(20) msymbol(smplus)) ///
       (dot differ order_differ if usa == 1, horizontal ndots(20) msymbol(smx)) ///
       (rspike differ_loci_p differ_hici order_differ, horizontal) ///
    , ///
    legend(label(1 "Europe") label(2 "US") label(3 "95 % CI") pos(1) ring(0)) ///
    xscale(range(0 30)) xlabel(0(10)30)  /// 
    ylabel(1/97, valuelabels ang(h) labsize(*.65)) ///
    ytitle("") ///
    xtitle("Absolute" "inequalities") ///
    xline(6.8754, lpattern(dash)) ///
    name(absolute, replace) ///
    xsize(2.7) ysize(10)

twoway (dot prev order_prev if usa == 0, horizontal ndots(20) msymbol(smplus)) ///
       (dot prev order_prev if usa == 1, horizontal ndots(20) msymbol(smx)) ///
    (rspike prev_loci prev_hici order_prev, horizontal) ///
    , ///
    legend(label(1 "Europe") label(2 "US") label(3 "95% CI") row(1) pos(12) size(vsmall) region(lwidth(vthin) lcolor(black)) bmargin(tiny) colgap(*.3)) ///
    ylabel(1/97, valuelabels ang(h) labsize(*.65)) ///
    ytitle("") ///
    xtitle("Prevalence" "poor health") ///
    xline(11.26938, lpattern(dash)) ///
    name(prevalence_leg, replace) ///
    xsize(2.7) ysize(10)    
    
// 8.27 × 11.69
*graph combine prevalence absolute relative, xsize(8.27) ysize(10) col(3) name(comb, replace)

grc1leg prevalence_leg absolute relative, xsize(8.27) ysize(10) col(3) imargin(small) legendfrom(prevalence_leg) name(oneleg, replace) pos(6) span // ignores size-commands

graph display oneleg, xsize(6.27) ysize(9.69)  // Thus redraw so that size commands take effect

May 24, 2015

Random graphs (48): Interaction plot with overlaid data points


// Simulate data
clear
set seed 1
set obs 500

generate  e  = 0 + (500 - 0) * runiform()     // To generate random variates over the
generate  x1 = 0 + (800 - 0) * runiform()     // interval [a,b), a+(b-a)*runiform()

generate  x2 = round(0 + (1 - 0) * runiform()) // Binary variable
generate  y  = (x1 + x2 + (x1*x2) + e) / 1000

// Calculations for Aiken & West (1991) style plot
regress y c.x1##x2

qui sum x1
local x1_minsd = r(mean) - (2*r(sd))
local x1_mean  = r(mean)
local x1_plusd = r(mean) + (2*r(sd))


// Calculate predicted values for plot
margins, at(c.x1 = (`x1_minsd' `x1_mean' `x1_plusd') ///
            c.x2 = (0 1)) vsquish
// Plot
qui marginsplot, recastci(rarea) ciopts(color(gs10)) ///
   title("Interaction plot with overlaid data points") ///
   ytitle("y") ///  
   ylabel(, format(%6.1f)) ///
   xtitle("") ///
   plotopts(msymbol(none)) ///        // Turn off markers
   plot1opts(lpattern(longdash)) ///  // Define line types here
   plot2opts(lpattern(solid)) ///
   addplot(scatter y x1 ///
         , symbol(o) ///
           xlabel(`x1_minsd' "-2 SD" ///       // The addplot seems to override
                  `x1_mean' "Average x1" ///   // the regular axis label 
                  `x1_plusd' "+2 SD") ///      // command
           legend(subtitle(x2) ///             // And the legend command
                  order(3 "x2 = 0" 4 "x2 = 1" 2 "95 % CI"))) ///    
   legend(pos(5) ring(0)) ///                  // But not completely
   name(plot, replace)

May 6, 2015

Random graphs (46): Plotting the contours of interactions

Next to the classical layout for interaction plots, it is also possible to make use of contour plots for visualizing interactions.
// Simulate data
clear
set seed 1
set obs 100

generate  e  = 0 + (100000 - 0) * runiform()  // To generate random variates over the
generate  x1 = 0 + (600 - 0) * runiform()     // interval [a,b), a+(b-a)*runiform()
generate  x2 = 0 + (600 - 0) * runiform()  
generate  y  = (x1 + x2 + (x1*x2) + e) / 1000

// Calculations for Aiken & West (1991) style plot
qui regress y c.x1##c.x2

qui sum x1
local x1_minsd = r(mean) - r(sd)
local x1_mean  = r(mean)
local x1_plusd = r(mean) + r(sd)

qui sum x2
local x2_minsd = r(mean) - r(sd)
local x2_mean  = r(mean)
local x2_plusd = r(mean) + r(sd)

// Calculate predicted values for plot
margins, at(c.x2 = (`x2_minsd' `x2_mean' `x2_plusd') ///
            c.x1 = (`x1_minsd' `x1_mean' `x1_plusd')) vsquish
// Plot
qui marginsplot, recastci(rarea) ciopts(color(gs10)) ///
   xlabel(`x2_minsd' "-1 SD" `x2_mean' "Average x2" `x2_plusd' "+1 SD") ///
   title("Aiken & West (1991)-style plot for interactions") ///
   ytitle("Predicted y") ///
   xtitle("") ///
   plotopts(msymbol(none)) ///        // Turn off markers
   plot1opts(lpattern(longdash)) ///  // Define line types here
   plot2opts(lpattern(solid)) ///
   plot3opts(lpattern(shortdash)) ///
   legend(subtitle(x1) ///
   order(4 "- 1 SD" 5 "Average x1" 6 "+ 1 SD" 2 "95% CI")) ///
   name(lines, replace)
        // For some reason, marginsplot ignores the -label option-, therefore
        // -order- is used here

// Calculations for contour plot
qui regress y c.x1##c.x2
predict y_hat
// Plot
qui twoway contour y_hat x1 x2, ///
       title("Contour plot for interactions") ///
    ztitle("Predicted y") ///
    name(contours, replace)
// Combine figures
graph combine lines contours, col(1) ysize(8)

Jul 9, 2013

Morgan and Winship's (2007) example for bias due to conditioning on a collider in Stata

Morgan and Winship (2007: p. 66) illustrate Pearl's (2009) concern about conditioning on a collider variable using a simple example.

The general problem of conditioning on a collider is as follows. Consider three variables A, B, and C, with both A and B being causes of C: A → C ← B. (Formally, any variable C that has two arrows pointing to it along a given path is a collider.) Unlike a confounder (an uncontrolled common cause of A and B), a collider does not induce a zero-order correlation between A and B. However, when handled inappropriately, a collider can induce a conditional correlation between A and B.

Morgan and Winship's example shows just that: A college admits applicants based on their SAT scores and ratings of their motivation based on an interview. Those in the top 15 per cent of the sum of SAT and motivation ratings are being admitted. SAT scores and motivation ratings are largely uncorrelated.

// Generate and label two variables
drawnorm sat motivation, ///
         n(250) ///
         means(.007, -.053) ///
         sds(1.01, 1.02) ///
         corr(1, .035, 1) cstorage(lower) ///
         clear seed(1)
label var sat        "SAT"
label var motivation "Motivation"

// Only the 15 per cent at the top are admitted
gen admission_sc = sat + motivation
_pctile admission_sc, percentiles(85)

gen admission = (admission_sc > r(r1))
label var admission "Admission status"
label define admission 1 "Admitted applicant" ///
                       0 "Rejected applicant"
label val admission admission
drop admission_sc

// Plot as in Morgan and Winship, p. 67:
twoway (scatter motivation sat if admission == 1) ///
       (scatter motivation sat if admission == 0) ///
       , ///
       legend(label(1 "Admitted applicants") ///
              label(2 "Rejected applicants") ///
              pos(5) ring(0)) ///
       ylabel(none) xlabel(none) ///
       name(collider1, replace)
// Enhanced plot with fitted lines and correlations
quietly cor motivation sat
local r_overall = round(r(rho), .01)
quietly cor motivation sat if admission == 1
local r_admitted = round(r(rho), .01)
quietly cor motivation sat if admission == 0
local r_rejected = round(r(rho), .01)
    
twoway (scatter motivation sat if admission == 1) ///
       (scatter motivation sat if admission == 0) ///
       (lfit motivation sat) ///
       (lfit motivation sat if admission == 1) ///
       (lfit motivation sat if admission == 0) ///
       , ///
       legend(label(1 "Admitted applicants") ///
              label(2 "Rejected applicants") ///
              label(3 "Overall fit, {it:r} = `r_overall'") ///
              label(4 "Fit for admitted, {it:r} = `r_admitted'") ///
              label(5 "Fit for rejected, {it:r} = `r_rejected'") ///
              pos(5) ring(0)) ///
       ylabel(none) xlabel(none) ytitle("Motivation")

What the Figures show is that the very small correlation between motivation and SAT score for the overall group turns out to be much larger when conditioning for admission status.

This also shows in an OLS regression:

regress motivation sat, beta
estimates store m1
regress motivation sat admission, beta
estimates store m2

estimates table m1 m2, b(%7.2f) se(%7.2f) stats(N) label

----------------------------------------------
                Variable |   m1        m2     
-------------------------+--------------------
                     SAT |    0.10     -0.21  
                         |    0.06      0.06  
        Admission status |              1.60  
                         |              0.17  
                Constant |   -0.12     -0.36  
                         |    0.06      0.06  
-------------------------+--------------------
                       N |     250       250  
----------------------------------------------
                                  legend: b/se

Cole et al. (2010) present additional illustrations for this problem.

References


Cole, Stephen R., Robert W. Platt, Enrique F. Schisterman, Haitao Chu, Daniel Westreich, David Richardson, and Charles Poole. 2010. "Illustrating Bias Due to Conditioning on a Collider." International Journal of Epidemiology 39(2):417-420. doi: 10.1093/ije/dyp334

Morgan, Stephen L., and Christopher Winship. 2007. Counterfactuals and Causal Inference. Methods and Principles for Social Research. Cambridge University Press.

Pearl, Judea. 2009. Causality. Models, Reasoning, and Inference, 2nd ed. Cambridge University Press.