Gradient and hessian of outer optimization

[simonbeyer1]

Hello,

I am using the nlmixr package for model fitting, specifically employing the foceiControl function with outerOpt = "nlminb". After fitting a model using nlmixr, I am wondering if the gradient and Hessian of the outer optimization problem are stored within the final fit object.

Example:

## Define the model

one.compartment <- function() {
  ini({
    tka <- log(2); label("Ka")
    tcl <- log(2.72); label("Cl")
    tv <- log(31.5); label("V")
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    logrel.sd <- log(0.1)
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    rel.sd <- exp(logrel.sd)
    d/dt(depot) <- -ka * depot
    d/dt(center) <- ka * depot - cl / v * center
    cp <- center / v
    cp ~ prop(rel.sd)
  })
}

## Perform model fitting
fit <- nlmixr(one.compartment, theo_sd, est = foceiControl(outerOpt = "nlminb", print = 1))

# final_outer_grad <- ?
# final_outer_hess <- ?

Are the gradient and Hessian of the outer optimization problem stored within the final fit object? (I couldn't find it). Or is there another easy way to get these?

Kind regards,
Simon

[mattfidler]

Hi Simon,

I am unsure exactly what you are looking for, but you can get some of these quantities as follows:

library(nlmixr2)
#> Warning: package 'nlmixr2' was built under R version 4.3.3
#> Loading required package: nlmixr2data
#> Warning: package 'nlmixr2data' was built under R version 4.3.3

one.compartment <- function() {
  ini({
    tka <- log(2); label("Ka")
    tcl <- log(2.72); label("Cl")
    tv <- log(31.5); label("V")
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    logrel.sd <- log(0.1)
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    rel.sd <- exp(logrel.sd)
    d/dt(depot) <- -ka * depot
    d/dt(center) <- ka * depot - cl / v * center
    cp <- center / v
    cp ~ prop(rel.sd)
  })
}

## Perform model fitting
fit <- nlmixr(one.compartment, theo_sd, est = foceiControl(outerOpt = "nlminb", print = 0))
#> ℹ infer estimation `focei` from control
#> rxode2 2.1.3 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> calculating covariance matrix
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 5952
#> → compress parHistData in nlmixr2 object, save 6088

# Hessian, though on the normal scale
fit$R*2
#>                  tka        tcl        tv   logrel.sd
#> tka        55.995808  -0.131367 -30.55977    3.923016
#> tcl        -0.131367 392.526362 -13.50022  -14.814516
#> tv        -30.559768 -13.500225 826.29329  -39.564042
#> logrel.sd   3.923016 -14.814516 -39.56404 1023.689801

# you can get the last gradient calculated
# In this case it is on the scaled problem, (unlike the Hessian above)
tail(fit$parHistData[fit$parHistData$type %in% c("Gill83 Gradient", "Mixed Gradient", "Forward Difference", "Central Difference"),])
#>     iter               type objf        tka        tcl         tv   logrel.sd
#> 177   37 Forward Difference   NA  0.5229459  0.9179864 -0.1828402 -1.19324012
#> 178   38 Forward Difference   NA  0.2299080 -0.2587128 -0.3280831 -0.08695498
#> 179   40 Forward Difference   NA  0.3816128  0.2721467  1.5885789  0.57061120
#> 180   41 Forward Difference   NA -0.1114214  0.2206542 -0.2527001 -0.05671918
#> 181   45 Forward Difference   NA  0.2399665  1.7626491 -0.6369708  0.00676509
#> 182   47 Forward Difference   NA -0.4125014 -0.7064363  0.8090009 -1.29373595
#>              o1            o2          o3
#> 177 -0.07643536 -0.2727113583 -0.59336620
#> 178  0.01603052 -0.2551059798 -0.57315364
#> 179 -0.01400599 -0.1937701768 -0.54774903
#> 180 -0.12720286 -0.2155226160 -0.55133902
#> 181 -0.86937944 -0.4221508811 -0.10350763
#> 182 -0.21256984  0.0003226439 -0.01438607

^{Created on 2024-07-10 with reprex v2.1.0}

This takes advantage of the fact R = Hessian / 2 . This is on the unscaled parameter space

Note that the gradients that you are seeing are the last observed gradient and on the scaled parameter space. nlminb does not always request a gradient at the last observation.

[simonbeyer1]

Hi Matt,

Thank you very much for the quick answer. And sorry for posting in 'issues'!

Thanks for the hessian on the normal scale. I was not aware of fit$R!

I am currently studying non-linear mixed-effects models. So, for clarification:

In my understanding, the objective function that is minimized is $- 2 \log \mathcal{L} (\theta, \Omega, \sigma)$, where

$$ \mathcal L (\theta,\Omega,\sigma) = \prod_i p(y^{\text{data}}_i \vert \theta, \Omega, \sigma) = \prod_i \int p(y^{\text{data}}_i \vert \eta_i , \theta, \sigma) p(\eta_i \vert \Omega) d \eta_i $$

$\theta$ are population parameters, $\eta$ are random effects, $\Omega$ is the covariance matrix, and $\sigma$ are error model parameters. I guess the first order approximation is than applied to $p(y^{\text{data}}_i \vert \eta_i , \theta, \sigma)$. But anyway, in this case, the objective function should have 7 free parameters. The gradient aligns with this (I assume o1, o2, o3 represent the diagonal elements of $\Omega$). The objective function value can be obtained with fit$objf. Is it possible to compute gradient and Hessian (which should be a 7x7 matrix in this example) at the bestfit $(\hat \theta, \hat \Omega, \hat \sigma)$?

[mattfidler]

Is it possible to compute gradient and Hessian (which should be a 7x7 matrix in this example) at the bestfit $(\hat \theta, \hat \Omega, \hat \sigma)$?

Yes it is possible but has not be implemented. The full hessian (and covariance in general) is on the todo list and is a frequent request.

There has been no request to get the final gradient, though it could be used as a metric to see how well the minimum is performing.