如何从nlme调用中获取Hessian
如何从nlme调用中获取Hessian
提问于 2017-08-09 06:15:53
library(nlme)
fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1,
start = c(Asym = -10311111, R0 = 8.5^4, lrc = 0.01),
verbose = TRUE)
**Iteration 1
LME step: Loglik: -312.2787, nlminb iterations: 23
reStruct parameters:
Seed
10.41021
Error in nlme.formula(height ~ SSasymp(age, Asym, R0, lrc), data = Loblolly, :
Singularity in backsolve at level 0, block 1我正试图研究为什么一些nlme模型不能成功地通过观察恒河。有办法提取这个矩阵吗？
我还研究了fdHess函数(也来自同一个空间)，该函数“使用有限差分计算标量函数的近似Hessian和梯度”，这是否等同于目前在函数nlme中实现的方法？
回答 2
Stack Overflow用户
回答已采纳
发布于 2019-12-03 10:50:09
我相信你的问题是由于起点选择不当造成的。向量c(Asym = 103, R0 = -8.5, lrc = -3.3)收敛，没有任何复杂情况：
nlme(height ~ SSasymp(age, Asym, R0, lrc),
data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1,
start = c(Asym = 103, R0 = -8.5, lrc = -3.3))
#> Nonlinear mixed-effects model fit by maximum likelihood
#> Model: height ~ SSasymp(age, Asym, R0, lrc)
#> Data: Loblolly
#> Log-likelihood: -114.7428
#> Fixed: Asym + R0 + lrc ~ 1
#> Asym R0 lrc
#> 101.449600 -8.627331 -3.233751
#>
#> Random effects:
#> Formula: Asym ~ 1 | Seed
#> Asym Residual
#> StdDev: 3.650642 0.7188625
#>
#> Number of Observations: 84
#> Number of Groups: 14 最终，模型拟合可以理解为一个优化问题。当您的模型是非线性的(例如混合效应模型)时，必须使用迭代优化算法来解决这个问题。因此，起始值的选择可能是相当关键的。这里有一篇很好的科学文章讨论这个主题：
鲍萨-坎托，E.，阿隆索，A.A. & Banga，J.R.生化网络动态建模的迭代辨识过程。BMC 4，11 (2010年) doi:10.1186/1752-0509-4-11
Stack Overflow用户
发布于 2018-01-29 21:50:08
我认为这里要做的第一件事就是看看这些控件：
nlmeControl()
# this gives in my case the following settings
$maxIter
[1] 50
$pnlsMaxIter
[1] 7
$msMaxIter
[1] 50
$minScale
[1] 0.001
$tolerance
[1] 1e-05
$niterEM
[1] 25
$pnlsTol
[1] 0.001
$msTol
[1] 1e-06
$returnObject
[1] FALSE
$msVerbose
[1] FALSE
$gradHess
[1] TRUE
$apVar
[1] TRUE
$.relStep
[1] 6.055454e-06
$minAbsParApVar
[1] 0.05
$opt
[1] "nlminb"
$natural
[1] TRUE
$sigma
[1] 0我还没有读过任何文档，但您至少可以运行以下代码来更好地了解情况：
# using additional controls list argument
fm1 <- nlme(
height ~ SSasymp(age, Asym, R0, lrc),
data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1,
start = c(Asym = -10311111, R0 = 8.5^4, lrc = 0.01),
control = list(opt = "nlm", msVerbose = 2, msTol = 1e-06),
verbose = TRUE
)这将为每次迭代打印输出，最后：
# ............
iteration = 9
Parameter:
[1] 7.180326
Function Value
[1] 379.1821
Gradient:
[1] -4.212256e-05
Relative gradient close to zero.
Current iterate is probably solution.
**Iteration 1
LME step: Loglik: -312.2787, nlm iterations: 9
reStruct parameters:
Seed
7.180326
Error in nlme.formula(height ~ SSasymp(age, Asym, R0, lrc), data = Loblolly, :
Singularity in backsolve at level 0, block 1也许您也想修改msTol或其他控件。请注意，如果我打印的话，nlm也允许返回恒河语：
$hessian
[,1]
[1,] 8.478483e-05如果您想知道我是如何打印的，就我而言，我正在编辑nlme.formula，并将我的新版本分配给一个名为nlme.formula_new的函数，然后将其插入nlme。
godmode:::assignAnywhere("nlme.formula", nlme.formula_new)godmode是在Github上的，但是肯定还有其他的方法来实现这一点。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/45582905
复制相关文章
相似问题
腾讯云开发者
