Scaling predictors for models with interactions

The scaling of predictors is an important issue for the interpretation of models with interaction terms. The Table below presents a simple OLS regression model of happiness regressed on age, sex, and their interaction.

When none of the predictors are standardized and 0 is not a meaningful value on a predictor, interpretation of the regression coefficients can become awkward. For instance, in the model with none of the predictors standardized, the coefficient for female – .348 – states that women score .348 points higher in terms of happiness than men among those who are 0 years old. The coefficient for age – -.004 – indicates that with every additional year, men's happiness decreases by .004 points. The coefficient for the interaction term – -.008 – denotes that women's happiness decrease more than men's with each additional year of age. The intercept is the average happiness of men who are 0 years old.

A first step towards giving the coefficients from an interaction model an easier interpretation is to standardize continuous variables. For the second model shown in the Table, age was standardized by subtracting its mean and dividing it by its standard deviation. Thus, for the standardized variable, the value of 0 refers to the average age, and a change of 1 refers to a change of one standard deviation. The age coefficient (-.07) now indicates that for every additional standard deviation of age, men's happiness decreases by .07. For men and women of average age, the happiness differential is .02. Women's happiness decreases by .15 per standard deviation of age stronger than men's. The intercept shows the average happiness of average-age men.

While it is a rather controversial technique, it is in principle also possible to standardize the dummy variable indicating the difference between men and women. Then, the intercept is the overall mean of happiness. However, an increase of 1 in the gender coefficient becomes difficult to interpret, as there is no such thing as a one-standard deviation increase in gender. When the dummy variable is approximately evenly distributed, then the coefficient reflects half of the difference in the outcome between the two groups.

Gelman and Hill (2007, p. 56) suggest to standardize variables by two rather than one standard deviation to keep binary and continuous variables on a roughly common scale. An increase of one in age then refers to the difference between one standard deviation below the mean and one standard deviation above the mean; an increase of one in age reflects roughly the difference between men and women.

As the plot above shows, none of these different scaling decisions affect the form of the interaction.

Reference

Gelman, Andrew, and Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.

use agea gndr happy using ESS1e06.3_F1.dta, clear

// Set up variables
rename agea age
recode gndr (1 = 0) (2 = 1) (.a = .), gen(female)
label var female "Female sex (ref. male)"
label var age    "Age"
label var happy  "Happiness"

// Listwise deletion
drop if missing(female, age, happy)

// Standardize variables
qui sum female
gen female_1z = (female - r(mean)) /   r(sd)
gen female_2z = (female - r(mean)) / (2*r(sd))
qui sum age 
gen age_1z = (age - r(mean)) /   r(sd)
gen age_2z = (age - r(mean)) / (2*r(sd))

// Fit models
reg happy c.age##female
estimates store unstandardized

reg happy c.age_1z##female
estimates store cont_standardized

reg happy c.age_1z##c.female_1z
estimates store standardized_1z

reg happy c.age_2z##c.female_2z
estimates store standardized_2z

// Create table
esttab unstandardized cont_standardized standardized_1z standardized_2z, ///
       rename(age_1z age ///
              age_2z age ///
              female_1z 1.female ///
              female_2z 1.female ///
              c.age_1z#c.female_1z 1.female#c.age ///
              c.age_2z#c.female_2z 1.female#c.age ///
              1.female#c.age_1z  1.female#c.age) ///
       drop(0b.female 0b.female#co.age 0b.female#co.age_1z) ///
       coeflabel(age "Age" ///
                 1.female "Female (ref. male)" ///
                 1.female#c.age "Age X Female" ///
                 _cons "Intercept") ///
       mtitles("Unstandardized" "Age standardized (1 SD)" ///
               "Standardized (1 SD)" "Standardized (2 SD)") ///
       nonumbers varwidth(18)  modelwidth(24) ///
       se r2 obslast

// Plot models
estimates restore unstandardized

qui sum age 
local age_minsd = r(mean) - r(sd)
local age_mean  = r(mean)
local age_plusd = r(mean) + r(sd)

qui margins, at(c.age = (`age_minsd' `age_mean' `age_plusd') ///
                female = (0 1)) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
             xlabel(`age_minsd' "-1 SD" `age_mean' "Average age" `age_plusd' "+1 SD") ///
             title("Unstandardized variables", size(*.9)) ///
             ytitle("Predicted happiness") ///
             xtitle("") ///
             plotopts(msymbol(none)) ///        // Turn off markers
             plot1opts(lpattern(dash)) ///  // Define line types here
             plot2opts(lpattern(solid)) ///
             legend(subtitle(Gender, size(*.8)) ///
                    order(3 "Males" 4 "Females") ring(0) pos(7) ///
                    size(*.8) symysize(*.7) symxsize(*.7)) ///
             name(unstandardized, replace)

estimates restore cont_standardized

qui margins, at(c.age_1z = (-1 0 1) ///
                female = (0 1)) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
             xlabel(-1 "-1 SD" 0 "Average age" 1 "+1 SD") ///
             title("Age standardized at 1 SD", size(*.9)) ///
             ytitle("Predicted happiness") ///
             xtitle("") ///
             plotopts(msymbol(none)) ///        // Turn off markers
             plot1opts(lpattern(dash)) ///  // Define line types here
             plot2opts(lpattern(solid)) ///
             legend(subtitle(Gender, size(*.8)) ///
                    order(3 "Males" 4 "Females") ring(0) pos(7) ///
                    size(*.8) symysize(*.7) symxsize(*.7)) ///
             name(cont_standardized, replace)			 
			 
estimates restore standardized_1z

qui sum female
local fem_hi = (1 - r(mean)) /   r(sd)
local fem_lo = (0 - r(mean)) /   r(sd)

qui margins, at(c.age_1z = (-1 0 1) ///
                c.female_1z = (`fem_lo' `fem_hi')) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
             xlabel(-1 "-1 SD" 0 "Average age" 1 "+1 SD") ///
             title("All variables standardized at 1 SD", size(*.9)) ///
             ytitle("Predicted happiness") ///
             xtitle("") ///
             plotopts(msymbol(none)) ///        // Turn off markers
             plot1opts(lpattern(dash)) ///  // Define line types here
             plot2opts(lpattern(solid)) ///
             legend(subtitle(Gender, size(*.8)) ///
                    order(3 "Males" 4 "Females") ring(0) pos(7) ///
                    size(*.8) symysize(*.7) symxsize(*.7)) ///
             name(standardized_1z, replace)
			 
estimates restore standardized_2z

qui sum female
local fem_hi = (1 - r(mean)) /  (2*r(sd))
local fem_lo = (0 - r(mean)) /  (2*r(sd))

qui margins, at(c.age_2z = (-.5 0 .5) ///
                c.female_2z = (`fem_lo' `fem_hi')) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
             xlabel(-.5 "-1 SD" 0 "Average age" .5 "+1 SD") ///
             title("All variables standardized at 2 SD", size(*.9)) ///
             ytitle("Predicted happiness") ///
             xtitle("") ///
             plotopts(msymbol(none)) ///        // Turn off markers
             plot1opts(lpattern(dash)) ///  // Define line types here
             plot2opts(lpattern(solid)) ///
             legend(subtitle(Gender, size(*.8)) ///
                    order(3 "Males" 4 "Females") ring(0) pos(7) ///
                    size(*.8) symysize(*.7) symxsize(*.7)) ///
             name(standardized_2z, replace)			 

graph combine unstandardized cont_standardized standardized_1z standardized_2z, ///
      row(2) col(2) ysize(7) name(combined, replace)

Random graphs (23): Between and within regression

clear

input j i xij x_j yij y_j
// Data from Snijders and Bosker (1999, Table 3.2):
1 1 1 2 5 6
1 2 3 2 7 6
2 1 2 3 4 5
2 2 4 3 6 5 
3 1 3 4 3 4
3 2 5 4 5 4
4 1 4 5 2 3
4 2 6 5 4 3
5 1 5 6 1 2
5 2 7 6 3 2
end


// Plot as Figure 3.4 of Snijders ans Bosker (1999):
twoway (scatter yij xij) ///
       (lfit yij xij) ///
       (lfit y_j x_j) /// 

       (lfit yij xij if j == 1, lpattern(dash_dot)) ///

       (lfit yij xij if j == 2, lpattern(dash_dot)) ///
       (lfit yij xij if j == 3, lpattern(dash_dot)) ///
       (lfit yij xij if j == 4, lpattern(dash_dot)) ///
       (lfit yij xij if j == 5, lpattern(dash_dot)) ///
     , legend(label(2 "Total regression") ///
       label(3 "Between regression") ///
       label(4 "Within regression") ///
       order(2 3 4) ///
       pos(1) ring(0)) ///
       xtitle(X) ytitle(Y) ///
       title("Within, between, and total relations") /// 

       xlabel(none) ylabel(none) ///
       name(wbt_legend, replace)
    

// Lines directly labeled
twoway (scatter yij xij) ///
       (lfit yij xij) ///
       (lfit y_j x_j) ///
       (lfit yij xij if j == 1, lpattern(dash_dot)) ///
       (lfit yij xij if j == 2, lpattern(dash_dot)) ///
       (lfit yij xij if j == 3, lpattern(dash_dot)) ///
       (lfit yij xij if j == 4, lpattern(dash_dot)) ///
       (lfit yij xij if j == 5, lpattern(dash_dot)) ///
     , legend(off) ///
       text(1.8 6.5 "Between") ///
       text(3.0 7.5 "Total") ///
       text(7.0 3.5 "Within") ///
       text(5.0 5.5 "Within") ///
       xtitle(X) ytitle(Y) ///
       title("Within, between, and total relations") ///
       xlabel(none) ylabel(none) ///
       name(wbt_nolegend, replace)

graph combine wbt_legend wbt_nolegend, col(1) ysize(8) xcommon ycommon
// ysize(4) and xsize(5.5) are defaults

Reference

 Snijders, Tom, and Roel Boskers. 1999. Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling. Sage.

Random graphs (19): Simple interaction plots with -marginsplot- for continuous variables

sysuse auto, clear

regress mpg c.weight##c.gear_ratio

// Obtain means and standard deviation for plots
qui sum weight
local w_minsd = r(mean) - r(sd)
local w_mean  = r(mean)
local w_plusd = r(mean) + r(sd)

qui sum gear_ratio
local g_minsd = r(mean) - r(sd)
local g_mean  = r(mean)
local g_plusd = r(mean) + r(sd)

// Calculate predicted values for plot
margins, at(c.gear_ratio = (`g_minsd' `g_mean' `g_plusd') ///
            c.weight = (`w_minsd' `w_mean' `w_plusd')) vsquish
// Plot
marginsplot, recastci(rarea) ciopts(color(gs10)) ///
   title("Weight moderates the gear ratio{char 0150}mileage association") ///
   xlabel(`g_minsd' "-1 SD" `g_mean' "Average gear ratio" `g_plusd' "+1 SD") ///
   ytitle("Predicted mileage") ///
   xtitle("Gear ratio") ///
   plotopts(msymbol(none)) ///        // Turn off markers
   plot1opts(lpattern(longdash)) ///  // Define line types here
   plot2opts(lpattern(solid)) ///
   plot3opts(lpattern(shortdash)) ///
   legend(subtitle(Weight) ///
   order(4 "- 1 SD" 5 "Mean" 6 "+ 1 SD" 2 "95% CI"))
        // For some reason, marginsplot ignores the -label option-, therefore
        // -order- is used here
gr export plot.png, as(png) replace

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.