如何在 ggplot2 中绘制模型的 mle2 拟合以及数据？

Question

我为模型创建了一个对数似然函数，并将其与 mle2() 中的起始值一起使用以拟合模型（见下文），但无法弄清楚如何在ggplot2 中的数据。我以前从未在这个网站上发过帖子，所以我不确定将数据文件放在哪里，但如果需要的话，我有一个可以重现的文件。

我花了几天时间试图找到一个具体说明我需要做什么的例子，但找不到任何相关的东西。显然 stat_smooth 具有最合适的选项，除了 mle，none 我可以将其用于此模型。这是一个渔业 Ricker 种群补充模型，它与假设对数正态误差的 mle 相吻合。

LL函数：

Ricker.LL <- function(a,b) {
  wf<-read.csv("wf_SR data.csv",sep=",",header=T)
  s <- wf$Adult.CPUE.t.1
  r <- wf$YOY.CPUE
  model.pred <- a*s*exp(-(b)*s)
  ndata <- length(s)
  NLL <- -sum(dlnorm(x=s,meanlog=model.pred,sdlog=1,log=TRUE))
  return(NLL)
}

mle2 适合：

mle2(minuslogl=Ricker.LL,start=list(a=0.4515,b=0.2665),method="Nelder-Mead",lower=-Inf,upper=Inf)

然后，我尝试将预测值分配给新的 df，以便用 geom_line 绘制这些值，但出现错误：

dat <- predict(fit)

Error : object of type 'symbol' is not subsettable
Error in gfun(object, newdata = newdata, location = location, op = "predict") : 
  can only use predict() if formula specified

因此，我尝试在调用 predict() 之前将公式包含在 mle2() 中：

fit<-mle2(YOY.CPUE~a*Adult.CPUE.t.1*exp(-(b)*Adult.CPUE.t.1),data=wf,start=list(a=0.4515,b=0.2665))

并得到错误：'Error in '*' (x=c(....):operator needs one or two arguments.

我只想要数据图 (s & r)，并覆盖相关拟合。我使用 nls() 和 stat_smooth() 没有问题，但必须使用 mle 来适应这个。

Answer 1

预赛：

library(bbmle)
library(ggplot2); theme_set(theme_bw())
rickerfun <- function(x,a,b) {
    a*x*exp(-b*x)
}

有多种方法可以做到这一点。预测的主要区别在于我们是预测响应分布的 中位数 （相当于对数正态情况下的几何平均值）还是它的平均值 ...

作为直接最大似然估计，对数均值等于 Ricker(A(t),a,b) [sdlog 参数在对数尺度上估计；在 a 和 b 上使用下限以避免令人讨厌的警告消息]

m1 <- mle2(YOY.CPUE ~ dlnorm(meanlog=log(rickerfun(Adult.CPUE.t.1,a,b)),
                         sdlog=exp(logsdlog)),
       method="L-BFGS-B",
       lower=c(a=1e-2,b=1e-2,logsdlog=-10),
       start=list(a=1,b=1,logsdlog=0),
       data=dat)

还需要从 mle2() 获得预测：

slnorm <- function(meanlog, sdlog) {
   list(title="Log-normal",
        median=exp(meanlog),
        mean=exp(meanlog+sdlog^2/2))
}

作为对数线性模型（如果log(Y) = loga + log(A) - b*A，两边取幂显示Y = a*A*exp(-b*A)；对数尺度上的正态误差对应于原始尺度上的对数正态误差）

 m2 <- lm(log(YOY.CPUE) ~ Adult.CPUE.t.1+ offset(log(Adult.CPUE.t.1)),
          data=dat)

作为具有对数 link 和 Gamma 响应的广义线性模型[Gamma 分布与对数正态分布具有相同的均值-方差关系，并且通常是充分的近似值]

 m3 <- glm(YOY.CPUE ~ Adult.CPUE.t.1+ offset(log(Adult.CPUE.t.1)),
      data=dat,
      family=Gamma(link="log"))

为了比较，nls() 适合（这也应该等同于 m3 和 family=gaussian(link="log"):

m4 <- nls(YOY.CPUE ~ a*Adult.CPUE.t.1*exp(-b*Adult.CPUE.t.1),
      start=list(a=0.45, b=0.27),
      data=dat)

计算所有模型的预测：

predframe <- data.frame(Adult.CPUE.t.1=seq(0,5.5,length=51))
predframe$mle2 <- predict(m1,newdata=predframe)
predframe$mle_med <- predict(m1,newdata=predframe,location="median")
predframe$loglm <- exp(predict(m2,newdata=predframe))
predframe$glm <- predict(m3,newdata=predframe,type="response")
predframe$nls <- predict(m4,newdata=predframe)

为了 ggplot 方便而融化：

predframe_m <- reshape2::melt(predframe,id.var="Adult.CPUE.t.1",
                              variable.name="model",
                              value.name="YOY.CPUE")

剧情：

library(ggplot2)
ggplot(dat,aes(Adult.CPUE.t.1,YOY.CPUE))+ geom_point() +
    geom_smooth(method="glm",
                formula=y~x + offset(log(x)),
                method.args=list(family=Gamma(link="log")))+
    geom_point(data=predframe_m,aes(colour=model,shape=model))

带回家的留言：

如上所述，预测中最大的差异在于均值 ("mle2") 和 median/geometric 均值（"mle2_med" 和 "loglm"）的预测之间。差异的大小让我感到惊讶：从中位数到均值相当于将所有预测乘以 exp(sdlog^2/2)
"mle2_med" 和 "loglm" 的预测是相同的（很难看出来，但是绿色方块恰好在黄色三角形上面）
GLM 预测和内置 geom_smooth() 预测是相同的（它们应该是！）
mle2 和 loglm 模型的系数相同，直到变换：

all.equal(coef(m1)[["a"]],exp(coef(m2)[["(Intercept)"]]), tol=1e-5)
all.equal(coef(m1)[["b"]],-coef(m2)[["Adult.CPUE.t.1"]], tol=1e-4)

Answer 2

这是考虑了年龄多样性的模型。我不得不删除下限以使其适合，并且出现 'longer object length not multiple of shorter object length' 错误，但它似乎仍然有效。再次感谢您的帮助。

rickerH <- function(x,z,a,b,f) {
  a*x*exp(-b*x)*exp(f*z)
}

fit<-mle2(YOY.CPUE~dlnorm(meanlog=log(rickerH(Adult.CPUE.t.1,H.t.1,a,b,f)),
sdlog=exp(logsdlog)),
          # method="L-BFGS-B",lower=c(a=0.01,b=0.01,f=0.01,logsdlog=-10),
          start=list(a=1,b=1,f=1,logsdlog=0),
          data=wf)
coef(fit)
slnorm <- function(meanlog,sdlog) {
  list(title="Log-normal",
       median=exp(meanlog),
       mean=exp(meanlog+sdlog^2/2))
}

# predframe <- data.frame(Adult.CPUE.t.1=seq(0,5.5,length=51))
predframe$mle2_f <- predict(fit,newdata=predframe)
# predframe$mle_med <- predict(fit,newdata=predframe,location="median")

predframe_melt <- reshape2::melt(predframe,id.var="Adult.CPUE.t.1",variable.name="model",
value.name="YOY.CPUE")

ggplot(wf,aes(Adult.CPUE.t.1,YOY.CPUE))+
  geom_point()+
  # geom_smooth(method="glm",formula=y~x+offset(log(x)),
method.args=list(family=Gamma(link="log")))+
  geom_line(data=predframe_melt,aes(color=model,shape=model),size=2)+
  theme_bw()

如何在 ggplot2 中绘制模型的 mle2 拟合以及数据？

How do I plot a mle2 fit of a model in ggplot2, along with the data?

r

ggplot2

mle