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)<- 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) fixThe [derived] block then uses them to compute:
EFFEMAX * IPRED / (EC50 + IPRED)Emax response at each time point
CMAX_INTERVALmax(IPRED, TAD < 24)Peak in current dosing interval
TAM_TAUintegral(1.0, IPRED > MEC, window=24, anchor=0, step=0.1)Hours above MEC per 24-h interval
TAM_D1integral(1.0, IPRED > MEC, from=0, to=24, step=0.1)Hours above MEC on day 1
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.
============================================================
NONLINEAR MIXED EFFECTS MODEL ESTIMATION
============================================================
Model: warfarin_derived_pkpd Dataset: warfarin
Method: FOCE | Gradient: ANALYTIC | Subjects: 10 | Obs: 110
STATUS: CONVERGED 67 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.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.4.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:
param transform estimate se rse_pct lower_95
TVCL TVCL identity 0.132969032 0.0066273862 4.984158 0.119979355
TVV TVV identity 7.730699911 0.2340548168 3.027602 7.271952470
TVKA TVKA identity 0.725207008 0.1245295553 17.171587 0.481129080
MEC MEC identity 0.500000000 0.0000000000 0.000000 0.500000000
EMAX EMAX identity 80.000000000 0.0000000000 0.000000 80.000000000
EC50 EC50 identity 3.000000000 0.0000000000 0.000000 3.000000000
ETA_CL ETA_CL variance 0.028594955 0.0127980762 44.756414 0.003510725
ETA_V ETA_V variance 0.009576924 0.0042973116 44.871523 0.001154193
ETA_KA ETA_KA variance 0.348964242 0.1608188259 46.084615 0.033759343
PROP_ERR PROP_ERR proportional 0.010748528 0.0009417291 8.761470 0.008902739
upper_95 estimate_natural lower_95_natural upper_95_natural
TVCL 0.14595871 NA NA NA
TVV 8.18944735 NA NA NA
TVKA 0.96928494 NA NA NA
MEC 0.50000000 NA NA NA
EMAX 80.00000000 NA NA NA
EC50 3.00000000 NA NA NA
ETA_CL 0.05367918 NA NA NA
ETA_V 0.01799965 NA NA NA
ETA_KA 0.66416914 NA NA NA
PROP_ERR 0.01259432 NA NA NA
init_as_sd weight
TVCL FALSE <NA>
TVV FALSE <NA>
TVKA FALSE <NA>
MEC FALSE <NA>
EMAX FALSE <NA>
EC50 FALSE <NA>
ETA_CL FALSE <NA>
ETA_V FALSE <NA>
ETA_KA FALSE <NA>
PROP_ERR TRUE <NA>
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)