[R] Use predict() after lmer() (library lme4)

[Marc Girondot]

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"]
)