【问题标题】:Using mle2 for parameter estimates with error and predict使用 mle2 进行有误差的参数估计并进行预测
【发布时间】:2020-02-26 06:24:01
【问题描述】:

我正在使用mle2 来估计非线性模型的参数,并且我想要估计参数估计周围的误差(标准误差)。同样,我想使用该模型然后使用 newdata 进行预测,并且在此过程中的几个步骤中我遇到了问题(错误)。

这是数据:

table<- "ac.temp performance
1      2.17   47.357923
4      2.17  234.255317
7      2.17  138.002633
10     2.17  227.545902
13     2.17   28.118072
16     9.95  175.638448
19     9.95  167.392218
22     9.95  118.162747
25     9.95  102.770622
28     9.95  191.874867
31    16.07  206.897159
34    16.07   74.741619
37    16.07  127.219884
40    16.07  208.231559
42    16.07   89.544612
44    20.14  314.946107
47    20.14  290.994063
50    20.14  243.322497
53    20.14  192.335133
56    20.14  133.841776
58    23.83  139.746673
61    23.83  224.135993
64    23.83  126.726493
67    23.83  246.443386
70    23.83  163.019896
83    28.04    4.614154
84    28.04    2.851866
85    28.04    2.935584
86    28.04  153.868415
87    28.04  103.884295
88    30.60    0.000000
89    29.60    0.000000
90    30.30    0.000000
91    29.90    0.000000
92    30.80    0.000000
93    28.90    0.000000
94    30.00    0.000000
95    30.20    0.000000
96    30.40    0.000000
97    30.70    0.000000
98    27.90    0.000000
99    28.60    0.000000
100   28.60    0.000000
101   30.40    0.000000
102   30.60    0.000000
103   29.70    0.000000
104   30.50    0.000000
105   29.30    0.000000
106   28.90    0.000000
107   30.40    0.000000
108   30.20    0.000000
109   30.10    0.000000
110   29.50    0.000000
111   31.00    0.000000
112   30.60    0.000000
113   30.90    0.000000
114   31.10    0.000000"

perfdat<- read.table(text=table, header=T)

首先,我必须为我的非线性模型设置几个固定参数,以了解动物在温度方面的表现

pi = mean(subset(perfdat, ac.temp<5)$performance)
ti = min(perfdat$ac.temp)

定义x 变量(温度)

t = perfdat$ac.temp

为非线性模型创建函数

tpc = function(tm, topt, pmax) {
    perfmu = pi+(pmax-pi)*(1+((topt-t)/(topt-tm)))*(((t-ti)/(topt-ti))^((tm-ti)/(topt-tm)))
    perfmu[perfmu<0] = 0.00001
    return(perfmu)
}

创建负对数似然函数

LL1 = function(tm, topt, pmax, performance=perfdat$performance) {
    perfmu = tpc(tm=tm, topt=topt, pmax=pmax)
    loglike = -sum(dnorm(x=performance, mean=perfmu, log=TRUE))
    return(loglike)
}

使用mle2 - maximum likelihood estimator 的模型性能

m1<- mle2(minuslogl = LL1, start=list(tm=15, topt=20, pmax=150), control=list(maxit=5000))

summary(m1)

这给了我参数估计值,但不是warning message: In sqrt(diag(object@vcov)) : NaNs produced 的误差估计值(标准错误)。但是,参数估计很好,可以让我做出有意义的预测。

Plot of non-linear curve using parameter estimates

我尝试了许多不同的优化器和方法,并得到了关于无法计算标准的相同错误。错误,通常带有关于无法反转粗麻布的警告。或者我对我的参数的估计非常不稳定,而且没有任何意义。

如果我使用:

confint(m1) 

我的每个参数都有 95% 的间隔,但我无法将它们合并到可以用来制作如下图的预测方法中,我使用 nls 模型和 predict() 制作了该图:

non-linear model with error graphed

如果我通过将模型公式嵌入mle2 函数来重新创建我的mle2() 模型:

tpcfun<- function(t, tm.f, topt.f, pmax.f) {
    perfmu = pi+(pmax.f-pi)*(1+((topt.f-t)/(topt.f-tm)))*(((t-ti)/(topt.f-ti))^((tm.f-ti)/(topt.f-tm.f)))
    perfmu[perfmu<0] = 0.00001
    return(perfmu)
}

m2<- mle2(performance ~ dnorm(mean=-sum(log(tpcfun(t=ac.temp, tm.f, topt.f, pmax.f))), sd=1),  method='L-BFGS-B', lower=c(tm.f=1, topt.f=1, pmax.f=1), control=list(maxit=5000, parscale=c(tm.f=1e1, topt.f=1e1, pmax.f=1e2)), start=list(tm.f=15, topt.f=20, pmax.f=150), data=perfdat)

summary(m2)

我的参数得到了不合理的估计值,但我仍然没有得到错误估计值。

我的问题是,任何人都可以看出我的模型(模型函数和似然函数)有什么问题,或者我做错了什么吗?我有一种感觉,我可能写错了似然函数,但我尝试了各种分布和不同的方式,但我可能完全搞砸了。

或者有没有一种方法可以让我对我的参数的误差进行估计,以便我可以使用它们在图表中可视化我的模型预测的误差?

谢谢,

赖利

PS。我决定制作一个只有点估计值的图表,然后是趋势线没有错误,但我想在每个点估计值上放置 95% CI 的条形图,但是 confint() 给了我非常小的 CI '甚至没有出现在图表上,因为它们比我使用的点字符小,哈。

【问题讨论】:

    标签: r predict mle hessian non-linear


    【解决方案1】:

    我认为问题出在“maxint”参数中。尝试使用良好的起始值并避免高交互。第二个问题是算法“L-BFGS-B”,除非默认。当我们使用 confint 函数时,通常会获得 mle2 优化是否收敛的区间。检查配置文件是否可以绘制为图(plotprofmle 函数)。它更安全。 如果您的数据在应用日志时包含零值,则通常会出现“NaN 生成”错误。我建议使用这个顺序:

    loglike = -sum(dnorm(x=performance, mean=perfmu, log=TRUE), na.rm=TRUE)

    验证结果是否合理。

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2012-07-08
      • 1970-01-01
      • 2021-10-21
      • 1970-01-01
      • 2017-04-27
      • 2013-03-29
      • 2020-06-16
      • 2013-01-23
      相关资源
      最近更新 更多