[R] Difficulty in calculating MLE through NLM
Ravi Varadhan
RVaradhan at jhmi.edu
Thu Jul 2 16:43:14 CEST 2009
Madan,
I "did" tell you in my previous email about what you should do. Did you try
my suggestions? Did any of them work for you? If not, please send me the
details of what you tried and what the results were.
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
----------------------------------------------------------------------------
--------
-----Original Message-----
From: Madan Gopal Kundu [mailto:Madan.Kundu at ranbaxy.com]
Sent: Thursday, July 02, 2009 1:32 AM
To: Ravi Varadhan
Cc: nashjc at uottawa.ca
Subject: RE: [R] Difficulty in calculating MLE through NLM
Dear Prof. Vardhan,
Thank you very much for your answer. I will be very much thankful to you if
you can give me some suggestion. In my earlier mail, I described the coding
part in R. Now I will detail some theoretical aspect of my problem.
I am in effort of Multi Gamma Poisson Shrinkage (MGPS) estimate of Relative
Risk (RR) or Relative Report Rate (RRR) [Ref: DuMouchel 1999, Bayesian Data
Mining in Large Frequency Tables, With an Application to the FDA Spontaneous
Reporting System, The American Statistician, 53(3): 177-202.]
As earlier I have shared, the data (dummy one) contains different observed
frequency (N) and expected frequency (E) of different Adverse Event (AE) and
Drug combinations. I want to find out the Maximum Likelihood Estimate (MLE)
of alpha1, beta1, alpha2, beta2 and p from the following marginal
distribution of N (eq. 5 of the cited ref.):
Pr[N=n] = P*f(n; alpha1, beta1, E) + (1-P)* f(n; alpha2, beta2, E) Where,
f(n; alpha, beta, E)= [(1+beta/E)^(-n)]* [(1+E/beta)^(-alpha)]
*gamma(alpha+n)/gamma(alpha)*n!
I am eagerly look forward to have your suggestion in finding MLEs.
Thanks & Regards
Madan Gopal Kundu
Biostatistician, CDM, MACR, Ranbaxy Labs. Ltd.
Tel(O): +91 (0) 1245194045 - Mobile: +91 (0) 9868788406
http://mgkundu.webs.com/
-----Original Message-----
From: Ravi Varadhan [mailto:rvaradhan at jhmi.edu]
Sent: Wednesday, July 01, 2009 8:22 PM
To: Madan Gopal Kundu
Cc: r-help; nashjc at uottawa.ca
Subject: Re: [R] Difficulty in calculating MLE through NLM
Hi Madan,
You are trying to find the MLE of a binary mixture distribution. So, there
are constraints on the parameters, as you have indicated. The nlm()
function cannot handle constraints. I would recommend one of the following
functions (not necessarily in any order), all of which can handle
box-constraints:
1. optim(), with method="L-BFGS-B")
2. nlminb()
3. spg() in package "BB"
4. Write an EM algorithm; this automatically imposes all the required
constraints, and is guaranteed to converge (to a local maximum).
John Nash and I are writing a package called "optimx", which contains a
function called optimx() that integrates all the different optimization
functions for smooth nonlinear problems. Thsi would make it extremely easy
for you to do steps (1) - (3) above, using a single call to optimx(). Let
me know if you are interested.
Hope this helps,
Ravi.
____________________________________________________________________
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology School of Medicine Johns
Hopkins University
Ph. (410) 502-2619
email: rvaradhan at jhmi.edu
----- Original Message -----
From: Madan Gopal Kundu <Madan.Kundu at ranbaxy.com>
Date: Wednesday, July 1, 2009 9:25 am
Subject: [R] Difficulty in calculating MLE through NLM
To: r-help <r-help at stat.math.ethz.ch>
> Hi R-friends,
>
> Attached is the SAS XPORT file that I have imported into R using
> following code
> library(foreign)
> mydata<-read.xport("C:\\ctf.xpt")
> print(mydata)
>
> I am trying to maximize logL in order to find Maximum Likelihood
> Estimate (MLE) of 5 parameters (alpha1, beta1, alpha2, beta2, p) using
> NLM function in R as follows.
>
> # Defining Log likelihood - In the function it is noted as logL >
> library(stats) > loglike1<- function(x) + { + alpha1<-x[1] +
> beta1<-x[2] + alpha2<-x[3] + beta2<-x[4] + p<-x[5] + n<- mydata[3]
> + e<-mydata[4] + f1<-
> ((1+beta1/e)^(-n))*((1+e/beta1)^(-alpha1))*gamma(alpha1+n)/(gamma(n+1)
> *gamma(alpha1)) + f2<-
> ((1+beta2/e)^(-n))*((1+e/beta2)^(-alpha2))*gamma(alpha2+n)/(gamma(n+1)
> *gamma(alpha2))
> + logL=sum(log(p*f1+(1-p)*f2))
> + logL<- -logL
> + }
>
> # Supplying starting parameter values > theta<-c(.2041,.0582,
> 1.4150, 1.8380,0.0969)
>
> # Calculating MLE using NLM function
> > result<- nlm(loglike1, theta, hessian=TRUE, print.level=1)
>
> Now the problem is, this is not working as there is no improvement in
> final parameter estimate over starting values and NLM just stops just
> after 1 iteration with gradient value of all the 5 parameters as zero.
> I have tried other set of starting values, but then also I am getting
> final parameter estimates similar to starting values and iteration
> stops just after one step. When I check for warnings, R displays
> following kind of warnings:
>
> * NA/Inf replaced by maximum positive value
> * value out of range in 'gammafn'
>
> Please suggest what I should do. I am expecting the final MLE of
> alpha1, alpha2, beta1 and beta2 greater than 0 and P should lie
> between 0 to 1.
>
> Thanks & Regards,
> Madan Gopal Kundu
> Biostatistician, CDM, MACR, Ranbaxy Labs. Ltd.
> Tel(O): +91 (0) 1245194045 - Mobile: +91 (0) 9868788406
>
>
>
> (i) The information contained in this e-mail message is intended only
>
> for the confidential use of the recipient(s) named above. This
> message is privileged and confidential. If the reader of this message
> is not
>
> the intended recipient or an agent responsible for delivering it to
> the intended recipient, you are hereby notified that you have
> received this document in error and that any review, dissemination,
> distribution, or copying of this message is strictly prohibited. If
> you have received this communication in error, please notify us
> immediately by e-mail, and delete the original message.
>
> (ii) Madan.Kundu at ranbaxy.com confirms that Ranbaxy shall not be
> responsible if this email message is used for any indecent,
> unsolicited or illegal purposes, which are in violation of any
> existing laws and the same shall solely be the responsibility of
> Madan.Kundu at ranbaxy.com and that Ranbaxy shall at all times be
> indemnified of any civil and/ or criminal liabilities or consequences
> there.
> ______________________________________________
> R-help at r-project.org mailing list
>
> PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.
(i) The information contained in this e-mail message is intended only for
the confidential use of the recipient(s) named above. This message is
privileged and confidential. If the reader of this message is not the
intended recipient or an agent responsible for delivering it to the intended
recipient, you are hereby notified that you have received this document in
error and that any review, dissemination, distribution, or copying of this
message is strictly prohibited. If you have received this communication in
error, please notify us immediately by e-mail, and delete the original
message.
(ii) Madan.Kundu at ranbaxy.com confirms that Ranbaxy shall not be responsible
if this email message is used for any indecent, unsolicited or illegal
purposes, which are in violation of any existing laws and the same shall
solely be the responsibility of Madan.Kundu at ranbaxy.com and that Ranbaxy
shall at all times be indemnified of any civil and/ or criminal liabilities
or consequences there.
More information about the R-help
mailing list