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]]
Sophia Kyriakou <sophia.kyriakou17 <at> 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? > Thankssnip 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 r-sig-mixed-models at r-project.org instead ... Ben Bolker
yes Ben, this works indeed! Thanks a million!! On Mon, Mar 9, 2015 at 7:17 PM, Ben Bolker <bbolker at gmail.com> wrote:> Sophia Kyriakou <sophia.kyriakou17 <at> 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 r-sig-mixed-models at r-project.org > instead ... > > Ben Bolker > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >[[alternative HTML version deleted]]