Random graphs (53): Visualizing multilevel models

use cntry happy mnactic agea gndr using ESS1e06_4.dta, clear

// Prepare variables
recode agea (999 = .), gen(age)
recode gndr (1 = 0 "Male") (2 = 1 "Female") (9 = .), gen(female)
recode mnactic (77 88 99 = .), gen(active)
label var female "Sex (ref. Male)"
label var age    "Age in decades"
label var happy  "Happiness"
label var active "Labor market status"
label val active mnactic

center happy age // Center outcome to make intercept smaller for plot
replace c_age = c_age / 10

// Listwise deletion
drop if missing(female, c_age, c_happy, active, cntry)
drop happy gndr agea mnactic // Unecessary variables can go

// Fit model
mixed c_happy c_age i.female i.active || cntry: c_age, ml cov(uns)

// Save random part
local var_age = round(exp(_b[lns1_1_1:_cons])^2, .0001)
local var_int = round(exp(_b[lns1_1_2:_cons])^2, .0001)
local covaria = round(tanh(_b[atr1_1_1_2:_cons]) * ///
                       exp(_b[lns1_1_1:_cons]) * ///
        exp(_b[lns1_1_2:_cons]), .0001)

// Plot fixed part
coefplot, xline(0) ///
xtitle(" " "Estimates and 95% CI's") ///
scheme(s1mono) ///
ciopts(recast(rcap)) ///
coeflabels(_cons = "{bf:Intercept}" ///
           c_age = "{bf:Age} in decades" ///
     8.active = "Housework") ///  // Shorten label
headings(c_age = " " ///
         1.female  = "{bf:Sex} ({it:ref.} Male)" ///
         2.active = `""{bf:Labor market status}" "({it:ref.} Paid work)""' ///
   _cons = " ") ///
mlabel format(%9.2f) mlabposition(12) mlabgap(*1.5) ///
xscale(alt) msymbol(x) ///
ysize(8) xsize(4) ///
title("Fixed part", span) ///
name(fixed_part, replace) nodraw

/* The default sizes of the available area are -ysize(4)- and -xsize(5.5)-,
   by the way. Letter size is -ysize(11)- and -xsize(8.5)-*/

// Plot random part   
  // Estimate residuals
capture drop u1* u0*
predict u1 u0, reffects 
predict u1se u0se, reses   

  // Plot intercept variation
preserve
egen pickone = tag(cntry)
keep if pickone
   
egen order_ = rank(-u0), unique
labmask order_, value(cntry)

gen high = u0 + (1.96 * u0se)
gen low  = u0 - (1.96 * u0se)

twoway (rcap u0 u0 order_, dsymbol(x)) ///
       (rspike high low order_) , ///
        yline(0) xlabel(1/22, val ang(v)) ///
        title("Intercept variance = `var_int'") ///
        xtitle("") ///
        ytitle("Random intercept residuals") ///
        legend(off) ///
        name(rand_int, replace) nodraw
restore 

  // Plot slope variation
preserve
egen pickone = tag(cntry)
keep if pickone

egen order_ = rank(-u1), unique
labmask order_, value(cntry)

gen high = u1 + (1.96 * u1se)
gen low  = u1 - (1.96 * u1se)

twoway (rcap u1 u1 order_, dsymbol(x)) ///
       (rspike high low order_) , ///
        yline(0) xlabel(1/22, val ang(v)) ///
        xtitle("") ///
        title("Slope variance = `var_age'") ///
        ytitle("Random slope residuals") ///
        ylabel(-1 (.5) 1) ///  // Make y-axis identical to other plot,
        legend(off) ///        // otherwise x-axis looks different, too
        name(rand_slope, replace) nodraw
restore 

  // Plot intercept-slope covariance
preserve

gen predRandomSlope= (_b[_cons] + u0) + ((_b[c_age] + u1) * c_age)
qui sum c_age
local hi1 =  1*r(sd)
local hi2 =  2*r(sd)
local hi3 =  3*r(sd)
local lo1 = -1*r(sd)
local lo2 = -2*r(sd)
      
sort cntry c_age
twoway (line predRandomSlope c_age, connect(ascending)), ///
        ytitle("Predicted Happiness") /// 
        xtitle("") ///
        title("Intercept{c 150}slope covariance = `covaria'") ///
        xlabel(`lo1' "-1 SD" 0 "Average age" `hi1' "+1 SD" `hi2' "+2 SD" `hi3' "+3 SD") ///
        ylabel(,format(%6.1f)) ///
        name(covariance, replace) nodraw
restore

// Combine Figures
graph combine rand_int rand_slope covariance, col(1) ysize(8) title("Random part") name(random_part, replace) nodraw
graph combine fixed_part random_part, col(2) ysize(8) xsize(8) altshrink title("Random coefficient model")

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 (34): Country scatterplots

use ART.dta, replace

// Fix missing values
recode avgemtrans (-1 -3 = .)
recode emhum      (-1 = .)

// Calculate r-squared
qui regress avgemtrans emhum
local r2 = round(e(r2), .01) // Round coeff.

// Create first panel
twoway (scatter avgemtrans emhum, mlab(countrycode) mlabpos(0) msymbol(i)) ///
       (lfit avgemtrans emhum) ///
    , xlabel(, format(%6.1f)) ylabel(, format(%6.1f)) ///
    ytitle("Avg. no. of fresh non-donor embryos" "per fresh embryo transfer cycle") ///
    xtitle("Avg. level disagreement that" "embryo is a human being") ///
       legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(8)) ///
    name(figure2a, replace)

// Fix missing values
recode afford clonehelpinf (-1 = .)

// Listwise deletion to control y-axis range
preserve
keep if !missing(afford, clonehelpinf)

// Calculate r-squared    
regress afford clonehelpinf
local r2 = round(e(r2), .01) // Round coeff.

// Create second panel
twoway (scatter afford clonehelpinf, mlab(countrycode) mlabpos(0) msymbol(i)) ///
       (lfit afford clonehelpinf) ///
    , xlabel(, format(%6.1f)) ylabel(0 (5) 25, format(%6.1f)) ///
    ytitle("Affordability (net cost of a fresh ART cycle" "as % of annual disposable income)") ///
    xtitle("Avg. level disagreement with" "cloning to help infertile couples") ///
       legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(8)) ///
    name(figure2b, replace)
restore

// Create final Figure
graph combine figure2a figure2b, row(1)