Skip to contents
nlmixr
nlmixr

Multiple endpoints

Joint PK/PD models, or PK/PD models where you fix certain components are common in pharmacometrics. A classic example, (provided by Tomoo Funaki and Nick Holford) is Warfarin.

In this example, we have a transit-compartment (from depot to gut to central volume) PK model and an effect compartment for the PCA measurement.

Below is an illustrated example of a model that can be applied to the data:

pk.turnover.emax <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)
    ##
    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    prop.err <- 0.1
    pkadd.err <- 0.1
    ##
    temax <- logit(0.8)
    #temax <- 7.5
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    ##
    eta.emax ~ .5
    eta.ec50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5
    ##
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    ##
    #poplogit = log(temax/(1-temax))
    emax=expit(temax+eta.emax)
    #logit=temax+eta.emax
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    effect(0) = e0
    kin = e0*kout
    ##
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center
    d/dt(effect) = kin*PD -kout*effect
    ##
    cp = center / v
    cp ~ prop(prop.err) + add(pkadd.err)
    effect ~ add(pdadd.err)
  })
}

Notice there are two endpoints in the model cp and effect. Both are modeled in nlmixr using the ~ “modeled by” specification.

To see more about how nlmixr will handle the multiple compartment model, it is quite informative to parse the model and print the information about that model. In this case an initial parsing would give:

ui <- nlmixr(pk.turnover.emax)
ui
#>  ── rxode2-based free-form 4-cmt ODE model ────────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>       tktr        tka        tcl         tv   prop.err  pkadd.err      temax 
#>  0.0000000  0.0000000 -2.3025851  2.3025851  0.1000000  0.1000000  1.3862944 
#>      tec50      tkout        te0  pdadd.err 
#> -0.6931472 -2.9957323  4.6051702 10.0000000 
#> 
#> Omega ($omega): 
#>          eta.ktr eta.ka eta.cl eta.v eta.emax eta.ec50 eta.kout eta.e0
#> eta.ktr        1      0      0     0      0.0      0.0      0.0    0.0
#> eta.ka         0      1      0     0      0.0      0.0      0.0    0.0
#> eta.cl         0      0      2     0      0.0      0.0      0.0    0.0
#> eta.v          0      0      0     1      0.0      0.0      0.0    0.0
#> eta.emax       0      0      0     0      0.5      0.0      0.0    0.0
#> eta.ec50       0      0      0     0      0.0      0.5      0.0    0.0
#> eta.kout       0      0      0     0      0.0      0.0      0.5    0.0
#> eta.e0         0      0      0     0      0.0      0.0      0.0    0.5
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1            depot
#> 2                  2              gut
#> 3                  3           center
#> 4                  4           effect
#>  ── Multiple Endpoint Model ($multipleEndpoint): ──  
#>     variable                   cmt                   dvid*
#> 1     cp ~ …     cmt='cp' or cmt=5     dvid='cp' or dvid=1
#> 2 effect ~ … cmt='effect' or cmt=4 dvid='effect' or dvid=2
#>   * If dvids are outside this range, all dvids are re-numered sequentially, ie 1,7, 10 becomes 1,2,3 etc
#> 
#>  ── μ-referencing ($muRefTable): ──  
#>   theta      eta level
#> 1  tktr  eta.ktr    id
#> 2   tka   eta.ka    id
#> 3   tcl   eta.cl    id
#> 4    tv    eta.v    id
#> 5 temax eta.emax    id
#> 6 tec50 eta.ec50    id
#> 7 tkout eta.kout    id
#> 8   te0   eta.e0    id
#> 
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     ini({
#>         tktr <- 0
#>         tka <- 0
#>         tcl <- -2.30258509299405
#>         tv <- 2.30258509299405
#>         prop.err <- c(0, 0.1)
#>         pkadd.err <- c(0, 0.1)
#>         temax <- 1.38629436111989
#>         tec50 <- -0.693147180559945
#>         tkout <- -2.99573227355399
#>         te0 <- 4.60517018598809
#>         pdadd.err <- c(0, 10)
#>         eta.ktr ~ 1
#>         eta.ka ~ 1
#>         eta.cl ~ 2
#>         eta.v ~ 1
#>         eta.emax ~ 0.5
#>         eta.ec50 ~ 0.5
#>         eta.kout ~ 0.5
#>         eta.e0 ~ 0.5
#>     })
#>     model({
#>         ktr <- exp(tktr + eta.ktr)
#>         ka <- exp(tka + eta.ka)
#>         cl <- exp(tcl + eta.cl)
#>         v <- exp(tv + eta.v)
#>         emax = expit(temax + eta.emax)
#>         ec50 = exp(tec50 + eta.ec50)
#>         kout = exp(tkout + eta.kout)
#>         e0 = exp(te0 + eta.e0)
#>         DCP = center/v
#>         PD = 1 - emax * DCP/(ec50 + DCP)
#>         effect(0) = e0
#>         kin = e0 * kout
#>         d/dt(depot) = -ktr * depot
#>         d/dt(gut) = ktr * depot - ka * gut
#>         d/dt(center) = ka * gut - cl/v * center
#>         d/dt(effect) = kin * PD - kout * effect
#>         cp = center/v
#>         cp ~ prop(prop.err) + add(pkadd.err)
#>         effect ~ add(pdadd.err)
#>     })
#> }

In the middle of the printout, it shows how the data must be formatted (using the cmt and dvid data items) to allow nlmixr to model the multiple endpoint appropriately.

Of course, if you are interested you can directly access the information in ui$multipleEndpoint.

ui$multipleEndpoint
#>     variable                   cmt                   dvid*
#> 1     cp ~ …     cmt='cp' or cmt=5     dvid='cp' or dvid=1
#> 2 effect ~ … cmt='effect' or cmt=4 dvid='effect' or dvid=2

Notice that the cmt and dvid items can use the named variables directly as either the cmt or dvid specification. This flexible notation makes it so you do not have to rename your compartments to run nlmixr model functions.

The other thing to note is that the cp is specified by an ODE compartment above the number of compartments defined in the rxode2 part of the nlmixr model. This is because cp is not a defined compartment, but a related variable cp.

The last thing to notice that the cmt items are numbered cmt=5 for cp or cmt=4 for effect even though they were specified in the model first by cp and cmt. This ordering is because effect is a compartment in the rxode2 system. Of course cp is related to the compartment center, and it may make more sense to pair cp with the center compartment.

If this is something you want to have you can specify the compartment to relate the effect to by the | operator. In this case you would change

cp ~ prop(prop.err) + add(pkadd.err)

to

cp ~ prop(prop.err) + add(pkadd.err) | center

With this change, the model could be updated to:

pk.turnover.emax2 <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)
    ##
    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    prop.err <- 0.1
    pkadd.err <- 0.1
    ##
    temax <- logit(0.8)
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    ##
    eta.emax ~ .5
    eta.ec50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5
    ##
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    ##
    emax=expit(temax+eta.emax)
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    effect(0) = e0
    kin = e0*kout
    ##
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center
    d/dt(effect) = kin*PD -kout*effect
    ##
    cp = center / v
    cp ~ prop(prop.err) + add(pkadd.err) | center
    effect ~ add(pdadd.err)
  })
}
ui2 <- nlmixr(pk.turnover.emax2)
ui2$multipleEndpoint
#>     variable                   cmt                   dvid*
#> 1     cp ~ … cmt='center' or cmt=3 dvid='center' or dvid=1
#> 2 effect ~ … cmt='effect' or cmt=4 dvid='effect' or dvid=2

Notice in this case the cmt variables are numbered sequentially and the cp variable matches the center compartment.

DVID vs CMT, which one is used

When dvid and cmt are combined in the same dataset, the cmt data item is always used on the event information and the dvid is used on the observations. nlmixr expects the cmt data item to match the dvid item for observations OR to be either zero or one for the dvid to replace the cmt information.

If you do not wish to use dvid items to define multiple endpoints in nlmixr, you can set the following option:

options(rxode2.combine.dvid=FALSE)
ui2$multipleEndpoint
#>     variable                   cmt
#> 1     cp ~ … cmt='center' or cmt=3
#> 2 effect ~ … cmt='effect' or cmt=4

Then only cmt items are used for the multiple endpoint models. Of course you can turn it on or off for different models if you wish:

options(rxode2.combine.dvid=TRUE)
ui2$multipleEndpoint
#>     variable                   cmt                   dvid*
#> 1     cp ~ … cmt='center' or cmt=3 dvid='center' or dvid=1
#> 2 effect ~ … cmt='effect' or cmt=4 dvid='effect' or dvid=2

Running a multiple endpoint model

With this information, we can use the built-in warfarin dataset in nlmixr2:

summary(warfarin)
#>        id             time             amt                dv          dvid    
#>  Min.   : 1.00   Min.   :  0.00   Min.   :  0.000   Min.   :  0.00   cp :283  
#>  1st Qu.: 8.00   1st Qu.: 24.00   1st Qu.:  0.000   1st Qu.:  4.50   pca:232  
#>  Median :15.00   Median : 48.00   Median :  0.000   Median : 11.40            
#>  Mean   :16.08   Mean   : 52.08   Mean   :  6.524   Mean   : 20.02            
#>  3rd Qu.:24.00   3rd Qu.: 96.00   3rd Qu.:  0.000   3rd Qu.: 26.00            
#>  Max.   :33.00   Max.   :144.00   Max.   :153.000   Max.   :100.00            
#>       evid               wt              age            sex     
#>  Min.   :0.00000   Min.   : 40.00   Min.   :21.00   female:101  
#>  1st Qu.:0.00000   1st Qu.: 60.00   1st Qu.:23.00   male  :414  
#>  Median :0.00000   Median : 70.00   Median :28.00               
#>  Mean   :0.06214   Mean   : 69.27   Mean   :31.85               
#>  3rd Qu.:0.00000   3rd Qu.: 78.00   3rd Qu.:36.00               
#>  Max.   :1.00000   Max.   :102.00   Max.   :63.00

Since dvid specifies pca as the effect endpoint, you can update the model to be more explicit making one last change:

cp ~ prop(prop.err) + add(pkadd.err)
effect ~ add(pdadd.err) 

to

cp ~ prop(prop.err) + add(pkadd.err)
effect ~ add(pdadd.err)  | pca
pk.turnover.emax3 <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)
    ##
    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    prop.err <- 0.1
    pkadd.err <- 0.1
    ##
    temax <- logit(0.8)
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    ##
    eta.emax ~ .5
    eta.ec50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5
    ##
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    emax = expit(temax+eta.emax)
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    effect(0) = e0
    kin = e0*kout
    ##
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center
    d/dt(effect) = kin*PD -kout*effect
    ##
    cp = center / v
    cp ~ prop(prop.err) + add(pkadd.err)
    effect ~ add(pdadd.err) | pca
  })
}

Run the models with SAEM

fit.TOS <- nlmixr(pk.turnover.emax3, warfarin, "saem", 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 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 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 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 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

print(fit.TOS)
#> ── nlmix SAEM OBJF by FOCEi approximation ──
#> 
#>          OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 1383.99 2309.685 2389.105      -1135.842        1711.869        19.29712
#> 
#> ── Time (sec $time): ──
#> 
#>            setup optimize covariance    saem  table compress
#> elapsed 0.002451    5e-06   0.060006 119.586 10.914    0.026
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>            Est.     SE  %RSE Back-transformed(95%CI) BSV(CV% or SD) Shrink(SD)%
#> tktr      0.426  0.443   104      1.53 (0.642, 3.65)           103.      52.3% 
#> tka       -0.19  0.225   118     0.827 (0.532, 1.28)           41.7      63.6% 
#> tcl       -1.97  0.051  2.58    0.139 (0.126, 0.153)           26.6      5.33% 
#> tv            2 0.0475  2.37       7.41 (6.75, 8.14)           21.4      17.9% 
#> prop.err  0.124                                0.124                           
#> pkadd.err 0.758                                0.758                           
#> temax      2.93  0.398  13.6     0.95 (0.896, 0.976)          0.107      83.9% 
#> tec50     -0.17  0.147  86.5     0.844 (0.633, 1.13)           49.2      8.67% 
#> tkout      -2.9  0.039  1.34 0.0548 (0.0507, 0.0591)           6.31      46.2% 
#> te0        4.57 0.0115 0.252       96.5 (94.3, 98.7)           5.15      17.8% 
#> pdadd.err  3.73                                 3.73                           
#>  
#>   Covariance Type ($covMethod): linFim
#>   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 
#>   Censoring ($censInformation): No censoring
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 483 × 44
#>   ID     TIME CMT      DV EPRED  ERES  NPDE    NPD    PDE     PD  PRED   RES
#>   <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl>
#> 1 1       0.5 cp      0    1.91 -1.91 -1.68 -1.50  0.0467 0.0667  1.45 -1.45
#> 2 1       1   cp      1.9  4.32 -2.42  1.83 -0.954 0.967  0.17    4.06 -2.16
#> 3 1       2   cp      3.3  8.12 -4.82 -2.05 -1.53  0.02   0.0633  8.47 -5.17
#> # ℹ 480 more rows
#> # ℹ 32 more variables: WRES <dbl>, IPRED <dbl>, IRES <dbl>, IWRES <dbl>,
#> #   CPRED <dbl>, CRES <dbl>, CWRES <dbl>, eta.ktr <dbl>, eta.ka <dbl>,
#> #   eta.cl <dbl>, eta.v <dbl>, eta.emax <dbl>, eta.ec50 <dbl>, eta.kout <dbl>,
#> #   eta.e0 <dbl>, depot <dbl>, gut <dbl>, center <dbl>, effect <dbl>,
#> #   ktr <dbl>, ka <dbl>, cl <dbl>, v <dbl>, emax <dbl>, ec50 <dbl>, kout <dbl>,
#> #   e0 <dbl>, DCP <dbl>, PD.1 <dbl>, kin <dbl>, tad <dbl>, dosenum <dbl>

SAEM Diagnostic plots

plot(fit.TOS)



v1s <- vpcPlot(fit.TOS, show=list(obs_dv=TRUE), scales="free_y") +
  ylab("Warfarin Cp [mg/L] or PCA") +
  xlab("Time [h]")

v2s <- vpcPlot(fit.TOS, show=list(obs_dv=TRUE), pred_corr = TRUE) +
  ylab("Prediction Corrected Warfarin Cp [mg/L] or PCA") +
  xlab("Time [h]")

v1s

v2s

FOCEi fits

## FOCEi fit/vpcs
fit.TOF <- nlmixr(pk.turnover.emax3, warfarin, "focei", control=list(print=0),
                  table=list(cwres=TRUE, npde=TRUE))
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:01:29 
#> done

FOCEi Diagnostic Plots

print(fit.TOF)
#> ── nlmix FOCEi (outer: nlminb) ──
#> 
#>           OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 1369.349 2295.043 2374.464      -1128.522         1007646        6234.877
#> 
#> ── Time (sec $time): ──
#> 
#>            setup optimize covariance table compress    other
#> elapsed 0.001961 89.30359   89.30359 1.852    0.008 184.7599
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>            Est.     SE  %RSE Back-transformed(95%CI) BSV(CV% or SD) Shrink(SD)%
#> tktr      0.589  0.303  51.4       1.8 (0.995, 3.26)           181.      63.8% 
#> tka       0.347  0.135  38.9       1.41 (1.09, 1.84)           151.      66.3% 
#> tcl       -2.03  0.239  11.8   0.132 (0.0827, 0.211)           31.4      8.61% 
#> tv         2.06 0.0192  0.93       7.89 (7.59, 8.19)           23.1      7.27% 
#> prop.err  0.117                                0.117                           
#> pkadd.err 0.249                                0.249                           
#> temax      5.45   2.18  40.1        0.996 (0.763, 1)          0.480      97.7% 
#> tec50     0.173 0.0853  49.4        1.19 (1.01, 1.4)           51.8      11.7% 
#> tkout     -2.94  0.014 0.476 0.0527 (0.0513, 0.0542)           10.4      29.6% 
#> te0        4.57 0.0156 0.341       96.4 (93.5, 99.4)           7.08      26.1% 
#> pdadd.err   3.8                                  3.8                           
#>  
#>   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: 483 × 44
#>   ID     TIME CMT      DV EPRED  ERES  NPDE    NPD    PDE    PD  PRED   RES
#>   <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
#> 1 1       0.5 cp      0    3.19 -3.19 -1.64 -2.05  0.05   0.02   2.40 -2.40
#> 2 1       1   cp      1.9  5.70 -3.80  1.99 -0.941 0.977  0.173  5.94 -4.04
#> 3 1       2   cp      3.3  8.52 -5.22 -2.13 -1.16  0.0167 0.123 10.3  -6.97
#> # ℹ 480 more rows
#> # ℹ 32 more variables: WRES <dbl>, IPRED <dbl>, IRES <dbl>, IWRES <dbl>,
#> #   CPRED <dbl>, CRES <dbl>, CWRES <dbl>, eta.ktr <dbl>, eta.ka <dbl>,
#> #   eta.cl <dbl>, eta.v <dbl>, eta.emax <dbl>, eta.ec50 <dbl>, eta.kout <dbl>,
#> #   eta.e0 <dbl>, depot <dbl>, gut <dbl>, center <dbl>, effect <dbl>,
#> #   ktr <dbl>, ka <dbl>, cl <dbl>, v <dbl>, emax <dbl>, ec50 <dbl>, kout <dbl>,
#> #   e0 <dbl>, DCP <dbl>, PD.1 <dbl>, kin <dbl>, tad <dbl>, dosenum <dbl>
plot(fit.TOF)


v1f <- vpcPlot(fit.TOF, show=list(obs_dv=TRUE), scales="free_y") +
  ylab("Warfarin Cp [mg/L] or PCA") +
  xlab("Time [h]")

v2f <- vpcPlot(fit.TOF, show=list(obs_dv=TRUE), pred_corr = TRUE) +
  ylab("Prediction Corrected Warfarin Cp [mg/L] or PCA") +
  xlab("Time [h]")

v1f

v2f