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]


Wayne (Yanwei) Zhang
Statistical Research
CNA
Phone: 312-822-6296
Email: yanwei.zh...@cna.com<mailto:yanwei.zh...@cna.com>




NOTICE:  This e-mail message, including any attachments and appended messages, 
is for the sole use of the intended recipients and may contain confidential and 
legally privileged information.
If you are not the intended recipient, any review, dissemination, distribution, 
copying, storage or other use of all or any portion of this message is strictly 
prohibited.
If you received this message in error, please immediately notify the sender by 
reply e-mail and delete this message in its entirety.

        [[alternative HTML version deleted]]

______________________________________________
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.

Reply via email to