Dear Marc,
This question is more suited for R-Sig-mixed models to which I'm forwarding
it.
Your manual predictions for out3 are wrong. Here are the correct manual
predictions.
fixed <- model.matrix(~effect, data = dataF) %*% fixef(out3)
random <- rowSums(model.matrix(~Sector, data = dataF) *
ranef(out3)$Beach[dataF$Beach, ])
range(fixed + random - predict(out3))
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
2016-02-03 8:08 GMT+01:00 Marc Girondot <marc_...@yahoo.fr>:
> 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.
>
[[alternative HTML version deleted]]
______________________________________________
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.
