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 (144): Automatically labeling figures with letters

// Open ESS round 5
use agea stflife cntry if inlist(cntry, "DE", "GB", "NL", "SE") using "ESS5e03_4.dta", clear

// Country variable
kountry cntry, from(iso2c)
rename NAMES_STD country

// Restrict to complete cases
keep if !missing(agea, stflife, cntry)

// Center age
qui sum agea, detail
replace agea = `r(p99)' if agea > `r(p99)' & !missing(agea)
bys country: center agea, gen(age)
drop agea
label var age "Age (centered)"

// Letters
tokenize "`c(ALPHA)'"
local i 1

levelsof country, local(country)

foreach z of local country {
  qui sum age if country == `"`z'"', detail 
  foreach x of numlist 10 25 50 75 90 {
    local x`x': di %9.1f r(p`x')
    di `x`x''
  }
  qui regress stflife c.age##c.age if country == `"`z'"'
  
  qui margins, at(age = (`x10' `x25' `x50' `x75' `x90'))
  marginsplot, title(`"{bf:``i''} `z'"') ///
               ytitle("Life satisfaction" "(predicted)") ///
               xtitle("Age (percentiles)") ///
               xlabel(`x10' "10th" ///
                      `x25' "25th" ///
                      `x50' "50th" ///
                      `x75' "75th" ///
                      `x90' "90th") ///
               ylabel(, format(%6.1f)) ///
               recastci(rarea) ciopts(color(gs10)) ///
               name("``i''", replace) nodraw
  local ++i
}

graph combine A B C D, col(2) ycommon

Modeling variance functions

version 12.1
set more off

// Open National Longitudinal Survey (Young women 14-26 years of age in 1968)
webuse nlswork, clear
// Level 1: annual measurement
// Level 2: young women

// Null model
xtmixed ln_w || id: , var
xtmrho // VPC in random intercept case w/out covariates equals ICC

// Prepare level-1 covariate
center tenure, replace
label var tenure "Job tenure (years)"

// Add level-1 covariate
xtmixed ln_w tenure || id: , var
xtmrho // VPC after adjusting for a covariate

// 1) Constant variance

// Create variables: constant variance function at level 2
gen lev2var = e(var_u1)
label var lev2var "Variance at level 2"

// Create variables: constant variance function at level 1
gen lev1var = e(var_e)
label var lev1var "Variance at level 1"

// Create variable: predicted values
predict fit, xb

// Create variables: constant bounds for between-neighborhood variation
gen lev2cbhi = fit + 1.96*(lev2var^.5)
gen lev2cblo = fit - 1.96*(lev2var^.5)
label var fit      "Fixed relationship"
label var lev2cbhi "Upper bound between-person variation"
label var lev2cblo "Lower bound between-person variation"

// Graphs: Heterogeneity as a constant function
qui twoway (line lev2var tenure, sort) ///
           (line lev1var tenure, sort) ///
           , legend(ring(0) pos(5)) ///
             ylabel(.085 (.01) .12,format(%6.3f)) ///
             ytitle("Variance") ///
             name(variancefunction1, replace)
  
qui twoway (line fit lev2cbhi lev2cblo tenure, sort) ///
           , legend(ring(0)) ///
             ylabel(1 (.5) 3,format(%6.1f)) ///
             ytitle("Predicted log wage") ///
             name(betweenperson1, replace)
    
graph combine variancefunction1 between1

            , col(1) ysize(8) xcommon ///
              title("Random intercept model") ///
              name(constantvariance, replace)

// Clean up
drop lev2var lev1var fit lev2cbhi lev2cblo

// Add random slope to model
xtmixed ln_w tenure || id: tenure, var cov(uns)
xtmrho // VPC after adjusting for a covariate and random slope

// 2) Quadratic variance
   // Getting at the covariance parameter 
   // and slope variation as
   // -xtmrho- hasn't stored it in a better format
local cov      = tanh([atr1_1_1_2]_b[_cons])* exp([lns1_1_1]_b[_cons])* exp([lns1_1_2]_b[_cons])
di `cov' 
local var_tenure = exp([lns1_1_1]_b[_cons])^2
   // In order to get the other variance components
   // without -xtmrho-, one needs:
   *local var_cons = exp([lns1_1_2]_b[_cons])^2    // Intercept variance
   *local var_res  = exp([lnsig_e]_b[_cons])^2     // Residual variance

// Create variables: quadratic variance function at level 2
gen lev2var = e(var_u1) ///
            + 2 * `cov' * tenure ///
   + `var_tenure' * tenure * tenure
label var lev2var "Variance function at level 2"

// Create variable: predicted values
predict fit, xb

// Create variables: quadratic bound bounds for between-neighborhood variation
gen lev2cbhi = fit + 1.96*(lev2var^.5)
gen lev2cblo = fit - 1.96*(lev2var^.5)
label var fit      "Fixed relationship"
label var lev2cbhi "Upper bound between-person variation"
label var lev2cblo "Lower bound between-person variation"

// Graphs: Heterogeneity as a quadratic function
qui twoway (line lev2var tenure, sort) ///
           , legend(ring(0)) ///
             ytitle("Variance") ///
             ylabel(, format(%6.3f)) ///
             legend(ring(0) pos(11)) ///
             name(variancefunction2, replace)

qui twoway (line fit lev2cbhi lev2cblo tenure, sort) ///

           , legend(ring(0) pos(11)) ///

             name(betweenperson2, replace)

qui graph combine variancefunction2 betweenperson2, ///
                  , col(1) xcommon name(quadraticvariance, replace) ///
                    title("Random coefficient model")
   
graph combine constantvariance quadraticvariance, row(1) xsize(12) ysize(8)

// Clean up
drop lev2var fit lev2cbhi lev2cblo

Grand-mean centering variables on a higher level

It is good practice to center predictor variables in multilevel modeling. The -center- command is extremely useful in these situations. Especially grand-mean centering is an approach that is often necessary. However, when the variables to be centered are located on the group level and groups at this higher level are of unequal sizes, centering can become inaccurate, as bigger groups get a higher weight when calculating the grand mean. A way to deal with this problem is as follows:

egen pickone = tag(group) 
  // Create a tag variable "pickone" that has a value of 1 for a single, randomly
  // chosen case per group and 0 for all other cases

center grouplevelpredictor if pickone
  // This then centers the variable "grouplevelpredictor" around the grand mean 
  // at the group level due to the if-condition and creates a new variable 
  // "c_grouplevelpredictor." The problem is now only that this variable has
  // missing values for all cases where "pickone" does not equal 1

bysort group (c_grouplevelpredictor): ///
 replace c_grouplevelpredictor = grouplevelpredictor[1]
  // This command then replaces the missing values per country with the centered,
  // non-missing value. This works because the data are sorted by the variable
  // "c_grouplevelpredictor." The centered, non-missing value will be the first
  // one as it is – at least in Stata – always smaller than the missing values.

Random graphs (4): Plotting intercept and slope variation

use V3 V51 ISCO88 using "C:\Users\User\work\data sets\issp - work orientations\work orientations iii (2005)\ZA4350_F1.dta", clear

// Generate variables
gen jobsatisf = 7 - V51
iskoisei isei, isko(ISCO88)
drop V51 ISCO88
preserve




// Calculating intercept and slope variation using OLS
statsby inter = _b[_cons] ///
        slope = _b[isei] ///
       , by(V3) ///
         saving(ols, replace): ///
         regress jobsatisf isei
merge m:1 V3 using ols
drop _merge 

// Visualizing intercept and slope variation
gen yhat_ols = inter + slope*isei
separate jobsatisf, by(V3)
separate yhat_ols, by(V3)

twoway (line yhat_ols1-yhat_ols43 isei, sort(V3 isei)) ///
       (lfit jobsatisf isei, clwidth(vvthick) clcolor(black)) ///
    , legend(off) ytitle("Job satisfaction") xtitle("ISEI") ///
      xlabel(16 25 50 75 90) ylabel(,format(%6.1f)) ///
      caption("{it:Source:} ISSP 2005 (Work Orientations III), own calculations", span) ///
         name(one, replace)
restore



// Calculating intercept and slope variation using a random effects model
center isei, inplace
mixed jobsatisf isei || V3: isei, var
predict u1 u0, reffects

// Visualizing intercept and slope variation
gen predRandomSlope = (_b[_cons] + u0) + ((_b[isei] + u1) * isei)
twoway (line predRandomSlope isei, connect(ascending) sort(V3 isei)), ///
        ytitle("Job satisfaction") /// 
        xtitle("ISEI (centered)") ///
        xlabel(-25 0 25 50) ///
        ylabel(,format(%6.1f)) ///
        caption("{it:Source:} ISSP 2005 (Work Orientations III), own calculations", span) ///
        name(two, replace)