[R] nested mixed-effect model: variance components

[Christoph Buser]

Dear John

I've put your mail back to the R list since I have no
explanation for the "lmer" result. Maybe someone else has an
idea. I adapted it to show some else that I do not understand. 

## Creation of the data (habitat is nested in lagoon):
set.seed(1)
dat <- data.frame(y = rnorm(100), lagoon = factor(rep(1:4,each = 25)),
                  habitat = factor(rep(1:20, each = 5)))

## I do not understand how the random intercepts for lagoon and
## lagoon:habitat can both be estimated. It seems a little bit
## strange to me that they are identical (0.46565).
library(lme4)
summary(reg3 <- lmer(y~habitat+(1|lagoon)+(1|lagoon:habitat), data = dat))

## Furthermore I do not understand why the standard errors for
## the fixed effect of habitat are 1.131 for some habitats
## and 1.487 for the others???


## If one removes (1|lagoon), the variance component
## (1|lagoon:habitat) does not change its value (still 0.46565)???
summary(reg3a <- lmer(y~habitat+(1|lagoon:habitat), data = dat))

## Now all standard errors for the fixed factor habitat are 1.131.

## Altogether it seems a little bit strange to me and with the
## warnings and errors of the lme and aov call, I'd be carefull
## by using the output of lmer in that case. In addition I do
## not understand the interpretation of the random effect lagoon
## in top of the nested FIXED factor habitat.

summary(aov(y~habitat + Error(lagoon/habitat), data = dat))

detach(package:Matrix)
detach(package:lme4)
library(nlme)
summary(reg2 <- lme(y~habitat, random = ~1|lagoon/habitat, data = dat))
anova(reg2)


Best regards,

Christoph Buser

--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C13
ETH Zurich	8092 Zurich	 SWITZERLAND
phone: x-41-44-632-4673		fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------



>> Christoph,
>> 
>> I am sending this off list bacause I tried 'lmer'and
>> it seems to work with your contrived data,but I don't know why.
>> Can you explain it ?
>> 
>> 
>> John
>> 
>> 
>> 
>> 
>> > detach(package:nlme)
>> > library(Matrix)
>> 
>> summary(lmer(y~habitat+(1|lagoon)+(1|lagoon:habitat), data = dat))
>> 
>> Linear mixed-effects model fit by REML
>> Formula: y ~ habitat + (1 | lagoon) + (1 | lagoon:habitat)
>>           Data: dat
>>       AIC      BIC    logLik MLdeviance REMLdeviance
>>  292.3627 349.6764 -124.1814   280.0245     248.3627
>> Random effects:
>>  Groups         Name        Variance Std.Dev.
>>  lagoon:habitat (Intercept) 0.46565  0.68239
>>  lagoon         (Intercept) 0.46565  0.68239
>>  Residual                   0.87310  0.93440
>> number of obs: 100, groups: lagoon:habitat, 20; lagoon, 4
>> 
>> Fixed effects:
>>               Estimate Std. Error  t value
>> (Intercept)  0.1292699  1.0516317  0.12292
>> habitat2     0.0058658  1.1316138  0.00518
>> habitat3    -0.0911469  1.1316138 -0.08055
>> habitat4     0.3302971  1.1316138  0.29188
>> habitat5    -0.0480394  1.1316138 -0.04245
>> habitat6    -0.4778469  1.4872319 -0.32130
>> habitat7    -0.0867301  1.4872319 -0.05832
>> habitat8     0.0696507  1.4872319  0.04683
>> habitat9    -0.0998728  1.4872319 -0.06715
>> habitat10    0.1096064  1.4872319  0.07370
>> habitat11   -0.0430979  1.4872319 -0.02898
>> habitat12    0.0714719  1.4872319  0.04806
>> habitat13    0.3380993  1.4872319  0.22733
>> habitat14    0.3057808  1.4872319  0.20560
>> habitat15   -0.4915582  1.4872319 -0.33052
>> habitat16   -0.2624539  1.4872319 -0.17647
>> habitat17   -0.2203461  1.4872319 -0.14816
>> habitat18    0.2165269  1.4872319  0.14559
>> habitat19    0.6932896  1.4872319  0.46616
>> habitat20   -0.7271468  1.4872319 -0.48893
>>  anova(summary(lmer(y~habitat+(1|lagoon)+(1|lagoon:habitat), data = dat)))
>> 
>> Analysis of Variance Table
>>         Df  Sum Sq Mean Sq
>> habitat 19 2.65706 0.13985
>> 
>> > VarCorr(summary(summary(lmer(y~habitat+(1|lagoon)+(1|lagoon:habitat), data
>> = dat)))
>> + )
>> $`lagoon:habitat`
>> 1 x 1 Matrix of class "dpoMatrix"
>>             (Intercept)
>> (Intercept)   0.4656545
>> 
>> $lagoon
>> 1 x 1 Matrix of class "dpoMatrix"
>>             (Intercept)
>> (Intercept)   0.4656545
>> 
>> attr(,"sc")
>> [1] 0.9343993
>> 
>> Christoph Buser wrote ---
>> 
>> Dear Eric
>> 
>> Do you really have habitats nested within lagoons or are they
>> partially crossed (meaning that you have the same habitats in
>> different lagoons)?
>> 
>> If you have them perfectly nested, I think that you cannot
>> calculate both a fixed effect for habitats and a random effect
>> for lagoon (see the example below, lme and aov).
>> 
>> You can compare e.g. two lagoons by defining a contrast of the
>> habitats of one lagoon against the habitats of the other (if you
>> think that this is a meaningful test to interpret), but you
>> cannot estimate a random effect lagoon in presence of a nested
>> FIXED effect habitat.
>> 
>> aov() will not return you the test and warn you about the
>> singular model.
>> 
>> lme() will estimate a variance component for lagoon, but does
>> not provide you a test for the fixed factor.
>> 
>> Regards,
>> 
>> Christoph Buser
>> 
>> 
>> 
>> set.seed(1)
>> dat <- data.frame(y = rnorm(100), lagoon = factor(rep(1:4,each = 25)),
>>                   habitat = factor(rep(1:20, each = 5)))
>> 
>> summary(aov(y~habitat + Error(lagoon/habitat), data = dat))
>> 
>> library(nlme)
>> summary(lme(y~habitat, random = ~1|lagoon/habitat, data = dat))
>> 
>> 
>>