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 ______________________________________________ R-help@r-project.org mailing list 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.
