help in maximum likelihood

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

help in maximum likelihood

Alaa Sindi-2

Hi all,

what is wrong with this code? I am trying to estimate the model parameters by maximizing the likelihood function and I am getting this warning


Warning message:
In nlm(fn, p = c(-50, 20), hessian = TRUE) :
  NA/Inf replaced by maximum positive value




x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
n <- c(59, 60, 62, 56, 63, 59, 62, 60)
fn <- function(p)
    sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))
            + log(choose(n, y)) ))
out <- nlm(fn, p = c(-50,20), hessian = TRUE)


Thanks
______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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 | Threaded
Open this post in threaded view
|

Re: help in maximum likelihood

J C Nash
It's useful to add "print.level=2" inside your call to find that there's
essentially nothing wrong.

Rvmmin doesn't give the msg and numDeriv gives a similar (BUT NOT
EXACTLY THE SAME!) hessian estimate.

It's almost always worthwhile turning on the diagnostic printing when
doing optimization, even if you throw away the output pretty well right
away, because it can often suggest whether the results are reasonable or
rubbish.

JN


On 16-03-01 03:09 PM, Alaa Sindi wrote:

>
> Hi all,
>
> what is wrong with this code? I am trying to estimate the model parameters by maximizing the likelihood function and I am getting this warning
>
>
> Warning message:
> In nlm(fn, p = c(-50, 20), hessian = TRUE) :
>   NA/Inf replaced by maximum positive value
>
>
>
>
> x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
> n <- c(59, 60, 62, 56, 63, 59, 62, 60)
> fn <- function(p)
>     sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))
>             + log(choose(n, y)) ))
> out <- nlm(fn, p = c(-50,20), hessian = TRUE)
>
>
> Thanks
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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 | Threaded
Open this post in threaded view
|

Re: help in maximum likelihood

Ravi Varadhan-2
In reply to this post by Alaa Sindi-2
There is nothing wrong with the optimization.  It is a warning message.  However, this is a good example to show that one should not simply dismiss a warning before understanding what it means.  The MLE parameters are also large, indicating that there is something funky about the model or the data or both.  In your case, there is one major problem with the data:  for the highest dose (value of x), you have all subjects responding, i.e. y = n.  Even for the next lower dose, there is almost complete response.  Where do these data come from? Are they real or fake (simulated) data?

Also, take a look at the eigenvalues of the hessian at the solution.  You will see that there is some ill-conditioning, as the eigenvalues are widely separated.

x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
n <- c(59, 60, 62, 56, 63, 59, 62, 60)

# note: there is no need to have the choose(n, y) term in the likelihood
fn <- function(p)
    sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))) )

out <- nlm(fn, p = c(-50,20), hessian = TRUE)

out

eigen(out$hessian)


Hope this is helpful,
Ravi



Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
Associate Professor,  Department of Oncology
Division of Biostatistics & Bionformatics
Sidney Kimmel Comprehensive Cancer Center
Johns Hopkins University
550 N. Broadway, Suite 1111-E
Baltimore, MD 21205
410-502-2619


        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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 | Threaded
Open this post in threaded view
|

Re: help in maximum likelihood

Alaa Sindi-2
Thank you very much prof. Ravi,

That was very helpful. Is there a way to get the t and p value for the coefficients?

Thanks
Alaa

> On Mar 2, 2016, at 10:05 AM, Ravi Varadhan <[hidden email]> wrote:
>
> There is nothing wrong with the optimization.  It is a warning message.  However, this is a good example to show that one should not simply dismiss a warning before understanding what it means.  The MLE parameters are also large, indicating that there is something funky about the model or the data or both.  In your case, there is one major problem with the data:  for the highest dose (value of x), you have all subjects responding, i.e. y = n.  Even for the next lower dose, there is almost complete response.  Where do these data come from? Are they real or fake (simulated) data?
>  
> Also, take a look at the eigenvalues of the hessian at the solution.  You will see that there is some ill-conditioning, as the eigenvalues are widely separated.
>  
> x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
> n <- c(59, 60, 62, 56, 63, 59, 62, 60)
>  
> # note: there is no need to have the choose(n, y) term in the likelihood
> fn <- function(p)
>     sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))) )
>  
> out <- nlm(fn, p = c(-50,20), hessian = TRUE)
>  
> out
>  
> eigen(out$hessian)
>  
>  
> Hope this is helpful,
> Ravi
>  
>  
>  
> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
> Associate Professor,  Department of Oncology
> Division of Biostatistics & Bionformatics
> Sidney Kimmel Comprehensive Cancer Center
> Johns Hopkins University
> 550 N. Broadway, Suite 1111-E
> Baltimore, MD 21205
> 410-502-2619
>  


        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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 | Threaded
Open this post in threaded view
|

Re: help in maximum likelihood

Ravi Varadhan-2
Standard error = sqrt(diag(solve(opt$hessian)))

Ravi

From: Alaa Sindi [mailto:[hidden email]]
Sent: Wednesday, March 02, 2016 3:22 PM
To: Ravi Varadhan <[hidden email]>
Cc: [hidden email]
Subject: Re: help in maximum likelihood

Thank you very much prof. Ravi,

That was very helpful. Is there a way to get the t and p value for the coefficients?

Thanks
Alaa
On Mar 2, 2016, at 10:05 AM, Ravi Varadhan <[hidden email]<mailto:[hidden email]>> wrote:

There is nothing wrong with the optimization.  It is a warning message.  However, this is a good example to show that one should not simply dismiss a warning before understanding what it means.  The MLE parameters are also large, indicating that there is something funky about the model or the data or both.  In your case, there is one major problem with the data:  for the highest dose (value of x), you have all subjects responding, i.e. y = n.  Even for the next lower dose, there is almost complete response.  Where do these data come from? Are they real or fake (simulated) data?

Also, take a look at the eigenvalues of the hessian at the solution.  You will see that there is some ill-conditioning, as the eigenvalues are widely separated.

x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
n <- c(59, 60, 62, 56, 63, 59, 62, 60)

# note: there is no need to have the choose(n, y) term in the likelihood
fn <- function(p)
    sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))) )

out <- nlm(fn, p = c(-50,20), hessian = TRUE)

out

eigen(out$hessian)


Hope this is helpful,
Ravi



Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
Associate Professor,  Department of Oncology
Division of Biostatistics & Bionformatics
Sidney Kimmel Comprehensive Cancer Center
Johns Hopkins University
550 N. Broadway, Suite 1111-E
Baltimore, MD 21205
410-502-2619



        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
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.