[R-sig-ME] LMER vs MLwiN

[W Robert Long]

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 
? 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