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(nlmixr)
one.compartment <- function() {
ini({
tka <- 0.45 # Log Ka
tcl <- 1 # Log Cl
## This works with interactive models
## You may also label the preceding line with label("label text")
tv <- 3.45; # log V
## the label("Label name") works with all models
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) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 116.8139 373.4137 393.5933 -179.7068 73.85496
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 5.489827 0.47566 0.475669 0.737 2.591844
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Log Ka 0.474 0.209 44.1 1.61 (1.07, 2.42) 68.9 0.384%
#> tcl Log Cl 1.01 0.0943 9.32 2.75 (2.29, 3.31) 26.8 3.87%
#> tv log V 3.46 0.0403 1.16 31.8 (29.4, 34.4) 13.9 10.3%
#> add.sd 0.696 0.696
#>
#> Covariance Type ($covMethod): r,s
#> Some strong fixed parameter correlations exist ($cor) :
#> cor:tcl,tka cor:tv,tka cor:tv,tcl
#> 0.131 0.396 0.710
#>
#>
#> 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):
#> 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:
#> nlmixr(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 132 × 26
#> ID TIME DV EPRED ERES NPDE NPD PRED RES WRES IPRED IRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.0967 0.643 1.24 1.24 0 0.74 1.06 0 0.74
#> 2 1 0.25 2.84 3.83 -0.986 -0.486 -0.486 3.29 -0.449 -0.242 3.85 -1.01
#> 3 1 0.57 6.57 6.15 0.422 -1.71 -1.71 5.87 0.705 0.287 6.79 -0.216
#> # … with 129 more rows, and 14 more variables: 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) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 116.8519 371.4517 388.7485 -179.7259 10.40787
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 4.365104 0.260299 0.260308 0.384 1.426289
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tka Log Ka 0.5 FIXED FIXED 1.65 72.8
#> tcl Log Cl 1.01 0.089 8.82 2.74 (2.3, 3.27) 26.8
#> tv log V 3.46 0.0488 1.41 31.8 (28.9, 35) 13.9
#> add.sd 0.694 0.694
#> Shrink(SD)%
#> tka 3.82%
#> tcl 4.04%
#> tv 10.1%
#> add.sd
#>
#> Covariance Type ($covMethod): r,s
#> Some strong fixed parameter correlations exist ($cor) :
#> cor:tv,tcl
#> 0.741
#>
#>
#> 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: 132 × 26
#> ID TIME DV EPRED ERES NPDE NPD PRED RES WRES IPRED IRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.0964 0.644 1.26 1.26 0 0.74 1.07 0 0.74
#> 2 1 0.25 2.84 3.93 -1.09 -0.341 -0.341 3.36 -0.520 -0.266 3.85 -1.01
#> 3 1 0.57 6.57 6.23 0.337 -1.79 -1.79 5.96 0.611 0.239 6.79 -0.220
#> # … with 129 more rows, and 14 more variables: 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))
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
print(one.ka.0.5)
#> ── nlmixr FOCEi (outer: nlminb) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 116.8095 371.4093 388.7061 -179.7047 7.965388
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 0.39741 0.27486 0.274868 0.375 1.967862
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tka Log Ka 0.45 FIXED FIXED 1.57 70.0
#> tcl Log Cl 1.01 0.0679 6.7 2.75 (2.41, 3.14) 26.8
#> tv log V 3.46 0.0523 1.51 31.8 (28.7, 35.2) 13.9
#> add.sd 0.695 0.695
#> Shrink(SD)%
#> tka 1.44%
#> tcl 4.02%
#> tv 10.3%
#> add.sd
#>
#> Covariance Type ($covMethod): r,s
#> Some strong fixed parameter correlations exist ($cor) :
#> cor:tv,tcl
#> 0.759
#>
#>
#> 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: 132 × 26
#> ID TIME DV EPRED ERES NPDE NPD PRED RES WRES IPRED IRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.0965 0.643 1.24 1.24 0 0.74 1.06 0 0.74
#> 2 1 0.25 2.84 3.77 -0.932 -0.496 -0.496 3.23 -0.387 -0.209 3.84 -1.00
#> 3 1 0.57 6.57 6.08 0.491 -1.71 -1.71 5.78 0.788 0.318 6.78 -0.213
#> # … with 129 more rows, and 14 more variables: 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) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 116.8133 373.4131 393.5927 -179.7065 49.0338
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 0.347456 0.481256 0.481264 0.373 1.604024
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Log Ka 0.471 0.202 42.8 1.6 (1.08, 2.38) 69.5 0.996%
#> tcl Log Cl 1.01 0.0757 7.49 2.75 (2.37, 3.19) 26.4 3.12%
#> tv log V 3.46 0.0543 1.57 31.9 (28.6, 35.4) 13.8 9.93%
#> add.sd 0.697 0.697
#>
#> Covariance Type ($covMethod): r,s
#> Some strong fixed parameter correlations exist ($cor) :
#> cor:tcl,tka cor:tv,tka cor:tv,tcl
#> 0.138 0.413 0.728
#>
#>
#> 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):
#> 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:
#> nlmixr(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 132 × 26
#> ID TIME DV EPRED ERES NPDE NPD PRED RES WRES IPRED IRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.0968 0.643 1.24 1.24 0 0.74 1.06 0 0.74
#> 2 1 0.25 2.84 3.82 -0.979 -0.468 -0.468 3.28 -0.438 -0.235 3.85 -1.01
#> 3 1 0.57 6.57 6.13 0.437 -1.75 -1.75 5.85 0.722 0.293 6.79 -0.215
#> # … with 129 more rows, and 14 more variables: 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 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) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 176.5764 431.1762 448.473 -209.5881 32.76372
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 4.1869 0.29998 0.299989 0.732 1.291131
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Log Ka 0.432 0.165 38.2 1.54 (1.11, 2.13)
#> tcl Log Cl 0.991 0.0762 7.69 2.69 (2.32, 3.13) 30.3 7.79%
#> tv log V 3.48 0.0476 1.37 32.5 (29.6, 35.6) 15.4 7.17%
#> add.sd 1.02 1.02
#>
#> 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: 132 × 25
#> ID TIME DV EPRED ERES NPDE NPD PRED RES WRES IPRED IRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.142 0.598 0.915 0.915 0 0.74 0.726 0 0.74
#> 2 1 0.25 2.84 3.23 -0.394 -0.573 -0.573 3.12 -0.276 -0.246 3.57 -0.730
#> 3 1 0.57 6.57 5.69 0.878 -0.449 -0.449 5.61 0.963 0.730 6.45 0.122
#> # … with 129 more rows, and 13 more variables: 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) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 116.8215 373.4213 393.6009 -179.7106 21.94615
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 4.572886 0.473139 0.473141 0.724 1.679834
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka Log Ka 0.473 0.385 81.4 1.6 (0.754, 3.41) 70.2 1.67%
#> tcl Log Cl 1.01 0.305 30.2 2.75 (1.51, 5.01) 26.3 3.10%
#> tv log V 3.46 0.179 5.17 31.8 (22.4, 45.2) 14.0 10.6%
#> add.sd 0.699 0.699
#>
#> Covariance Type ($covMethod): s
#> Some strong fixed parameter correlations exist ($cor) :
#> cor:tcl,tka cor:tv,tka cor:tv,tcl
#> 0.0188 -0.152 -0.843
#>
#>
#> 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):
#> 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:
#> nlmixr(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 132 × 26
#> ID TIME DV EPRED ERES NPDE NPD PRED RES WRES IPRED IRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.0971 0.643 1.45 1.45 0 0.74 1.06 0 0.74
#> 2 1 0.25 2.84 3.63 -0.792 -0.696 -0.696 3.28 -0.442 -0.235 3.85 -1.01
#> 3 1 0.57 6.57 5.90 0.668 -1.83 -1.83 5.85 0.716 0.288 6.79 -0.215
#> # … with 129 more rows, and 14 more variables: 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.06688996 0.00000000 0.00000
#> eta.v 0.00000000 0.01952671 0.00000
#> bsv.ka 0.00000000 0.00000000 0.40026Note 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 %>% ## FIXME could be based on if it finds the covarite in the last used nlmixr data.
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) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 116.8401 373.4398 393.6194 -179.7199 74527.99
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 4.45397 0.504872 0.504875 0.766 3.356283
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tka Log Ka 0.5 FIXED FIXED 1.65 70.3
#> tcl Log Cl 4.04 0.0793 1.96 57.1 (48.9, 66.7) 26.6
#> tv log V 3.46 0.0422 1.22 31.9 (29.3, 34.6) 14.0
#> add.sd 0.695 0.695
#> WT 1.23 0.171 14 1.23 (0.89, 1.56)
#> covWtPow 0.75 FIXED FIXED 0.75
#> Shrink(SD)%
#> tka 1.61%
#> tcl 3.68%
#> tv 10.4%
#> add.sd
#> WT
#> covWtPow
#>
#> Covariance Type ($covMethod): |r|,|s|
#> Some strong fixed parameter correlations exist ($cor) :
#> cor:tv,tcl cor:WT,tcl cor:WT,tv
#> -0.708 -0.721 1.00
#>
#>
#> 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):
#> 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:
#> nlmixr(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 132 × 26
#> ID TIME DV EPRED ERES NPDE NPD PRED RES WRES IPRED IRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.0965 0.643 1.26 1.26 0 0.74 1.06 0 0.74
#> 2 1 0.25 2.84 3.90 -1.06 -0.394 -0.394 3.35 -0.514 -0.269 3.85 -1.01
#> 3 1 0.57 6.57 6.22 0.349 -1.79 -1.79 5.95 0.621 0.249 6.79 -0.219
#> # … with 129 more rows, and 14 more variables: 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 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) fit ────────────────────────────────────────────
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#> FOCEi 104.3608 362.9606 386.023 -173.4803 23.74029
#>
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#>
#> setup optimize covariance table other
#> elapsed 4.691085 0.493411 0.493413 0.747 4.789091
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tka Log Ka 0.406 0.395 97.3 1.5 (0.692, 3.25) 69.8
#> tcl Log Cl 1.02 0.248 24.2 2.79 (1.71, 4.53) 26.0
#> tv log V 3.47 0.146 4.22 32 (24, 42.6) 12.8
#> add.err 0.28 0.28
#> prop.err 0.133 0.133
#> Shrink(SD)%
#> tka 2.63%
#> tcl 1.75%
#> tv 15.6%
#> add.err
#> prop.err
#>
#> Covariance Type ($covMethod): s
#> Some strong fixed parameter correlations exist ($cor) :
#> cor:tcl,tka cor:tv,tka cor:tv,tcl
#> 0.0832 -0.211 -0.772
#>
#>
#> 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):
#> 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:
#> nlmixr(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 132 × 26
#> ID TIME DV EPRED ERES NPDE NPD PRED RES WRES IPRED IRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0.0389 0.701 2.47 2.47 0 0.74 2.64 0 0.74
#> 2 1 0.25 2.84 3.60 -0.755 -1.07 -1.07 3.09 -0.252 -0.146 3.47 -0.630
#> 3 1 0.57 6.57 5.85 0.716 -0.515 -0.515 5.59 0.982 0.401 6.24 0.328
#> # … with 129 more rows, and 14 more variables: 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>