Example: Derived PK/PD quantities

This example uses [derived] to compute pharmacodynamic metrics from PK predictions after the fit: an Emax response, a TAD-gated Cmax, and time-above-MEC per dosing interval.

The PD parameters (MEC, EMAX, EC50) are fixed at assumed values — they do not enter the PK likelihood and cannot be estimated from PK-only data. Only TVCL, TVV, TVKA, and the random effects are estimated.

For a model where PD observations are also fitted jointly, see Example: simultaneous PK/PD (Emax).

Model

library(ferx)

ex  <- ferx_example("warfarin_derived_pkpd")
ferx_model_show(ex$model)
# model: warfarin_derived_pkpd.ferx 
# Warfarin with PD effect and time-above-MEC derived quantities

[parameters]
  theta TVCL(0.2, 0.001, 10.0)
  theta TVV(10.0, 0.1, 500.0)
  theta TVKA(1.5, 0.01, 50.0)
  theta MEC(0.5)    fix
  theta EMAX(80.0)  fix
  theta EC50(3.0)   fix

  omega ETA_CL ~ 0.09
  omega ETA_V  ~ 0.04
  omega ETA_KA ~ 0.30

  sigma PROP_ERR ~ 0.02 (sd)

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV  * exp(ETA_V)
  KA = TVKA * exp(ETA_KA)

[structural_model]
  pk one_cpt_oral(cl=CL, v=V, ka=KA)

[error_model]
  DV ~ proportional(PROP_ERR)

[derived]
  KE     = CL / V
  T_HALF = 0.6931472 / KE

  # Emax PD effect at current concentration
  EFF    = EMAX * IPRED / (EC50 + IPRED)

  # TAD-gated Cmax: peak within the current dosing interval
  CMAX_INTERVAL = max(IPRED, TAD < 24)

  # Time above MEC per dosing interval (hours)
  TAM_TAU = integral(1.0, IPRED > MEC, window=24, anchor=0, step=0.1)

  # Time above MEC on day 1 specifically
  TAM_D1  = integral(1.0, IPRED > MEC, from=0, to=24, step=0.1)

[output]
  CL V KA

[fit_options]
  method  = foce
  maxiter = 300

The fixed PD thetas appear with fix in the [parameters] block:

theta MEC(0.5)    fix
theta EMAX(80.0)  fix
theta EC50(3.0)   fix

The [derived] block then uses them to compute:

Column Expression Description
EFF EMAX * IPRED / (EC50 + IPRED) Emax response at each time point
CMAX_INTERVAL max(IPRED, TAD < 24) Peak in current dosing interval
TAM_TAU integral(1.0, IPRED > MEC, window=24, anchor=0, step=0.1) Hours above MEC per 24-h interval
TAM_D1 integral(1.0, IPRED > MEC, from=0, to=24, step=0.1) Hours above MEC on day 1

Fit

fit <- ferx_fit(ex$model, ex$data, verbose = FALSE)
print(fit)
============================================================
 NONLINEAR MIXED EFFECTS MODEL ESTIMATION
============================================================
 Model: warfarin_derived_pkpd    Dataset: warfarin
 Method: FOCE | Gradient: ANALYTIC (DUAL2) | Subjects: 10 | Obs: 110

 STATUS: CONVERGED   102 iterations   2.5s
 OFV: -280.3640    AIC: -266.3640    BIC: -247.4606

MODEL STRUCTURE (auto-derived)
------------------------------------------------------------
  Structural:  1-cpt oral  (TVCL, TVV, TVKA, MEC, EMAX, EC50)
  IIV:         ETA_CL, ETA_V, ETA_KA
  IOV:         none
  Residual:    proportional

THETA
------------------------------------------------------------
Parameter            Estimate           SE       %RSE
----------------------------------------------------
TVCL                 0.132970     0.006628        5.0
TVV                  7.730695     0.234062        3.0
TVKA                 0.725187     0.125303       17.3
MEC                  0.500000     0.000000        0.0
EMAX                80.000000     0.000000        0.0
EC50                 3.000000     0.000000        0.0

OMEGA  (between-subject variability)
------------------------------------------------------------
  ETA_CL                   [log-normal]  = 0.028585  CV% = 17.0  SE = 0.012790
  ETA_V                    [log-normal]  = 0.009575  CV% = 9.8  SE = 0.004294
  ETA_KA                   [log-normal]  = 0.348959  CV% = 64.6  SE = 0.160968

SIGMA  (residual error)
------------------------------------------------------------
  PROP_ERR         [proportional] = 0.010749  (var = 0.000116, CV% = 1.1)  SE = 0.000942  [initial specified as SD]

SHRINKAGE
------------------------------------------------------------
 ETA_CL: 0.0%   ETA_V: 0.1%   ETA_KA: 0.2%   EPS: 16.1%

DIAGNOSTICS
------------------------------------------------------------
 Covariance: computed   Cond: 2.7   DW: 2.61 [negative autocorrelation]   IWRES lag-1 r: -0.365

RUN INFO
------------------------------------------------------------
 Gradient (requested): auto   (used: analytic (Dual2))
 ferx v0.1.6 (core v0.1.6)

SETTINGS  (model file / call-time override)
------------------------------------------------------------
  method                       foce  [model only]
  maxiter                      300  [model only]

------------------------------------------------------------
 4 warning  --  call ferx_warnings(fit) for details
============================================================

The PK thetas have proper SEs; MEC, EMAX, and EC50 show SE = 0 because they are fixed:

ferx_estimates(fit)
      param    transform     estimate           se   rse_pct     lower_95
1      TVCL     identity  0.132969570 0.0066284653  4.984949  0.119977778
2       TVV     identity  7.730695262 0.2340617937  3.027694  7.271934147
3      TVKA     identity  0.725187222 0.1253027264 17.278673  0.479593879
4       MEC     identity  0.500000000 0.0000000000  0.000000  0.500000000
5      EMAX     identity 80.000000000 0.0000000000  0.000000 80.000000000
6      EC50     identity  3.000000000 0.0000000000  0.000000  3.000000000
7    ETA_CL     variance  0.028585310 0.0127900421 44.743410  0.003516827
8     ETA_V     variance  0.009575317 0.0042943833 44.848474  0.001158325
9    ETA_KA     variance  0.348959133 0.1609684369 46.128163  0.033460997
10 PROP_ERR proportional  0.010748634 0.0009421476  8.765277  0.008902025
      upper_95 estimate_natural lower_95_natural upper_95_natural init_as_sd
1   0.14596136               NA               NA               NA      FALSE
2   8.18945638               NA               NA               NA      FALSE
3   0.97078057               NA               NA               NA      FALSE
4   0.50000000               NA               NA               NA      FALSE
5  80.00000000               NA               NA               NA      FALSE
6   3.00000000               NA               NA               NA      FALSE
7   0.05365379               NA               NA               NA      FALSE
8   0.01799231               NA               NA               NA      FALSE
9   0.66445727               NA               NA               NA      FALSE
10  0.01259524               NA               NA               NA       TRUE

PD derived columns

head(fit$sdtab[, c("ID", "TIME", "IPRED", "EFF", "TAM_TAU", "TAM_D1",
                   "CMAX_INTERVAL")], 12)
   ID  TIME     IPRED      EFF TAM_TAU TAM_D1 CMAX_INTERVAL
1   1   0.5  5.351751 51.26351      24     24      11.38329
2   1   1.0  8.283413 58.72984      24     24      11.38329
3   1   2.0 10.708454 62.49256      24     24      11.38329
4   1   4.0 11.383286 63.31397      24     24      11.38329
5   1   8.0 10.757786 62.55533      24     24      11.38329
6   1  12.0 10.069303 61.63636      24     24      11.38329
7   1  24.0  8.254797 58.67576      24     24      11.38329
8   1  48.0  5.547784 51.92255      24     24      11.38329
9   1  72.0  3.728488 44.33077      24     24      11.38329
10  1  96.0  2.505797 36.40958      24     24      11.38329
11  1 120.0  1.684066 28.76246      24     24      11.38329
12  2   0.5  3.227396 41.46062      24     24      12.10572

EFF varies per time point (it depends on IPRED). CMAX_INTERVAL, TAM_TAU, and TAM_D1 are aggregate quantities — constant per 24-h window or per subject.

Why step= for time-above-threshold

TAM_TAU uses step=0.1 to force grid evaluation at 0.1-hour intervals. Without step=, only observation time points are evaluated. If the observation design is sparse (large gaps between samples), the trapezoidal rule can miss brief periods above threshold between those times. The step= argument fills in the gaps:

# Without step= (observation-time only -- can miss threshold crossings):
TAM_approx = integral(1.0, IPRED > MEC, window=24, anchor=0)

# With step=0.1 (fine grid -- catches brief excursions):
TAM_TAU    = integral(1.0, IPRED > MEC, window=24, anchor=0, step=0.1)

See also