[R-sig-ME] Predict for glmer
Ben Bolker
bbo|ker @end|ng |rom gm@||@com
Fri Oct 9 21:10:29 CEST 2020
I don't know if anyone ever answered this ... you need to construct
the RE model matrix without an intercept, i.e.
random <- rowSums(model.matrix(~ R1 -1 , data = d) * er$R2[d$R2, ])
FWIW you can also use getME(.,"X") and getME(.,"Z") to retrieve the
fixed- and random-effects model matrices (I know that's not the point
here ...)
cheers
Ben Bolker
On 9/27/20 4:52 AM, Marc Girondot via R-sig-mixed-models wrote:
> 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
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list