# Equivalence of Fixed Effects in Contextual Models with and without Random Slopes

When estimating "contextual models" (i.e., models that contain level-1 predictors as well as their cluster means on level-2), the estimation of the fixed effects should be unaffected by the centering scheme (grand mean centering, cgm, vs. group mean centering, cwc), according to the (psychological) literature (see, e.g., Kreft et al., 1995, or Enders & Tofighi, 2007).

However, if I try this in R with lme4, this is (numerically) true only if the models are specified without random slopes.

Moreover, the well-known relationship for the context effect (or its variants), $$\gamma_{10}^{cgm}+\gamma_{01}^{cgm}=\gamma_{01}^{cwc}$$, also only holds in this case (see, e.g., Enders & Tofighi, 2007).

Enders and Tofighi (2007) write: Because of a lack of precision in the estimation process, the algebraic identities outlined by Kreft et al. (1995) may not hold exactly when substituting the estimated parameter values into the equations outlined in this section.

My questions are:

1. Why is this so?
2. Why is the imprecision only there for models without random slopes?
3. Which of the two is (cgm or cwc) is the "better" estimate?

Also note that the problem remains when using ML and not REML for estimation.

Example code (Example comes from here)

library(dplyr)
library(lme4)

select(id, negaff, pss) %>%
mutate(stress = 4 - pss) %>%
mutate(stress_trait = mean(stress, na.rm = TRUE), .by = id) %>%
mutate(stress_trait_c = stress_trait - mean(stress_trait)) %>%
mutate(stress_cgm = stress - mean(stress, na.rm = TRUE),
stress_cwc = stress - stress_trait)

ri_cgm <- lmer(negaff ~ stress_cgm + stress_trait_c + (1 | id), data = AMIB)
ri_cwc <- lmer(negaff ~ stress_cwc + stress_trait_c + (1 | id), data = AMIB)

fixef(ri_cgm)
sum(fixef(ri_cgm)[-1])
fixef(ri_cwc)

rs_cgm <- lmer(negaff ~ stress_cgm + stress_trait_c + (stress_cgm | id), data = AMIB)
rs_cwc <- lmer(negaff ~ stress_cwc + stress_trait_c + (stress_cwc | id), data = AMIB)

fixef(rs_cgm)
sum(fixef(rs_cgm)[-1])
fixef(rs_cwc)


Output:

> fixef(ri_cgm)
(Intercept)     stress_cgm stress_trait_c
2.4567674      0.8429788      0.1996167
> sum(fixef(ri_cgm)[-1])
[1] 1.042596
> fixef(ri_cwc)
(Intercept)     stress_cwc stress_trait_c
2.4587394      0.8429788      1.0425955
>
> fixef(rs_cgm)
(Intercept)     stress_cgm stress_trait_c
2.4473444      0.7803662      0.2232565
> sum(fixef(rs_cgm)[-1])
[1] 1.003623
> fixef(rs_cwc)
(Intercept)     stress_cwc stress_trait_c
2.459474       0.775756       1.019996


Literature

Kreft, I. G., De Leeuw, J., & Aiken, L. S. (1995). The effect of different forms of centering in hierarchical linear models. Multivariate behavioral research, 30(1), 1-21.

Enders, C. K., & Tofighi, D. (2007). Centering predictor variables in cross-sectional multilevel models: a new look at an old issue. Psychological methods, 12(2), 121.

• One immediate problem is that your stress_trait calculation is producing a single value applied to all ids. Instead, you need to group_by id and calculate a person mean: AMIB <- AMIB %>% group_by(id) %>% mutate(stress_pmn = mean(stress, na.rm = TRUE)) %>% ungroup() Commented Jun 3 at 19:31
• Thank you for your comment and answer. In fact, in the line mutate(stress_trait = mean(stress, na.rm = TRUE), .by = id) the argument .by = id should take care of that, it's a shortcut such that we don't have to group/ungroup. Commented Jun 4 at 13:57
• Ah, that must be a feature in a newer version of dplyr that I didn't have on my old laptop. Sorry about that! Commented Jun 4 at 14:15

• Thank you for your answer, although I don't think it addresses the question, which refers to why the two estimates ($\gamma_{10}^{cwc}$ and $\gamma_{10}^{cgm}$) which should be the same are numerically different when there are random slopes. Commented Jun 5 at 7:40