Dear all,
In order to be sure I understand glmm correctly, I try to reproduce by
hand a simple result. Here is a reproducible code. The questions are in
_________________
Of course I have tried to find the solution using internet but I was not
able to find a solution. I have also tried to follow glmer but it is
very complicated code!
Thanks for any help.
Marc
# Generate set of df with nb successes and failures
# and ID being A, B or C (the random effect)
# and x being the fixed effect
set.seed(1)
df <- rbind(matrix(data = c(sample(x=5:30, size=40, replace = TRUE),
rep(10, 40)), ncol=2),
matrix(data = c(sample(x=10:30, size=40, replace = TRUE),
rep(10, 40)), ncol=2),
matrix(data = c(sample(x=20:30, size=40, replace = TRUE),
rep(10, 40)), ncol=2))
ID <- as.factor(c(rep("A", 40), rep("B", 40), rep("C", 40)))
x <- c(runif(40, min=10, max=30), runif(40, min=20, max=30), runif(40,
min=40, max=60))
x <- (x-min(x))/(max(x)-min(x))
# In g0, I have the results of the glmm
library(lme4)
g0 <- glmer(formula = df ~ x + (1 | ID), family =
binomial(link="logit"), nAGQ=1)
-logLik(g0) # 'log Lik.' 268.0188 (df=3)
# I get the fitted parameters
fixep <- fixef(g0)
par <- getME(g0, c("theta","beta"))
# _______________________________________________________________________
# Question 1: how theta is converted into the specific effect on
(intercept) for the random effect ?
# Then how a theta parameter is converted into intercepts?
# _______________________________________________________________________
intercepts <- ranef(g0)$ID
# This part is ok, the predict is correct
pfit <- 1-c(1/(1+exp(fixep["(Intercept)"]+intercepts["A",
1]+x[ID=="A"]*fixep["x"])),
1/(1+exp(fixep["(Intercept)"]+intercepts["B",
1]+x[ID=="B"]*fixep["x"])),
1/(1+exp(fixep["(Intercept)"]+intercepts["C", 1]+x[ID=="C"]*fixep["x"])))
predict(g0, type = "response")
# _______________________________________________________________________
# Why I obtain 266.4874 and not 268.0188 as in -logLik(g0)?
# _______________________________________________________________________
-sum(dbinom(x=df[, 1], size=df[, 1]+df[, 2], prob=pfit, log=TRUE)) #
266.4874
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
