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: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml ---------------------------------------------------------------------------- -------- -----Original Message----- From: Madan Gopal Kundu [mailto:madan.ku...@ranbaxy.com] Sent: Thursday, July 02, 2009 1:32 AM To: Ravi Varadhan Cc: nas...@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:rvarad...@jhmi.edu] Sent: Wednesday, July 01, 2009 8:22 PM To: Madan Gopal Kundu Cc: r-help; nas...@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: rvarad...@jhmi.edu ----- Original Message ----- From: Madan Gopal Kundu <madan.ku...@ranbaxy.com> Date: Wednesday, July 1, 2009 9:25 am Subject: [R] Difficulty in calculating MLE through NLM To: r-help <r-h...@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.ku...@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.ku...@ranbaxy.com and that Ranbaxy shall at all times be > indemnified of any civil and/ or criminal liabilities or consequences > there. > ______________________________________________ > R-help@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.ku...@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.ku...@ranbaxy.com and that Ranbaxy shall at all times be indemnified of any civil and/ or criminal liabilities or consequences there. ______________________________________________ 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.