[R] Try reproduce glmm by hand

[Marc Girondot]

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