如何在ggplot2中绘制模型的mle2拟合以及数据? [英] How do I plot a mle2 fit of a model in ggplot2, along with the data?

问题描述

我为该模型创建了一个对数似然函数，并将其与mle2()中的起始值配合使用以拟合该模型(请参见打击)，但无法弄清楚如何将该模型拟合到数据的顶部ggplot2.我以前从未在此站点上发布过，所以我不确定将数据文件放在哪里，但是如果需要的话，我有一个文件可供复制.

I created a log likelihood function for the model, and use this with starting values in mle2() to fit the model (see blow), but can't figure out how to plot this model fit over top of the data in ggplot2. I've never posted on this site before, so I'm not sure where to put the data file, but I have one for reproducibility if needed.

我花了几天时间试图找到一个具体说明我需要做的事的例子，但找不到任何相关的东西.显然stat_smooth除了mle之外，具有最合适的选项，我无法在该模型中使用这些选项.这是一个渔业里克鱼库存补充模型，假设对数正态误差，适合于mle.

I have spent days trying to find an example of specifically what I need to do, and can't find anything relevant. Apparently stat_smooth has most fitting options except mle, none of which I can use for this model. This is a fisheries Ricker stock-recruitment model which is fit with mle assuming log-normal errors.

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绘制这些预测值，但出现错误:

Then, I tried to assign predicted values to a new df in order to plot these with geom_line, but got the error:

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

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

So, I tried to include the formula in mle2() before calling predict():

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

并收到错误:'*'错误(x = c(....):操作员需要一个或两个参数.

and got the error: 'Error in '*' (x=c(....):operator needs one or two arguments.

我只想要一个数据图(s& r)，并关联相关拟合.使用nls()和stat_smooth()时我没有问题，但必须使用mle来解决这个问题.

I just want a plot of the data (s & r), with the associated fit overlain. I have had no problem using nls() and stat_smooth() but must use mle to fit this.

推荐答案

初步:

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

有多种方法可以做到这一点.预测的主要区别在于我们是预测响应分布的中位数(等于对数正态时的几何均值)还是其均值...

There are a variety of ways to do this. The main difference in the predictions turns out to be whether we predict the median of the response distribution (equivalent to the geometric mean in the case of the log-Normal) or its mean ...

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

1. As a direct maximum likelihod estimate, with the log-mean equal to Ricker(A(t),a,b) [the sdlog parameter is estimated on the log scale; use lower bounds on a and b to avoid nuisance warning messages]

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

1. 作为对数线性模型(如果 log(Y)= loga + log(A)-b * A ，对两边求幂表示 Y = a * A * exp(-b * A);对数刻度上的常规误差与原始刻度上的对数常规误差相对应)

1. As a log-linear model (if log(Y) = loga + log(A) - b*A, exponentiating both sides shows Y = a*A*exp(-b*A); Normal errors on the log scale correspond to log-Normal errors on the original scale)

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

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

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

1. 为了进行比较，将 nls()拟合(这也应与 family = gaussian(link ="log") m3 等效)>:

1. For comparison, the nls() fit (this should also be equivalent to m3 with 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")和中位数/几何均值("mle2_med"和"loglm")的预测之间.差异的大小让我感到惊讶:从中位数到均值等于将所有预测乘以 exp(sdlog ^ 2/2)
• "mle2_med"和"loglm"预测是相同的(很难看到，但是绿色方块恰好在黄色三角形的顶部)
• GLM预测和内置的 geom_smooth()预测是相同的(它们应该是！)
• mle2和loglm模型的系数相同，直到变换为止:

• as stated above, the biggest difference in predictions is between predictions of the mean ("mle2") and of the median/geometric mean ("mle2_med" and "loglm"). The size of the difference surprised me: going from median to mean is equivalent to multiplying all predictions by exp(sdlog^2/2)
• the "mle2_med" and "loglm" predictions are identical (it's hard to see, but the green squares are exactly on top of the yellow triangles)
• the GLM prediction and the built-in geom_smooth() prediction are identical (they should be!)
• the coefficients of the mle2 and loglm models are identical, up to transformations:

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)   

这篇关于如何在ggplot2中绘制模型的mle2拟合以及数据?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！