search

nlmixrdevelopment / nlmixr

nlmixr: an R package for population PKPD modeling
https://nlmixrdevelopment.github.io/nlmixr/
GNU General Public License v2.0
116 stars 45 forks source link

Estimate posthoc predictions and/or likelihood without doing a full estimation #32

Closed mattfidler closed 6 years ago

mattfidler commented 6 years ago

Based on a model and data, many people want to:

This is disccussed in Issue #28.

JSC67 commented 6 years ago

Hi Matt (or whoever takes up this issue),

Since the ability to calculate post-hoc predictions given only a fully specified model and some data is essential to the ability to do cross-validation (using the testing set), I sure hope that this issue will be bumped up the priority list. I tried to do so using nlme's predict() function as so:

predict(object=trainfits$`BS_iteration: 1`,
        newdata=ans4[ans4$BsIteration==1,])
# returns this:
 Compiling model...done
 Error in rxSolve.default(object = <environment>, newdata = ..1, events = list( :
   ..1 used in an incorrect context, no ... to look in

This error message is opaque to me. Also unclear to me is nlmixr's prediction() function for gnlmm objects since the example in the documentation has very little in the way of explanation. What would the fit and pred arguments be in my case, just wanting to do cross-validation? I'd like to be able to do a fitting—I don't much care if it's a gnlmm fitting or one more like in the vignette tutorial—and then use the results of that fitting on new data, i.e. to do cross-validation. I'm not sure if it would be better to open a new issue for predict() or prediction() or just comment on this open issue, but for now I opt for the latter.

Thanks for hearing me out.

James

mattfidler commented 6 years ago

Hi James,

The FOCEi branch includes the new posthoc feature you are requesting. However, it is a major rewrite of the nlmixr output object, so it isn't compatible with shinyMixR or xpose.nlmixr right now.

I don't know how nlme's predict function works. If I get time, I will look into it.

Higher on the priority list for me is this posthoc and a few other issues like multiple end-points.

JSC67 commented 6 years ago

Hi Matt,

Just to clarify, the FOCEi branch is still "under construction", correct? If it's working, could you please let us know the syntax it needs? As a major rewrite of the nlmixr output object, my guess is that it will require a new version number to download, no?

Thanks for you help! -James

mattfidler commented 6 years ago

While it is working, the syntax isn't finalized. I would like to review suggestions with the nlmixr team first.

mattfidler commented 6 years ago

Yes. It will require a new version download.

mattfidler commented 6 years ago

With the new version you woud do the following:

post <- nlmixr(model,newData,"posthoc")

We are working toward allowing you to pip the results of a model to a posthoc prediction possibly with R's predict function

JSC67 commented 6 years ago

I wonder if you could elaborate a bit on this syntax, in particular the word "model". The structural model is given by the UUI, but how does one specify exactly the parameters? For example after a bootstrap run the parameters would be given as fit$theta, fit$omega, etc. and will vary from one bootstrap sample to the next.

Using the theo_sd data set included in nlmixr, I fit the data as in the vignette and found this behavior to be strange when trying to update the value of theta:

> fit <- nlmixr(cf1cmt.adderr, theo_sd, est="nlme")  # for demo and testing

**Iteration 1
LME step: Loglik: -182.2318, nlminb iterations: 1
reStruct  parameters:
      ID1       ID2       ID3 
1.0022617 0.2783009 2.0758984 
 Beginning PNLS step: ..  completed fit_nlme() step.
PNLS step: RSS =  63.1492 
 fixed effects: -3.211935  0.447998  -0.7859319  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
    fixed  reStruct 
0.2723749 3.2673493 

**Iteration 2
LME step: Loglik: -179.291, nlminb iterations: 9
reStruct  parameters:
       ID1        ID2        ID3 
0.96385384 0.08333465 1.63552496 
 Beginning PNLS step: ..  completed fit_nlme() step.
PNLS step: RSS =  63.28224 
 fixed effects: -3.211529  0.4441499  -0.7863391  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
      fixed    reStruct 
0.008664056 0.172140362 

**Iteration 3
LME step: Loglik: -179.337, nlminb iterations: 8
reStruct  parameters:
       ID1        ID2        ID3 
0.96127446 0.07085101 1.63887826 
 Beginning PNLS step: ..  completed fit_nlme() step.
PNLS step: RSS =  63.21956 
 fixed effects: -3.211695  0.445849  -0.7861748  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
      fixed    reStruct 
0.003810953 0.083041222 

**Iteration 4
LME step: Loglik: -179.3201, nlminb iterations: 7
reStruct  parameters:
       ID1        ID2        ID3 
0.96243389 0.07619035 1.63736920 
 Beginning PNLS step: ..  completed fit_nlme() step.
PNLS step: RSS =  63.24586 
 fixed effects: -3.211622  0.4451181  -0.7862475  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
      fixed    reStruct 
0.001642014 0.034519536 

**Iteration 5
LME step: Loglik: -179.328, nlminb iterations: 6
reStruct  parameters:
       ID1        ID2        ID3 
0.96192227 0.07381955 1.63801887 
 Beginning PNLS step: ..  completed fit_nlme() step.
PNLS step: RSS =  63.23423 
 fixed effects: -3.211654  0.4454379  -0.7862161  
 iterations: 7 
Convergence crit. (must all become <= tolerance = 1e-05):
       fixed     reStruct 
0.0007179909 0.0153297789 

**Iteration 6
LME step: Loglik: -179.3247, nlminb iterations: 5
reStruct  parameters:
       ID1        ID2        ID3 
0.96215590 0.07486979 1.63772980 
 Beginning PNLS step: ..  completed fit_nlme() step.
PNLS step: RSS =  63.24087 
 fixed effects: -3.211654  0.4454379  -0.7862161  
 iterations: 1 
Convergence crit. (must all become <= tolerance = 1e-05):
      fixed    reStruct 
0.000000000 0.008689733 

**Iteration 7
LME step: Loglik: -179.3247, nlminb iterations: 1
reStruct  parameters:
       ID1        ID2        ID3 
0.96214642 0.07485047 1.63773353 
 Beginning PNLS step: ..  completed fit_nlme() step.
PNLS step: RSS =  63.24087 
 fixed effects: -3.211654  0.4454379  -0.7862161  
 iterations: 1 
Convergence crit. (must all become <= tolerance = 1e-05):
       fixed     reStruct 
0.000000e+00 1.297069e-11 
Calculating residuals/tables
done.
Warning message:
In nlmixrUI.nlme.var(obj) :
  Initial condition for additive error ignored with nlme
> fit$theta
    logtcl     logtka      logtv    add.err 
-3.2116536  0.4454379 -0.7862161  0.6921686 
> length(fit$theta)
[1] 4
> fit$theta[2] <- 0.4  # change the second value to something close
> length(fit$theta)
[1] 132
> fit$theta
  [1] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
  [7] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
 [13] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
 [19] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
 [25] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
 [31] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
 [37] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
 [43] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
 [49] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
 [55] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
 [61] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
 [67] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
 [73] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
 [79] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
 [85] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
 [91] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
 [97] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
[103] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
[109] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
[115] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
[121] -3.2116536  0.4000000 -0.7862161  0.6921686 -3.2116536  0.4000000
[127] -0.7862161  0.6921686 -3.2116536  0.4000000 -0.7862161  0.6921686
>

Why did the length of fit$theta change? What would be the syntactically best way to set values of theta and omega and their correlation, as e.g. using the results of one fitting?

Thanks for your help, again. James

mattfidler commented 6 years ago

Why did the length of fit$theta change?

Fundamentally, a fit object is a data.frame. Therefore assigning fit$theta changes the data frame, not the properties you are trying to change. While the theta can be changed in a different way, it will not change the fit parameters, and really shouldn't be done.

What would be the syntactically best way to set values of theta and omega and their correlation, as e.g. using the results of one fitting?

Currently the fit has a parsed UI object attached at fit$uif. Not only does this have the last fit information but that is what should be modified to run a new fit.

Note the parameter values and properties are stored in a bounds object fit$uif$ini. You can sort of see the structure by executing as.data.frame(fit$uif$ini) If you modify these parameters in a new ui object, you can then call nlmixr(newParsedUI, data, "estimation")

Eventually I'm hoping that I can get a general function to do this. In the last revision, I made ini and model nlmixr functions, so you could possibly use a pipeline in the future

fit %>% ini({cl=1}) %>% nlmixr(data,"posthoc")

or

fit %>% ini(list(cl=1)) %>% nlmixr(data,"posthoc")

However these details haven't been worked out yet. As a side note, now nlmixr functions now mostly evaluate to the parsed UI object instead of throwing an error. Some information (like parameter labels) are ignored by this approach, however.

mattfidler commented 6 years ago

Note what actually happend with your theta is the following:

  1. The theta vector was modified
  2. Since the data frame is a multiple of the number of thetas, it assigns and repeats this new theta in your modeling data frame.
JSC67 commented 6 years ago

You wrote, "If you modify these parameters in a new ui object, you can then call nlmixr(newParsedUI, data, "estimation")." I've tried various things (below), but nothing seems to work. How would I get newParsedUI, the first argument to nlmixr? Or is this something that I'll just have to wait on (which for me is okay for a bit, as I'll be unavailable to work on nlmixr stuff for the next 10 days)?

> fit$uif$ini
Fixed Effects ($theta):
    logtcl     logtka      logtv 
-3.2116536  0.4454379 -0.7862161 

Omega ($omega):
           eta.cl    eta.ka     eta.v
eta.cl 0.06993812 0.0000000 0.0000000
eta.ka 0.00000000 0.4124862 0.0000000
eta.v  0.00000000 0.0000000 0.0181095
> fit$uif$ini$theta
    logtcl     logtka      logtv 
-3.2116536  0.4454379 -0.7862161 
> fit$uif$ini$omega
           eta.cl    eta.ka     eta.v
eta.cl 0.06993812 0.0000000 0.0000000
eta.ka 0.00000000 0.4124862 0.0000000
eta.v  0.00000000 0.0000000 0.0181095
> is.matrix(fit$uif$ini$omega)
[1] TRUE
> is.vector(fit$uif$ini$theta)
[1] TRUE
> as.data.frame(fit$uif$ini)
  ntheta neta1 neta2    name lower         est upper   fix  err
1      1    NA    NA  logtcl  -Inf -3.21165355   Inf FALSE <NA>
2      2    NA    NA  logtka  -Inf  0.44543795   Inf FALSE <NA>
3      3    NA    NA   logtv  -Inf -0.78621611   Inf FALSE <NA>
4      4    NA    NA add.err  -Inf  0.69216856   Inf FALSE  add
5     NA     1     1  eta.cl  -Inf  0.06993812   Inf FALSE <NA>
6     NA     2     2  eta.ka  -Inf  0.41248621   Inf FALSE <NA>
7     NA     3     3   eta.v  -Inf  0.01810950   Inf FALSE <NA>
                        label condition
1 log typical value Cl (L/hr)      <NA>
2 log typical value Ka (1/hr)      <NA>
3    log typical value V  (L)      <NA>
4                        <NA>      <NA>
5                        <NA>        ID
6                        <NA>        ID
7                        <NA>        ID
> current.est <- as.data.frame(fit$uif$ini)$est
> current.est
[1] -3.21165355  0.44543795 -0.78621611  0.69216856  0.06993812  0.41248621
[7]  0.01810950
> current.est[2] <- 0.4 ; current.est
[1] -3.21165355  0.40000000 -0.78621611  0.69216856  0.06993812  0.41248621
[7]  0.01810950
> new.est <- current.est
> as.data.frame(fit$uif$ini)
  ntheta neta1 neta2    name lower         est upper   fix  err
1      1    NA    NA  logtcl  -Inf -3.21165355   Inf FALSE <NA>
2      2    NA    NA  logtka  -Inf  0.44543795   Inf FALSE <NA>
3      3    NA    NA   logtv  -Inf -0.78621611   Inf FALSE <NA>
4      4    NA    NA add.err  -Inf  0.69216856   Inf FALSE  add
5     NA     1     1  eta.cl  -Inf  0.06993812   Inf FALSE <NA>
6     NA     2     2  eta.ka  -Inf  0.41248621   Inf FALSE <NA>
7     NA     3     3   eta.v  -Inf  0.01810950   Inf FALSE <NA>
                        label condition
1 log typical value Cl (L/hr)      <NA>
2 log typical value Ka (1/hr)      <NA>
3    log typical value V  (L)      <NA>
4                        <NA>      <NA>
5                        <NA>        ID
6                        <NA>        ID
7                        <NA>        ID
> as.data.frame(fit$uif$ini)$est <- current.est
Error in as.data.frame(fit$uif$ini)$est <- current.est : 
  could not find function "as.data.frame<-"
> as.data.frame(fit$uif$ini)$est
[1] -3.21165355  0.44543795 -0.78621611  0.69216856  0.06993812  0.41248621
[7]  0.01810950
> as.data.frame(fit$uif$ini)
  ntheta neta1 neta2    name lower         est upper   fix  err
1      1    NA    NA  logtcl  -Inf -3.21165355   Inf FALSE <NA>
2      2    NA    NA  logtka  -Inf  0.44543795   Inf FALSE <NA>
3      3    NA    NA   logtv  -Inf -0.78621611   Inf FALSE <NA>
4      4    NA    NA add.err  -Inf  0.69216856   Inf FALSE  add
5     NA     1     1  eta.cl  -Inf  0.06993812   Inf FALSE <NA>
6     NA     2     2  eta.ka  -Inf  0.41248621   Inf FALSE <NA>
7     NA     3     3   eta.v  -Inf  0.01810950   Inf FALSE <NA>
                        label condition
1 log typical value Cl (L/hr)      <NA>
2 log typical value Ka (1/hr)      <NA>
3    log typical value V  (L)      <NA>
4                        <NA>      <NA>
5                        <NA>        ID
6                        <NA>        ID
7                        <NA>        ID
> (as.data.frame(fit$uif$ini)$est) <- current.est
Error in (as.data.frame(fit$uif$ini)$est) <- current.est : 
  could not find function "(<-"
> (as.data.frame(fit$uif$ini))$est <- current.est
Error in (as.data.frame(fit$uif$ini))$est <- current.est : 
  could not find function "(<-"
> df <- as.data.frame(fit$uif$ini)
> df
  ntheta neta1 neta2    name lower         est upper   fix  err
1      1    NA    NA  logtcl  -Inf -3.21165355   Inf FALSE <NA>
2      2    NA    NA  logtka  -Inf  0.44543795   Inf FALSE <NA>
3      3    NA    NA   logtv  -Inf -0.78621611   Inf FALSE <NA>
4      4    NA    NA add.err  -Inf  0.69216856   Inf FALSE  add
5     NA     1     1  eta.cl  -Inf  0.06993812   Inf FALSE <NA>
6     NA     2     2  eta.ka  -Inf  0.41248621   Inf FALSE <NA>
7     NA     3     3   eta.v  -Inf  0.01810950   Inf FALSE <NA>
                        label condition
1 log typical value Cl (L/hr)      <NA>
2 log typical value Ka (1/hr)      <NA>
3    log typical value V  (L)      <NA>
4                        <NA>      <NA>
5                        <NA>        ID
6                        <NA>        ID
7                        <NA>        ID
> df$est <- current.est
> df
  ntheta neta1 neta2    name lower         est upper   fix  err
1      1    NA    NA  logtcl  -Inf -3.21165355   Inf FALSE <NA>
2      2    NA    NA  logtka  -Inf  0.40000000   Inf FALSE <NA>
3      3    NA    NA   logtv  -Inf -0.78621611   Inf FALSE <NA>
4      4    NA    NA add.err  -Inf  0.69216856   Inf FALSE  add
5     NA     1     1  eta.cl  -Inf  0.06993812   Inf FALSE <NA>
6     NA     2     2  eta.ka  -Inf  0.41248621   Inf FALSE <NA>
7     NA     3     3   eta.v  -Inf  0.01810950   Inf FALSE <NA>
                        label condition
1 log typical value Cl (L/hr)      <NA>
2 log typical value Ka (1/hr)      <NA>
3    log typical value V  (L)      <NA>
4                        <NA>      <NA>
5                        <NA>        ID
6                        <NA>        ID
7                        <NA>        ID
> # Why did this work but 'as.data.frame(fit$uif$ini)$est <- current.est' failed?
> nlmixr(df, theo_sd, "estimation")
Error in UseMethod("nlmixr") : 
  no applicable method for 'nlmixr' applied to an object of class "data.frame"
> fit$uif$ini$est
[1] -3.21165355  0.44543795 -0.78621611  0.69216856  0.06993812  0.41248621
[7]  0.01810950
> fit$uif$ini$est <- current.est
Error in `$<-.data.frame`(`*tmp*`, uif, value = list(ini = list(ntheta = c(1,  : 
  replacement has 5 rows, data has 132
> length(current.est)
[1] 7
> ## Huh- "replacement has 5 rows"?
> model(df)
> model.df <- model(df)
> model.df
NULL
> model(df)
> df
  ntheta neta1 neta2    name lower         est upper   fix  err
1      1    NA    NA  logtcl  -Inf -3.21165355   Inf FALSE <NA>
2      2    NA    NA  logtka  -Inf  0.40000000   Inf FALSE <NA>
3      3    NA    NA   logtv  -Inf -0.78621611   Inf FALSE <NA>
4      4    NA    NA add.err  -Inf  0.69216856   Inf FALSE  add
5     NA     1     1  eta.cl  -Inf  0.06993812   Inf FALSE <NA>
6     NA     2     2  eta.ka  -Inf  0.41248621   Inf FALSE <NA>
7     NA     3     3   eta.v  -Inf  0.01810950   Inf FALSE <NA>
                        label condition
1 log typical value Cl (L/hr)      <NA>
2 log typical value Ka (1/hr)      <NA>
3    log typical value V  (L)      <NA>
4                        <NA>      <NA>
5                        <NA>        ID
6                        <NA>        ID
7                        <NA>        ID
> ini.df <- ini(df)
> ini.df
NULL
> # So how to reparse it?

Thanks again for all your help. I really look forward to this functionality! :)

James

JSC67 commented 6 years ago

Hi Matt,

You wrote, "Note the parameter values and properties are stored in a bounds object fit$uif$ini. You can sort of see the structure by executing as.data.frame(fit$uif$ini). If you modify these parameters in a new ui object, you can then call nlmixr(newParsedUI, data, "estimation")." While fit$uif$ini may indeed be a data frame as you wrote it is (in my case with 7 observations of 11 variables), it doesn't seem to be an ordinary data frame since it also has other fields associated with it:

> fit$uif$ini
Fixed Effects ($theta):
    logtcl     logtka      logtv 
-3.2116536  0.4454379 -0.7862161 

Omega ($omega):
           eta.cl    eta.ka     eta.v
eta.cl 0.06993812 0.0000000 0.0000000
eta.ka 0.00000000 0.4124862 0.0000000
eta.v  0.00000000 0.0000000 0.0181095
> str(fit$uif$ini)
'data.frame':   7 obs. of  11 variables:
 $ ntheta   : num  1 2 3 4 NA NA NA
 $ neta1    : num  NA NA NA NA 1 2 3
 $ neta2    : num  NA NA NA NA 1 2 3
 $ name     : Factor w/ 7 levels "logtcl","logtka",..: 1 2 3 4 5 6 7
 $ lower    : num  -Inf -Inf -Inf -Inf -Inf ...
 $ est      : num  -3.2117 0.4454 -0.7862 0.6922 0.0699 ...
 $ upper    : num  Inf Inf Inf Inf Inf ...
 $ fix      : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ err      : chr  NA NA NA "add" ...
 $ label    : chr  "log typical value Cl (L/hr)" "log typical value Ka (1/hr)" "log typical value V  (L)" NA ...
 $ condition: chr  NA NA NA "" ...
 $ theta     : num ... (theta estimates)
 $ theta.full: num ... (theta estimates, including error terms)
 $ omega     : matrix ... (omega matrix)
 $ random    : matrix class ... (Based on Between Subject Random effects)
 $ fixed.form: formula  ... (Fixed effect parameters based on theta.)
 $ focei.upper: Upper bounds for FOCEi
 $ focei.lower: Lower bounds for FOCEi
 $ focei.err.type: Residual Error type for FOCEi thetas
 $ eta.names: Eta names
 $ focei.names: Theta names for FOCEi
> print.default(fit$uif$ini)
$`ntheta`
[1]  1  2  3  4 NA NA NA

$neta1
[1] NA NA NA NA  1  2  3

$neta2
[1] NA NA NA NA  1  2  3

$name
[1] logtcl  logtka  logtv   add.err eta.cl  eta.ka  eta.v  
Levels: logtcl logtka logtv add.err eta.cl eta.ka eta.v

$lower
[1] -Inf -Inf -Inf -Inf -Inf -Inf -Inf

$est
[1] -3.21165355  0.44543795 -0.78621611  0.69216856  0.06993812  0.41248621
[7]  0.01810950

$upper
[1] Inf Inf Inf Inf Inf Inf Inf

$fix
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE

$err
[1] NA    NA    NA    "add" NA    NA    NA   

$label
[1] "log typical value Cl (L/hr)" "log typical value Ka (1/hr)"
[3] "log typical value V  (L)"    NA                           
[5] NA                            NA                           
[7] NA                           

$condition
[1] NA   NA   NA   ""   "ID" "ID" "ID"

attr(,"class")
[1] "nlmixrBounds" "data.frame"  
>

Can you please elaborate a little bit on what sort of data frame this is (how it's different from a simple textbook example of a data frame), or perhaps provide a reference where I could learn more?

You also wrote, "Currently the fit has a parsed UI object attached at fit$uif. Not only does this have the last fit information but that is what should be modified to run a new fit." How exactly should this be modified? Are there any accessor or mutator functions we can use for such an object? It would be really great if you could supply us with a preferred syntax for changing the fit$uif$ini parameters so one can then call ' nlmixr(newParsedUI, data, "estimation") ' as you wrote above, and also use something like ' post <- nlmixr(model, newData, "posthoc") ' for cross-validation while keeping the model and its parameters the same as for fitting on a training data set.

Much thanks!

James

mattfidler commented 6 years ago

I'll opten a new feature request for this. This feature didn't make 1.0

JSC67 commented 2 years ago

Hello Matt,

Almost 4 years ago I wrote in Issue # 28, "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." On 2018-09-22 you closed issue # 32 noting that this feature request didn't make it into version 1.0. It seems that a work-around for Bill Denney's issue was adequate for his needs (i.e. set EVID = 2), and that this feature request may have fallen by the wayside. Since it was essential for my interest in cross-validation, I had to drop nlmixr as a viable option for my work on using bootstrap cross-validation for model selection and used Monolix' R API instead. (This has since been published at https://www.tandfonline.com/doi/full/10.1080/19466315.2020.1828159 in 2020.) For cross-validation (and for using cross-validation for model selection), it is necessary to have a way to compare likelihoods for testing data. In NONMEM this would be with the MAXEVAL=0 option.

So now, nearly 4 years later, I would again like to take a serious interest in nlmixr and see its capabilities. I'm glad to see that it's matured in the interim. Can an equivalent to MAXEVAL=0 please be included?

Thank-you, James

mattfidler commented 2 years ago

Hi James,

This has been present in nlmixr for awhile.

The easiest way to use this is:

fit <- nlmixr(model, data, "posthoc")

However you can technically do the following too:

fit <- nlmixr(model, data, "focei", control=foceiControl(maxOuterIterations=0))

This will find the empirical Bayes estimates based on the prior model through optimizing the "inner" problem. This is what maxevals=0 does.

You can even take it a step further than NONMEM does by supplying the etas (in the order of the UI and in matrix form):

fit <- nlmixr(model, data, "focei", control=foceiControl(maxOuterIterations=0, maxInnerIterations=0, etaMat=etas))

which evaluates an objective function at the ETAs that were pre-specified.

This feature has been there for a very long time, near the time of 1.0. It is used to get the focei likelihood for saem

JSC67 commented 2 years ago

Matt, I appreciate your reply. I was unaware of this, so thank-you! -James

From: Matthew Fidler @.> Sent: Monday, January 17, 2022 5:41 AM To: nlmixrdevelopment/nlmixr @.> Cc: JSC67 @.>; Comment @.> Subject: Re: [nlmixrdevelopment/nlmixr] Estimate posthoc predictions and/or likelihood without doing a full estimation (#32)

Hi James,

This has been present in nlmixr for awhile.

The easiest way to use this is:

fit <- nlmixr(model, data, "posthoc")

However you can technically do the following too:

fit <- nlmixr(model, data, "focei", control=foceiControl(maxOuterIterations=0))

This will find the empirical Bayes estimates based on the prior model through optimizing the "inner" problem. This is what maxevals=0 does.

You can even take it a step further than NONMEM does by supplying the etas (in the order of the UI and in matrix form):

fit <- nlmixr(model, data, "focei", control=foceiControl(maxOuterIterations=0, maxInnerIterations=0, etaMat=etas))

which evaluates an objective function at the ETAs that were pre-specified.

This feature has been there for a very long time, near the time of 1.0. It is used to get the focei likelihood for saem

— Reply to this email directly, view it on GitHubhttps://github.com/nlmixrdevelopment/nlmixr/issues/32#issuecomment-1014561363, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AHIJWUNYLEFROQBNHMSUM23UWQMA5ANCNFSM4ESE2TCA. Triage notifications on the go with GitHub Mobile for iOShttps://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Androidhttps://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you commented.Message ID: @.***>