Re: [R] nested mixed-effect model: variance components

From: Christoph Buser <buser_at_stat.math.ethz.ch>
Date: Tue 20 Jun 2006 - 01:33:21 EST

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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Jun 20 01:38:29 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 20 Jun 2006 - 04:11:42 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.