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

Question

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)

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

标签： rstatisticsdata-modeling

Answer 1

为下面的脑残代码道歉......

我做了一个类似于你的可重复的例子，似乎给出了合理的结果。

• 您是否收到有关无法收敛/“达到最大迭代次数”的任何警告？
• 您的代码中似乎有一些关于处理的未使用/剩余的东西；这是个好主意，但仅适用于公式界面（见下文）

一些辅助函数：

## 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)
}

模拟数据：

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

拟合（使用公式界面）：

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(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