[R] MaxLik estimation issues

[Arne Henningsen]

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