[R] MaxLik estimation issues

[Arne Henningsen]

Dear Filipe

On 14 October 2013 17:42, Filipe Ribeiro <flipjribeiro at hotmail.com> wrote:
> Dear Arne,
>
> First of all, thank you very much for your reply. Secondly, please accept my
> apologies for only give you an answer now, but I was moving from one country
> to another and only today I was able to get back to work.
> I did not write the model specification because it was just a trial, but my
> main idea was to apply the "Nelder-Mead" model.
> Since the beginning my guess was that something is wrong with the
> log-likelihood function, however, I don't find any problem with the same
> exact function applying the "optim" function:
>
> loglike.GGompi <- function(theta,age,deaths,exposures) {
>   alpha <- exp(theta[1])
>   beta <- exp(theta[2])
>   gamma <- exp(theta[3])
>   first <- alpha*exp(beta*age)
>   second <- 1+(((alpha*gamma)/beta)*(exp(beta*age)-1))
>   mu <- first/second
>   llk <- -sum((deaths * log(mu)) + (- mu*exposures))
>   return(llk)
> }
>
>
> fit1 <- optim(par=c(-4.1402817, -0.6375773, -1.6945914),
>               loglike.GGompi,
>               age=0:6,
>               deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
>               exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
>                           231951.64, 69502.65, 15798.72))
>
> exp(fit1$par)
>
>> exp(fit1$par)
> [1] 0.01585683 0.53471945 0.25368426
>
> Due to this results I can't understand why...

Please note that optim() is (by default) *minimizing* the provided
function, while maxLik() is *maximizing* the provided log-likelihood
function.

Cheers,
Arne



> A 26/09/2013, às 05:28, Arne Henningsen escreveu:
>
> Dear Filipe
>
> On 25 September 2013 14:23, Filipe Ribeiro <flipjribeiro at hotmail.com> wrote:
>
> Hello everybody!
>
>
> I'm having some trouble to compute maximum likelihood
>
> estimations using maxLik package and I hope that you
>
> could give me a hint.
>
> The main problem is that I'm not able to get a result not
>
> even close to the ones given by glm() directly, and the
>
> second one is: "Error in maxNRCompute(fn = logLikAttr,
>
> fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in
>
> gradient".
>
>
> The codes:
>
> loglike.GGompiMaxLik <- function(theta,age,deaths,exposures) {
>
> alpha <- exp(theta[1])
>
> beta <- exp(theta[2])
>
> gamma <- exp(theta[3])
>
> first <- alpha*exp(beta*age)
>
> second <- 1+(((alpha*gamma)/beta)*(exp(beta*age)-1))
>
> mu <- first/second
>
> llk <- -sum((deaths * log(mu)) + (- mu*exposures))
>
> return(llk)
>
> }
>
>
>
> fit1 <- maxLik(loglike.GGompiMaxLik,
>
>              age=0:6,
>
>              deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
>
>              exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
>
> 231951.64, 69502.65, 15798.72),
>
>              start=c(-4.1402817, -0.6375773, -1.6945914))
>
>
> Do you know how I can solve this problem?
>
>
> You did not write which model specification you want to estimate but I
> am pretty sure that something in your log-likelihood function is
> incorrect. The log-likelihood value at the starting values of the
> parameters is so large that R even cannot calculate the likelihood
> value:
>
> a <- loglike.GGompiMaxLik(c(-4.1402817, -0.6375773, -1.6945914), age=0:6,
>
> +     deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
> +     exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
> +      231951.64, 69502.65, 15798.72))
>
> a
>
> [1] 580365.2
>
> exp(a)
>
> [1] Inf
>
> In the second iteration, the first parameter gets so small (large in
> absolute terms, -5e+10) that the log-likelihood value become extremely
> (numerically infinitely) large and the gradients cannot be computed
> (by the finite-difference method):
>
> fit1 <- maxLik(loglike.GGompiMaxLik,
>
> +     age=0:6,
> +     deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
> +     exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
> +         231951.64, 69502.65, 15798.72),
> +     start=c(-4.1402817, -0.6375773, -1.6945914))
> Iteration 2
> Parameter:
> [1] -5.174233e+10 -3.839076e+02  5.988668e+00
> Gradient:
>     [,1] [,2] [,3]
> [1,]  NaN  NaN  NaN
> Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad,
> hessOrig = hess,  :
>  NA in gradient
>
> b <- loglike.GGompiMaxLik(c(-5.174233e+10, -3.839076e+02, 5.988668e+00),
> age=0:6,
>
> +     deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
> +     exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
> +         231951.64, 69502.65, 15798.72))
>
> b
>
> [1] Inf
>
> Please note that you can also find hints and ask questions about the
> maxLik package in the forums at maxLik's R-Forge site:
>
> https://r-forge.r-project.org/projects/maxlik/
>
> ...and please do not forget to cite the maxLik package in your publications:
>
> http://cran.r-project.org/web/packages/maxLik/citation.html
>
> Best wishes,
> Arne



-- 
Arne Henningsen
http://www.arne-henningsen.name