I'm fitting some hierarchical models in R using lmer, and am trying to understand why the results change as they do when I either include or exclude a fixed-effects intercept term. I have already looked elsewhere on CrossValidated and the broader web, to no avail.
Some quick background to help interpretation of the below snippets: N=149 participants each viewed the same set of M=40 videos, and responded to each video in various ways; in the models below, the dependent measure is the Arousal subscale of the Self Assessment Manikin. The reason for using a random intercept for participants is, I think, self-explanatory; the reason for using a random intercept for videos as well is that I want to ensure the results can be interpreted as generalizing to videos outside of these specific 40.
When I fit the model with a fixed-effects intercept, these are the results; in both of the following snippets, the key output is the reported random effects variance for the Video intercept, along with the calculated values reported at the bottom:
> summary(lmer(SAM_Aro ~ (1 | Participant) + (1 | Video), data=dat))
Linear mixed model fit by REML ['lmerMod']
Formula: SAM_Aro ~ (1 | Participant) + (1 | Video)
Data: dat
REML criterion at convergence: 22263
Scaled residuals:
Min 1Q Median 3Q Max
-5.2073 -0.5432 0.0458 0.5935 3.5653
Random effects:
Groups Name Variance Std.Dev.
Participant (Intercept) 2.4008 1.5494
Video (Intercept) 0.2045 0.4522
Residual 2.1914 1.4804
Number of obs: 5960, groups: Participant, 149; Video, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.9903 0.1469 40.77
> var(ranef(lmer(SAM_Aro ~ (1 | Participant) + (1 | Video), data=dat))$Video[["(Intercept)"]])
[1] 0.1907348
Of course, the variance of the random effects terms won't exactly match the modeled variance of the parametric distribution from which those terms are assumed to arise by the model, but intuitively, they should be close, as they are here.
Excluding the fixed-effects intercept, the picture changes drastically:
> summary(lmer(SAM_Aro ~ (1 | Participant) + (1 | Video) + 0, data=dat))
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: SAM_Aro ~ (1 | Participant) + (1 | Video) + 0
Data: dat
AIC BIC logLik deviance df.resid
22470.1 22490.1 -11232.0 22464.1 5957
Scaled residuals:
Min 1Q Median 3Q Max
-5.2162 -0.5444 0.0471 0.5939 3.6062
Random effects:
Groups Name Variance Std.Dev.
Participant (Intercept) 2.413 1.553
Video (Intercept) 34.799 5.899
Residual 2.191 1.480
Number of obs: 5960, groups: Participant, 149; Video, 40
> var(ranef(lmer(SAM_Aro ~ (1 | Participant) + (1 | Video) + 0, data=dat))$Video[["(Intercept)"]])
[1] 0.2189779
> 1/40*sum(ranef(lmer(SAM_Aro ~ (1 | Participant) + (1 | Video) + 0, data=dat))$Video[["(Intercept)"]]^2)
[1] 34.76858
In this case, the reported variance for the Video random intercept term increases by a factor of 170 compared with the first model. The actual variance of the terms themselves remains essentially unchanged, as expected. On a hunch, I computed the second raw moment of the terms, and it is clear that this value corresponds very closely with what was reported in the lmer summary output.
My questions are: 1) is this intended behavior on the part of lmer, or some sort of bug? And 2) if I want to report the (modeled) variances associated with each grouping variable, should I just resign myself to using the estimates from the model with a fixed-effects intercept? There's no reason in principle to disprefer it, except that the model without seems somehow cleaner in terms of interpreting the parcelling of variance components to each of the two groups plus residual.
Any clarifications on this specific issue, as well as more philosophical thoughts as to the appropriateness of using either model, are greatly appreciated.
