CODEWITHSUNDEEP
rstatisticsdata-modeling
Francine van Neerbos
Francine van Neerbos

Reputation: 21

Function with optimized parameters does not come close to data using mle2 in R

So I've been trying to optimize a Michaelis-Menten relationship with a gamma error distribution to model the average of some data that I have collected. However, no matter how I optimize the function, the lowest AIC I get is for parameters that don't even come close to the data. Is there any way I could solve this?

Here's my code:

I first create a maximum likelihood function:

MicNLL <- function(a,b){
  #a=150.6727
  #b=319.7007 optim val
  top <- a*x
  bot <- b+x
  Mic <- top/bot
  nll <- -sum(dgamma(y, shape=(Mic^2/var(x)), scale=(var(x)/Mic), log=TRUE))
  return(nll)
}

Then I have written the optimization function using the mle2() function in the bbmle package:

MN <- mle2(minuslogl = MicNLL, parameters=list(a~Treatment, b~Treatment), start=list(a=100,b=260), data=list(x=NSug3$VolpulT, y=NSug3$SugarpugT), control=list(maxit=1e4), method="SANN", hessian=T)

MN 
AICMN <- (2*2)-(2*logLik(MN))
AICMN

While the eyeballed parameters of a=100 & b=260 would fit nicely with my data, it usually optimizes the parameters to a=242 & b=182, which results 

Michealis <- function(a, b, x){
  top <- a*x
  bot <- b+x
  Mic <- top/bot
  return(Mic)
}
ggplot(NSug3, aes(x=VolpulT, y=SugarpugT))+
  geom_point(stat="identity", size=0.8)+
  theme_classic()+
  ggtitle("help")+
  ylab("Sugar concentration")+
  xlab("Volume per Extra floral nectary")+
  stat_function(fun= Michealis, args=c(a=100, b=260), colour="Orange", size=0.725)+
  stat_function(fun= Michealis, args=c(a=MN@coef[[1]], b=MN@coef[[2]]), colour="Red", size=0.725)

data model

So long story short, how can I make sure my optimized model actually runs through my data?

Upvotes: 0

Views: 180

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226057

With apologies for brain-dumped code below ...

I made a reproducible example similar to yours that seems to give reasonable results.

  • are you getting any warnings about failure to converge/"max iterations reached"?
  • there seems to be some unused/leftover stuff about Treatment in your code; that's a good idea, but will only work with the formula interface (see below)

Some helper functions:

## Gamma parameterized by mean and variance
## m = a*s, v = a*s^2 -> s=v/m; a=m^2/v
rgamma2 <- function(n, m, v) {
    rgamma(n, shape=m^2/v, scale=v/m)
}
dgamma2 <- function(x, m, v, log=FALSE) {
    dgamma(x, shape=m^2/v, scale=v/m, log=log)
}
sgamma2<- function(m, v) {   ## for predict()
    list(title="Gamma", mean=m, sd=sqrt(v))
}
mm <- function(x, a=100, b=260) {
    a*x/(b+x)
}

Simulate data:

set.seed(101)
x <- rlnorm(100,meanlog=4,sdlog=1)
dd <- data.frame(x,y=rgamma2(100,m=mm(x), v= 100))

Fit (using formula interface):

library(bbmle)
m1 <-mle2(y~dgamma2(m=mm(x,a,b),v=exp(logv)),
     start=list(a=50,b=200,logv=0),
     data=dd,
     control=list(maxit=1000))

Plot results:

plot(y~x,data=dd)
lines(sort(dd$x),mm(sort(dd$x)),col=2)     ## true
lines(sort(dd$x),sort(predict(m1)),col=3)  ## predicted

Upvotes: 1

Related Questions