具有优化参数的函数无法使用R中的mle2接近数据

Question

因此，我一直在尝试使用伽马误差分布来优化Michaelis-Menten关系，以对我收集的某些数据的平均值进行建模。但是，无论我如何优化功能，获得的最低AIC都是针对那些甚至不接近数据的参数。有什么办法可以解决这个问题？

这是我的代码：

我首先创建一个最大似然函数：

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

然后，我使用mle2()包中的bbmle函数编写了优化函数：

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

虽然a = 100和b = 260的眼球参数非常适合我的数据，但通常会将参数优化为a = 242和b = 182，结果是

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)

长话短说，如何确定优化的模型实际上遍历我的数据？

Answer 1

[为下面的转储代码道歉...

我做了一个与您相似的可复制示例，似乎给出了合理的结果。

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

一些帮助功能：

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