Skip to content

very basic HLM question

3 messages · Paul Johnson, Sebastián Daza

Original
Sebastián Daza
# 
Hi everyone,

I need to get a between-component variance (e.g. random effects Anova), 
but using lmer I don't get the same results (variance component) than 
using random effects Anova. I am using a database of students, clustered 
on schools (there is not the same number of students by school).

According to the ICC1 command, the interclass correlation is .44

 > ICC1(anova1)
[1] 0.4414491

However, I cannot get the same ICC from the lmer output:

 > anova2 <- lmer(math ~ 1 + (1|schoolid), data=nels88)
 > summary(anova2 <- 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 intercept random effect is 34.011. If I compute the ICC manually I get:

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

According to my Anova analysis, the between-component variance is 59.004.
Does anyone know what the problem is? How can I get the 59.004 figure 
using R?
Paul Johnson
# 
2011/2/5 Sebasti?n Daza <sebastian.daza at gmail.com>:
If you don't tell us exactly what model you are calculating in
"anova1", how would we guess if there is something wrong?

Similarly, I get this
Error: object 'ICC1' not found

so it must mean you've loaded a package or written a function, which
you've not shown us.

I googled my way to a package called "multilevel" that has ICC1, and
its code for ICC1 shows a formula that does not match the one you used
to calculate ICC from lmer.

function (object)
{
    MOD <- summary(object)
    MSB <- MOD[[1]][1, 3]
    MSW <- MOD[[1]][2, 3]
    GSIZE <- (MOD[[1]][2, 1] + (MOD[[1]][1, 1] + 1))/(MOD[[1]][1,
        1] + 1)
    OUT <- (MSB - MSW)/(MSB + ((GSIZE - 1) * MSW))
    return(OUT)
}

I'm not saying that's right or wrong, just not obviously identical to
the formula you proposed.
Instead, do this (same thing, fits model only once):
Note that lmer is going to estimate a normally distributed random
effect for each school, as well as an individual observation random
effect (usual error term) that is assumed independent of the
school-level effect.  What is "anova1" estimating?
Sebastián Daza
# 
Thank you for your reply and sorry for my ambiguity.

I computed:

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

ICC1 comes from the multilevel package. I have found an article where it 
is pointed out that with this formula:

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

I should get practically an identical result than using ICC1. I have 
used the database of that article and I have got an identical result. 
But with the dataset I am using (schools) it 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)

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


I have attached the small database where this "equality" doesn't exist.
Thank you in advance.
On 2/5/2011 11:07 PM, Paul Johnson wrote: