[R] Help on glm and optim

[Thomas Lumley]

On Thu, 2 Sep 2010, Zhang,Yanwei wrote:

> Dear all,
>
> I'm trying to use the "optim" function to replicate the results  from the "glm" using an example from the help page of "glm", but I could not get the "optim" function to work. Would you please point out where I did wrong? Thanks a lot.
>
> The following is the code:
>
> # Step 1: fit the glm
> clotting <- data.frame(
>    u = c(5,10,15,20,30,40,60,80,100),
>    lot1 = c(118,58,42,35,27,25,21,19,18),
>    lot2 = c(69,35,26,21,18,16,13,12,12))
> fit1 <- glm(lot1 ~ log(u), data=clotting, family=Gamma)
>
> # Step 2: use optim
> # define loglikelihood function to be maximized over
> # theta is a vector of three parameters: intercept, cofficient for log(u) and dispersion parameter
> loglik <- function(theta,data){
>        E <- 1/(theta[1]+theta[2]*log(data$u))
>        V <- theta[3]*E^2
>        loglik <- sum(dgamma(data$lot1,shape=1/theta[3],rate=1/(E*theta[3]),log=T))
>        return(loglik)
> }
>
> # use the glm result as initial values
> theta <- c(as.vector(coef(fit1)),0.002446059)
> fit2 <- optim(theta, loglik,  clotting, gr = NULL, hessian = TRUE,
>        control = list(fnscale = -1))
>
> # Then I got the following error message:
> Error in optim(theta, loglik, clotting, gr = NULL, hessian = TRUE, control = list(fnscale = -1)) :
>  non-finite finite-difference value [3]
>

If you use trace(loglik, tracer=quote(print(theta))) to trace the inputs to loglik() you will find that it is being called with negative values of theta[3] to get finite differences.   One fix is to reparametrize and use the log scale rather than the scale as a parameter.

    -thomas


Thomas Lumley
Professor of Biostatistics
University of Washington, Seattle