

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)) :
nonfinite finitedifference value [3]
Wayne (Yanwei) Zhang
Statistical Research
CNA
Phone: 3128226296
Email: [hidden email]<mailto: [hidden email]>
NOTICE: This email 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 email and delete this message in its entirety.
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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)) :
> nonfinite finitedifference 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
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Thomas,
Thanks a lot. This solves my problem.
Wayne (Yanwei) Zhang
Statistical Research
>CNA
Original Message
From: Thomas Lumley [mailto: [hidden email]]
Sent: Thursday, September 02, 2010 11:24 AM
To: Zhang,Yanwei
Cc: [hidden email]
Subject: Re: [R] Help on glm and optim
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)) :
> nonfinite finitedifference 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
NOTICE: This email 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 email and delete this message in its entirety.
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.

