[R] Help with optim() to maximize log-likelihood

[Sophia Kyriakou]

hello, I am using the optim function to maximize the log likelihood of a
generalized linear mixed model and I am trying to replicate glmer's
estimated components. If I set both the sample and subject size to q=m=100
I replicate glmer's results for the random intercept model with parameters
 beta=-1 and sigma^2=1. But if I change beta to 2 glmer works and optim
gives me the error message "function cannot be evaluated at initial
parameters".

If anyone could please help?
Thanks

 # likelihood function
 ll <- function(x,Y,m){
 beta <- x[1]
 psi <- x[2]
 q <- length(Y)
  p <- 20
 rule20 <- gaussHermiteData(p)
 wStar <- exp(rule20$x * rule20$x + log(rule20$w))
 # Integrate over(-Inf, +Inf) using adaptive Gauss-Hermite quadrature
 g <- function(alpha, beta, psi, y, m) {-y+m*exp(alpha + beta)/(1 +
exp(alpha + beta)) + alpha/exp(psi)}
 DDfLik <- deriv(expression(-y+m*exp(alpha + beta)/(1 + exp(alpha + beta))
+ alpha/exp(psi)),
 namevec = "alpha", func = TRUE,function.arg = c("alpha", "beta", "psi",
"y", "m"))
   int0 <- rep(NA,q)
 piYc_ir <- matrix(NA,q,p)
 for (i in 1:q){
 muHat <- uniroot(g, c(-10, 10),extendInt ="yes", beta = beta, psi = psi, y
= Y[i], m = m)$root
 jHat <- attr(DDfLik(alpha = muHat, beta, psi, Y[i], m), "gradient")
 sigmaHat <- 1/sqrt(jHat)
 z <- muHat + sqrt(2) * sigmaHat * rule20$x
 piYc_ir[i,] <-
choose(m,Y[i])*exp(Y[i]*(z+beta))*exp(-z^2/(2*exp(psi)))/((1+exp(z+beta))^m*sqrt(2*pi*exp(psi)))
 int0[i] <- sqrt(2)*sigmaHat*sum(wStar*piYc_ir[i,])
 }
 ll <- -sum(log(int0))
 ll
 }

 beta <- 2
 sigma2 <- 1
 m <- 100
 q <- 100

 cl <- seq.int(q)
 tot <- rep(m,q)

 set.seed(123)
 alpha <- rnorm(q, 0, sqrt(sigma2))
 Y <- rbinom(q,m,plogis(alpha+beta))

 dat <- data.frame(y = Y, tot = tot, cl = cl)
 f1 <- glmer(cbind(y, tot - y) ~ 1 + (1 | cl), data = dat,family =
binomial(),nAGQ = 20)
 betaH <- summary(f1)$coefficients[1]
 sigma2H <- as.numeric(summary(f1)$varcor)
 thetaglmer <- c(betaH,sigma2H)

 logL <- function(x) ll(x,Y,m)
 thetaMLb <- optim(c(plogis(sum(Y/m)),log(sigma2H)),fn=logL)$par
 Error in optim(c(plogis(sum(Y/m)), log(sigma2H)), fn = logL) :  function
cannot be evaluated at initial parameters

thetaglmer
[1] 2.1128529 0.8311484
 (thetaML <- c(thetaMLb[1],exp(thetaMLb[2])))

	[[alternative HTML version deleted]]