Bonjour, (don't worry, after I will write in English [at least I will
try ;) ])
I try to understand better mixed models and then I have generated data
and I try to understand how the fixed and the random effects are used in
predict(). I understand when the random effect is of the form (1 | rf]
but I don't understand for the form (rf1 | rf2]. Let do an example. The
last formula does not give the same than predict().
Thanks if someone could explain what formula use to reproduce the
predict() results.
Marc
# 1/ Generate data in a data.frame, 1 response (number), one effect
(effect) and two factors (sector and beach) that I want use as random
effect. These two factors are hierarchical beach within sector
Sector <- c(rep("I", 100), rep("II", 100))
Beach <- c(
sample(c("A", "B", "C", "D", "E"), 100, replace=TRUE),
sample(c("F", "G"), 100, replace=TRUE)
)
number <- rnorm(200, 10, 1)
# Sector effect
number[1:100] <- number[1:100] +0.1
number[101:200] <- number[101:200] +0.5
# beach effect
beach.value <- 1:7
names(beach.value) <- LETTERS[1:7]
number <- number + unname(beach.value[Beach])
dataF <- data.frame(number=number, effect= number/10+runif(200, 0, 2),
Sector=Sector, Beach=Beach)
plot(dataF$number, dataF$effect)
##############
library("lme4")
##############
##############
# 2/ Random effect is (1 | Beach). I can reproduce the predict()
##############
out1 <- lmer(formula = number ~ effect + (1 | Sector) , data=dataF)
head(predict(out1))
ef <- fixef(out1)
er <- ranef(out1)
head(ef["(Intercept)"]+dataF$effect*ef["effect"]+
er$Sector[dataF$Sector, "(Intercept)"]
)
##############
# 3/ Random effect is (1 | Beach) + (1 | Sector). I can reproduce the
predict()
##############
out2 <- lmer(formula = number ~ effect + (1 | Sector) + (1 | Beach),
data=dataF)
head(predict(out2))
ef <- fixef(out2)
er <- ranef(out2)
head(ef["(Intercept)"]+dataF$effect*ef["effect"]+
er$Sector[dataF$Sector, "(Intercept)"] +
er$Beach[dataF$Beach, "(Intercept)"]
)
##############
# 4/ Random effect if (Sector | Beach). I don't understand how to
reproduce the predict()
##############
out3 <- lmer(formula = number ~ effect + (Sector | Beach), data=dataF)
head(predict(out3))
ef <- fixef(out3)
er <- ranef(out3)
head(ef["(Intercept)"]+dataF$effect*ef["effect"]+
er$Beach[dataF$Beach, "(Intercept)"]+
er$Beach[dataF$Beach, "SectorII"]
)
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
