Let's say I have measures of a response and covariate from each of several populations. Each population has its own mean for $y_i$ and slope for the $y_i \sim x_i $ relationship. The model looks like this:
$$y_{i,j} = \alpha_j + \beta_j \times x_{i,j} + \epsilon_{i,j}$$ $$ \alpha_j\sim N(\mu_\alpha,\sigma^2_\alpha) $$ $$ \beta_j\sim N(\mu_\beta,\sigma^2_\beta) $$ $$ \epsilon_i \sim N(0,\sigma^2) $$
Now imagine I am simply interested in estimating the average effect of the covariate, $\mu_\beta$, its SE, and its among-population variance, $\sigma^2_\beta$. I am uninterested in the differing means of $ y $ among populations (imagine, for example, they are artifacts of study design). I could center the response by population such that within-population mean=0, and thereby eliminate the need to model this as a random intercept (?). If I'm only interested in the effect of the covariate, is there negative consequences to standardizing the response by population? A simple simulation in R seems to indicate No (see below). Are there other reasons I've not thought of why group-centering could be problematic in more complicated situations?
Would this also extend to differing variances by group? i.e., I could standardize the response by group to mean=0, SD=1, and thereby not need to model the differing variances (this could be helpful for users of lme4 package which doesn't currently allow variance-by-group.)
Here is R code to simulate the centering problem:
require(ggplot2)
require(lme4)
require(dplyr)
require(broom)
#set up experimental design
n.groups <- 15
n.sample <- 30
n <- n.groups * n.sample
pop <- gl(n = n.groups, k = n.sample) # Indicator for population
covariate <- rnorm(n)
#model matrix
Xmat <- model.matrix(~pop*covariate-1-covariate)
intercept.mean <- 100 # mu_alpha
intercept.sd <- 25 # sigma_alpha
slope.mean <- 20 # mu_beta
slope.sd <- 10 # sigma_beta
intercept.effects<-rnorm(n = n.groups, mean = intercept.mean, sd = intercept.sd)
slope.effects <- rnorm(n = n.groups, mean = slope.mean, sd = slope.sd)
all.effects <- c(intercept.effects, slope.effects) # Put them all together
lin.pred <- Xmat[,] %*% all.effects # Value of linear predictor
eps <- rnorm(n = n, mean = 0, sd = 15) # residuals
y <- lin.pred + eps # response = linear predictor + residual
df<- data.frame(y,covariate,pop)
#visualize data
ggplot(df, aes(x=covariate,y=y)) +
geom_point() +
facet_wrap(~pop)
#fit data using lmer(), random intercept and slope (no correlation)
fit1<- lmer(y~covariate + (covariate||pop), data=df)
#center the response by group, refit without random intercept
df2<- df %>% group_by(pop) %>% mutate(centered_y= scale(y, scale=FALSE)[,1]) %>% ungroup()
fit2<- lmer(centered_y~covariate + (0+covariate||pop), data=df2)
#compare coefficients
tidy(fit1)
tidy(fit2)
