nlmixr is an R package for fitting population pharmacokinetic (PK) and pharmacokinetic-pharmacodynamic (PKPD) models.
nlmixr(
object,
data,
est = NULL,
control = list(),
table = tableControl(),
...,
save = NULL,
envir = parent.frame()
)
# S3 method for `function`
nlmixr(
object,
data,
est = NULL,
control = list(),
table = tableControl(),
...,
save = NULL,
envir = parent.frame()
)
# S3 method for nlmixrFitCore
nlmixr(
object,
data,
est = NULL,
control = list(),
table = tableControl(),
...,
save = NULL,
envir = parent.frame()
)
# S3 method for nlmixrUI
nlmixr(
object,
data,
est = NULL,
control = list(),
...,
save = NULL,
envir = parent.frame()
)Fitted object or function specifying the model.
Dataset to estimate. Needs to be RxODE compatible (see https://nlmixrdevelopment.github.io/RxODE/articles/RxODE-event-types.html for detailed dataset requirements).
Estimation method
Estimation control options. They could be
nlmeControl, saemControl or
foceiControl
A list controlling the table options (i.e. CWRES, NPDE etc).
See tableControl.
Other parameters
Boolean to save a nlmixr object in a rds file in the
working directory. If NULL, uses option "nlmixr.save"
Environment that nlmixr is evaluated in.
Either a nlmixr model or a nlmixr fit object
The nlmixr generalized function allows common access to the nlmixr estimation routines.
Rationale
nlmixr estimation routines each have their own way of specifying
models. Often the models are specified in ways that are most
intuitive for one estimation routine, but do not make sense for
another estimation routine. Sometimes, legacy estimation
routines like nlme have their own syntax that
is outside of the control of the nlmixr package.
The unique syntax of each routine makes the routines themselves
easier to maintain and expand, and allows interfacing with
existing packages that are outside of nlmixr (like
nlme). However, a model definition language
that is common between estimation methods, and an output object
that is uniform, will make it easier to switch between estimation
routines and will facilitate interfacing output with external
packages like Xpose.
The nlmixr mini-modeling language, attempts to address this issue by incorporating a common language. This language is inspired by both R and NONMEM, since these languages are familiar to many pharmacometricians.
Initial Estimates and boundaries for population parameters
nlmixr models are contained in a R function with two blocks:
ini and model. This R function can be named
anything, but is not meant to be called directly from R. In fact
if you try you will likely get an error such as Error: could
not find function "ini".
The ini model block is meant to hold the initial estimates
for the model, and the boundaries of the parameters for estimation
routines that support boundaries (note nlmixr's saem
and nlme do not currently support parameter boundaries).
To explain how these initial estimates are specified we will start with an annotated example:
f <- function(){ ## Note the arguments to the function are currently
## ignored by nlmixr
ini({
## Initial conditions for population parameters (sometimes
## called theta parameters) are defined by either `<-` or '='
lCl <- 1.6 #log Cl (L/hr)
## Note that simple expressions that evaluate to a number are
## OK for defining initial conditions (like in R)
lVc = log(90) #log V (L)
## Also a comment on a parameter is captured as a parameter label
lKa <- 1 #log Ka (1/hr)
## Bounds may be specified by c(lower, est, upper), like NONMEM:
## Residuals errors are assumed to be population parameters
prop.err <- c(0, 0.2, 1)
})
## The model block will be discussed later
model({})
}As shown in the above examples:
Simple parameter values are specified as a R-compatible assignment
Boundaries my be specified by c(lower, est, upper).
Like NONMEM, c(lower,est) is equivalent to c(lower,est,Inf)
Also like NONMEM, c(est) does not specify a lower bound, and is equivalent
to specifying the parameter without R's `c` function.
The initial estimates are specified on the variance scale, and in analogy with NONMEM, the square roots of the diagonal elements correspond to coefficients of variation when used in the exponential IIV implementation
These parameters can be named almost any R compatible name. Please note that:
Residual error estimates should be coded as population estimates (i.e. using an '=' or '<-' statement, not a '~').
Naming variables that start with "_" are not supported. Note that R does not
allow variable starting with "_" to be assigned without quoting them.
Naming variables that start with "rx_" or "nlmixr_" is not supported since
RxODE and nlmixr use these prefixes internally for certain estimation
routines and calculating residuals.
Variable names are case sensitive, just like they are in R. "CL" is not the
same as "Cl".
Initial Estimates for between subject error distribution (NONMEM's $OMEGA)
In mixture models, multivariate normal individual deviations from
the population parameters are estimated (in NONMEM these are
called eta parameters). Additionally the
variance/covariance matrix of these deviations is also estimated
(in NONMEM this is the OMEGA matrix). These also have initial
estimates. In nlmixr these are specified by the `~` operator that
is typically used in R for "modeled by", and was chosen to
distinguish these estimates from the population and residual error
parameters.
Continuing the prior example, we can annotate the estimates for the between subject error distribution
f <- function(){
ini({
lCl <- 1.6 #log Cl (L/hr)
lVc = log(90) #log V (L)
lKa <- 1 #log Ka (1/hr)
prop.err <- c(0, 0.2, 1)
## Initial estimate for ka IIV variance
## Labels work for single parameters
eta.ka ~ 0.1 # BSV Ka
## For correlated parameters, you specify the names of each
## correlated parameter separated by a addition operator `+`
## and the left handed side specifies the lower triangular
## matrix initial of the covariance matrix.
eta.cl + eta.vc ~ c(0.1,
0.005, 0.1)
## Note that labels do not currently work for correlated
## parameters. Also do not put comments inside the lower
## triangular matrix as this will currently break the model.
})
## The model block will be discussed later
model({})
}As shown in the above examples:
Simple variances are specified by the variable name and the estimate separated by `~`.
Correlated parameters are specified by the sum of the variable labels and then the lower triangular matrix of the covariance is specified on the left handed side of the equation. This is also separated by `~`.
Currently the model syntax does not allow comments inside the lower triangular matrix.
Model Syntax for ODE based models (NONMEM's $PK, $PRED, $DES and $ERROR)
Once the initialization block has been defined, you can define a
model in terms of the defined variables in the ini block. You can
also mix in RxODE blocks into the model.
The current method of defining a nlmixr model is to specify the parameters, and then possibly the RxODE lines:
Continuing describing the syntax with an annotated example:
f <- function(){
ini({
lCl <- 1.6 #log Cl (L/hr)
lVc <- log(90) #log Vc (L)
lKA <- 0.1 #log Ka (1/hr)
prop.err <- c(0, 0.2, 1)
eta.Cl ~ 0.1 ## BSV Cl
eta.Vc ~ 0.1 ## BSV Vc
eta.KA ~ 0.1 ## BSV Ka
})
model({
## First parameters are defined in terms of the initial estimates
## parameter names.
Cl <- exp(lCl + eta.Cl)
Vc = exp(lVc + eta.Vc)
KA <- exp(lKA + eta.KA)
## After the differential equations are defined
kel <- Cl / Vc;
d/dt(depot) = -KA*depot;
d/dt(centr) = KA*depot-kel*centr;
## And the concentration is then calculated
cp = centr / Vc;
## Last, nlmixr is told that the plasma concentration follows
## a proportional error (estimated by the parameter prop.err)
cp ~ prop(prop.err)
})
}A few points to note:
Parameters are defined before the differential equations. Currently directly defining the differential equations in terms of the population parameters is not supported.
The differential equations, parameters and error terms are in a single block, instead of multiple sections.
State names, calculated variables cannot start with either "rx_"
or "nlmixr_" since these are used internally in some estimation routines.
Errors are specified using the `~`. Currently you can use either add(parameter)
for additive error, prop(parameter) for proportional error or add(parameter1) + prop(parameter2)
for additive plus proportional error. You can also specify norm(parameter) for the additive error,
since it follows a normal distribution.
Some routines, like saem require parameters in terms of Pop.Parameter + Individual.Deviation.Parameter + Covariate*Covariate.Parameter.
The order of these parameters do not matter. This is similar to NONMEM's mu-referencing, though
not quite so restrictive.
The type of parameter in the model is determined by the initial block; Covariates used in the
model are missing in the ini block. These variables need to be present in the modeling
dataset for the model to run.
Model Syntax for solved PK systems
Solved PK systems are also currently supported by nlmixr with the `linCmt()` pseudo-function. An annotated example of a solved system is below:
##'
f <- function(){
ini({
lCl <- 1.6 #log Cl (L/hr)
lVc <- log(90) #log Vc (L)
lKA <- 0.1 #log Ka (1/hr)
prop.err <- c(0, 0.2, 1)
eta.Cl ~ 0.1 ## BSV Cl
eta.Vc ~ 0.1 ## BSV Vc
eta.KA ~ 0.1 ## BSV Ka
})
model({
Cl <- exp(lCl + eta.Cl)
Vc = exp(lVc + eta.Vc)
KA <- exp(lKA + eta.KA)
## Instead of specifying the ODEs, you can use
## the linCmt() function to use the solved system.
##
## This function determines the type of PK solved system
## to use by the parameters that are defined. In this case
## it knows that this is a one-compartment model with first-order
## absorption.
linCmt() ~ prop(prop.err)
})
}A few things to keep in mind:
Currently the solved systems support either oral dosing, IV dosing or IV infusion dosing and does not allow mixing the dosing types.
While RxODE allows mixing of solved systems and ODEs, this has not been implemented in nlmixr yet.
The solved systems implemented are the one, two and three compartment models with or without first-order absorption. Each of the models support a lag time with a tlag parameter.
In general the linear compartment model figures out the model by the parameter names. nlmixr currently knows about numbered volumes, Vc/Vp, Clearances in terms of both Cl and Q/CLD. Additionally nlmixr knows about elimination micro-constants (ie K12). Mixing of these parameters for these models is currently not supported.
Checking model syntax
After specifying the model syntax you can check that nlmixr is
interpreting it correctly by using the nlmixr function on
it.
Using the above function we can get:
> nlmixr(f)
## 1-compartment model with first-order absorption in terms of Cl
## Initialization:
################################################################################
Fixed Effects ($theta):
lCl lVc lKA
1.60000 4.49981 0.10000
Omega ($omega):
[,1] [,2] [,3]
[1,] 0.1 0.0 0.0
[2,] 0.0 0.1 0.0
[3,] 0.0 0.0 0.1
## Model:
################################################################################
Cl <- exp(lCl + eta.Cl)
Vc = exp(lVc + eta.Vc)
KA <- exp(lKA + eta.KA)
## Instead of specifying the ODEs, you can use
## the linCmt() function to use the solved system.
##
## This function determines the type of PK solved system
## to use by the parameters that are defined. In this case
## it knows that this is a one-compartment model with first-order
## absorption.
linCmt() ~ prop(prop.err)
In general this gives you information about the model (what type of solved system/RxODE), initial estimates as well as the code for the model block.
Using the model syntax for estimating a model
Once the model function has been created, you can use it and a dataset to estimate the parameters for a model given a dataset.
This dataset has to have RxODE compatible events IDs. Both
Monolix and NONMEM use a different dataset description. You may
convert these datasets to RxODE-compatible datasets with the
nmDataConvert function. Note that steady state
doses are not supported by RxODE, and therefore not supported by
the conversion function.
As an example, you can use a simulated rich 1-compartment dataset.
d <- Oral_1CPT
d <- d[,names(d) != "SS"];
d <- nmDataConvert(d);
Once the data has been converted to the appropriate format, you
can use the nlmixr function to run the appropriate code.
The method to estimate the model is:
fit <- nlmixr(model.function, rxode.dataset, est="est",control=estControl(options))Currently nlme and saem are implemented. For example, to run the
above model with saem, we could have the following:
> f <- function(){
ini({
lCl <- 1.6 #log Cl (L/hr)
lVc <- log(90) #log Vc (L)
lKA <- 0.1 #log Ka (1/hr)
prop.err <- c(0, 0.2, 1)
eta.Cl ~ 0.1 ## BSV Cl
eta.Vc ~ 0.1 ## BSV Vc
eta.KA ~ 0.1 ## BSV Ka
})
model({
## First parameters are defined in terms of the initial estimates
## parameter names.
Cl <- exp(lCl + eta.Cl)
Vc = exp(lVc + eta.Vc)
KA <- exp(lKA + eta.KA)
## After the differential equations are defined
kel <- Cl / Vc;
d/dt(depot) = -KA*depot;
d/dt(centr) = KA*depot-kel*centr;
## And the concentration is then calculated
cp = centr / Vc;
## Last, nlmixr is told that the plasma concentration follows
## a proportional error (estimated by the parameter prop.err)
cp ~ prop(prop.err)
})
}
> fit.s <- nlmixr(f,d,est="saem",control=saemControl(n.burn=50,n.em=100,print=50));
Compiling RxODE differential equations...done.
c:/Rtools/mingw_64/bin/g++ -I"c:/R/R-34~1.1/include" -DNDEBUG -I"d:/Compiler/gcc-4.9.3/local330/include" -Ic:/nlmixr/inst/include -Ic:/R/R-34~1.1/library/STANHE~1/include -Ic:/R/R-34~1.1/library/Rcpp/include -Ic:/R/R-34~1.1/library/RCPPAR~1/include -Ic:/R/R-34~1.1/library/RCPPEI~1/include -Ic:/R/R-34~1.1/library/BH/include -O2 -Wall -mtune=core2 -c saem3090757b4bd1x64.cpp -o saem3090757b4bd1x64.o
In file included from c:/R/R-34~1.1/library/RCPPAR~1/include/armadillo:52:0,
from c:/R/R-34~1.1/library/RCPPAR~1/include/RcppArmadilloForward.h:46,
from c:/R/R-34~1.1/library/RCPPAR~1/include/RcppArmadillo.h:31,
from saem3090757b4bd1x64.cpp:1:
c:/R/R-34~1.1/library/RCPPAR~1/include/armadillo_bits/compiler_setup.hpp:474:96: note: #pragma message: WARNING: use of OpenMP disabled; this compiler doesn't support OpenMP 3.0+
#pragma message ("WARNING: use of OpenMP disabled; this compiler doesn't support OpenMP 3.0+")
^
c:/Rtools/mingw_64/bin/g++ -shared -s -static-libgcc -o saem3090757b4bd1x64.dll tmp.def saem3090757b4bd1x64.o c:/nlmixr/R/rx_855815def56a50f0e7a80e48811d947c_x64.dll -Lc:/R/R-34~1.1/bin/x64 -lRblas -Lc:/R/R-34~1.1/bin/x64 -lRlapack -lgfortran -lm -lquadmath -Ld:/Compiler/gcc-4.9.3/local330/lib/x64 -Ld:/Compiler/gcc-4.9.3/local330/lib -Lc:/R/R-34~1.1/bin/x64 -lR
done.
1: 1.8174 4.6328 0.0553 0.0950 0.0950 0.0950 0.6357
50: 1.3900 4.2039 0.0001 0.0679 0.0784 0.1082 0.1992
100: 1.3894 4.2054 0.0107 0.0686 0.0777 0.1111 0.1981
150: 1.3885 4.2041 0.0089 0.0683 0.0778 0.1117 0.1980
Using sympy via SnakeCharmR
## Calculate ETA-based prediction and error derivatives:
Calculate Jacobian...................done.
Calculate sensitivities.......
done.
## Calculate d(f)/d(eta)
## ...
## done
## ...
## done
The model-based sensitivities have been calculated
Calculating Table Variables...
done
The options for saem are controlled by saemControl.
You may wish to make sure the minimization is complete in the case
of saem. You can do that with traceplot which shows the
iteration history with the divided by burn-in and EM phases. In
this case, the burn in seems reasonable; you may wish to increase
the number of iterations in the EM phase of the estimation.
Overall it is probably a semi-reasonable solution.
nlmixr output objects
In addition to unifying the modeling language sent to each of the estimation routines, the outputs currently have a unified structure.
You can see the fit object by typing the object name:
> fit.s
-- nlmixr SAEM fit (ODE); OBJF calculated from FOCEi approximation -------------
OBJF AIC BIC Log-likelihood Condition Number
62337.09 62351.09 62399.01 -31168.55 82.6086
-- Time (sec; fit.s$time): -----------------------------------------------------
saem setup Likelihood Calculation covariance table
elapsed 430.25 31.64 1.19 0 3.44
-- Parameters (fit.s$par.fixed): -----------------------------------------------
Parameter Estimate SE
lCl log Cl (L/hr) 1.39 0.0240 1.73 4.01 (3.83, 4.20) 26.6
lVc log Vc (L) 4.20 0.0256 0.608 67.0 (63.7, 70.4) 28.5
lKA log Ka (1/hr) 0.00924 0.0323 349. 1.01 (0.947, 1.08) 34.3
prop.err prop.err 0.198 19.8
Shrink(SD)
lCl 0.248
lVc 1.09
lKA 4.19
prop.err 1.81
No correlations in between subject variability (BSV) matrix
Full BSV covariance (fit.s$omega) or correlation (fit.s$omega.R; diagonals=SDs)
Distribution stats (mean/skewness/kurtosis/p-value) available in fit.s$shrink
-- Fit Data (object fit.s is a modified data.frame): ---------------------------
# A tibble: 6,947 x 22
ID TIME DV PRED RES WRES IPRED IRES IWRES CPRED CRES
* <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.25 205. 198. 6.60 0.0741 189. 16.2 0.434 198. 6.78
2 1 0.5 311. 349. -38.7 -0.261 330. -19.0 -0.291 349. -38.3
3 1 0.75 389. 464. -74.5 -0.398 434. -45.2 -0.526 463. -73.9
# ... with 6,944 more rows, and 11 more variables: CWRES <dbl>, eta.Cl <dbl>,
# eta.Vc <dbl>, eta.KA <dbl>, depot <dbl>, centr <dbl>, Cl <dbl>, Vc <dbl>,
This example shows what is typical printout of a nlmixr fit object. The elements of the fit are:
Metrics of goodness of fit (AIC, BIC,
and logLik).
To align the comparison between methods, the FOCEi likelihood objective is calculated regardless of the method used and used for goodness of fit metrics.
This FOCEi likelihood has been compared to NONMEM's objective function and gives the same values (based on the data in Wang 2007)
Also note that saem does not calculate an objective function,
and the FOCEi is used as the only objective function for the fit.
Even though the objective functions are calculated in the same manner, caution should be used when comparing fits from various estimation routines.
The next item is the timing of each of the steps of the fit.
These can be also accessed by (fit.s$time).
As a mnemonic, the access for this item is shown in the printout. This is true for almost all of the other items in the printout.
After the timing of the fit, the parameter estimates are displayed (can be accessed by
fit.s$par.fixed)
While the items are rounded for R printing, each estimate without rounding is still accessible by the `$` syntax. For example, the `$Untransformed` gives the untransformed parameter values.
The Untransformed parameter takes log-space parameters and back-transforms them to normal parameters. Not the CIs are listed on the back-transformed parameter space.
Proportional Errors are converted to
Omega block (accessed by fit.s$omega)
The table of fit data. Please note:
A nlmixr fit object is actually a data frame. Saving it as a Rdata object and then loading it without nlmixr will just show the data by itself. Don't worry; the fit information has not vanished, you can bring it back by simply loading nlmixr, and then accessing the data.
Special access to fit information (like the $omega) needs nlmixr to extract the information.
If you use the $ to access information, the order of precedence is:
Fit data from the overall data.frame
Information about the parsed nlmixr model (via $uif)
Parameter history if available (via $par.hist and $par.hist.stacked)
Fixed effects table (via $par.fixed)
Individual differences from the typical population parameters (via $eta)
Fit information from the list of information generated during the post-hoc residual calculation.
Fit information from the environment where the post-hoc residual were calculated
Fit information about how the data and options interacted with the specified model (such as estimation options or if the solved system is for an infusion or an IV bolus).
While the printout may displays the data as a data.table object or tbl
object, the data is NOT any of these objects, but rather a derived data frame.
Since the object is a data.frame, you can treat it like one.
In addition to the above properties of the fit object, there are a few additional that may be helpful for the modeler:
$theta gives the fixed effects parameter estimates (in NONMEM the
thetas). This can also be accessed in fixed.effects
function. Note that the residual variability is treated as a fixed effect parameter
and is included in this list.
$eta gives the random effects parameter estimates, or in NONMEM the
etas. This can also be accessed in using the random.effects
function.
# \donttest{
f_ode <- function(){
ini({
lCl <- 1.6 #log Cl (L/hr)
lVc <- log(80) #log Vc (L)
lKA <- 0.3 #log Ka (1/hr)
prop.err <- c(0, 0.2, 1)
eta.Cl ~ 0.3 ## BSV Cl
eta.Vc ~ 0.2 ## BSV Vc
eta.KA ~ 0.1 ## BSV Ka
})
model({
## First parameters are defined in terms of the initial estimates
## parameter names.
Cl <- exp(lCl + eta.Cl)
Vc = exp(lVc + eta.Vc)
KA <- exp(lKA + eta.KA)
## After the differential equations are defined
kel <- Cl / Vc;
d/dt(depot) = -KA*depot;
d/dt(centr) = KA*depot-kel*centr;
## And the concentration is then calculated
cp = centr / Vc;
## Last, nlmixr is told that the plasma concentration follows
## a proportional error (estimated by the parameter prop.err)
cp ~ prop(prop.err)
})
}
f_linCmt <- function(){
ini({
lCl <- 1.6 #log Cl (L/hr)
lVc <- log(90) #log Vc (L)
lKA <- 0.1 #log Ka (1/hr)
prop.err <- c(0, 0.2, 1)
add.err <- c(0, 0.01)
eta.Cl ~ 0.1 ## BSV Cl
eta.Vc ~ 0.1 ## BSV Vc
eta.KA ~ 0.1 ## BSV Ka
})
model({
Cl <- exp(lCl + eta.Cl)
Vc = exp(lVc + eta.Vc)
KA <- exp(lKA + eta.KA)
## Instead of specifying the ODEs, you can use
## the linCmt() function to use the solved system.
##
## This function determines the type of PK solved system
## to use by the parameters that are defined. In this case
## it knows that this is a one-compartment model with first-order
## absorption.
linCmt() ~ add(add.err) + prop(prop.err)
})
}
# Use nlme algorithm
fit_linCmt_nlme <- try(nlmixr(f_ode, Oral_1CPT, est="nlme",
control=nlmeControl(maxstepsOde = 50000, pnlsTol=0.4)))
#> ℹ parameter labels from comments will be replaced by 'label()'
#>
#>
#> **Iteration 1
#> LME step: Loglik: -40817.93, nlminb iterations: 25
#> reStruct parameters:
#> ID1 ID2 ID3
#> 0.1221898 -0.1682207 -0.3301838
#> Beginning PNLS step: .. completed fit_nlme() step.
#> PNLS step: RSS = 640.6597
#> fixed effects: 1.385724 4.204827 0.02666615
#> iterations: 4
#> Convergence crit. (must all become <= tolerance = 1e-05):
#> fixed reStruct
#> 10.25022 102.28910
#>
#> **Iteration 2
#> LME step: Loglik: -37497.72, nlminb iterations: 13
#> reStruct parameters:
#> ID1 ID2 ID3
#> -0.2767254 -0.3391281 -0.4783564
#> Beginning PNLS step: .. completed fit_nlme() step.
#> PNLS step: RSS = 269.6296
#> fixed effects: 1.388601 4.201673 -0.006723139
#> iterations: 2
#> Convergence crit. (must all become <= tolerance = 1e-05):
#> fixed reStruct
#> 4.966324 1.571577
#>
#> **Iteration 3
#> LME step: Loglik: -37492.25, nlminb iterations: 4
#> reStruct parameters:
#> ID1 ID2 ID3
#> -0.2680315 -0.3305057 -0.4791520
#> Beginning PNLS step: .. completed fit_nlme() step.
#> PNLS step: RSS = 274.0337
#> fixed effects: 1.388601 4.201673 -0.006723139
#> iterations: 1
#> Convergence crit. (must all become <= tolerance = 1e-05):
#> fixed reStruct
#> 0.00000000 0.03338791
#>
#> **Iteration 4
#> LME step: Loglik: -37492.25, nlminb iterations: 1
#> reStruct parameters:
#> ID1 ID2 ID3
#> -0.2680392 -0.3305202 -0.4791585
#> Beginning PNLS step: .. completed fit_nlme() step.
#> PNLS step: RSS = 274.0337
#> fixed effects: 1.388601 4.201673 -0.006723139
#> iterations: 1
#> Convergence crit. (must all become <= tolerance = 1e-05):
#> fixed reStruct
#> 0 0
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#>
#> → calculate sensitivities
#>
#> → calculate ∂(f)/∂(η)
#>
#> → calculate ∂(R²)/∂(η)
#>
#> → finding duplicate expressions in inner model...
#>
#> → optimizing duplicate expressions in inner model...
#>
#> → finding duplicate expressions in EBE model...
#>
#> → optimizing duplicate expressions in EBE model...
#>
#> → compiling inner model...
#>
#> ✔ done
#> → finding duplicate expressions in FD model...
#>
#> → optimizing duplicate expressions in FD model...
#>
#> → compiling EBE model...
#>
#> ✔ done
#> → compiling events FD model...
#>
#> ✔ done
#> Calculating residuals/tables
#> done
if (!inherits(fit_linCmt_nlme, "try-error")) print(fit_linCmt_nlme)
#> ── nlmixr nlme by maximum likelihood (ODE; μ-ref & covs) fit ───────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 62237.06 75018.79 75066.72 -37502.40 1.727187
#> nlme 62216.76 74998.49 75046.42 -37492.25 1.727187
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> nlme setup table cwres other
#> elapsed 41.728 4.473892 1.164 5.81 0.721108
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> lCl log Cl (L/hr) 1.39 0.0238 1.72 4.01 (3.83, 4.2) 26.4
#> lVc log Vc (L) 4.2 0.0256 0.61 66.8 (63.5, 70.2) 28.1
#> lKA log Ka (1/hr) -0.00672 0.0313 466 0.993 (0.934, 1.06) 32.9
#> prop.err 0.198 0.198
#> Shrink(SD)%
#> lCl 0.0750%
#> lVc 0.932%
#> lKA 7.35%
#> prop.err
#>
#> Covariance Type ($covMethod): nlme
#> Fixed parameter correlations in $cor
#> 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
#>
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 6,960 × 24
#> ID EVID TIME DV PRED RES WRES IPRED IRES IWRES CPRED CRES
#> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.25 205. 196. 8.82 0.102 185. 19.8 0.541 196. 8.89
#> 2 1 0 0.5 311. 346. -35.3 -0.245 321. -10.6 -0.166 346. -35.2
#> 3 1 0 0.75 389. 460. -70.9 -0.388 421. -31.7 -0.380 460. -70.5
#> # … with 6,957 more rows, and 12 more variables: CWRES <dbl>, eta.Cl <dbl>,
#> # eta.Vc <dbl>, eta.KA <dbl>, depot <dbl>, centr <dbl>, Cl <dbl>, Vc <dbl>,
#> # KA <dbl>, kel <dbl>, tad <dbl>, dosenum <dbl>
# Use Focei algorithm
fit_linCmt_focei <- try(nlmixr(f_linCmt, Oral_1CPT, est="focei"))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#>
#> → calculate ∂(f)/∂(η)
#>
#> → calculate ∂(R²)/∂(η)
#>
#> → finding duplicate expressions in inner model...
#>
#> → optimizing duplicate expressions in inner model...
#>
#> → finding duplicate expressions in EBE model...
#>
#> → optimizing duplicate expressions in EBE model...
#>
#> → compiling inner model...
#>
#> ✔ done
#> → finding duplicate expressions in FD model...
#>
#> → optimizing duplicate expressions in FD model...
#>
#> → compiling EBE model...
#>
#> ✔ done
#> → compiling events FD model...
#>
#> ✔ done
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximation
#> F: Forward difference gradient approximation
#> C: Central difference gradient approximation
#> M: Mixed forward and central difference gradient approximation
#> Unscaled parameters for Omegas=chol(solve(omega));
#> Diagonals are transformed, as specified by foceiControl(diagXform=)
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> | #| Objective Fun | lCl | lVc | lKA | prop.err |
#> |.....................| add.err | o1 | o2 | o3 |
#> | 1| 62378.086 | -0.2917 | 1.000 | -0.9599 | -0.9154 |
#> |.....................| -1.000 | -0.2123 | -0.2123 | -0.2123 |
#> | U| 62378.086 | 1.600 | 4.500 | 0.1000 | 0.2000 |
#> |.....................| 0.01000 | 1.778 | 1.778 | 1.778 |
#> | X| 62378.086 | 4.953 | 90.00 | 1.105 | 0.2000 |
#> |.....................| 0.01000 | 1.778 | 1.778 | 1.778 |
#> | G| Gill Diff. | 487.4 | 681.2 | 187.2 | 70.06 |
#> |.....................| -0.5594 | 46.70 | 288.5 | 69.87 |
#> | 2| 62504.532 | -0.8261 | 0.2531 | -1.165 | -0.9922 |
#> |.....................| -0.9994 | -0.2635 | -0.5286 | -0.2889 |
#> | U| 62504.532 | 1.066 | 3.753 | -0.1053 | 0.1923 |
#> |.....................| 0.01000 | 1.687 | 1.216 | 1.642 |
#> | X| 62504.532 | 2.903 | 42.65 | 0.9001 | 0.1923 |
#> |.....................| 0.01000 | 1.687 | 1.216 | 1.642 |
#> | 3| 62250.138 | -0.5264 | 0.6720 | -1.050 | -0.9491 |
#> |.....................| -0.9997 | -0.2348 | -0.3512 | -0.2460 |
#> | U| 62250.138 | 1.365 | 4.172 | 0.009862 | 0.1966 |
#> |.....................| 0.01000 | 1.738 | 1.531 | 1.718 |
#> | X| 62250.138 | 3.917 | 64.83 | 1.010 | 0.1966 |
#> |.....................| 0.01000 | 1.738 | 1.531 | 1.718 |
#> | F| Forward Diff. | -51.16 | -40.52 | 27.63 | -82.65 |
#> |.....................| -0.6161 | -181.1 | -306.6 | -11.29 |
#> | 4| 62272.286 | -0.5342 | 0.6210 | -1.114 | -0.8526 |
#> |.....................| -0.9988 | -0.007246 | 0.002245 | -0.2419 |
#> | U| 62272.286 | 1.358 | 4.121 | -0.05440 | 0.2063 |
#> |.....................| 0.01001 | 2.143 | 2.160 | 1.726 |
#> | X| 62272.286 | 3.887 | 61.61 | 0.9470 | 0.2063 |
#> |.....................| 0.01001 | 2.143 | 2.160 | 1.726 |
#> | 5| 62209.348 | -0.5000 | 0.6929 | -1.064 | -0.9064 |
#> |.....................| -0.9994 | -0.1413 | -0.1929 | -0.2401 |
#> | U| 62209.348 | 1.392 | 4.193 | -0.004409 | 0.2009 |
#> |.....................| 0.01000 | 1.905 | 1.813 | 1.729 |
#> | X| 62209.348 | 4.022 | 66.21 | 0.9956 | 0.2009 |
#> |.....................| 0.01000 | 1.905 | 1.813 | 1.729 |
#> | F| Forward Diff. | 4.259 | -27.99 | 2.216 | 176.2 |
#> |.....................| -0.5135 | -47.31 | -69.76 | -4.942 |
#> | 6| 62278.664 | -0.5041 | 0.7202 | -1.066 | -1.078 |
#> |.....................| -0.9989 | -0.09516 | -0.1249 | -0.2353 |
#> | U| 62278.664 | 1.388 | 4.220 | -0.006568 | 0.1837 |
#> |.....................| 0.01001 | 1.987 | 1.934 | 1.737 |
#> | X| 62278.664 | 4.005 | 68.04 | 0.9935 | 0.1837 |
#> |.....................| 0.01001 | 1.987 | 1.934 | 1.737 |
#> | 7| 62205.639 | -0.5007 | 0.6978 | -1.065 | -0.9368 |
#> |.....................| -0.9993 | -0.1331 | -0.1808 | -0.2393 |
#> | U| 62205.639 | 1.391 | 4.198 | -0.004791 | 0.1979 |
#> |.....................| 0.01000 | 1.919 | 1.834 | 1.730 |
#> | X| 62205.639 | 4.019 | 66.53 | 0.9952 | 0.1979 |
#> |.....................| 0.01000 | 1.919 | 1.834 | 1.730 |
#> | F| Forward Diff. | 5.482 | -13.33 | 1.146 | -8.050 |
#> |.....................| -0.5783 | -34.90 | -51.22 | -2.300 |
#> | 8| 62204.377 | -0.5036 | 0.7048 | -1.065 | -0.9325 |
#> |.....................| -0.9990 | -0.1146 | -0.1537 | -0.2381 |
#> | U| 62204.377 | 1.388 | 4.205 | -0.005399 | 0.1983 |
#> |.....................| 0.01000 | 1.952 | 1.883 | 1.733 |
#> | X| 62204.377 | 4.007 | 67.00 | 0.9946 | 0.1983 |
#> |.....................| 0.01000 | 1.952 | 1.883 | 1.733 |
#> | F| Forward Diff. | -5.005 | 6.076 | 0.02355 | 17.86 |
#> |.....................| -0.5667 | -6.648 | -8.466 | -0.6091 |
#> | 9| 62205.622 | -0.4836 | 0.6824 | -1.063 | -0.9374 |
#> |.....................| -0.9972 | -0.1046 | -0.1427 | -0.2368 |
#> | U| 62205.622 | 1.408 | 4.182 | -0.002890 | 0.1978 |
#> |.....................| 0.01001 | 1.970 | 1.902 | 1.735 |
#> | X| 62205.622 | 4.088 | 65.51 | 0.9971 | 0.1978 |
#> |.....................| 0.01001 | 1.970 | 1.902 | 1.735 |
#> | 10| 62204.331 | -0.5025 | 0.7036 | -1.065 | -0.9363 |
#> |.....................| -0.9989 | -0.1132 | -0.1518 | -0.2379 |
#> | U| 62204.331 | 1.389 | 4.203 | -0.005404 | 0.1979 |
#> |.....................| 0.01001 | 1.955 | 1.886 | 1.733 |
#> | X| 62204.331 | 4.012 | 66.91 | 0.9946 | 0.1979 |
#> |.....................| 0.01001 | 1.955 | 1.886 | 1.733 |
#> | F| Forward Diff. | -0.7833 | 2.682 | 0.08084 | -6.082 |
#> |.....................| -0.5753 | -4.469 | -5.489 | -0.1852 |
#> | 11| 62204.314 | -0.5020 | 0.7021 | -1.065 | -0.9334 |
#> |.....................| -0.9986 | -0.1110 | -0.1492 | -0.2378 |
#> | U| 62204.314 | 1.390 | 4.202 | -0.005419 | 0.1982 |
#> |.....................| 0.01001 | 1.958 | 1.891 | 1.733 |
#> | X| 62204.314 | 4.014 | 66.81 | 0.9946 | 0.1982 |
#> |.....................| 0.01001 | 1.958 | 1.891 | 1.733 |
#> | F| Forward Diff. | 0.6437 | -1.884 | 0.2261 | 12.41 |
#> |.....................| -0.5679 | -1.207 | -1.444 | -0.1867 |
#> | 12| 62204.318 | -0.5039 | 0.7017 | -1.066 | -0.9356 |
#> |.....................| -0.9953 | -0.1092 | -0.1499 | -0.2378 |
#> | U| 62204.318 | 1.388 | 4.202 | -0.005695 | 0.1980 |
#> |.....................| 0.01002 | 1.962 | 1.889 | 1.733 |
#> | X| 62204.318 | 4.006 | 66.79 | 0.9943 | 0.1980 |
#> |.....................| 0.01002 | 1.962 | 1.889 | 1.733 |
#> | 13| 62204.314 | -0.5020 | 0.7021 | -1.065 | -0.9334 |
#> |.....................| -0.9986 | -0.1110 | -0.1492 | -0.2378 |
#> | U| 62204.314 | 1.390 | 4.202 | -0.005419 | 0.1982 |
#> |.....................| 0.01001 | 1.958 | 1.891 | 1.733 |
#> | X| 62204.314 | 4.014 | 66.81 | 0.9946 | 0.1982 |
#> |.....................| 0.01001 | 1.958 | 1.891 | 1.733 |
#> calculating covariance matrix
#> done
#> Calculating residuals/tables
#> done
#> Warning: initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Warning: ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> Warning: gradient problems with initial estimate and covariance; see $scaleInfo
if (!inherits(fit_linCmt_focei, "try-error")) print(fit_linCmt_focei)
#> ── nlmixr FOCEi (outer: nlminb) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 62204.31 74988.05 75042.81 -37486.02 1.826486
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 5.212046 6.688262 6.688269 0.798 7.979423
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> lCl log Cl (L/hr) 1.39 0.0242 1.74 4.01 (3.83, 4.21) 26.5
#> lVc log Vc (L) 4.2 0.0258 0.615 66.8 (63.5, 70.3) 28.5
#> lKA log Ka (1/hr) -0.00542 0.0327 603 0.995 (0.933, 1.06) 34.2
#> prop.err 0.198 0.198
#> add.err 0.01 0.01
#> Shrink(SD)%
#> lCl 0.396%
#> lVc 1.42%
#> lKA 6.27%
#> prop.err
#> add.err
#>
#> Covariance Type ($covMethod): r,s
#> Fixed parameter correlations in $cor
#> 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
#> Minimization message ($message):
#> relative convergence (4)
#>
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 6,960 × 21
#> ID EVID TIME DV PRED RES WRES IPRED IRES IWRES CPRED CRES
#> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.25 205. 196. 8.64 0.0980 184. 20.6 0.566 196. 8.86
#> 2 1 0 0.5 311. 346. -35.6 -0.242 321. -10.2 -0.161 346. -35.1
#> 3 1 0 0.75 389. 460. -71.1 -0.382 422. -32.4 -0.387 460. -70.3
#> # … with 6,957 more rows, and 9 more variables: CWRES <dbl>, eta.Cl <dbl>,
#> # eta.Vc <dbl>, eta.KA <dbl>, Cl <dbl>, Vc <dbl>, KA <dbl>, tad <dbl>,
#> # dosenum <dbl>
# The ODE model can be fitted using the saem algorithm, more
# iterations should be used for real applications
fit_ode_saem <- try(nlmixr(f_ode, Oral_1CPT, est = "saem",
control = saemControl(n.burn = 50, n.em = 100, print = 50)))
#> ℹ parameter labels from comments will be replaced by 'label()'
#>
#> → generate SAEM model
#> ✔ done
#> 1: 1.8298 4.6775 0.1602 0.2850 0.1900 0.0950 0.6393
#> 50: 1.3912 4.2090 0.0133 0.0695 0.0794 0.1075 0.1982
#> 100: 1.3877 4.2020 0.0017 0.0681 0.0790 0.1126 0.1978
#> 150: 1.3877 4.2022 0.0023 0.0680 0.0787 0.1124 0.1978
#> Calculating covariance matrix
#>
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → compiling EBE model...
#>
#> ✔ done
#> Calculating residuals/tables
#> done
if (!inherits(fit_ode_saem, "try-error")) print(fit_ode_saem)
#> ── nlmixr SAEM(ODE); OBJF not calculated fit ───────────────────────────────────
#>
#> Gaussian/Laplacian Likelihoods: AIC() or $objf etc.
#> FOCEi CWRES & Likelihoods: addCwres()
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> saem setup table covariance other
#> elapsed 339.898 1.667443 0.284 0.768 1.001557
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> lCl log Cl (L/hr) 1.39 0.024 1.73 4.01 (3.82, 4.2) 26.5
#> lVc log Vc (L) 4.2 0.026 0.619 66.8 (63.5, 70.3) 28.6
#> lKA log Ka (1/hr) 0.0023 0.0325 1.41e+03 1 (0.94, 1.07) 34.5
#> prop.err 0.198 0.198
#> Shrink(SD)%
#> lCl -0.393%
#> lVc 1.29%
#> lKA 4.36%
#> prop.err
#>
#> Covariance Type ($covMethod): linFim
#> Fixed parameter correlations in $cor
#> 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
#>
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 6,960 × 21
#> ID EVID TIME DV PRED RES IPRED IRES IWRES eta.Cl eta.Vc eta.KA
#> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.25 205. 197. 7.36 173. 31.5 0.919 0.0862 0.141 0.0121
#> 2 1 0 0.5 311. 348. -37.6 305. 5.25 0.0870 0.0862 0.141 0.0121
#> 3 1 0 0.75 389. 463. -73.4 405. -16.2 -0.202 0.0862 0.141 0.0121
#> # … with 6,957 more rows, and 9 more variables: cp <dbl>, depot <dbl>,
#> # centr <dbl>, Cl <dbl>, Vc <dbl>, KA <dbl>, kel <dbl>, tad <dbl>,
#> # dosenum <dbl>
# }