我正在使用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)
}
mle
2-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=真(
验证结果是否合理。