Hi,
I try to understand exactly how mixed model is working. Let do this
simple model with fake data:
library(lme4)
invlogit <- function(n) {1/(1 + exp(-n))}
d <- data.frame(A=c(10, 20, 10, 20, 30, 15, 17, 19, 20, 30, 20, 30, 21,
23),
??????????????? B=c(2, 3, 7, 9, 10, 8, 5, 7, 9, 10, 2, 3, 7, 8),
??????????????? D=c(10, 11, 12, 13, 14, 10, 11, 12, 13, 14, 13, 12, 12,
13),
??????????????? R1=c("A", "A", "A", "A", "A", "A", "A", "B", "B", "B",
"B", "B", "B", "B"),
??????????????? R2=c("C", "C", "C", "D", "D", "D", "D", "E", "E", "E",
"F", "F", "F", "F"))
m <- glmer(formula = cbind(A, B) ~ D + (1 | R1 / R2),
????? data=d,
????? family = binomial(link = "logit"))
It seems to work well (except the warning for boundary because probably
I don't introduce enough data).
Let do a predict:
predict(m, type="response")
??????? 1???????? 2???????? 3???????? 4???????? 5 6???????? 7????????
8???????? 9??????? 10??????? 11
0.7706127 0.7665437 0.7624248 0.7502806 0.7459699 0.7629174 0.7587547
0.7572462 0.7530161 0.7487368 0.7696229
?????? 12??????? 13??????? 14
0.7736541 0.7736541 0.7696229
It is ok. Now I try to do the predict by hand to be sure that I
understand how it works:
ef <- fixef(m)
er <- ranef(m)
fixed <- model.matrix(~ D, data = d) %*% ef
random <- rowSums(model.matrix(~ R1, data = d) * er$R2[d$R2, ])
invlogit(fixed[, 1] + random)
??????? 1???????? 2???????? 3???????? 4???????? 5???????? 6 7????????
8???????? 9??????? 10??????? 11
0.7706127 0.7665437 0.7624248 0.7502806 0.7459699 0.7629174 0.7587547
0.7523144 0.7480270 0.7436906 0.7809063
?????? 12??????? 13??????? 14
0.7847953 0.7847953 0.7809063
The first 7 estimates are ok but not the last 7. So it is related to R1
factor... but I don't understand why it is badly estimated.
If someone can help me...
Thanks a lot
Marc