[R-sig-ME] LMER vs MLwiN

[Douglas Bates]

On Fri, Feb 17, 2012 at 10:06 AM, W Robert Long <longrob604 at gmail.com> wrote:
> Hi Ben
>
> Thanks very much.
>
> I have removed PUPIL as random effect (makes 100% sense to me - it's
> variance is fixed - I should have realised that myself earlier !).
>
> glmer(LOSS~(1|COMMUNE/SCHOOL/CLASS),data=dt,family=binomial(link = "logit"))
>
> The output is:
>
> Generalized linear mixed model fit by the Laplace approximation
> Formula: LOSS ~ (1 | COMMUNE/SCHOOL/CLASS)
>   Data: dt
>   AIC   BIC logLik deviance
>  10380 10408  -5186    10372
> Random effects:
>  Groups                 Name        Variance Std.Dev.
>  CLASS:(SCHOOL:COMMUNE) (Intercept) 0.173259 0.41624
>  SCHOOL:COMMUNE         (Intercept) 0.736581 0.85824
>  COMMUNE                (Intercept) 0.024233 0.15567
> Number of obs: 9162, groups: CLASS:(SCHOOL:COMMUNE), 310; SCHOOL:COMMUNE,
> 98; COMMUNE, 20
>
> Fixed effects:
>            Estimate Std. Error z value Pr(>|z|)
> (Intercept)   0.9874     0.1030    9.59   <2e-16 ***
>
> The estimates look pretty close, but the standard errors for the REs are
> quite different - I seem to remember the sampling variance of REs has a
> skewed distribution, but I don't know if this has anything to do with it ?

Those are not standard errors in the glmer output.  They are simply
the variance estimates on the standard deviation scale (i.e. 0.15567 =
sqrt(0.024233)).  The reason that glmer does not provide a standard
error for an estimate of a variance component is because they don't
make sense in most cases.  The distribution of the estimator is highly
skewed.

> Was this what you were getting at when you said "Note that it is a bit
> harder to get uncertainty estimates on the variance parameters in lme4" ?
>
> BTW, MLwiN is using penalised quasi-likelihood (order=2).
>
> I was very interested to back-test using your simulated R code - would I
> need to change it, bearing in mind that pupil isn't a random effect ? I did
> export the data into MLwiN but it gave zero estimates for all the random
> effects (and SEs) and a negative estimate for the fixed effect :(
>
> Thanks again
> RL
>
>
>
> On 17/02/2012 3:15 PM, Ben Bolker wrote:
>>
>> W Robert Long<longrob604 at ...>  writes:
>>
>>>
>>> Hi Federico
>>>
>>> Thank you very much.
>>>
>>> However, the model I'm trying to fit has no fixed effects. The MLwiN
>>> output gives me:
>>> beta_0_hat: 1.012(0.107)  conditional log odds
>>> var_f0_hat: 0.031(0.066)  variance between communues
>>> var_v0_hat: 0.760(0.141)  variance between schools within commune
>>> var_u0_hat: 0.186(0.38)   variance between classes within schools
>>>
>>> This is what I am trying to duplicate in R with the same data.
>>>
>>> BTW I am using 64 bit R
>>
>>
>>   OK, based on this output you shouldn't include 'pupil' in your
>> random effects specification (now that I think of it, you probably
>> shouldn't anyway, because it's unidentifiable for a Bernoulli outcome).
>>
>>   If you wanted you could redo my example with your observed effects
>> (e.g. u.commune = rnorm(n.comm,sd=sqrt(0.031)) ...)
>>
>>   Note that it is a bit harder to get uncertainty estimates on the
>> variance parameters in lme4.
>>
>>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models