nlmeControl documentation issues

[JSC67]

I am interested in seeing the likelihood (or log-likelihood) for a combination of model, known parameter values, and a data set; this would be the value prior to any parameter optimization. I've read the nlmixr documentation fairly carefully but it's not yet clear to me how best to do this. I assume that this should be controlled by an option in nlmeControl or saemControl (but only the latter has an entry in the documentation [let's call this sub-issue 1]). In the Package 'nlmixr' documentation (2017-11-09) under gnlmm, in the Examples section, there is this line:

fit = gnlmm(llik, pump, inits,
control=list(
reltol.outer=1e-4,
optim.outer="nmsimplex",
nAQD=5
)
)

What does nAQD control? [sub-issue 2] It would be nice to have some discussion or reference to the control parameters.

Would something like this be suitable for getting only the initial objective function value or log-likelihood?

fit <- nlmixr(one.compartment.model, data, est="nlme",
              control = nlmeControl(maxIter = 1, pnlsMaxIter = 1, msMaxIter = 1))

Is this overkill or the best way to get this information? [sub-issue 3]

Thank-you,
James

[mattfidler]

My comments:

• nlme is controlled by a different package and its documentation is controlled by them. Therefore nlmeControl cannot be updated by us.
• gnlmm uses a Laplacian approximation to the likelihood function, whereas nlme and saem can output FOCEi objective functions (See very different SAEM results and other tutorial issues #29 for more disucssion). Therefore these will not be comparable.
• While it is possible to have the objective function calculated based on a set of parameter estimates and individual ETAs, I haven't figured out a clean way to do this from the nlmixr function. This will take a bit more thought.

[mattfidler]

I believe that nADQ referes to the number of points in the adapative gaussian quadrature used to help the laplacian minimization (I forgot to mention this).

Currently gnlmm and gnlmm2 have not be expanded to the single function model, but there are plans for this in the future. It will take some time.

[RikSchoemaker]

Just to add on Matt's comment: nlme control options are in the nlme package, but running

?nlmeControl

brings the documentation up perfectly!

[mattfidler]

Also if you are interested,

?saemControl

brings up the documenation of saem as well 😄

[JSC67]

Thank-you both for your comments. Matt, you wrote, "While it is possible to have the objective function calculated based on a set of parameter estimates and individual ETAs, I haven't figured out a clean way to do this from the nlmixr function. This will take a bit more thought." Getting such values is critical to what I want to do with them, so I hope that you'll be able to find time to work on exactly this. Any comments on control = nlmeControl(maxIter = 1, pnlsMaxIter = 1, msMaxIter = 1) ? Thanks!

[mattfidler]

Yes. control=nlmeControl(maxIter=1, pnlsMaxIter=1, nsMaxIter=1) will still estimate the ETAs, and it unclear if these esimates are going to be great since you are interrupting them early.

[JSC67]

Although this issue is marked as closed, to me it's still open. I tried R control=nlmeControl(maxIter=1, pnlsMaxIter=1, nsMaxIter=1) as you wrote would work, but I got an error and no object was created. I'm sorry if I'm overlooking something obvious, but can you please help me in this? It may seem pointless now, but it's actually critical for what I want to do with it. Thanks.

R version 3.4.3 (2017-11-30) -- "Kite-Eating Tree"
Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Preliminaries
> setwd("C:/Users/James/Documents/pharmacometrics/testing")
> library(nlmixr)
Loading required package: nlme
Loading required package: RxODE
> # rxVersion(echo=TRUE)
> rxClean()
[1] TRUE
> rxVersion(echo=TRUE)
  _____         ____  _____  ______
 |  __ \       / __ \|  __ \|  ____| 0.7-0
 | |__) |__  _| |  | | |  | | |__
 |  _  / \ \/ / |  | | |  | |  __|
 | | \ \  >  <| |__| | |__| | |____
 |_|  \_\/_/\_\\____/|_____/|______|
> # Use some pre-loaded data.
> str(theo_sd)  # like Theoph but with ID instead of Subject, also other NONMEM-like
'data.frame':	144 obs. of  7 variables:
 $ ID  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ TIME: num  0 0 0.25 0.57 1.12 2.02 3.82 5.1 7.03 9.05 ...
 $ DV  : num  0 0.74 2.84 6.57 10.5 9.66 8.58 8.36 7.47 6.89 ...
 $ AMT : num  4.02 0 0 0 0 0 0 0 0 0 ...
 $ EVID: int  101 0 0 0 0 0 0 0 0 0 ...
 $ CMT : int  1 2 2 2 2 2 2 2 2 2 ...
 $ WT  : num  79.6 79.6 79.6 79.6 79.6 79.6 79.6 79.6 79.6 79.6 ...
> # 'uif' for Unified Interface Function, 1 compartment, lin for linCmt()
> uif.1cm.lin <- function() {
+   ini({
+     logtcl <- -3.2      # log typical value Cl (L/hr)
+     logtka <-  0.5      # log typical value Ka (1/hr)
+     logtv  <- -1        # log typical value V  (L)
+     
+     # error model
+     add.err <- 0.1
+     
+     # Initial estimates for IIV variances
+     # Labels work for single parameters.
+     eta.cl ~ 2
+     eta.ka ~ 1    ## BSV Ka
+     eta.v  ~ 1
+   })
+   model({
+     cl <- exp(logtcl + eta.cl)
+     ka <- exp(logtka + eta.ka)
+     v  <- exp(logtv  + eta.v)
+     linCmt() ~ add(add.err)
+   })
+ }
> # NLME fittings
> 
> fit.nlme.1cm.lin.theo_sd     <- nlmixr(uif.1cm.lin, theo_sd, est="nlme")

**Iteration 1
LME step: Loglik: -182.2318, nlminb iterations: 1
reStruct  parameters:
      ID1       ID2       ID3 
1.0022617 0.2783009 2.0758984 
PNLS step: RSS =  63.14922 
 fixed effects: -3.211937  0.4479979  -0.7859318  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
   fixed reStruct 
0.272375 3.267309 

**Iteration 2
LME step: Loglik: -179.291, nlminb iterations: 9
reStruct  parameters:
       ID1        ID2        ID3 
0.96385527 0.08333435 1.63552295 
PNLS step: RSS =  63.28221 
 fixed effects: -3.211535  0.4441505  -0.7863391  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
      fixed    reStruct 
0.008662381 0.172135330 

**Iteration 3
LME step: Loglik: -179.3381, nlminb iterations: 8
reStruct  parameters:
       ID1        ID2        ID3 
0.96129928 0.07085834 1.63885629 
PNLS step: RSS =  63.2196 
 fixed effects: -3.211732  0.4458486  -0.7861743  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
      fixed    reStruct 
0.003808629 0.082919683 

**Iteration 4
LME step: Loglik: -179.3201, nlminb iterations: 7
reStruct  parameters:
       ID1        ID2        ID3 
0.96243507 0.07618717 1.63736945 
PNLS step: RSS =  63.22991 
 fixed effects: -3.211823  0.4451494  -0.7862432  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
      fixed    reStruct 
0.001570583 0.013884126 

**Iteration 5
LME step: Loglik: -179.3277, nlminb iterations: 4
reStruct  parameters:
       ID1        ID2        ID3 
0.96197111 0.07396744 1.63796900 
PNLS step: RSS =  63.24112 
 fixed effects: -3.211722  0.4454457  -0.7862145  
 iterations: 6 
Convergence crit. (must all become <= tolerance = 1e-05):
       fixed     reStruct 
0.0006651255 0.0147544693 

**Iteration 6
LME step: Loglik: -179.3246, nlminb iterations: 1
reStruct  parameters:
       ID1        ID2        ID3 
0.96215433 0.07487426 1.63772430 
PNLS step: RSS =  63.24112 
 fixed effects: -3.211722  0.4454457  -0.7862145  
 iterations: 1 
Convergence crit. (must all become <= tolerance = 1e-05):
       fixed     reStruct 
0.000000e+00 3.910435e-10 
Using sympy via SnakeCharmR
## Calculate ETA-based prediction and error derivatives:
Calculate Jacobian...................done.
Calculate sensitivities.......
done.
## Calculate d(f)/d(eta) 
## ...
## done
## ...
## done
C:/RTools/3.4/mingw_64/bin/gcc  -I"C:/R/R-34~1.3/include" -DNDEBUG          -O2 -Wall  -std=gnu99 -mtune=generic -c rx_ebf0e5096f80d0203de27cbc48be0f26_x64.c -o rx_ebf0e5096f80d0203de27cbc48be0f26_x64.o
C:/RTools/3.4/mingw_64/bin/gcc -shared -s -static-libgcc -o rx_ebf0e5096f80d0203de27cbc48be0f26_x64.dll tmp.def rx_ebf0e5096f80d0203de27cbc48be0f26_x64.o -LC:/R/R-34~1.3/bin/x64 -lRblas -LC:/R/R-34~1.3/bin/x64 -lRlapack -lgfortran -lm -lquadmath -LC:/R/R-34~1.3/bin/x64 -lR
C:/RTools/3.4/mingw_64/bin/gcc  -I"C:/R/R-34~1.3/include" -DNDEBUG          -O2 -Wall  -std=gnu99 -mtune=generic -c rx_40eae078bd4f9e0651450c3b310d5e9a_x64.c -o rx_40eae078bd4f9e0651450c3b310d5e9a_x64.o
C:/RTools/3.4/mingw_64/bin/gcc -shared -s -static-libgcc -o rx_40eae078bd4f9e0651450c3b310d5e9a_x64.dll tmp.def rx_40eae078bd4f9e0651450c3b310d5e9a_x64.o -LC:/R/R-34~1.3/bin/x64 -lRblas -LC:/R/R-34~1.3/bin/x64 -lRlapack -lgfortran -lm -lquadmath -LC:/R/R-34~1.3/bin/x64 -lR
The model-based sensitivities have been calculated.
It will be cached for future runs.
Diagonal form: sqrt
     [,1]
[1,]    1
Calculate symbolic inverse...done
Calculate symbolic determinant of inverse...done
Calculate d(Omega)/d(Est) and d(Omega^-1)/d(Est)...
.0
...done.
Calculating Table Variables...
done
Warning message:
In nlmixrUI.nlme.var(obj) :
  Initial condition for additive error ignored with nlme
nlmixr UI combined dataset and properties
 $ par.hist         : Parameter history (if available)
 $ par.hist.stacked : Parameter history in stacked form for easy plotting (if available)
 $ par.fixed        : Fixed Effect Parameter Table
 $ eta              : Individual Parameter Estimates
 $ seed             : Seed (if applicable)
 $ model.name       : Model name (from R function)
 $ data.name        : Name of R data input
> # Only 1 iteration for DFN.
> fit.nlme.1cm.lin.theo_sd.init <- nlmixr(uif.1cm.lin, theo_sd, est="nlme",
+                                         control = nlmeControl(maxIter = 1,
+                                                               pnlsMaxIter = 1,
+                                                               msMaxIter = 1)
+                                         )

**Iteration 1
LME step: Loglik: -182.2318, nlminb iterations: 1
reStruct  parameters:
      ID1       ID2       ID3 
1.0022617 0.2783009 2.0758984 
PNLS step: RSS =  87.52968 
 fixed effects: -3.230002  0.5346704  -0.7523971  
 iterations: 1 
Convergence crit. (must all become <= tolerance = 1e-05):
       fixed     reStruct 
3.290853e-01 4.201928e-06 
Error in nlme.formula(model = DV ~ (nlmixr::nlmeModList("user_fn"))(logtcl,  : 
  maximum number of iterations (maxIter = 1) reached without convergence
In addition: Warning message:
In nlmixrUI.nlme.var(obj) :
  Initial condition for additive error ignored with nlme
> # No object created.
> # So try to expand the tolerance:
> fit.nlme.lin.theo_sd.init <- nlmixr(uif.1cm.lin, theo_sd, est="nlme",
+                                     control = nlmeControl(maxIter = 1,
+                                                           pnlsMaxIter = 1,
+                                                           msMaxIter = 1),
+                                     tolerance = 1 # deliberately a huge tolerance
+                                     )
Error in (function (model, data = sys.frame(sys.parent()), fixed, random = fixed,  : 
  unused argument (tolerance = 1)
In addition: Warning message:
In nlmixrUI.nlme.var(obj) :
  Initial condition for additive error ignored with nlme
> # No object created.

[mattfidler]

I was mistaken; nlme doesn't return the estimates if the maximum iterations were reached.

I have opened a feature request bug for you #32

[JSC67]

Thank-you, Matt. Having something like NONMEM's MAXEVAL=0 would indeed be helpful to get just likelihoods and related stats (AIC, BIC, ...) without doing minimizations.