very basic HLM question

2011/2/7 Sebasti?n Daza<sebastian.daza at gmail.com>:

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).

I think the issue is that the estimates of the variance components are
different, which would then lead to different estimates of the ICC.

You mentioned in your original posting that the data are unbalanced
(i.e. there are different numbers of students in different schools).
The estimates returned from lmer are the REML estimates in this case.
Do you know how the other estimates are being calculated?  I'm not
sure that the naive calculation of within- and between- sums of
squares and equating expected mean squares with observed mean squares
works with unbalanced data.

First, I get the intraclass using ICC1 in R:

As Paul indicated, it is more informative to say, "using ICC1 from the
multilevel package for R".  ICC1 is not part of the base R nor the
required packages.

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

[1] 0.4414491

That is using the raw mean squares, not the estimates of the variance
components.  I have, thankfully, forgotten most of what I know about
expected and observed mean squares but I am pretty sure that those
don't correspond to the REML estimates, even in the balanced case.

Consider the results from the Dyestuff data, which is part of the lme4
package and is balanced

print(fm1<- lmer(Yield ~ 1|Batch, Dyestuff))

Linear mixed model fit by REML ['merMod']
Formula: Yield ~ 1 | Batch
    Data: Dyestuff
REML criterion at convergence: 319.6543

Random effects:
  Groups   Name        Variance Std.Dev.
  Batch    (Intercept) 1764     42.00
  Residual             2451     49.51
Number of obs: 30, groups: Batch, 6

Fixed effects:
             Estimate Std. Error t value
(Intercept)  1527.50      19.38    78.8

VarCorr(fm1)

$Batch
             (Intercept)
(Intercept)     1764.05
attr(,"stddev")
(Intercept)
     42.0006
attr(,"correlation")
             (Intercept)
(Intercept)           1

attr(,"sc")
[1] 49.5101

vc<- VarCorr(fm1)
var1hat<- as.vector(vc$Batch)
var1hat

[1] 1764.05

varhat<- attr(vc, "sc")^2
varhat

[1] 2451.25

var1hat/(var1hat + varhat)

[1] 0.4184874

anova1<- aov(Yield ~ 1|Batch, Dyestuff)

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
   0 (non-NA) cases
In addition: Warning message:
In Ops.factor(1, Batch) : | not meaningful for factors

anova1<- aov(Yield ~ Batch, Dyestuff)
summary(anova1)

             Df Sum Sq Mean Sq F value   Pr(>F)
Batch        5  56358 11271.5  4.5983 0.004398
Residuals   24  58830  2451.2

You can see that the estimates of the residual variance are the same
for the two model fits but the mean square for batch is not the
estimate of the between-batch variance.

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