Search code examples
rcurve-fittinglog-likelihood

How to get coefficients from a model in R?


I am trying to store coefficients from a mle model just like a normally do for other models. Surprisingly, it throws an error. Any other methods to store it?

set.seed(111)
y.size <- rnorm(100,4,1)
p <- rnorm(100,5,1)

df <- data.frame(p, y.size)

glae <- function(A,K,Ka,q, sigma, c) {
  lnqp <- if(q == 0) log(p) else ((p^q - 1)) / q
  y.pred <- ((A *((p*K/Ka)-1))-1)* lnqp + c
  ll <-   -sum(dnorm(y.size, mean = y.pred, sd = sigma, log=TRUE ))
  ll
}

mle2.model <- mle(glae, start = list(A = -2, K = 100, Ka = 200, q= 0.1, sigma = 0.1, c=3.8))
summary(mle2.model)

Maximum likelihood estimation

Call:

mle(minuslogl = glae, start = list(A = -2, K = 100, Ka = 200, q = 0.1, sigma = 0.1, c = 3.8))

Coefficients:
        Estimate Std. Error
A      12.511258        NaN
K      99.537509        NaN
Ka    200.231236 1684.30918
q     -32.959393        NaN
sigma 363.906677        NaN
c       8.135369   35.47982

-2 log L: 1363.185 


A.final <- summary(mle2.model)$Coefficients[1]

`Error in summary(mle2.model)$Coefficients : $ operator not defined for this S4 class

In addition: Warning message: In sqrt(diag(object@vcov)) : NaNs produced`


Solution

  • You can access this "mle" S4 object contents using the @ operator, e.g.

    library(bbmle)
    #> Loading required package: stats4
    
    set.seed(111)
    y.size <- rnorm(100,4,1)
    p <- rnorm(100,5,1)
    
    df <- data.frame(p, y.size)
    
    glae <- function(A,K,Ka,q, sigma, c) {
      lnqp <- if(q == 0) log(p) else ((p^q - 1)) / q
      y.pred <- ((A *((p*K/Ka)-1))-1)* lnqp + c
      ll <-   -sum(dnorm(y.size, mean = y.pred, sd = sigma, log=TRUE ))
      ll
    }
    
    mle2.model <- mle(glae, start = list(A = -2, K = 100, Ka = 200, q= 0.1, sigma = 0.1, c=3.8))
    coefs <- mle2.model@coef
    coefs
    #>          A          K         Ka          q      sigma          c 
    #>  12.511258  99.537509 200.231236 -32.959393 363.906677   8.135369
    

    Created on 2022-05-27 by the reprex package (v2.0.1)