Objective function value using SAEM

karthikl6 commented 5 years ago

Hi,

When I run the example from here, the output of fit does not show any OBJF using SAEM estimation method. But says OBJF not calculated. In the example, it shows OBJF calculated from FOCEi approximation.

Am I missing anything?

Thanks

mattfidler commented 5 years ago

Hi @karthikl6,

SAEM does not calculate the objective function during optimization. At the last release we were working with SAEM logLik based on adaptive Gaussian quadrature, so the idea was you could choose what method you wanted when calculating likelihood.

The github version is more explicit about this and accepts things like AIC(fit) and calculates the likelihood as needed. However, you can still get the FOCEi objective function.

Missing objective function (github version)

In this case no options are added, and the objective function calculation is calculated on demand.

fit <- nlmixr(one.compartment.saem,theo_sd,"saem",control=saemControl(print=500));

``` r library(nlmixr) one.compartment.saem <- function() { ini({ tka <- .5 # Log Ka tcl <- -3.2 # Log Cl tv <- -1 # Log V eta.ka ~ 1 eta.cl ~ 2 eta.v ~ 1 add.err <- 0.1 }) model({ ka <- exp(tka + eta.ka) cl <- exp(tcl + eta.cl) v <- exp(tv + eta.v) d/dt(depot) = -ka * depot d/dt(center) = ka * depot - cl / v * center cp = center / v cp ~ add(add.err) }) } fit <- nlmixr(one.compartment.saem,theo_sd,"saem",control=saemControl(print=500)); #> Compiling RxODE equations...done. #> 1: 0.1503 -3.3314 -0.6447 0.9500 1.9000 0.9500 9.6167 #> 500: 0.4546 -3.2153 -0.7848 0.4413 0.0713 0.0180 0.4774 #> Calculating residuals/tables #> done. fit #> ── nlmixr SAEM(ODE); OBJF not calculated fit ────────────────────────────── #> Gaussian/Laplacian Likelihoods: AIC(fit) or fit$objf etc. #> FOCEi CWRES & Likelihoods: addCwres(fit) #> #> ── Time (sec; fit$time): ────────────────────────────────────────────────── #> saem setup optimize covariance table #> 4.074 0.771227 4e-06 5e-06 0.006 #> #> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──────────────── #> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) #> tka Log Ka 0.455 0.198 43.4 1.58 (1.07, 2.32) 74.5% #> tcl Log Cl -3.22 0.0814 2.53 0.0401 (0.0342, 0.0471) 27.2% #> tv Log V -0.785 0.0426 5.43 0.456 (0.42, 0.496) 13.5% #> add.err 0.691 0.691 #> Shrink(SD)% #> tka -0.520% #> tcl 3.76% #> tv 11.0% #> add.err #> #> Covariance Type (fit$covMethod): fim #> No correlations in between subject variability (BSV) matrix #> Full BSV covariance (fit$omega) #> or correlation (fit$omegaR; diagonals=SDs) #> Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink #> #> ── Fit Data (object fit is a modified tibble): ──────────────────────────── #> # A tibble: 132 x 18 #> ID TIME DV EVID PRED RES IPRED IRES IWRES eta.ka eta.cl #> #> 1 1 0 0.74 0 0 0.74 0 0.74 1.07 0.117 -0.625 #> 2 1 0.25 2.84 0 2.84 0.00441 3.86 -1.02 -1.48 0.117 -0.625 #> 3 1 0.570 6.57 0 5.07 1.50 6.79 -0.219 -0.316 0.117 -0.625 #> # … with 129 more rows, and 7 more variables: eta.v , ka , #> # cl , v , cp , depot , center ## Now you can access AIC and other sorts of metrics, even though it wasn't calculated. AIC(fit) #> Calculating Laplace -2LL (nsd=3) #> [1] 283.8811 ## Afterward the fit is updated with the objective function method used fit #> ── nlmixr SAEM(ODE); OBJF by Lapalcian (n.sd=3) fit ─────────────────────── #> OBJF AIC BIC Log-likelihood Condition Number #> laplace3 269.8811 283.8811 286.6664 -134.9406 21.72835 #> #> ── Time (sec; fit$time): ────────────────────────────────────────────────── #> saem setup optimize covariance table #> 4.074 0.771227 4e-06 5e-06 0.006 #> #> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──────────────── #> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) #> tka Log Ka 0.455 0.198 43.4 1.58 (1.07, 2.32) 74.5% #> tcl Log Cl -3.22 0.0814 2.53 0.0401 (0.0342, 0.0471) 27.2% #> tv Log V -0.785 0.0426 5.43 0.456 (0.42, 0.496) 13.5% #> add.err 0.691 0.691 #> Shrink(SD)% #> tka -0.520% #> tcl 3.76% #> tv 11.0% #> add.err #> #> Covariance Type (fit$covMethod): fim #> No correlations in between subject variability (BSV) matrix #> Full BSV covariance (fit$omega) #> or correlation (fit$omegaR; diagonals=SDs) #> Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink #> #> ── Fit Data (object fit is a modified tibble): ──────────────────────────── #> # A tibble: 132 x 18 #> ID TIME DV EVID PRED RES IPRED IRES IWRES eta.ka eta.cl #> #> 1 1 0 0.74 0 0 0.74 0 0.74 1.07 0.117 -0.625 #> 2 1 0.25 2.84 0 2.84 0.00441 3.86 -1.02 -1.48 0.117 -0.625 #> 3 1 0.570 6.57 0 5.07 1.50 6.79 -0.219 -0.316 0.117 -0.625 #> # … with 129 more rows, and 7 more variables: eta.v , ka , #> # cl , v , cp , depot , center ``` ^{Created on 2019-03-27 by the [reprex package](https://reprex.tidyverse.org) (v0.2.1)}

Requesting CWRES (and FOCEi objective function directly)

This is simply by adding the table options:

fit <- nlmixr(one.compartment.saem,theo_sd,"saem",control=saemControl(print=500),
        table=tableControl(cwres=TRUE));

Also note you can add npde residuals by this option too.

``` r library(nlmixr) one.compartment.saem <- function() { ini({ tka <- .5 # Log Ka tcl <- -3.2 # Log Cl tv <- -1 # Log V eta.ka ~ 1 eta.cl ~ 2 eta.v ~ 1 add.err <- 0.1 }) model({ ka <- exp(tka + eta.ka) cl <- exp(tcl + eta.cl) v <- exp(tv + eta.v) d/dt(depot) = -ka * depot d/dt(center) = ka * depot - cl / v * center cp = center / v cp ~ add(add.err) }) } fit <- nlmixr(one.compartment.saem,theo_sd,"saem",control=saemControl(print=500), table=tableControl(cwres=TRUE)); #> Compiling RxODE equations...done. #> 1: 0.1503 -3.3314 -0.6447 0.9500 1.9000 0.9500 9.6167 #> 500: 0.4546 -3.2153 -0.7848 0.4413 0.0713 0.0180 0.4774 #> Calculate ETA-based prediction and error derivatives: #> Using sympy via reticulate #> Calculate Jacobian ``` #> Calculate Sensitivities #> Load into sympy...done #> Calculate ∂(f)/∂(η) #> Calculate ∂(R²)/∂(η) #> Freeing Python/SymPy memory...done #> ################################################################################ #> Optimizing expressions in Predictions/EBE model...done #> Compiling Predictions/EBE model...done #> Optimizing expressions in Inner model...done #> Compiling Inner model...done #> The model-based sensitivities have been calculated. #> Calculating residuals/tables #> done. fit #> ── nlmixr SAEM(ODE); OBJF not calculated fit ────────────────────────────── #> OBJF AIC BIC Log-likelihood Condition Number #> FOCEi 116.1543 150.3709 153.1562 -68.18546 21.72835 #> #> ── Time (sec; fit$time): ────────────────────────────────────────────────── #> saem setup optimize covariance table #> 4.126 10.94404 3e-06 5e-06 0.013 #> #> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──────────────── #> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) #> tka Log Ka 0.455 0.198 43.4 1.58 (1.07, 2.32) 74.5% #> tcl Log Cl -3.22 0.0814 2.53 0.0401 (0.0342, 0.0471) 27.2% #> tv Log V -0.785 0.0426 5.43 0.456 (0.42, 0.496) 13.5% #> add.err 0.691 0.691 #> Shrink(SD)% #> tka -0.520% #> tcl 3.76% #> tv 11.0% #> add.err #> #> Covariance Type (fit$covMethod): fim #> No correlations in between subject variability (BSV) matrix #> Full BSV covariance (fit$omega) #> or correlation (fit$omegaR; diagonals=SDs) #> Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink #> #> ── Fit Data (object fit is a modified tibble): ──────────────────────────── #> # A tibble: 132 x 22 #> ID TIME DV EVID PRED RES WRES IPRED IRES IWRES CPRED #> #> 1 1 0 0.74 0 0 0.74 1.06 0 0.74 1.07 0 #> 2 1 0.25 2.84 0 2.84 0.00441 0.00634 3.86 -1.02 -1.48 2.68 #> 3 1 0.570 6.57 0 5.07 1.50 2.15 6.79 -0.219 -0.316 4.84 #> # … with 129 more rows, and 11 more variables: CRES , CWRES , #> # eta.ka , eta.cl , eta.v , ka , cl , v , #> # cp , depot , center ^{Created on 2019-03-27 by the [reprex package](https://reprex.tidyverse.org) (v0.2.1)}

Requesting FOCEi objective function after fit

In this case, you have a fit without an objective function, but want to add the FOCEi objective function to it; It can be added with the addCwres() function

fit <- nlmixr(one.compartment.saem,theo_sd,"saem",control=saemControl(print=500));
fit <- addCwres(fit) # Add CWRES & FOCEi objective function

Note you can also add npdes to a fit with addNpde(fit)

``` r library(nlmixr) one.compartment.saem <- function() { ini({ tka <- .5 # Log Ka tcl <- -3.2 # Log Cl tv <- -1 # Log V eta.ka ~ 1 eta.cl ~ 2 eta.v ~ 1 add.err <- 0.1 }) model({ ka <- exp(tka + eta.ka) cl <- exp(tcl + eta.cl) v <- exp(tv + eta.v) d/dt(depot) = -ka * depot d/dt(center) = ka * depot - cl / v * center cp = center / v cp ~ add(add.err) }) } fit <- nlmixr(one.compartment.saem,theo_sd,"saem",control=saemControl(print=500)); #> Compiling RxODE equations...done. #> 1: 0.1503 -3.3314 -0.6447 0.9500 1.9000 0.9500 9.6167 #> 500: 0.4546 -3.2153 -0.7848 0.4413 0.0713 0.0180 0.4774 #> Calculating residuals/tables #> done. ## No FOCEi objective function fit #> ── nlmixr SAEM(ODE); OBJF not calculated fit ────────────────────────────── #> Gaussian/Laplacian Likelihoods: AIC(fit) or fit$objf etc. #> FOCEi CWRES & Likelihoods: addCwres(fit) #> #> ── Time (sec; fit$time): ────────────────────────────────────────────────── #> saem setup optimize covariance table #> 4.091 0.828354 2e-06 3e-06 0.007 #> #> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──────────────── #> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) #> tka Log Ka 0.455 0.198 43.4 1.58 (1.07, 2.32) 74.5% #> tcl Log Cl -3.22 0.0814 2.53 0.0401 (0.0342, 0.0471) 27.2% #> tv Log V -0.785 0.0426 5.43 0.456 (0.42, 0.496) 13.5% #> add.err 0.691 0.691 #> Shrink(SD)% #> tka -0.520% #> tcl 3.76% #> tv 11.0% #> add.err #> #> Covariance Type (fit$covMethod): fim #> No correlations in between subject variability (BSV) matrix #> Full BSV covariance (fit$omega) #> or correlation (fit$omegaR; diagonals=SDs) #> Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink #> #> ── Fit Data (object fit is a modified tibble): ──────────────────────────── #> # A tibble: 132 x 18 #> ID TIME DV EVID PRED RES IPRED IRES IWRES eta.ka eta.cl #> #> 1 1 0 0.74 0 0 0.74 0 0.74 1.07 0.117 -0.625 #> 2 1 0.25 2.84 0 2.84 0.00441 3.86 -1.02 -1.48 0.117 -0.625 #> 3 1 0.570 6.57 0 5.07 1.50 6.79 -0.219 -0.316 0.117 -0.625 #> # … with 129 more rows, and 7 more variables: eta.v , ka , #> # cl , v , cp , depot , center fit <- addCwres(fit) #> Calculating residuals/tables #> done. ## Now fit has the FOCEi objective function fit #> ── nlmixr SAEM(ODE); OBJF not calculated fit ────────────────────────────── #> OBJF AIC BIC Log-likelihood Condition Number #> FOCEi 116.1543 150.3709 153.1562 -68.18546 21.72835 #> #> ── Time (sec; fit$time): ────────────────────────────────────────────────── #> saem setup optimize covariance table #> 0 0.205693 2e-06 3e-06 0.01 #> #> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──────────────── #> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) #> tka Log Ka 0.455 0.198 43.4 1.58 (1.07, 2.32) 74.5% #> tcl Log Cl -3.22 0.0814 2.53 0.0401 (0.0342, 0.0471) 27.2% #> tv Log V -0.785 0.0426 5.43 0.456 (0.42, 0.496) 13.5% #> add.err 0.691 0.691 #> Shrink(SD)% #> tka -0.520% #> tcl 3.76% #> tv 11.0% #> add.err #> #> Covariance Type (fit$covMethod): fim #> No correlations in between subject variability (BSV) matrix #> Full BSV covariance (fit$omega) #> or correlation (fit$omegaR; diagonals=SDs) #> Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink #> #> ── Fit Data (object fit is a modified tibble): ──────────────────────────── #> # A tibble: 132 x 22 #> ID TIME DV EVID PRED RES IPRED IRES IWRES eta.ka eta.cl #> #> 1 1 0 0.74 0 0 0.74 0 0.74 1.07 0.117 -0.625 #> 2 1 0.25 2.84 0 2.84 0.00441 3.86 -1.02 -1.48 0.117 -0.625 #> 3 1 0.570 6.57 0 5.07 1.50 6.79 -0.219 -0.316 0.117 -0.625 #> # … with 129 more rows, and 11 more variables: eta.v , ka , #> # cl , v , cp , depot , center , WRES , #> # CRES , CWRES , CPRED ``` ^{Created on 2019-03-27 by the [reprex package](https://reprex.tidyverse.org) (v0.2.1)}

karthikl6 commented 5 years ago

Hi @mattfidler ,

Thank you very much for the detailed explanation. This helps a lot.

I have another question: Is there a way that we can do change point models like using MTIME in NONMEM?

mattfidler commented 5 years ago

There is a new (undocumented) feature that allows change-point models but it is only available on the github release. But I haven't tested it in a model yet. The MTIME is not the same as what NONMEM had either;

The syntax for mtime can be found here:

https://github.com/nlmixrdevelopment/RxODE/blob/master/tests/testthat/test-mtime.R

This is also different than how NONMEM handles mtimes.

karthikl6 commented 5 years ago

Hi @mattfidler ,

Great!! I will give it a try. Thank you very much!!

mattfidler commented 5 years ago

Great; Let me know if it works for you.

karthikl6 commented 5 years ago

Hi @mattfidler ,

Another quick question: Can we use an explicit equation (Eg: Eff = baseline - slope*Cp) for PD modeling without ODE? I want to do only PD modeling

mattfidler commented 5 years ago

You can

nlmixrdevelopment / nlmixr