[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