Centering predictor variables in OLS regression

clear 

clear 

set seed 1
set obs 1000

// Generate variables
generate female = (runiform() > .48)
generate health = 1 - female * rnormal(0.5, 1) + 0.1 * rnormal(3, 1)

// Recode predictor variable
center female
recode female  (0 = -1), gen(e_female)

eststo clear
eststo: regress health female
eststo: regress health c_female
eststo: regress health e_female

esttab, cells(b(fmt(2))) rename(c_female female e_female female) ///
        mtitle("Dummy" "Centering" "Effect") nonumber ///
		coeflabel(female "Female (ref. male)" _cons "Intercept") ///
		varwidth(18) collabel("") stats(N, fmt(%9.0gc)) ///
		title(Comparing different types of centering)

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 (135): Scatterplot with OLS 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

// 1) Basic scatterplot
twoway (scatter questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Number of questions answered correctly") ///
        name(figure1, replace) ylabel(0 (5) 15) xsize(4) ysize(4)

// 2) Scatterplot with regression line
twoway (scatter questions hours) ///
       (lfit questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 (5) 15) legend(off) name(figure2, replace) xsize(4) ysize(4)

// 3) Scatterplot with regression line and equation
regress questions hours
local intercept = round(_b[_cons], .1)
local x = round(_b[hours], .1)

twoway (scatter questions hours) ///
       (lfit questions hours), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(14 2 "y = `intercept' + `x' + e", size(large)) ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 (5) 15) legend(off) name(figure3, replace) xsize(4) ysize(4)

// 4) Scatterplot with regression line and equation and labels for components
twoway (scatter questions hours) ///
       (lfit questions hours, lpattern(solid) range(0 5)) ///
       (function y = 1.6 + 2.5, range(1 2)) ///
       (function y = 2, range(4.1 6.6) horizontal lpattern(solid) lcolor(red)) ///
       (scatteri 1.6 0 (3) "Intercept", msymbol(o) mlabcolor(red)) ///
       (scatteri 5 2 "Slope", msymbol(i) mlabcolor(red)), ///
        xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
        text(14 2 "y = `intercept' + `x' x + e", size(large)) ///
        ytitle("Number of questions answered correctly") ///
        ylabel(0 1.6 5 10 15) legend(off) name(figure4, replace) xsize(4) ysize(4)
    
graph combine figure1 figure2 figure3 figure4, col(2) xsize(8) ysize(8) altshrink name(figures14, replace)

Dummy variable adjustment for missing values in Stata

This piece of code replicates Table 3.1 in Allison (2002).

clear
set seed 1

// Generate data
set obs 10000
drawnorm x z, ///
         corr(1, .5, 1) cstorage(lower) 
generate e = rnormal()
generate y = x + z + e

// Drop 1/2 of values from z
generate d = (runiform() > . 5)
generate zstar1 = z if d
replace  zstar1 = . if !d

// Substitute missing values
qui sum zstar1
generate zstar2 = zstar1
replace  zstar2 = r(mean) if !d

eststo clear
eststo: regress y x z
eststo: regress y x zstar1
eststo: regress y x zstar2 d

esttab using test.tex, b(2) not nostar nocons rename(zstar1 z zstar2 z) ///
        mtitles("Full data" "Listwise deletion" "Dummy variable adjustment") ///
        title(Regression in Simulated Data for Three Methods) replace booktabs 

Reference

Allison, Paul D. 2002. Missing Data. Sage. doi: 10.4135/9781412985079

Omitted variable bias

clear 
set obs 10000
set seed 1

// Simulate data
generate x  = rnormal(0,1)                     // Exogenous variable
generate w  = rnormal(0,1)                     // Instrumental variable
generate u  = rnormal(0,1)                     // Omitted variable
generate e1 = rnormal(0,1)                     // Outcome eqation error
generate e2 = rnormal(0,1)                     // Endogenous regressor equation error
generate y2 =      x  + .2 * w + .5 * u + e2   // Endogenous regressor equation
generate y1 = .5 * y2 + .5 * x + .5 * u + e1   // Outcome equation

// Fit models
eststo clear
eststo: regress        y1 y2 x u
eststo: regress        y1 y2 x
eststo: ivregress 2sls y1 x (y2 = w)

coefplot est1 est2 est3, xscale(alt) xtitle(Coefficient) ///
                         xline(.5) ///
                         coeflabel(_cons = "Intercept" ///
                                   u     = "Omitted variable" ///
                                   x     = "Exogenous predictor" ///
                                   y2    = "Endogenous predictor") ///
                         legend(order(2 "Without omitted variable bias" ///
                                      4 "With omitted variable bias" ///
                                      6 "2SLS estimate") pos(11) ring(0) col(1))

Random graphs (122): Specification checks of OLS regression

// Simulate data set
clear
set obs 1000
set seed 1
generate x = runiform()
generate z = runiform()
generate u = rnormal()
generate y1 =      0 +      x   + 0 * z + u
generate y2 =      0 + 10 * x^2 + 0 * z + u
generate y3 = exp(-1 +      x   + 0 * z + u)

// Correctly specified model
regress y1 x z
rvfplot,   name(rvf1, replace)  title(Fitted values vs. residuals) nodraw
rvpplot x, name(rvpx1, replace) title(Predictor variable x vs. residuals) nodraw
rvpplot z, name(rvpz1, replace) title(Predictor variable z vs. residuals) nodraw
graph combine rvf1 rvpx1 rvpz1, col(3) xsize(7) ysize(2) /// 
                                title(Residual plots of correctly specified model) ///
                                name(correct, replace) nodraw

// Missing out on quadratic relationship
regress y2 x z
rvfplot,   name(rvf2, replace)  title(Fitted values vs. residuals) nodraw
rvpplot x, name(rvpx2, replace) title(Predictor variable x vs. residuals) nodraw
rvpplot z, name(rvpz2, replace) title(Predictor variable z vs. residuals) nodraw
graph combine rvf2 rvpx2 rvpz2, col(3) xsize(7) ysize(2) ///
                                title(Residual plots of model missing quadratic term) ///
                                name(quadractic, replace) nodraw

// Missing out on exponential relationship
regress y3 x z
rvfplot,   name(rvf3, replace)  title(Fitted values vs. residuals) nodraw
rvpplot x, name(rvpx3, replace) title(Predictor variable x vs. residuals) nodraw
rvpplot z, name(rvpz3, replace) title(Predictor variable z vs. residuals) nodraw
graph combine rvf3 rvpx3 rvpz3, col(3) xsize(7) ysize(2) ///
                                title(Residual plots of model missing exponential shape) ///
                                name(exponential, replace) nodraw

graph combine correct quadractic exponential, col(1) xsize(7) ysize(6) altshrink

Multiple imputation in Stata

clear

// Generate data
set obs 10000
generate e  = rnormal()
drawnorm x1 x2, ///
         corr(1, .4, 1) cstorage(lower) 
generate y = 1 + 1 * x1 + 1 * x2 + 25 * e
drop e

// Model with data-generating process
eststo clear
eststo: regress y x1 x2

// Drop cases conditional on other covariate
replace x1 = . if rnormal() > .1 & x2 > -1

// Model with missing data
eststo: regress y x1 x2

// Impute data
mi set mlong                  // Declare data format
mi register imputed x1        // Declare variable to be imputed
mi impute chained (regress) x1 = y x2, add(10) // Declare imputation model

// Model with imputed data
eststo: mi estimate, post:  regress y x1 x2 // post option is important for -esttab-

// Table
esttab, r2 se mtitles("DGP" "Missing data" "Imputed")

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

Measurement error in the outcome vs. in the predictor

clear
set seed 1
set obs 5000

// Create variables
gen e = rnormal()
gen x_true = rnormal()
gen y_true = 1 + 2 * x_true + 3 * e
gen y_observed = y_true + e
gen x_observed = x_true + e

label var y_true "Y"
label var x_true "X"
label var y_observed "Y w/ error"
label var x_observed "X w/ error"

// Loop over variables
foreach y of varlist y_observed y_true {
 foreach x of varlist x_observed x_true {
  regress `y' `x'
  local  b = round( _b[`x'], .01)
  local se = round(_se[`x'], .01)  
  twoway (scatter `y' `x', msymbol(p)) ///
         (lfit `y' `x') ///
        , name(`y'`x', replace) ///
          ytitle(`: variable label `y'') ///
          xtitle(`: variable label `x'') ///
          title("{it:B} = `b', {it:SE} = `se'") ///
          legend(off) nodraw
 }
}

// Combine plots
graph combine y_truex_true ///
              y_truex_observed ///
              y_observedx_true ///
              y_observedx_observed, ///
              xcommon ycommon col(2) row(2) ///
              title("Measurement error in the outcome vs. in the predictor")

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

Analyzing censored data with the Tobit model

This replicates Table 2.2 in Breen (1996). 

clear

// Simulate data
set obs 2000
generate ui = rnormal(0, 2)
generate xi = rnormal()
generate yi_star = 1 + 2*xi + ui
drop ui

// Fit OLS model
regress yi_star xi
eststo ols1
estadd scalar sigma = e(rmse)

// Truncate variable
generate yi = yi_star if yi_star > 0
replace  yi = 0       if yi_star <= 0
recode   yi (0 = .), gen(yi_h) 

// (A) OLS (using all observations on y
//     including y1 = 0)
regress yi xi
eststo ols2
estadd scalar sigma = e(rmse)

// (B) OLS (yi > 0) 
regress yi xi if yi >0
eststo ols3
estadd scalar sigma = e(rmse)

// (C) Heckman 2-step
heckman yi_h xi, select(xi) twostep 
eststo heckman

// (D) Tobit
tobit yi xi, ll(0)
eststo tobit
estadd scalar sigma = _b[sigma:_cons]

// Table 2.2


esttab ols2 ols3 heckman tobit ols1, b(3) se(3) nostar stat(sigma) ///
       mtitles("(A) OLS incl. yi = 0" ///
               "(B) OLS yi > 0" ///
               "(C) Heckman 2-step" ///
               "(D) Tobit" ///
               "OLS yi_star") ///
       coeflabel(_cons "alpha" xi "beta") ///
       collabels() ///
       drop(mills:lambda sigma:_cons) ///
       order(_cons xi) ///
       unstack ///
       modelwidth(20) nonumber

Reference

Breen, Richard. 1996. Regression Models. Censored, Sample Selected, or Truncated Data. Sage. doi: 10.4135/9781412985611

Simulating multilevel data

clear

set obs 200              // Number of level-2 units
gen j = _n              // ID for level-2 units
gen c_j = rnormal(0,1)  // Level-2 covariate
gen u_j = rnormal(0,1)  // Level-2 error term
expand 100              // Number of level-1 units per level-2 unit
bysort j: gen i = _n    // ID for level-1 units
gen x_ij = rnormal(0,1) // Level-1 covariate
gen e_ij = rnormal(0,1) // Level-1 error term
gen y_ij = 1 + 1 * x_ij + 1 * c_j + u_j + e_ij // Regression equation

mixed y_ij x_ij c_j

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

Interactions with and without main effects

set scheme s1mono
clear

set obs 10000
set seed 1

gen a = rnormal()
gen b = rnormal()
gen e = 2*rnormal()

gen y = 2 + 3*a + 4*b + 5*a*b + e

summarize y a b

estimates clear
regress y c.a##c.b
estimates store M1
regress y c.a  c.b
estimates store M2
regress y      c.b c.a#c.b
estimates store M3
regress y c.a      c.a#c.b
estimates store M4
regress y          c.a#c.b
estimates store M5

estimates restore M1
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M1, replace) noci ///
             title("M1: Main terms and interaction") legend(col(1))

estimates restore M2
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M2, replace) noci ///
             title("M2: Main terms, no interaction") legend(off)

estimates restore M3
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M3, replace) noci ///
             title("M3: Interaction, one main term missing") legend(off)

estimates restore M4
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M4, replace) noci ///
             title("M4: Interaction, other main term missing") legend(off)

estimates restore M5
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M5, replace) noci ///
             title("M5: Interaction term, both main terms missing") legend(off)

grc1leg M1 M2 M3 M4 M5, col(2) ysize(11) xsize(8) ring(0) pos(5) ycommon
estimates table M1 M2 M3 M4 M5, b(%9.2f) stat(F r2_a)

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)

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.