Prof J C Nash (U30A)
2015-Mar-10 13:39 UTC
[R] Help with optim() to maximize log-likelihood
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]] >