2
$\begingroup$

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)
Improve this question
$\endgroup$
4
  • $\begingroup$ Should the predictors in the first equation be $x_i$ instead of $x_i,_j$ ? $\endgroup$
    – Moss Murderer
    Oct 25, 2017 at 20:03
  • $\begingroup$ @MossMurderer I think $x_{i,j}$ is right, but $\epsilon_i$ should be $\epsilon_{i,j} $. Based on this, or are you getting at something else? $\endgroup$
    – Dave M
    Oct 26, 2017 at 17:19
  • $\begingroup$ No. I just wasn't sure I understood the question. Now I do. $\endgroup$
    – Moss Murderer
    Oct 27, 2017 at 4:44
  • 1
    $\begingroup$ For very simple random effects structures -- e.g. one grouping variable, zero correlations -- this seems like it might work, but it won't generalize well to more complicated structures. It also won't be trivial to use this to simplify complicated structures because it amounts to marginalizing out a particular grouping term with respect to the others (i.e. ignores the interrelationship between random effect terms) -- not a problem for certain types of nested random effects and classical experimental designs, but losing many of the advantages of more recent mixed model approaches. $\endgroup$
    – Livius
    Oct 28, 2017 at 11:53

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.