# How to center in multilevel models
# R code from https://philippmasur.de/blog/2018-05-23-centering-in-multilevel-models/

library(tidyverse) # for data wrangling
library(multilevel) # for simulating multilevel data
library(lme4) # for estimating multilevel models
library(texreg) # for better outputs
select <- dplyr::select # to avoid clashes with the MASS package

set.seed(1000)
x <- rnorm(100, 2, 0.5)
y <- 3 + 3*x + rnorm(100, 2, .65)
d <- cbind(x,y) %>%
  as_tibble
summary(d)

# predicting y with the raw variable x
lm1 <- lm(y ~ x, d)

# Results
lm1$coef %>%
  round(2)

# centering by subtracting mean
d$x.c <- d$x - mean(d$x)
d[, c("x", "x.c")]

# plot distributions
d %>% 
  gather(Variable, val, c("x", "x.c")) %>%
  ggplot(aes(x = val, group = Variable, fill = Variable)) + 
  geom_density(alpha = .4) + 
  labs(x = "Variable X" , y = "Count") +
  geom_vline(xintercept = mean(d$x),
             linetype = "dashed") +
  geom_vline(xintercept = mean(d$x.c),
             linetype = "dashed") +
  xlim(-1.5, 4) +
  labs(x = "" , y = "Count") +
  theme_classic()

# ggsave("plot.png", height = 3.4, width = 7.5)

# predicting y with the mean-centered variable x.c
lm2 <- lm(y ~ x.c, d)

# Comparing result from lm1 and lm2
cbind(uncentered = lm1$coef, 
      centered = lm2$coef) %>%
  round(2)

# Simulation of multilevel data with two level 1 variables
set.seed(15324)
d <- sim.icc(gsize = 10,
             ngrp = 100,
             icc1 = .30,
             nitems = 2, 
             item.cor = .50)

# Renaming variables
d <- d %>%
  rename(id = GRP,
         dv = VAR1, 
         iv = VAR2) %>%
  as_tibble
d

# estimating the null model using lme4
m0 <- lmer(dv ~ 1 + (1|id), d)

# Function to compute the intraclass correlation coefficient (ICC)
compute_icc <- function(lmer_object){
  var_dat <- lmer_object %>% VarCorr %>% as.data.frame
  icc <- var_dat$vcov[1]/(var_dat$vcov[1]+var_dat$vcov[2])
  return(icc)
}
compute_icc(m0) %>%
  round(2)

d <- d %>% 
  # Grand mean centering (CMC)
  mutate(iv.gmc = iv-mean(iv)) %>%
  # Person mean centering (more generally, centering within cluster)
  group_by(id) %>% 
  mutate(iv.cm = mean(iv),
         iv.cwc = iv-iv.cm) %>%
  ungroup %>%
  # Grand mean centering of the aggregated variable
  mutate(iv.cmc = iv.cm-mean(iv.cm))

# Comparing the results of the centering approaches
d %>% 
  select(iv, iv.gmc, iv.cm, 
         iv.cwc, iv.cmc) %>%
  summary

m1.nc <- lmer(dv ~ 1 + iv + (1|id), d) # no centering
m1.gmc <- lmer(dv ~ 1 + iv.gmc + (1|id), d) # grand mean centering
m1.cwc <- lmer(dv ~ 1 + iv.cwc + (1|id), d) # centering within cluster
m1.cmc <- lmer(dv ~ 1 + iv.cmc + iv.cwc + (1|id), d) # including the person's mean (cluster's mean)

# Create a nice output (use screenreg() to get a readable output in the console)
screenreg(list(m1.nc, m1.gmc, m1.cwc, m1.cmc), 
          single.row = T, 
          stars = numeric(0),
          caption = "",
          custom.note = "Model 1 = Uncentered predictor, 
                         Model 2 = grand-mean centered predictor, 
                         Model 3 = person-mean centered predictor, 
                         Model 4 = person-mean centered predictor and centered person mean")

# Variance reduction approach according to Hox, 2010, pp. 69-78
var_reduction = function(m0, m1){
  library(tidyverse)
  VarCorr(m0) %>% 
    as.data.frame %>% 
    select(grp, var_m0 = vcov) %>% 
    left_join(VarCorr(m1) %>% 
                as.data.frame %>% 
                select(grp, var_m1 = vcov)) %>% 
    mutate(var_red = 1 - var_m1 / var_m0) 
}

# Create comparison table
cbind(M1 = var_reduction(m0, m1.nc)[,4],
      M2 = var_reduction(m0, m1.gmc)[,4],
      M3 = var_reduction(m0, m1.cwc)[,4],
      M4 = var_reduction(m0, m1.cmc)[,4]) %>%
  round(2)