[R-sig-ME] Predict for glmer
Thierry Onkelinx
th|erry@onke||nx @end|ng |rom |nbo@be
Mon Sep 28 13:25:04 CEST 2020
Dear Marc,
You probably want random <- er$R1[d$R1, ] + er$R2[d$R2, ]
Best regards,
ir. Thierry Onkelinx
Statisticus / Statistician
Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx using inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be
///////////////////////////////////////////////////////////////////////////////////////////
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
///////////////////////////////////////////////////////////////////////////////////////////
<https://www.inbo.be>
Op zo 27 sep. 2020 om 10:53 schreef Marc Girondot via R-sig-mixed-models <
r-sig-mixed-models using r-project.org>:
> 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
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list