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?

Answer 1

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