我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)
}
模型性能使用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
。然而,参数估计很好,让我做出有意义的预测。
我尝试了许多不同的优化器和方法,并得到了关于无法计算标准的相同错误。错误,通常带有关于无法反转粗麻布的警告。或者我对我的参数的估计非常不稳定,这没有意义。
如果我使用:
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)
我对我的参数进行了无意义的估计,但我仍然没有得到错误估计。
我的问题是,任何人都可以看到我的模型(模型函数和似然函数)有什么问题,或者我做错了什么吗?我有一种感觉,我可能写错了似然函数,但我尝试了各种分布和不同的方式,但我可能完全搞砸了。
或者有没有一种方法可以让我对我的参数的误差进行估计,以便我可以使用它们在图表中可视化我的模型预测的误差?
谢谢,
赖利
PS。我决定制作一个仅包含点估计值的图表,然后是没有错误的趋势线,但我想在每个点估计值上放置 95% CI 的条形图,但是 confint() 给了我非常小的 CI '甚至没有出现在图表上,因为它们比我使用的点字符小,哈。