[R] Try reproduce glmm by hand

[Viechtbauer, Wolfgang (NP)]

Dear Marc,

Plugging the BLUPs into the likelihood isn't going to give you the ll of the model. You need to integrate over the random effect.

intfun <- function(nu, xi, mi, pred, theta)
   dbinom(xi, xi+mi, plogis(nu)) * dnorm(nu, pred, theta)

lli <- rep(NA, nrow(df))

for (i in 1:nrow(df)) {
   lli[i] <- log(integrate(intfun, xi=df[i,1], mi=df[i,2], pred=fixep[1] + fixep[2]*x[i], theta=par$theta, lower=-Inf, upper=Inf)$value)
}

-sum(lli)

Best,
Wolfgang

> -----Original Message-----
> From: R-help <r-help-bounces using r-project.org> On Behalf Of Marc Girondot via R-
> help
> Sent: Saturday, December 2, 2023 17:36
> To: Marc Girondot via R-help <r-help using r-project.org>
> Subject: [R] Try reproduce glmm by hand
>
> 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