Help with optim() to maximize log-likelihood

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

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

bbolker
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 [hidden email]
instead ...

  Ben Bolker

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

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

Sophia Kyriakou
yes Ben, this works indeed! Thanks a million!!

On Mon, Mar 9, 2015 at 7:17 PM, Ben Bolker <[hidden email]> 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 [hidden email]
> instead ...
>
>   Ben Bolker
>
> ______________________________________________
> [hidden email] 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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

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

Prof J C Nash (U30A)
In reply to this post by Sophia Kyriakou
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, [hidden email] wrote:

> Message: 12
> Date: Mon, 9 Mar 2015 16:18:06 +0200
> From: Sophia Kyriakou <[hidden email]>
> To: [hidden email]
> Subject: [R] Help with optim() to maximize log-likelihood
> Message-ID:
> <[hidden email]>
> 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]]
>

______________________________________________
[hidden email] 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.