R: error with optim and nlme

Question

library(nlme)
mydat <- 
  structure(list(class = c(91L, 91L, 91L, 91L, 91L, 91L, 92L, 92L, 
                           92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 
                           92L, 93L, 93L, 93L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 
                           94L, 94L, 94L, 94L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 
                           95L, 95L, 95L), days = c(5.7, 11.1, 13.9, 15.3, 18.3, 18.9, 1.9, 
                                                    2.1, 2.9, 3.4, 4.4, 5, 6.9, 10.4, 11.6, 13, 13.4, 15.7, 15.9, 
                                                    17.3, 17.7, 19.4, 2.3, 12.6, 15.4, 2, 2.4, 4.9, 5.4, 5.7, 7, 
                                                    7.1, 7.7, 9.1, 9.6, 12.6, 16.6, 17.1, 1, 2, 4.4, 5.6, 5.7, 10.4, 
                                                    12.1, 12.9, 13, 15.6, 16.1, 18.6), outcome = c(3.31586988521325, 
                                                                                                   3.26226473964254, 3.26098046236747, 3.26086828734987, 3.26081582971007, 
                                                                                                   3.2608136610639, 1.54175217273846, 1.74336277564818, 2.48010804039677, 
                                                                                                   2.73455940271066, 2.86602619542132, 2.85753365781511, 2.81667209739959, 
                                                                                                   2.80430238913247, 2.80395479988755, 2.80383291979961, 2.8038189189449, 
                                                                                                   2.80379174103878, 2.80379113213262, 2.80378896928776, 2.80378871890839, 
                                                                                                   2.80378827755366, 1.96537220180574, 3.00124636046136, 3.00096700482166, 
                                                                                                   2.05608815148142, 2.44248026102198, 4.03918455971327, 4.08570704450138, 
                                                                                                   4.09781416453829, 4.09869791687544, 4.09759058744364, 4.09084045815843, 
                                                                                                   4.07921200433542, 4.07668896637335, 4.07081047825795, 4.06993724272757, 
                                                                                                   4.06991925387706, 1.07225026715462, 3.72090308875724, 4.93988448353623, 
                                                                                                   4.7209971449984, 4.70681751285687, 4.49435510591282, 4.48824648431355, 
                                                                                                   4.4870870783591, 4.48698191392076, 4.48574152736339, 4.48566703246688, 
                                                                                                   4.48551396890595), s = c(0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                                                                                                                            0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                                                                                                                            0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                                                                                                                            0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                                                                                                                            0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                                                                                                                            0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                                                                                                                            0.0693, 0.0693, 0.0693, 0.0693, 0.0693)), .Names = c("class", 
                                                                                                                                                                                 "days", "outcome", "s"), row.names = 1001:1050, class = "data.frame")

library(nlme)
obj_NM <- function(arg){
  model = nlme(outcome ~ exp(beta1) * s * days / (s * days + exp(1 - beta3 * s * days)) +
                 exp(beta2) * exp(1 - beta4 * s * days) / (s * days + exp(1 - beta4 * s * days)),
               data = mydat,
               fixed = list(beta1 ~ 1, beta2 ~ 1, beta3 ~ 1, beta4 ~ 1),
               random = list(class = pdDiag(list(beta1 ~ 1, beta2 ~ 1, beta4 ~ 1))),
               start = list(fixed = c(beta1 = arg[1], beta2 = arg[2], beta3 = arg[3], beta4 = arg[4])), verbose = FALSE)
  return(model$logLik)
}

control = list(fnscale = -1)
optim(par = c(1.310239, -4.668217, 17.01345, 2.402943), fn = obj_NM, hessian = TRUE, control = control)

Running the above code gives me the error and warning:

Error in chol.default((value + t(value))/2) : 
  the leading minor of order 2 is not positive definite In addition: Warning messages:
1: In nlme.formula(outcome ~ exp(beta1) * s * days/(s * days + exp(1 -  :
  Singular precision matrix in level -1, block 1
2: In nlme.formula(outcome ~ exp(beta1) * s * days/(s * days + exp(1 -  :
  Singular precision matrix in level -1, block 1

My goal here is to obtain the hessian matrix and investigate why my nlme model may not be fitting. I am trying to maximize my objective function, therefore I set fnscale = -1 (documentation says that it should be negative in order for optim to perform maximization). However, I am not sure what to make of the error message. Is there a way for optim to output the hessian matrix? It seems that an error from nlme has stopped it from doing so.

You are trying to optimize starting values for an optimizer? I don't get it. — Roland, Aug 17 '17 at 09:40
I agree with @Roland - this is highly unconventional, and does not seem to make sense. Before trying to fix your code, we need to know if what you are trying to do is really what you want. Or is this an [XY problem](https://meta.stackexchange.com/questions/66377/what-is-the-xy-problem). Could you concisely explain in the question why you think optimising starting values is required. — dww, Aug 18 '17 at 14:49
... If what you are *actually* trying to achieve is to find the value of the Hessian at the point arrived at by nlme, then wrapping inside optim is completely the wrong approach. Instead you should either 1) just use a numerical Hessian method such as e.g. `pracma::hessian` (with parameter values from nlmeObject providing the input values to the hessian calculation. Or 2) since your formula looks pretty simple, just solve the hessian matrix analytically by symbolic differentiation. — dww, Aug 18 '17 at 15:27

score 0 · Answer 1 · answered Oct 12 '22 at 01:30

I have been able to optimize this function with an initial "pre-optimization step" with DEoptim :

mydat <- structure(list(class = c(91L, 91L, 91L, 91L, 91L, 91L, 92L, 92L, 
                                  92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 
                                  92L, 93L, 93L, 93L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 
                                  94L, 94L, 94L, 94L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 
                                  95L, 95L, 95L),
                        
                        days = c(5.7, 11.1, 13.9, 15.3, 18.3, 18.9, 1.9, 
                                 2.1, 2.9, 3.4, 4.4, 5, 6.9, 10.4, 11.6, 13, 13.4, 15.7, 15.9, 
                                 17.3, 17.7, 19.4, 2.3, 12.6, 15.4, 2, 2.4, 4.9, 5.4, 5.7, 7, 
                                 7.1, 7.7, 9.1, 9.6, 12.6, 16.6, 17.1, 1, 2, 4.4, 5.6, 5.7, 10.4, 
                                 12.1, 12.9, 13, 15.6, 16.1, 18.6), 
                        
                        outcome = c(3.31586988521325, 3.26226473964254, 3.26098046236747, 3.26086828734987, 3.26081582971007, 
                                    3.2608136610639, 1.54175217273846, 1.74336277564818, 2.48010804039677, 
                                    2.73455940271066, 2.86602619542132, 2.85753365781511, 2.81667209739959, 
                                    2.80430238913247, 2.80395479988755, 2.80383291979961, 2.8038189189449, 
                                    2.80379174103878, 2.80379113213262, 2.80378896928776, 2.80378871890839, 
                                    2.80378827755366, 1.96537220180574, 3.00124636046136, 3.00096700482166, 
                                    2.05608815148142, 2.44248026102198, 4.03918455971327, 4.08570704450138, 
                                    4.09781416453829, 4.09869791687544, 4.09759058744364, 4.09084045815843, 
                                    4.07921200433542, 4.07668896637335, 4.07081047825795, 4.06993724272757, 
                                    4.06991925387706, 1.07225026715462, 3.72090308875724, 4.93988448353623, 
                                    4.7209971449984, 4.70681751285687, 4.49435510591282, 4.48824648431355, 
                                    4.4870870783591, 4.48698191392076, 4.48574152736339, 4.48566703246688, 
                                                                                         4.48551396890595), 
                        s = c(0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                              0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                              0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                              0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                              0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                              0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 
                              0.0693, 0.0693, 0.0693, 0.0693, 0.0693)), 
                   .Names = c("class", "days", "outcome", "s"),
                   row.names = 1001 : 1050,
                   class = "data.frame")

library(nlme)
library(DEoptim)

obj_NM <- function(arg, bool_Print = FALSE)
{
  logLik <- tryCatch(nlme(outcome ~ exp(beta1) * s * days / (s * days + exp(1 - beta3 * s * days)) +
                         exp(beta2) * exp(1 - beta4 * s * days) / (s * days + exp(1 - beta4 * s * days)), data = mydat,
                         fixed = list(beta1 ~ 1, beta2 ~ 1, beta3 ~ 1, beta4 ~ 1),
                         random = list(class = pdDiag(list(beta1 ~ 1, beta2 ~ 1, beta4 ~ 1))),
                         start = list(fixed = c(beta1 = arg[1], beta2 = arg[2], beta3 = arg[3], beta4 = arg[4])), verbose = FALSE)$logLik,
                    error = function(e) NA)
  
  if(bool_Print == TRUE)
  {
    print(logLik)
  }
  
  if(is.na(logLik))
  {
    return(10 ^ 30)
    
  }else
  {
    return(-logLik) 
  }
}

obj_Iter <- DEoptim(fn = obj_NM, lower = rep(-1000, 4), upper = rep(1000, 4))

optim(par = obj_Iter$optim$bestmem, fn = obj_NM, hessian = TRUE, method = "BFGS")

$par
par1      par2      par3      par4 
184.59178 178.77184 179.57188  90.01865 

$value
[1] 2.010262

$counts
function gradient 
68        1 

$convergence
[1] 0

$message
NULL

R: error with optim and nlme

1 Answers1