How to obtain Hessian from nlme call

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

I am trying to investigate why some nlme models do not fit successfully by looking at the hessian. Is there a way to extract this matrix somehow?

I am also looking into the fdHess function (also from the same pacakge), which "Evaluate an approximate Hessian and gradient of a scalar function using finite differences" would this be equivalent to what is currently implemented in the function nlme?

Have you seen varFix from here https://stat.ethz.ch/R-manual/R-devel/library/nlme/html/nlmeObject.html ? A method that I used in the past to get around numerical problems, is to use different starting values, here using argument start. — papgeo, Jul 30 '18 at 07:13

score 1 · Accepted Answer · answered Dec 03 '19 at 10:50

I believe your problem is caused by a poor selection of starting points. The vector c(Asym = 103, R0 = -8.5, lrc = -3.3) converges without any complication:

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

At the end of the day, model fitting can be understood as an optimization problem. When your model is non-linear (such as a mixed-effects model), that problem has to be solved using iterative optimization algorithm. Hence, the selection of starting value can be quite critical. Here is a nice scientific article discussing this topic:

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

score 0 · Answer 2 · answered Jan 29 '18 at 21:50

I think the first thing to do here is to have a look at the controls:

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

I haven't read any documentation, but you could at least run the following to get a better clue of the situation:

# 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
)

This would print you output for every iteration, finally:

# ............
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

Maybe you want to modify msTol or other controls as well. Note that nlm allows to return the hessian too, if I print that I get:

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

If you wonder how I got the hessian printed, for my part I am editing nlme.formula and assigning my new version to a function called nlme.formula_new, which I then plug back into nlme:

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

godmode is on Github, but for sure there are other ways to achieve this.

How to obtain Hessian from nlme call

2 Answers2