Функция с оптимизированными параметрами не приближается к данным, использующим mle2 в R

Question

Итак, я пытался оптимизировать отношения Михаэлиса-Ментена с помощью распределения гамма-ошибок, чтобы смоделировать среднее значение некоторых данных, которые я собрал. Однако, как бы я ни оптимизировал функцию, самый низкий AI C, который я получаю, - это параметры, которые даже близко не подходят к данным. Можно ли как-нибудь решить эту проблему?

Вот мой код:

Сначала я создаю функцию максимального правдоподобия:

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

С извинениями за приведенный ниже код ...

Я сделал воспроизводимый пример, похожий на ваш, который, кажется, дает разумные результаты.

• Вы получаете какие-либо предупреждения о не удалось объединить / «достигнуто максимальное количество итераций»?
• в вашем коде, похоже, есть неиспользуемые / оставшиеся вещи о обработке; Это хорошая идея, но она будет работать только с интерфейсом формул (см. ниже)

Некоторые вспомогательные функции:

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

Fit (используя интерфейс формулы):

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