Let's say you have a response variable and an independent variable. Your data is measured across several levels of a categorical independent variable. One approach in analysing these data would be to use linear regression to estimate a slope at each level of the categorical independent variable. This is the approach I've used here, using
sleepstudy dataset from the
lme4 package (I've stored the betas from each model in
library(lme4); library(plyr); library(ggplot2) lmBetas <- daply(sleepstudy, .(Subject), function(x) coef(lm(Reaction ~ Days, data=x))["Days"])
Another approach in analysing these data would be to use a mixed effects model to estimate slopes for each level of the categorical independent variable, which in this case is
Subject. This is the approach I've taken here (I've stored the betas from the model in
lmerBetas <- coef(lmer(Reaction ~ Days + (Days | Subject), data=sleepstudy))$Subject[,"Days"]
I have learned that a single mixed effects model, as implemented through the
lmer function in R, is more accurate at estimating slopes than a multiple linear regression model applied to multilevel data. This can be demonstrated with this plot of betas from the above models.
betas <- data.frame(method.betas = c(lmerBetas, lmBetas)) betas$method <- c(rep("lmer", 18), rep("lm", 18)) ggplot(betas, aes(method.betas)) + geom_histogram() + facet_grid(method ~ .)
The top histogram shows betas estimated using linear regression, and the bottom histogram shows betas estimated using mixed effects. You can see betas estimated using linear regression are more widely spread than those estimated through the mixed effects model.
So finally, my questions:
Is a mixed effects model's higher accuracy in betas estimation connected with the fact that it models intercepts and slopes for each level of the categorical independent variable under a joint probability distribution?
Generally speaking, why is a mixed effects model more accurate in its betas estimation?