r-使用mle2进行有误差的参数估计和预测



我正在使用mle2来估计非线性模型的参数,我希望估计参数估计的误差(std.error(。此外,我想使用该模型来预测新数据,在这个过程中的几个步骤我遇到了问题(错误(。

数据如下:

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的误差(std.error(进行估计。然而,参数估计是好的,让我的预测是有意义的。

使用参数估计绘制非线性曲线

我尝试了许多不同的优化器和方法,但都得到了相同的错误,即无法计算std.error,通常会出现无法反转hessian的警告。或者,我对自己的参数做出了非常不稳定的估计,这些估计毫无意义。

如果我使用:

confint(m1) 

我的每个参数都有95%的区间,但我无法将这些区间纳入一种预测方法中,我可以使用这种方法来制作下面这样的图,我使用nls模型和predict():

具有误差图的非线性模型

如果我通过将模型公式嵌入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)

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

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

或者,有没有一种方法可以让我得到参数周围的误差估计,这样我就可以用它们在图中可视化我的模型预测周围的误差?

谢谢,

Rylee

PS。我决定画一张只有点估计的图,然后画一条没有误差的趋势线,但我想在每个点估计上加上95%置信区间的条形图,但confint((给了我一个非常小的置信区间,它们甚至没有出现在图上,因为它们比我使用的点字符小,哈。

我认为问题出在"maxint"arg中。尽量使用良好的起始值,避免高迭代。第二个问题是算法"L-BFGS-B",除非默认。当我们使用confint函数时,获得mle2优化是否收敛的区间是正常的。检查纵断面是否可以作为绘图(plotprofmle函数(。它更安全。如果您的数据在应用日志时包含零值,则通常会出现"NaNs生成"错误。我建议使用这个序列:

loglike=-sum(dnorm(x=性能,平均值=性能,log=真(,na.rm=真(

验证结果是否合理。

最新更新