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
Predict for glmer
2 messages · Marc Girondot, Thierry Onkelinx
1 day later
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 at 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 at 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 at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models