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

[Prof J C Nash (U30A)]

1) It helps to include the require statements for those of us who work
outside your particular box.
   lme4 and (as far as I can guess) fastGHQuad
are needed.

2) Most nonlinear functions have domains where they cannot be
evaluated. I'd be richer than Warren Buffett if I got $5 for
each time someone said "your optimizer doesn't work" and I
found   f(start, ...) was NaN or Inf, as in this case, i.e.,

 start <- c(plogis(sum(Y/m)),log(sigma2H))
 cat("starting params:")
 print(start)
 tryf0 <- ll(start,Y,m)
 print(tryf0)


It really is worthwhile actually computing your function at the initial
parameters EVERY time. (Or turn on the trace etc.)

JN

On 15-03-10 07:00 AM, r-help-request at r-project.org wrote:
> Message: 12
> Date: Mon, 9 Mar 2015 16:18:06 +0200
> From: Sophia Kyriakou <sophia.kyriakou17 at gmail.com>
> To: r-help at r-project.org
> Subject: [R] Help with optim() to maximize log-likelihood
> Message-ID:
> 	<CAO4gA+qokumHoZwvbU7EY3xaBBo2LnQjRcWxQkzcHm3U9OZ6kw at mail.gmail.com>
> Content-Type: text/plain; charset="UTF-8"
> 
> 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]]
>