mceta -- Monte-Carlo focei?

nlmixr2 4.0 mceta

Two of popular algorithms for fitting nonlinear fixed effects models in pharmacometrics are first order conditional estimation (with interaction), sometimes called focei. This is often considered a classical estimation method

As time progressed, Monte-Carlo methods started being used, such as stochastic approximation estimation-maximization (saem). This uses Mote-Carlo Methods to try to find the maximum expectation of the nonlinear-mixed effects model.

The new feature I am discussing is mceta for focei, which combines the Monte-Carlo sampling with focei.

What is mceta

One of the most important parts of getting a good model in non-linear mixed effects models is to choose reasonable initial estimates for the model. This is something pharmacometricians do to get the best and most reasonable model for the data. This can take several iterations to get right (depending on the data, model and algorithm).

Usually we set things like the population estimates of the model, and the variances/covariance of between subject variability as well as the residual errors.

One thing you can do in both nlmixr2 and NONMEM is to set the initial between subject variability estimates (which we will call etas) to get an even better initial estimate for the problem.

While true, this is infrequently done by pharmacometricians. Typically they accept the defaults of the algorithm and hope for a good final fit.

However with mceta, we can change the algorithm’s default. In nlmixr2 these options are:

• mceta=-1 – Use the best eta estimate from the last iteration (nlmixr2’s default, and not available in NONMEM at the time of this writing)

• mceta=0 – Start at all etas equal 0 (NONMEM’s default)

• mceta=1 – Try the eta at the last best fit estimate as well as the etas set to zero, and start the optimization at the eta with the lowest objective function.

• mceta=N – Try the etas as described in mceta=1 and then take N-1 random multivariate normal samples (assuming the specified between subject variability covariance matrix, omega, is true). Start the optimization at the eta value with the smallest individual objective function.

At values of mceta > 1, the focei algorithm becomes a Monte-Carlo method (though not a mcmc method).

This means that some of the strengths and robustness of the new algorithms can be achieved by increasing the mceta value; of course, this also means that the estimation method is non-deterministic too.

This extra feature was added to help include focei linearization of models (which I helped Omar Elashkar work on during his summer internship).

In general, this helps with the focei linearization (in both nlmixr2 and in NONMEM) but it can also help find better models with ODEs (at the cost of computation time).

I have not seen much discussion about this feature, though I have heard it mentioned a couple of times by some analysts. I think it can help in complex analyses where you want to use focei.

Simple implementation example

To implement mceta, you simply need to put it in the focei-related methods like focei, foce, laplace, and agq. For simplicity sake, I will add it to the thepohylline example:

library(nlmixr2)

## ── Attaching packages ───────────────────────────────────────────────────────── nlmixr2 4.0.1 ──

## ✔ lotri         1.0.2          ✔ monolix2rx    0.0.6     
## ✔ nlmixr2data   2.0.9          ✔ nlmixr2lib    0.3.0.9000
## ✔ nlmixr2est    4.1.0.9000     ✔ nlmixr2rpt    0.2.1     
## ✔ nlmixr2extra  3.0.3          ✔ nonmem2rx     0.1.8     
## ✔ nlmixr2plot   3.0.3          ✔ posologyr     1.2.8     
## ✔ rxode2        4.1.0.9000     ✔ shinyMixR     0.5.2     
## ✔ babelmixr2    0.1.9          ✔ xpose.nlmixr2 0.4.1     
## ✔ ggPMX         1.3.2

## ── Optional Packages Loaded/Ignored ─────────────────────────────────────────── nlmixr2 4.0.1 ──

## ✔ babelmixr2        ✔ nonmem2rx    
## ✔ ggPMX             ✔ posologyr    
## ✔ monolix2rx        ✔ shinyMixR    
## ✔ nlmixr2lib        ✔ xpose.nlmixr2
## ✔ nlmixr2rpt

## ── Conflicts ───────────────────────────────────────────────────────────── nlmixr2conflicts() ──
## ✖ rxode2::boxCox()     masks nlmixr2est::boxCox()
## ✖ rxode2::yeoJohnson() masks nlmixr2est::yeoJohnson()

# For reproducibility, set the seed for R and rxode2
set.seed(42)
rxSetSeed(42)

one.cmt <- function() {
  ini({
    tka <- 0.45 # Log Ka
    tcl <- log(c(0, 2.7, 100)) # Log Cl
    tv <- 3.45; label("log V")
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.sd <- 0.7
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    linCmt() ~ add(add.sd)
  })
}

fit <- nlmixr(one.cmt, theo_sd, est="focei",
              control=list(print=0, mceta=5))

## ℹ parameter labels from comments are typically ignored in non-interactive mode

## ℹ Need to run with the source intact to parse comments

## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====] 0:00:00 
## done

## → Calculating residuals/tables

## ✔ done

## → compress origData in nlmixr2 object, save 5952

## → compress parHistData in nlmixr2 object, save 5088

print(fit)

## ── nlmixr² FOCEi (outer: nlminb) ──
## 
##          OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
## FOCEi 116.806 373.4058 393.5854      -179.7029        68.43734        9.370183
## 
## ── Time (sec $time): ──
## 
##            setup optimize covariance table compress    other
## elapsed 0.001671 0.358515   0.358516 0.074    0.007 2.543298
## 
## ── Population Parameters ($parFixed or $parFixedDf): ──
## 
##        Parameter  Est.     SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
## tka              0.469  0.195 41.6        1.6 (1.09, 2.34)     70.7      2.04% 
## tcl               1.01 0.0751 7.41       2.75 (2.38, 3.19)     26.7      3.90% 
## tv         log V  3.46 0.0436 1.26       31.8 (29.2, 34.6)     14.1      10.7% 
## add.sd           0.694                               0.694                     
##  
##   Covariance Type ($covMethod): r,s
##   No correlations in between subject variability (BSV) matrix
##   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
##   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
##   Information about run found ($runInfo):
##    • gradient problems with initial estimate and covariance; see $scaleInfo 
##    • last objective function was not at minimum, possible problems in optimization 
##    • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.)) 
##    • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=)) 
##   Censoring ($censInformation): No censoring
##   Minimization message ($message):  
##     relative convergence (4) 
## 
## ── Fit Data (object is a modified tibble): ──
## # A tibble: 132 × 22
##   ID     TIME    DV  PRED    RES   WRES IPRED   IRES  IWRES CPRED   CRES  CWRES
##   <fct> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>
## 1 1      0     0.74  0     0.74   1.07   0     0.74   1.07   0     0.74   1.07 
## 2 1      0.25  2.84  3.28 -0.437 -0.525  3.85 -1.01  -1.45   3.49 -0.654 -0.744
## 3 1      0.57  6.57  5.85  0.721  0.672  6.79 -0.216 -0.311  6.16  0.405  0.344
## # ℹ 129 more rows
## # ℹ 10 more variables: eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>, depot <dbl>,
## #   central <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>, dosenum <dbl>

In this very simplistic case, there was not a huge difference, however with complex models this could make the difference between a good solution and a poor solution.

Happy modeling!

Roulette wheel originally designed by macrovector / Freepik with modifications for this post.