如何从nlme调用获取Hessian

Question

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 函数（也是来自同一个包），"Evaluate an approximate Hessian and gradient of a scalar function using finite differences" 这是否等同于当前在函数 nlme 中实现的内容？

Answer 1

我想在这里要做的第一件事是看一下控件：

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 也允许 return 粗麻布，如果我打印我得到：

$hessian
             [,1]
[1,] 8.478483e-05

如果您想知道我是如何打印 hessian 的，就我而言，我正在编辑 nlme.formula 并将我的新版本分配给一个名为 nlme.formula_new 的函数，然后我将其重新插入 [=19] =]:

godmode:::assignAnywhere("nlme.formula", nlme.formula_new)

godmode 在 Github 上，但肯定还有其他方法可以实现。

Answer 2

我认为你的问题是起点选择不当造成的。向量 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

归根结底，模型拟合可以理解为一个优化问题。当您的模型是非线性的（例如混合效应模型）时，必须使用迭代优化算法来解决该问题。因此，起始值的选择可能非常关键。这是一篇讨论该主题的不错的科学文章：

Balsa-Canto, E., Alonso, A.A. & Banga, J.R. An iterative identification procedure for dynamic modeling of biochemical networks. BMC Syst Biol 4, 11 (2010) doi:10.1186/1752-0509-4-11

如何从nlme调用获取Hessian

How to obtain Hessian from nlme call

r

nlme