Help with optim() to maximize log-likelihood

4 messages
Open this post in threaded view
|

Help with optim() to maximize log-likelihood

 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]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

Re: Help with optim() to maximize log-likelihood

 Sophia Kyriakou gmail.com> writes: > > 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  snip to make gmane happy. It looks like you're getting floating-point under/overflow.  If you do all the computations on the log scale first and then exponentiate, it seems to work, i.e.:         piYc_ir[i,] <- lchoose(m,Y[i]) + Y[i]*(z+beta) + (-z^2/(2*exp(psi))) -             m*(log1p(exp(z+beta))) - 0.5*(log(2*pi)+psi)         piYc_ir[i,] <- exp(piYc_ir[i,]) follow-ups should probably go to [hidden email] instead ...   Ben Bolker ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.