[R] Lmer BLUPS: was(lmer multilevel)

[Douglas Bates]

Sorry to be late coming to this discussion.   I was just talking with
Harold on the telephone and we decided why the inconsistency exists.

The lme function is designed to fit models with strictly nested random
effects.  (It can be used to fit models with crossed random effects
but the process is awkward and slow - it is best to think of it as
just fitting models with nested random effects.)  In that case it
makes sense to talk about coefficients at different levels.

The lmer function allows for models with crossed or partially crossed
random effects.  As soon as you have the slightest bit of crossing you
can no longer specify coefficients at various levels.  That is why you
must specify exactly which grouping factor you want the coefficient
for.
In Harold's example the call

coef(fm2)[[1]]

is asking for coefficients constructed from the fixed effects and the
random effects for the student only.  The call

coef(fm1, level = 2)

is asking for coefficients constructed from the fixed effects and the
random effects for the student and the random effects for the school
that the student attends because it is assumed that the student is
nested within schools.


On 3/29/06, Doran, Harold <HDoran at air.org> wrote:
>
>
> Paul:
>
> I may have found the issue (which is similar to your conclusion).  I checked
> using egsingle in the mlmRev package as these individuals are strictly
> nested in this case:
>
> library(mlmRev)
> library(nlme)
>
> fm1 <- lme(math ~ year, random=~1|schoolid/childid, egsingle)
> fm2 <- lmer(math ~ year +(1|schoolid:childid) + (1|schoolid), egsingle)
>
> Checking the summary of both models, the output is exactly the same w.r.t.
> the fixed effects, variance components, standard errors etc. The prior two
> lines of code fit the same models. However, the following does not generate
> similar output:
>
> head(coef(fm2)[[1]])
>                (Intercept)      year
> 2020:273026452 -0.30213394 0.7461233
> 2020:273030991  0.41469885 0.7461233
> 2020:273059461 -0.07443003 0.7461233
> 2020:278058841  0.61676579 0.7461233
> 2020:292017571  0.29529524 0.7461233
> 2020:292020281 -1.03811716 0.7461233
>
> > head(coef(fm1, level=2))
>                (Intercept)      year
> 2020/273026452   0.3248913 0.7461233
> 2020/273030991   1.0417241 0.7461233
> 2020/273059461   0.5525952 0.7461233
> 2020/278058841   1.2437910 0.7461233
> 2020/292017571   0.9223205 0.7461233
> 2020/292020281  -0.4110920 0.7461233
>
> Although they are similar for schools:
>
> > head(coef(fm1, level=1))
>      (Intercept)      year
> 2020  -0.1534565 0.7461233
> 2040  -0.6985680 0.7461233
> 2180  -1.0621073 0.7461233
> 2330  -0.6262567 0.7461233
> 2340  -1.0090541 0.7461233
> 2380  -0.6095037 0.7461233
>
> > head(coef(fm2)[[2]])
>      (Intercept)      year
> 2020  -0.1534566 0.7461233
> 2040  -0.6985681 0.7461233
> 2180  -1.0621072 0.7461233
> 2330  -0.6262567 0.7461233
> 2340  -1.0090541 0.7461233
> 2380  -0.6095037 0.7461233
>
> I checked a bit further to see if there was any pattern to the differences:
>
> tt <- cbind(coef(fm1, level=2)[1], coef(fm2)[[1]][1])
> tt$dif <- tt[,1] - tt[,2]
>
> There is clearly a systematic pattern to the difference between the two.
> Now, when you go and look at the random effects
>
>  head(ranef(fm2)[[2]])
>
> You can see that this is the value that is not being added in to the
> coefficients in lmer() and this accounts for the difference. So, these need
> to be added in for now.
>
> I hope this helps,
> Harold
>
>
>
> My question relates to problems that I'm having matching lme and lmer
> examples in P&B.
> using Matix 0.995
>
> In the Oxide example in p167-170 I can't get the level 2 coefficient
> estimates to match the fm1Oxide model in lme is
>
> data(Oxide,package="nlme")
> lme(Thickness~1,Oxide)
>
> which I translate in Lmer syntax to
>
> fm3Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)
> #or alternatively which gives the same result
> Oxide$LW<-with(Oxide,Lot:Wafer)[drop=TRUE]
>
> fm4Oxide<-lmer(Thickness~ (1|Lot)+(1|LW),data=Oxide)
>
> however if you look at say Lot 8, lme gives
>
> 8    1993.767
>
> 8/1    1993.677
> 8/2    1995.170
> 8/3    1990.693
>
>
> while  lmer gives
> coef(fm3Oxide)
> ....
> 8    1993.767
> 8:1 2000.062
>
> 8:2 2001.555
>
> 8:3 1997.078
>
> To me this looks like lmer in not including the lot random effect (8=
> -6.385129, Intercept 2000.153 ). Is this because I'm not specifying the
> model correctly?
>
> Thanks Paul