nlmixr
Changing models via piping
As in the running nlmixr vignette, Let’s start with a very simple PK example, using the single-dose theophylline dataset generously provided by Dr. Robert A. Upton of the University of California, San Francisco:
library(nlmixr2)
one.compartment <- function() {
ini({
tka <- 0.45; label("Ka")
tcl <- 1; label("Cl")
tv <- 3.45; label("V")
eta.ka ~ 0.6
eta.cl ~ 0.3
eta.v ~ 0.1
add.sd <- 0.7
})
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.sd)
})
}We can try the First-Order Conditional Estimation with Interaction (FOCEi) method to find a good solution:
fit <- nlmixr(one.compartment, theo_sd, est="focei",
control=list(print=0),
table=list(npde=TRUE, cwres=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(fit)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 116.8116 373.4114 393.591 -179.7057 115.6797 9.731882
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.020364 0.498618 0.498619 0.838 0.012 5.673399
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tka Ka 0.474 0.199 42.1 1.61 (1.09, 2.38) 69.2
#> tcl Cl 1.01 0.0742 7.34 2.75 (2.38, 3.18) 26.8
#> tv V 3.46 0.0337 0.973 31.8 (29.8, 34) 14.0
#> add.sd 0.695 0.695
#> Shrink(SD)%
#> tka 0.659%
#> tcl 4.19%
#> tv 10.6%
#> add.sd
#>
#> Covariance Type ($covMethod): r,s
#> 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
#> Information about run found ($runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> relative convergence (4)
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 28
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.109 0.631 0.440 0.878 0.67 0.81 0 0.74 1.06
#> 2 1 0.25 2.84 3.64 -0.798 -0.505 -0.394 0.307 0.347 3.29 -0.451 -0.242
#> 3 1 0.57 6.57 5.95 0.617 -1.75 0.332 0.04 0.63 5.87 0.702 0.284
#> # ℹ 129 more rows
#> # ℹ 16 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>,
#> # depot <dbl>, center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> # dosenum <dbl>Changing and fixing parameter values in models
Something that you may want to do is change initial estimates with a model. It is simple to modify the model definition and change them yourself, but you may also want to change them in a specific way; For example try a range of starting values to see how the system behaves (either by full estimation or by a posthoc estimation). In these situations it can be come tedious to modify the models by hand.
nlmixr provides the ability to:
- Change parameter estimates before or after running a model. (ie
ini(tka=0.5)) - Fix parameters to arbitrary values, or estimated values (ie
ini(tka=fix(0.5))orini(tka=fix))
The easiest way to illustrate this is by showing a few examples of piping changes to the model:
## Example 1 -- Set inital estimate to 0.5 (shown w/posthoc)
one.ka.0.5 <- fit %>%
ini(tka=0.5) %>%
nlmixr(est="posthoc", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
print(one.ka.0.5)
## Example 2 -- Fix tka to 0.5 and re-estimate.
one.ka.0.5 <- fit %>%
ini(tka=fix(0.5)) %>%
nlmixr(est="focei", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(one.ka.0.5)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 116.8471 371.4469 388.7437 -179.7234 17.52287 12.09477
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.002326 0.290734 0.290736 0.799 0.012 3.011204
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Ka 0.5 FIXED FIXED 0.5 69.2 0.600%
#> tcl Cl 1.01 0.312 30.9 2.75 (1.49, 5.07) 26.8 4.21%
#> tv V 3.46 0.174 5.04 31.8 (22.6, 44.7) 14.0 10.6%
#> add.sd 0.695 0.695
#>
#> Covariance Type ($covMethod): s
#> 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
#> Information about run found ($runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • using S matrix to calculate covariance, can check sandwich or R matrix with $covRS and $covR
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> false convergence (8)
#> In an ODE system, false convergence may mean "useless" evaluations were performed.
#> See https://tinyurl.com/yyrrwkce
#> It could also mean the convergence is poor, check results before accepting fit
#> You may also try a good derivative free optimization:
#> nlmixr2(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 28
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.109 0.631 0.440 0.878 0.67 0.81 0 0.74 1.06
#> 2 1 0.25 2.84 3.71 -0.867 -0.486 -0.422 0.313 0.337 3.36 -0.520 -0.275
#> 3 1 0.57 6.57 6.03 0.539 -1.75 0.253 0.04 0.6 5.96 0.610 0.246
#> # ℹ 129 more rows
#> # ℹ 16 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>,
#> # depot <dbl>, center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> # dosenum <dbl>
## Example 3 -- Fix tka to model estimated value and re-estimate.
one.ka.0.5 <- fit %>%
ini(tka=fix) %>%
nlmixr(est="focei", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(one.ka.0.5)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 116.8131 371.4128 388.7096 -179.7064 17.54673 12.12971
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.002228 0.299125 0.299127 0.867 0.022 3.02552
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Ka 0.474 FIXED FIXED 0.474 69.2 0.670%
#> tcl Cl 1.01 0.313 30.9 2.75 (1.49, 5.08) 26.8 4.19%
#> tv V 3.46 0.175 5.05 31.8 (22.6, 44.8) 14.0 10.6%
#> add.sd 0.695 0.695
#>
#> Covariance Type ($covMethod): s
#> 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
#> Information about run found ($runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • using S matrix to calculate covariance, can check sandwich or R matrix with $covRS and $covR
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> false convergence (8)
#> In an ODE system, false convergence may mean "useless" evaluations were performed.
#> See https://tinyurl.com/yyrrwkce
#> It could also mean the convergence is poor, check results before accepting fit
#> You may also try a good derivative free optimization:
#> nlmixr2(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 28
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.109 0.631 0.440 0.878 0.67 0.81 0 0.74 1.06
#> 2 1 0.25 2.84 3.64 -0.798 -0.505 -0.394 0.307 0.347 3.29 -0.451 -0.242
#> 3 1 0.57 6.57 5.95 0.617 -1.75 0.332 0.04 0.63 5.87 0.702 0.284
#> # ℹ 129 more rows
#> # ℹ 16 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>,
#> # depot <dbl>, center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> # dosenum <dbl>
## Example 4 -- Change tka to 0.7 in orginal model function and then estimate
one.ka.0.7 <- one.compartment %>%
ini(tka=0.7) %>%
nlmixr(theo_sd, est="focei", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(one.ka.0.7)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 116.9872 373.587 393.7666 -179.7935 77.33962 11.86369
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.001873 0.50367 0.503672 0.393 0.01 2.441785
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Ka 0.481 0.213 44.2 1.62 (1.07, 2.46) 70.5 1.61%
#> tcl Cl 1.02 0.0737 7.22 2.77 (2.4, 3.2) 28.0 7.96%
#> tv V 3.46 0.0516 1.49 31.8 (28.8, 35.2) 15.3 15.1%
#> add.sd 0.693 0.693
#>
#> Covariance Type ($covMethod): r,s
#> 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
#> Information about run found ($runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> false convergence (8)
#> In an ODE system, false convergence may mean "useless" evaluations were performed.
#> See https://tinyurl.com/yyrrwkce
#> It could also mean the convergence is poor, check results before accepting fit
#> You may also try a good derivative free optimization:
#> nlmixr2(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 28
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.109 0.631 0.422 0.890 0.663 0.813 0 0.74 1.07
#> 2 1 0.25 2.84 3.67 -0.826 -0.477 -0.394 0.317 0.347 3.31 -0.468 -0.246
#> 3 1 0.57 6.57 5.97 0.595 -1.79 0.297 0.0367 0.617 5.89 0.680 0.270
#> # ℹ 129 more rows
#> # ℹ 16 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>,
#> # depot <dbl>, center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> # dosenum <dbl>Changing parameter labels and order
For aesthetic reasons, sometimes it is preferred to update parameter labels and the order of parameters. These changes do not affect the estimation of the parameters. They only affect the output tables and order of parameters.
By using these, you can modify a model with model piping and still have the desired output table format ready to use in a report.
For example, you can change the label from "Ka" to
"Absorption rate" as follows:
fit %>%
ini(
tka <- label("Absorption rate")
)
#> ── rxode2-based free-form 2-cmt ODE model ──────────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tka tcl tv add.sd
#> 0.4741772 1.0112664 3.4593293 0.6952432
#>
#> Omega ($omega):
#> eta.ka eta.cl eta.v
#> eta.ka 0.3913202 0.0000000 0.00000000
#> eta.cl 0.0000000 0.0695591 0.00000000
#> eta.v 0.0000000 0.0000000 0.01935147
#>
#> States ($state or $stateDf):
#> Compartment Number Compartment Name
#> 1 1 depot
#> 2 2 center
#> ── μ-referencing ($muRefTable): ──
#> theta eta level
#> 1 tka eta.ka id
#> 2 tcl eta.cl id
#> 3 tv eta.v id
#>
#> ── Model (Normalized Syntax): ──
#> function() {
#> ini({
#> tka <- 0.474177171359262
#> label("Absorption rate")
#> tcl <- 1.01126635373108
#> label("Cl")
#> tv <- 3.45932933156061
#> label("V")
#> add.sd <- c(0, 0.695243197595218)
#> eta.ka ~ 0.391320225986854
#> eta.cl ~ 0.0695590964689484
#> eta.v ~ 0.019351468331187
#> })
#> 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.sd)
#> })
#> }And, if you’d prefer for volume to come before clearance in the
parameter table (fit$parFixed), you can change that,
too.
fit %>%
ini(
tv <- label("Central volume"),
append = "tcl"
)
#> ── rxode2-based free-form 2-cmt ODE model ──────────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tka tcl tv add.sd
#> 0.4741772 1.0112664 3.4593293 0.6952432
#>
#> Omega ($omega):
#> eta.ka eta.cl eta.v
#> eta.ka 0.3913202 0.0000000 0.00000000
#> eta.cl 0.0000000 0.0695591 0.00000000
#> eta.v 0.0000000 0.0000000 0.01935147
#>
#> States ($state or $stateDf):
#> Compartment Number Compartment Name
#> 1 1 depot
#> 2 2 center
#> ── μ-referencing ($muRefTable): ──
#> theta eta level
#> 1 tka eta.ka id
#> 2 tcl eta.cl id
#> 3 tv eta.v id
#>
#> ── Model (Normalized Syntax): ──
#> function() {
#> ini({
#> tka <- 0.474177171359262
#> label("Ka")
#> tcl <- 1.01126635373108
#> label("Cl")
#> tv <- 3.45932933156061
#> label("Central volume")
#> add.sd <- c(0, 0.695243197595218)
#> eta.ka ~ 0.391320225986854
#> eta.cl ~ 0.0695590964689484
#> eta.v ~ 0.019351468331187
#> })
#> 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.sd)
#> })
#> }See the documentation for ini
for more about how you can modify parameters with model piping.
Changing model features
When developing models, often you add and remove between subject
variability to parameters, add covariates to the effects, and/or change
the residual errors. You can change lines in the model by piping the fit
or the nlmixr model specification function to a model
Adding or Removing between subject variability
Often in developing a model you add and remove between subject variability to certain model parameters. For example, you could remove the between subject variability in the ka parameter by changing that line in the model;
For example to remove a eta from a prior fit or prior model
specification function, simply pipe it to the model function. You can
then re-estimate by piping it to the nlmixr function
again.
## Remove eta.ka on ka
noEta <- fit %>%
model(ka <- exp(tka)) %>%
nlmixr(est="focei", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(noEta)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 176.5777 431.1775 448.4743 -209.5887 33.12954 7.050454
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.002348 0.313753 0.313754 0.838 0.011 3.935145
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Ka 0.433 0.167 38.6 1.54 (1.11, 2.14)
#> tcl Cl 0.988 0.076 7.7 2.69 (2.31, 3.12) 30.4 7.95%
#> tv V 3.48 0.0481 1.38 32.4 (29.5, 35.6) 15.4 7.24%
#> add.sd 1.02 1.02
#>
#> Covariance Type ($covMethod): r,s
#> 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
#> Information about run found ($runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> relative convergence (4)
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 27
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.160 0.580 0.0753 0.486 0.53 0.687 0 0.74 0.724
#> 2 1 0.25 2.84 3.18 -0.339 -1.38 -0.245 0.0833 0.403 3.12 -0.282 -0.250
#> 3 1 0.57 6.57 5.68 0.892 -0.515 0.674 0.303 0.75 5.62 0.954 0.722
#> # ℹ 129 more rows
#> # ℹ 15 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.cl <dbl>, eta.v <dbl>, depot <dbl>,
#> # center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>, dosenum <dbl>Of course you could also add an eta on a parameter in the same way;
addBackKa <- noEta %>%
model({ka <- exp(tka + bsv.ka)}) %>%
ini(bsv.ka=0.1) %>%
nlmixr(est="focei", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(addBackKa)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 116.8978 373.4976 393.6772 -179.7488 50.95108 7.394634
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.002214 0.501822 0.501823 0.839 0.011 5.620141
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Ka 0.462 0.207 44.9 1.59 (1.06, 2.38) 69.1 0.451%
#> tcl Cl 1.01 0.0715 7.08 2.74 (2.39, 3.16) 27.4 6.47%
#> tv V 3.46 0.0505 1.46 31.9 (28.9, 35.2) 15.0 13.8%
#> add.sd 0.693 0.693
#>
#> Covariance Type ($covMethod): r,s
#> 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
#> Information about run found ($runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> false convergence (8)
#> In an ODE system, false convergence may mean "useless" evaluations were performed.
#> See https://tinyurl.com/yyrrwkce
#> It could also mean the convergence is poor, check results before accepting fit
#> You may also try a good derivative free optimization:
#> nlmixr2(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 28
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.109 0.631 0.431 0.890 0.667 0.813 0 0.74 1.07
#> 2 1 0.25 2.84 3.67 -0.826 -0.903 -0.271 0.183 0.393 3.25 -0.411 -0.222
#> 3 1 0.57 6.57 5.96 0.608 -1.41 0.314 0.08 0.623 5.81 0.758 0.306
#> # ℹ 129 more rows
#> # ℹ 16 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.cl <dbl>, eta.v <dbl>, bsv.ka <dbl>,
#> # depot <dbl>, center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> # dosenum <dbl>You can see the name change by examining the omega
matrix:
addBackKa$omega
#> eta.cl eta.v bsv.ka
#> eta.cl 0.07240395 0.00000000 0.0000000
#> eta.v 0.00000000 0.02217241 0.0000000
#> bsv.ka 0.00000000 0.00000000 0.3907193Note that new between subject variability parameters are distinguished from other types of parameters (ie population parameters, and individual covariates) by their name. Parameters starting or ending with the following names are assumed to be between subject variability parameters:
- eta (from NONMEM convention)
- ppv (per patient variability)
- psv (per subject variability)
- iiv (inter-individual variability)
- bsv (between subject variability)
- bpv (between patient variability)
Adding Covariate effects
## Note currently cov is needed as a prefix so nlmixr knows this is an
## estimated parameter not a parameter
wt70 <- fit %>%
model({cl <- exp(tcl + eta.cl)*(WT/70)^covWtPow}) %>%
ini(covWtPow=fix(0.75)) %>%
ini(tka=fix(0.5)) %>%
nlmixr(est="focei", control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(wt70)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 116.2103 370.8101 388.1069 -179.4051 3.142512 2.615431
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.002386 0.310883 0.310885 0.811 0.011 3.585846
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tka Ka 0.5 FIXED FIXED 0.5 69.4
#> tcl Cl 1.02 0.0777 7.62 2.77 (2.38, 3.23) 26.8
#> tv V 3.46 0.0576 1.66 31.8 (28.4, 35.6) 14.0
#> add.sd 0.695 0.695
#> covWtPow 0.75 FIXED FIXED 0.75
#> Shrink(SD)%
#> tka 1.26%
#> tcl 6.79%
#> tv 12.6%
#> add.sd
#> covWtPow
#>
#> Covariance Type ($covMethod): r,s
#> 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
#> Information about run found ($runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> false convergence (8)
#> In an ODE system, false convergence may mean "useless" evaluations were performed.
#> See https://tinyurl.com/yyrrwkce
#> It could also mean the convergence is poor, check results before accepting fit
#> You may also try a good derivative free optimization:
#> nlmixr2(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 29
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.109 0.631 0.394 0.878 0.653 0.81 0 0.74 1.06
#> 2 1 0.25 2.84 3.71 -0.867 -0.431 -0.422 0.333 0.337 3.36 -0.520 -0.274
#> 3 1 0.57 6.57 6.02 0.553 -1.79 0.271 0.0367 0.607 5.95 0.622 0.252
#> # ℹ 129 more rows
#> # ℹ 17 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>,
#> # depot <dbl>, center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> # dosenum <dbl>, WT <dbl>Changing residual errors
Changing the residual errors is also just as easy, by simply specifying the error you wish to change:
## Since there are 0 predictions in the data, these are changed to
## 0.0150 to show proportional error change.
d <- theo_sd
d$DV[d$EVID == 0 & d$DV == 0] <- 0.0150
addPropModel <- fit %>%
model({cp ~ add(add.err)+prop(prop.err)}) %>%
ini(prop.err=0.1) %>%
nlmixr(d,est="focei",
control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(addPropModel)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 104.9613 363.561 386.6235 -173.7805 29.38653 16.08647
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.002125 0.521314 0.521316 0.858 0.011 8.315245
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tka Ka 0.424 0.326 76.9 1.53 (0.806, 2.9) 72.2
#> tcl Cl 1.02 0.111 10.9 2.77 (2.23, 3.44) 26.4
#> tv V 3.46 0.176 5.08 31.8 (22.6, 44.9) 13.4
#> add.err 0.315 0.315
#> prop.err 0.12 0.12
#> Shrink(SD)%
#> tka 3.73%
#> tcl 1.93%
#> tv 14.3%
#> add.err
#> prop.err
#>
#> Covariance Type ($covMethod): r,s
#> 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
#> Information about run found ($runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> relative convergence (4)
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 28
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.0494 0.691 0.750 2.22 0.773 0.987 0 0.74 2.35
#> 2 1 0.25 2.84 3.50 -0.660 -0.603 -0.305 0.273 0.38 3.16 -0.317 -0.177
#> 3 1 0.57 6.57 5.80 0.765 -0.486 0.358 0.313 0.64 5.69 0.885 0.350
#> # ℹ 129 more rows
#> # ℹ 16 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>,
#> # depot <dbl>, center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> # dosenum <dbl>There is much more you can do with piping. For a more complete
discussion see see
rxode2 piping documentation. Since rxode2 and
nlmixr2 models can share the same functional form the
piping applies to fits as well as model definitions.