param (eta) sometimes do not change with L-BFGS-B

[FelicienLL]

####
library(tidyverse)
library(mapbayr)
library(mrgsolve)
#> 
#> Attachement du package : 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
packageVersion("mapbayr")
#> [1] '0.2.0.9003'
setwd("C:/Documents and Settings/le_louedec/Documents/REPREX/002/")
#### Dealing with the parameters "stuck" with L-BFGS-B

model <- mread("mrg_001.cpp")
#> Building mrg_001_cpp ...
#> done.
mrgsolve::see(model)
#> 
#> Model file:  mrg_001.cpp 
#> $PROB Reference model with IIV on 7 parameters
#> 
#> $PARAM @annotated
#> TVCL   :  0.2 : Clearance (L/h)
#> TVKA   :  1.0 : Absorption rate (h-1)
#> TVVC   :  5.0 : Central volume of distribution (L)
#> TVVP   : 10.0 : Peripheral volume of distribution (L)
#> TVF    :  0.8 : Bioavailability ()
#> TVQ    :  2.0 : Intercompartimental Clearance (L/h)
#> TVLAG  :  1.5 : Lagtime (h) 
#> 
#> ETA1 : 0 : Random effect on CL
#> ETA2 : 0 : Random effect on VC
#> ETA3 : 0 : Random effect on KA
#> ETA4 : 0 : Random effect on VP
#> 
#> $OMEGA
#> 0.4 // CL
#> 0.4 // VC
#> 0.4 // KA 
#> 0.4 // VP
#> 
#> $SIGMA 
#> 0.05 // err prop
#> 0   //  err additive 
#> 
#> 
#> $CMT @annotated
#> DEPOT   : Depot ()       [ADM]
#> CENTRAL : Central (mg/L) [OBS]
#> PERIPH  : Peripheral ()  []
#> 
#> $TABLE
#> double DV = (CENTRAL / VC) * (1 + EPS(1)) + EPS(2) ;
#> 
#> $MAIN
#> double CL  = TVCL  * exp(ETA(1) + ETA1 ) ; 
#> double VC  = TVVC  * exp(ETA(2) + ETA2 ) ;
#> double KA  = TVKA  * exp(ETA(3) + ETA3 ) ; 
#> double VP  = TVVP  * exp(ETA(4) + ETA4 ) ;
#> double F   = TVF   ;
#> double Q   = TVQ   ;
#> double LAG = TVLAG ;
#> 
#> double K20 = CL / VC ;
#> double K23 = Q  / VC ;
#> double K32 = Q  / VP ;
#> 
#> F_DEPOT = F ; 
#> ALAG_DEPOT = LAG ;
#> 
#> $ODE
#> dxdt_DEPOT   = - KA * DEPOT ;
#> dxdt_CENTRAL = - (K20 + K23) * CENTRAL + K32 * PERIPH + KA * DEPOT ;
#> dxdt_PERIPH  = K23 * CENTRAL - K32 * PERIPH ;
#> 
#> 
#> $CAPTURE @annotated
#> DV : Concentration central
#> ETA1 : 
#> ETA2 :
#> ETA3 : 
#> ETA4 :
data739 <- model %>% 
  adm_lines(amt = 1000, addl = 20, ii = 24, realize_addl = T) %>% 
  obs_lines(time = 40, DV = 76) %>% 
  obs_lines(time = 381, DV = 151) %>% 
  obs_lines(time = 528, DV = 94) %>% 
  see_data()
data739
#> # A tibble: 24 x 10
#>       ID  time  evid  addl    ii   amt   mdv   cmt  rate    DV
#>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
#>  1     1     0     1     0     0  1000     1     1     0    NA
#>  2     1    24     1     0     0  1000     1     1     0    NA
#>  3     1    40     0     0     0     0     0     2     0    76
#>  4     1    48     1     0     0  1000     1     1     0    NA
#>  5     1    72     1     0     0  1000     1     1     0    NA
#>  6     1    96     1     0     0  1000     1     1     0    NA
#>  7     1   120     1     0     0  1000     1     1     0    NA
#>  8     1   144     1     0     0  1000     1     1     0    NA
#>  9     1   168     1     0     0  1000     1     1     0    NA
#> 10     1   192     1     0     0  1000     1     1     0    NA
#> # ... with 14 more rows
#NONMEM founds (with DV values not rounded): 
#ETA1 = 0.0487
#ETA2 = -0.0317
#ETA3 = -0.000690
#ETA4 = -0.0743

#MBR with NEWUOA : 
model %>% 
  mbrest(data739)
#> 
#> ID 1... done.
#> Mrg 001_cpp 
#> ID : 1  individual(s).
#> OBS: 3  observation(s).
#> ETA: 4  parameter(s) to estimate.
#> 
#> Estimates: 
#>   ID       ETA1        ETA2          ETA3        ETA4
#> 1  1 0.04774575 -0.03142508 -0.0006947108 -0.07382542
#> 
#> Output (24 lines): 
#>   ID time evid addl ii  amt mdv cmt rate DV    IPRED     PRED       ETA1
#> 1  1    0    1    0  0 1000   1   1    0 NA  0.00000  0.00000 0.04774575
#> 2  1   24    1    0  0 1000   1   1    0 NA 38.07479 37.06504 0.04774575
#> 3  1   40    0    0  0    0   0   2    0 76 73.01766 71.30260 0.04774575
#> 4  1   48    1    0  0 1000   1   1    0 NA 65.16316 64.35634 0.04774575
#> 5  1   72    1    0  0 1000   1   1    0 NA 84.43529 84.45128 0.04774575
#> 6  1   96    1    0  0 1000   1   1    0 NA 98.14652 99.24745 0.04774575
#>          ETA2          ETA3        ETA4
#> 1 -0.03142508 -0.0006947108 -0.07382542
#> 2 -0.03142508 -0.0006947108 -0.07382542
#> 3 -0.03142508 -0.0006947108 -0.07382542
#> 4 -0.03142508 -0.0006947108 -0.07382542
#> 5 -0.03142508 -0.0006947108 -0.07382542
#> 6 -0.03142508 -0.0006947108 -0.07382542
#good

#MBR with L-BFGS-B : 
model %>% 
  mbrest(data739, method = "L-BFGS-B")
#> 
#> ID 1... done.
#> Mrg 001_cpp 
#> ID : 1  individual(s).
#> OBS: 3  observation(s).
#> ETA: 4  parameter(s) to estimate.
#> 
#> Estimates: 
#>   ID ETA1 ETA2 ETA3 ETA4
#> 1  1    0    0    0    0
#> 
#> Output (24 lines): 
#>   ID time evid addl ii  amt mdv cmt rate DV    IPRED     PRED ETA1 ETA2 ETA3
#> 1  1    0    1    0  0 1000   1   1    0 NA  0.00000  0.00000    0    0    0
#> 2  1   24    1    0  0 1000   1   1    0 NA 37.06504 37.06504    0    0    0
#> 3  1   40    0    0  0    0   0   2    0 76 71.30260 71.30260    0    0    0
#> 4  1   48    1    0  0 1000   1   1    0 NA 64.35634 64.35634    0    0    0
#> 5  1   72    1    0  0 1000   1   1    0 NA 84.45128 84.45128    0    0    0
#> 6  1   96    1    0  0 1000   1   1    0 NA 99.24745 99.24745    0    0    0
#>   ETA4
#> 1    0
#> 2    0
#> 3    0
#> 4    0
#> 5    0
#> 6    0
#wrong. stick to 0.

# Perform sensitivity analysis
sensi739 <- map(1:50, function(x){
  cat(x, " ")
  ini <- set_names(runif(4,-1,1), paste0("ETA",1:4))
  mbrest(x = model, data = data739, method = "L-BFGS-B", force_initial_eta = ini, verbose = F)
})
#> 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50

list(
  T0 = sensi739 %>% map("arg.optim") %>% map('par') %>% bind_rows() %>% mutate(ID = row_number()),
  T1 = sensi739 %>% map("final_eta") %>% map("1") %>% bind_rows() %>% mutate(ID = row_number()) 
) %>% 
  bind_rows(.id = "TYPE") %>% 
  pivot_longer(-c(ID, TYPE)) %>% 
  ggplot(aes(TYPE, value, group = ID))+
  facet_wrap("name",ncol = 4)+
  geom_line()

tab739 <- list(
  T0 = sensi739 %>% map("arg.optim") %>% map('par') %>% bind_rows() %>% mutate(ID = row_number()),
  T1 = sensi739 %>% map("final_eta") %>% map("1") %>% bind_rows() %>% mutate(ID = row_number()) 
) %>% 
  bind_rows(.id = "TYPE") %>% 
  pivot_longer(starts_with("ETA")) %>% 
  pivot_wider(names_from = TYPE) %>% 
  mutate(OK = T0==T1) %>% 
  select(-T1) %>% 
  pivot_wider(names_from = name, values_from = T0)

list(.x = c("ETA1", "ETA1", "ETA1", "ETA2", "ETA2", "ETA3"), 
     .y = c("ETA2", "ETA3", "ETA4", "ETA3", "ETA4", "ETA4")) %>% 
  pmap(function(.x, .y){
    tab739 %>% 
      ggplot(aes(.data[[.x]], .data[[.y]]))+
      geom_point(aes(col = OK))
  })  %>% 
  ggpubr::ggarrange(plotlist = ., common.legend = T)

^{Created on 2020-12-22 by the reprex package (v0.3.0)}

[FelicienLL]

To sum it up, it seems like L-BFGS-B sometimes thinks the optimal parameter values are the initial values. Convergence achieved without changing the parameter values, nor the OFV. No error thrown in the optimx output (convcode = 0).

When we run sensitivity analysis - I mean changing the initial values (instead of 0, the default with L-BFGS-B) - the problem is solved in most of cases. Except for a "space of parameter", for instance negative values of ETA1&3 +positive value for ETA2&4. I am pretty sure that a 4-dimension representation would show that the blue "TRUE" point are all all stuck together. It means that initial values must not be taken from this zone (Bermuda triangle...).

Possible solution : decrease the tolerance of the optimizer ? I am not sure of what I am doing here, or if my vocabulary is correct. I tried to tweak the L-BFGS-B arguments (ptol, lmm etc...) very empirically, but nothing happens (not shown).

Conclusion: For now, I think I will just add a "while" loop that re-run the optimization if there is no change in OFV. New initial value sampled from a rather "large interval" (i.e. runif(n, -1,1)) to be sure to get out of the "Bermuda triangle".