Integrate external optimisation results into msm object

[marc-vaisband]

I have a model specified in msm, where the built-in optimization terminates, but fails, using all available methods, to converge to a satisfactory result. By accessing the likelihood directly via

msm:::lik.msm(x, msmdata, msm_res$qmodel, msm_res$qcmodel, msm_res$cmodel, msm_res$hmodel, msm_res$paramdata),

I am doing the optimization using a different package, which yields a likelihood that is 10000 points better than what BFGS is able to obtain (even after extensive fiddling with fnscale etc.). How can I now use that external result with the msm object? I called

msm::updatepars.msm(msm_res, better_candidate),

but that seems to only set the estimate field, without updating anything else (not even the minus2loglik attribute). I would like to, for instance, check the Hessian at the new estimate, and obtain standard errors etc.

[chjackson]

I've attempted to write some documentation here for how to add a new optimisation function to msm. If any of these steps go wrong for you, I'll refine these instructions. It is easier than you might think, because you don't need to download or edit the msm source code.

Suppose your external optimisation function is called NEWOPTFUNCTION. If you call msm(..., opt.method="NEWOPTFUNCTION"), then msm will look for a function called msm.optim.NEWOPTFUNCTION in your search path which connects this optimiser to msm. Then everything is handled automatically from there.

So you just need to write this function, based on the following template:

msm.optim.NEWOPTFUNCTION <- function(p, gr, hessian, msmdata, qmodel, qcmodel, cmodel, hmodel, ...) {
  optim.args <- list(...)
  optim.args <- c(optim.args,
                  list(# Supply the required arguments to NEWOPTFUNCTION here,
                       # The function to be optimised is lik.msm,
                       # The corresponding gradient function is gr (can ignore if the optimiser doesn't use gradients)
                       # The initial values are in p$inits
                       msmdata=msmdata, qmodel=qmodel, qcmodel=qcmodel,
                       cmodel=cmodel, hmodel=hmodel, paramdata=p))
  p$opt <- do.call(NEWOPTFUNCTION, optim.args)

  ## Fill in the following with the appropriate outputs of NEWOPTFUNCTION:
  # minimised -2 log likelihood
  p$lik <-                    

  # parameter estimates, including only those being optimised, ignoring those fixed at their initial values
  p$params[p$optpars] <- 

  p$opt$par <- p$params[p$optpars]
  p$opt$hessian <-        # the Hessian.  
  #  if NEWOPTFUNCTION does not return the Hessian, could calculate it numerically with: 
  #  p$opt$hessian <- stats::optimHess(par=p$opt$par, fn=lik.msm, gr=gr,
                                 msmdata=msmdata, qmodel=qmodel, qcmodel=qcmodel,
                                 cmodel=cmodel, hmodel=hmodel, paramdata=p)

  p
}

To see examples of how this function is written in practice, look at msm.optim.bobyqa (which uses minqa::bobyqa as the optimiser) and msm.optim.nlm (which uses stats::nlm) in optim.R.

You could also include error handling, e.g. to catch any error code returned by NEWOPTFUNCTION if it doesn't converge, and give an informative message.

Out of interest, what is the optimiser you are using? It may be worth making your msm.optim.NEWOPTFUNCTION available more widely if it is more reliable. On the other hand, if you are having to fiddle with optimisers to get results, then consider whether a simpler model will suit your purpose - beware of trying to squeeze information out of your data that isn't there.

[marc-vaisband]

Thank you for your reply! The optimisation routine I'm using is actually in Julia - an adaptive differential evolution algorithm provided from the BlackBoxOptim.jl package -, where I'm executing all R code in an R environment provided via the RCall.jl package. (The reverse would presumably also be possible by grabbing this optimisation package through JuliaCall). For the moment, I'll try to just write a msm.optim.dummy routine which returns the hard-coded best parameters / objective function value.

For using stats::optimHess for Hessian estimation, what would the call to msm:::lik.msm look like? Can the gradient as passed via gr be used immediately inside the optimHess call?

[chjackson]

Thanks - I've updated the comment above for optimHess. See also the examples I referred to in optim.R in the msm source.

I am not experienced with Julia, but I've had "proof of concept" success with the JuliaConnectoR package for calling a Julia routine from R while touching Julia as little as possible.