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)
Warning in .ferx_compute_cor_matrix(result$cov_matrix): One or more diagonal
elements are non-positive; correlation matrix may not be meaningful.
print(fit)
============================================================
 NONLINEAR MIXED EFFECTS MODEL ESTIMATION
============================================================
 Model: warfarin_derived_pkpd    Dataset: warfarin
 Method: FOCE | Gradient: ANALYTIC | Subjects: 10 | Obs: 110

 STATUS: CONVERGED   50 iterations   2.3s
 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.132969     0.006627        5.0
TVV                  7.730700     0.234055        3.0
TVKA                 0.725207     0.124530       17.2
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.028595  CV% = 17.0  SE = 0.012798
  ETA_V                    [log-normal]  = 0.009577  CV% = 9.8  SE = 0.004297
  ETA_KA                   [log-normal]  = 0.348964  CV% = 64.6  SE = 0.160819

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.6   DW: 2.61 [negative autocorrelation]   IWRES lag-1 r: -0.365

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

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

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

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

fit$estimates
      param    transform     estimate           se   rse_pct     lower_95
1      TVCL     identity  0.132969032 0.0066273862  4.984158  0.119979355
2       TVV     identity  7.730699900 0.2340548162  3.027602  7.271952460
3      TVKA     identity  0.725207009 0.1245295550 17.171587  0.481129081
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.028594955 0.0127980762 44.756414  0.003510725
8     ETA_V     variance  0.009576924 0.0042973116 44.871522  0.001154193
9    ETA_KA     variance  0.348964241 0.1608188257 46.084615  0.033759343
10 PROP_ERR proportional  0.010748528 0.0009417291  8.761470  0.008902739
      upper_95 estimate_natural lower_95_natural upper_95_natural init_as_sd
1   0.14595871               NA               NA               NA      FALSE
2   8.18944734               NA               NA               NA      FALSE
3   0.96928494               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.05367918               NA               NA               NA      FALSE
8   0.01799965               NA               NA               NA      FALSE
9   0.66416914               NA               NA               NA      FALSE
10  0.01259432               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.383285 63.31396      24     24      11.38329
5   1   8.0 10.757785 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