Working with multiple endpoints
2023-05-30
Source:vignettes/multiple-endpoints.Rmd
multiple-endpoints.Rmd
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=2Notice 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) | centerWith 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=2Notice 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=4Then 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=2Running 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.00Since 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)
#> ── nlmixr² SAEM OBJF by FOCEi approximation ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 1383.99 2309.685 2389.105 -1135.842 1675.519 20.41992
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance saem table compress
#> elapsed 0.00248 5e-06 0.074007 186.716 14.98 0.034
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Est. SE %RSE Back-transformed(95%CI) BSV(CV% or SD) Shrink(SD)%
#> tktr 0.426 0.441 104 1.53 (0.644, 3.64) 103. 52.3%
#> tka -0.19 0.222 117 0.827 (0.535, 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.0476 2.38 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.412 14 0.95 (0.894, 0.977) 3.00 83.9%
#> tec50 -0.17 0.151 89 0.844 (0.627, 1.13) 49.2 8.67%
#> tkout -2.9 0.0394 1.36 0.0548 (0.0507, 0.0592) 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.95 -1.95 0.534 -1.48 0.703 0.07 1.45 -1.45
#> 2 1 1 cp 1.9 4.52 -2.62 0.515 -0.967 0.697 0.167 4.06 -2.16
#> 3 1 2 cp 3.3 8.24 -4.94 -2.94 -1.53 0.00167 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
#> doneFOCEi Diagnostic Plots
print(fit.TOF)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 1345.599 2271.293 2350.714 -1116.647 18451.62 74.84733
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.00279 89.75984 89.75984 2.497 0.012 472.1575
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Est. SE %RSE Back-transformed(95%CI) BSV(CV% or SD)
#> tktr -0.0547 0.0619 113 0.947 (0.839, 1.07) 68.2
#> tka 0.388 0.0646 16.7 1.47 (1.3, 1.67) 120.
#> tcl -2.02 0.0195 0.965 0.132 (0.128, 0.138) 28.1
#> tv 2.08 0.0151 0.728 7.99 (7.76, 8.23) 24.1
#> prop.err 0.109 0.109
#> pkadd.err 0.327 0.327
#> temax 5.7 0.71 12.5 0.997 (0.987, 0.999) 3.00
#> tec50 0.175 0.0477 27.3 1.19 (1.08, 1.31) 44.8
#> tkout -2.94 0.0265 0.901 0.0527 (0.05, 0.0555) 10.4
#> te0 4.57 0.00692 0.152 96.4 (95.1, 97.7) 7.08
#> pdadd.err 3.8 3.8
#> Shrink(SD)%
#> tktr 59.0%
#> tka 57.2%
#> tcl -1.32%
#> tv 9.16%
#> prop.err
#> pkadd.err
#> temax 95.7%
#> tec50 3.75%
#> tkout 31.3%
#> te0 27.3%
#> pdadd.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):
#> 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 2.04 -2.04 -0.297 -2.05 0.383 0.02 1.47 -1.47
#> 2 1 1 cp 1.9 4.56 -2.66 0.928 -0.854 0.823 0.197 4.06 -2.16
#> 3 1 2 cp 3.3 7.87 -4.57 -2.94 -1.32 0.00167 0.0933 8.29 -4.99
#> # ℹ 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