I realise this question reflects my ignorance more than anyone else, but if anyone can give me an answer I'd really appreciate it. I'm trying to understand mixed effects models. Let's say you've got a model with a random effect(group) with 50 levels and a fixed effect that you're fitting as a covariate. If you fitted a fixed-effects model with a group*covariate interaction you'd get a separate estimate for the slope and intercept for each group, but you'd use up a shedload of df. If you fit a random intercept and slopes model I get the bit about estimating the components of the variance arising from the groups represented by a random effect, but I'm finding it hard to work out how the model then produces group-level estimates for the intercept and slope without using up the same number of degrees of freedom that you'd use for the straight fixed-effects equivalent. Here's a toy example:
library(lme4)
library(dplyr)
# Fake explanatory variables
group <- rep(1:50, each = 5)
covar <- runif(250)
# Generate response variable
response <- 0.01*group + 2 * covar + rnorm(250)
# Fit random intercepts and slopes model
mod1 <- lmer(response ~ covar + (1+ covar|group))
# Extract coefficients
coef(mod1)$group %>% head(5)
Which gives you output something like this:
(Intercept) covar
1 -0.34799353 2.519426
2 0.24860259 2.161321
3 -0.01285849 2.301451
4 -0.30915518 2.424565
5 -0.30823363 2.635467
So the model is producing an estimate for the intercept and slope separately for each group: how does this work and why is it different from a fixed effects model? I guess the answer is somewhere in the black magic of how the random effect is fitted but I can't find a clear explanation of why this is different from doing it with a fixed effects model anywhere.
Any help much appreciated.