# Calculating level 2 variable based on mean of level 1 variable with multilevel model in R

I have 299 surveys collected from 299 individuals working at 26 different locations. I want to understand how the location specific features relate to the individual survey scores. The only inference I have as to location features is gathered from the individual survey scores. Is it a valid strategy to calculate means for each location based on the individual scores, and include this as a level 2 variable? Further, does it also make sense to include the same variable but as the level 1 variable, with slope varying freely between locations, if I want to compare the relative usefulness of the mean (best estimate of 'reality') to a persons individual score? (their perception of reality and response biases).

I feel like I may have some circularity in the logic. My implementation in R for one of the variables of interest follows, any feedback is welcome!

lmer(X21~X25+meanX25+(X25|X1),data=datai)
Linear mixed model fit by REML
Formula: X21 ~ X25 + meanX25 + (X25 | X1)
Data: datai
AIC  BIC logLik deviance REMLdev
1079 1105 -532.7     1056    1065
Random effects:
Groups   Name        Variance Std.Dev. Corr
X1       (Intercept) 0.384983 0.62047
X25         0.012382 0.11127  -1.000
Residual             1.936068 1.39143
Number of obs: 299, groups: X1, 26
Fixed effects:
Estimate Std. Error t value
(Intercept)  1.13616    0.38013   2.989
X25          0.56683    0.05265  10.766
meanX25      0.33897    0.12213   2.775

Correlation of Fixed Effects:
(Intr) X25
X25     -0.119
meanX25 -0.838 -0.389

-
I don't know why my R output posts so badly, it look ok in the editor, I can't seem to fix it... – Charlie Nov 15 '12 at 8:37
You have to indent it with 4 spaces manually or by selecting and clicking {} button. – mbq Nov 15 '12 at 9:43

If the predictor variable that you are interested in ("X25" in this case) varies both between and within units (that is, it varies "at level 1" in multilevel model parlance), then yes, the procedure you described can make sense. This is in fact exactly the procedure that Bafumi and Gelman (LINK) describe when discussing the problem of fitting multilevel models where the group effects are correlated with one or more of the predictors. Another paper I recently found that I think is very helpful in explaining these issues can be found HERE.

Although both of these papers are primarily concerned with solving a particular statistical problem that you did not specifically bring up here, the point in bringing them to your attention is both to point out that fitting models like the one you described can indeed be a sensible thing to do, and also to direct you to some information on interpreting the resulting parameter estimates (see in particular the second paper on this point).

-
lmer(X21~X25+meanX25+(X25|X1),data=datai)


There are two implicit intercepts provided in lmer unless you override it, so this is equivalent to:

lmer (X21 ~ 1 + X25 + meanX25 + (1 + X25 | X1), data=datai)


Where the "fixed" X25 represents the overall mean across locations, and the "random" 1 represents the average difference (intercept) of each location from the overall mean. (I'm assuming X1 is the location.)

So I'm a bit confused as to why you'd also include a manually-calculated mean for each group (meanX25), instead of simply:

lmer (X21 ~ 1 + X25 + (1 + X25 | X1), data=datai)


where you can take the overall (fixed effect) intercept and add it to each location's difference (random effect) to come up with each location's average. It's already in there, by virtue of a mixed-effects approach -- as I understand it.

If it does make sense to include meanX25, you should explicitly remove the intercept from the random part:

lmer (X21 ~ X25 + meanX25 + (0 + X25 | X1), data=datai)

-