[R-sig-ME] very basic HLM question

Sebastián Daza sebastian.daza at gmail.com
Mon Feb 7 17:51:51 CET 2011


Hi everyone,
I need to get a between-component variance (e.g. random effects Anova) 
in order to compute by hand the intraclass correlation. However, using 
lmer I don't get the same results (variance component) than using random 
effects Anova (with SPSS). I am using a database of students, clustered 
on schools (there is not the same number of students by school).

First, I get the intraclass using ICC1 in R:

summary(anova1 <- aov(math ~ as.factor(schoolid), data=nels88))
ICC1(anova1)

[1] 0.4414491

ICC1 comes from the multilevel package. Using lmer I get:

mod1 <- lmer(math ~ 1 + (1|schoolid), data=nels88)

Linear mixed model fit by REML
Formula: math ~ 1 + (1 | schoolid)
    Data: nels88
   AIC  BIC logLik deviance REMLdev
  1878 1888 -935.8     1875    1872
Random effects:
  Groups   Name        Variance Std.Dev.
  schoolid (Intercept) 34.011   5.8319
  Residual                72.256   8.5003
Number of obs: 260, groups: schoolid, 10

Fixed effects:
             Estimate Std. Error t value
(Intercept)   48.861      1.927   25.36

The random intercept effect is 34.011. If I compute the ICC manually 
(according to some articles I have read) I get:

 > 34.011/(34.011+72.256)
[1] 0.3200523

According to my Anova analysis, the between-component variance should be 
59.004. If I use 59.004 the formula works relatively well:

  59.004/(59.004+72.256)
[1] 0.449520037

Not equal, but at least similar in comparison with ICC from ICC1 command.
.
Does anyone know why I have got those differences? How can I get the 
59.004 figure using R? I should get practically an identical result than 
using ICC1 but I don't. I have used the database of that article and I 
have got an identical result between ICCs. But with the dataset I am 
using (nels88) doesn't work.

This is the example of the article mentioned above:

library(multilevel)
base(bh1996)

summary(lmer(WBEING ~ 1 + (1|GRP), data=bh1996))

0.035801/(0.035801+0.789497)
0.043379482

summary(test <- aov(WBEING~as.factor(GRP),data=bh1996))
ICC1(test)
[1] 0.04336905

Any suggestions?

-- 
Sebastián Daza
sebastian.daza at gmail.com




More information about the R-sig-mixed-models mailing list