Post-fit diagnostics

After fitting, ferx returns a rich ferx_fit object with all the diagnostics needed for model evaluation. This article covers the key diagnostic functions and fields.

Tip

For a complete diagnostics walkthrough with plots, see Chapter 5 (Diagnostics) in the ferx book.

Parameter estimates — `ferx_estimates()`

ferx_estimates() collects all estimated parameters (theta, omega, sigma) into a tidy data frame. Run it after any fit that used covariance = TRUE:

ex  <- ferx_example("warfarin")
fit <- ferx_fit(ex$model, ex$data)

ferx_estimates(fit)
#   param           transform    estimate      se  rse_pct  lower_95  upper_95
#   TVCL            identity       0.134  0.0108     8.1     0.113     0.155
#   TVV             identity       7.94   0.420      5.3     7.12      8.76
#   TVKA            identity       1.07   0.160     14.9     0.756     1.38
#   ETA_CL          variance       0.042  0.012     28.6    0.019     0.066
#   ETA_V           variance       0.020  0.009     44.9    0.002     0.037
#   ETA_KA          variance       0.157  0.048     30.6    0.063     0.251
#   PROP_ERR        proportional   0.011  0.001      9.1    0.009     0.013

Columns:

Column	Description
`param`	Parameter name
`transform`	Scale used by the estimator (`"identity"`, `"log"`, `"logit"`, `"variance"` for omega)
`estimate`	Point estimate on the estimation scale
`se`	Standard error (requires `covariance = TRUE`)
`rse_pct`	Relative SE: \|SE / estimate\| × 100
`lower_95`, `upper_95`	95% CI on the estimation scale
`estimate_natural`	Back-transformed estimate (logit theta only)
`lower_95_natural`, `upper_95_natural`	Back-transformed CI (logit theta only)

See Parameter transforms for how the transform column relates to log-normal, logit-normal, and additive ETAs.

Parameter correlation matrix — `ferx_cor_matrix()`

ferx_cor_matrix() converts the estimated parameter covariance matrix to correlations. A value near ±1 between two parameters signals a structural identifiability problem worth investigating:

ferx_cor_matrix(fit)
#                TVCL   TVV  TVKA  ...
#   TVCL         1.00  0.12  0.03
#   TVV          0.12  1.00 -0.07
#   TVKA         0.03 -0.07  1.00

High correlations between CL and V often arise in poorly-sampled datasets where absorption and elimination are not separately identifiable.

ETA-covariate correlations — `ferx_eta_cov()`

ferx_eta_cov() computes Pearson correlations between individual ETA estimates (EBEs) and numeric covariate columns. Identifies covariates worth including in a formal covariate search.

The warfarin dataset has no covariate columns. The following example uses two_cpt_oral_cov, which includes body weight (WT) and creatinine clearance (CRCL):

ex_cov   <- ferx_example("two_cpt_oral_cov")
fit_cov  <- ferx_fit(ex_cov$model, ex_cov$data)
data_cov <- read.csv(ex_cov$data)
ferx_eta_cov(fit_cov, data_cov)
#   eta      covariate    r   p_val  flag
#   ETA_CL   WT        0.38  0.008   [!]
#   ETA_CL   CRCL      0.41  0.003   [!]
#   ETA_V1   WT        0.29  0.047

[!] flags correlations with |r| ≥ 0.3 as candidates for formal covariate testing. Pass the data frame (not the path); any numeric column that is constant within a subject and not a standard NONMEM column (TIME, DV, AMT, EVID, MDV, CMT, RATE) is treated as a covariate.

IWRES autocorrelation — `check_diagnostics()`

check_diagnostics() returns a structured list with an $autocorrelation data frame (Durbin-Watson + lag-1 Pearson r) and a tidy $shrinkage data frame covering all ETA and EPS components:

diag <- check_diagnostics(fit)

diag$autocorrelation
#   dw_statistic  lag1_r               flag
#           1.84  0.0821  no autocorrelation

diag$shrinkage
#    param type shrinkage shrinkage_pct
#   ETA_CL  eta     0.142          14.2
#    ETA_V  eta     0.283          28.3
#   ETA_KA  eta     0.452          45.2
#      EPS  eps     0.094           9.4

print(fit) shows the same values in a --- Diagnostics --- block. A message() with ordered remedies is emitted automatically when DW < 1.5 (positive autocorrelation) or DW > 2.5 (negative autocorrelation).

DW range	Interpretation	Suggested action
< 1.5	Positive autocorrelation	Check time ordering; add transit compartment; consider IOV on ka/F; SDE
1.5 – 2.5	No concern	—
> 2.5	Negative autocorrelation	Over-parameterisation or misspecified error model

For a full walkthrough with plots see Chapter 5 (Diagnostics) in the ferx book.

Goodness-of-fit plots from `fit$sdtab`

fit$sdtab is a data frame with one row per observation containing PRED, IPRED, CWRES, IWRES, and additional columns depending on the model:

Column	Always present	When present
`ID`, `TIME`, `DV`	Yes	—
`PRED`, `IPRED`, `CWRES`, `IWRES`, `EBE_OFV`	Yes	—
`CMT`	No	Multi-endpoint models (any `CMT != 1` in the data)
`OCC`	No	IOV models
`CENS`	No	When the dataset has a `CENS` column (BLOQ)

For multi-endpoint models, split GOF plots by CMT:

library(dplyr)
sdtab <- fit$sdtab |>
  mutate(endpoint = factor(CMT, levels = c(2, 3), labels = c("PK", "PD")))

sdtab <- fit$sdtab

Standard GOF plots with ggplot2:

library(ggplot2)

# DV vs PRED
ggplot(sdtab, aes(PRED, DV)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, colour = "steelblue") +
  labs(title = "DV vs PRED", x = "Population prediction", y = "Observed")

# DV vs IPRED
ggplot(sdtab, aes(IPRED, DV)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, colour = "steelblue") +
  labs(title = "DV vs IPRED", x = "Individual prediction", y = "Observed")

# CWRES vs TIME
ggplot(sdtab, aes(TIME, CWRES)) +
  geom_point(alpha = 0.4) +
  geom_hline(yintercept = c(-2, 0, 2),
             linetype = c("dashed", "solid", "dashed"), colour = "steelblue") +
  labs(title = "CWRES vs TIME")

# CWRES vs PRED
ggplot(sdtab, aes(PRED, CWRES)) +
  geom_point(alpha = 0.4) +
  geom_hline(yintercept = c(-2, 0, 2),
             linetype = c("dashed", "solid", "dashed"), colour = "steelblue") +
  labs(title = "CWRES vs PRED")

Structured warnings — `ferx_warnings()`

ferx classifies non-fatal issues from a run into severity levels and categories and makes them available via ferx_warnings(). Critical issues are also surfaced inline in the STATUS line of print(fit).

ferx_warnings(fit)
#   severity  category              message
#   warning   shrinkage             ETA_KA shrinkage = 0.58 (> 0.40): EBEs
#                                   unreliable for this parameter
#   info      mu_referencing        ETA_CL detected as log-normal mu-referenced
#   info      mu_referencing        ETA_V detected as log-normal mu-referenced

Severity levels:

Severity	Meaning
`critical`	Likely convergence failure or data issue; inspect before using results
`warning`	Non-fatal issue worth investigating (high shrinkage, near-boundary estimate)
`info`	Informational (mu-referencing detection, unused parameter notice)

Common categories:

Category	Examples
`shrinkage`	ETA shrinkage > 40%; EPS shrinkage > 30%
`convergence`	Outer optimizer hit maxiter; OFV gradient not converged
`mu_referencing`	Which ETAs were detected as mu-referenced
`unused_parameter`	Theta/omega declared in `[parameters]` but not used in the model
`boundary`	Estimated parameter is at or near its declared bound

The raw warning strings are also available on the fit object:

fit$warnings             # character vector of all warning messages
fit$warnings_structured  # data frame: severity, category, message

Run log — `ferx_runlog()`

ferx_runlog() produces a NONMEM-style .lst summary in the console — useful for archiving results, quick convergence checks, and sharing with collaborators who expect NONMEM-style output.

ex  <- ferx_example("warfarin")
fit <- ferx_fit(ex$model, ex$data, method = "foce", covariance = TRUE)
ferx_runlog(fit)

The output for the warfarin example looks like:

============================================================
  ferx run log
  Model:   warfarin
  Method:  FOCE
  ferx:    0.1.5
  Date:    2026-06-04 00:54:58
============================================================
============================================================
MODEL FILE
------------------------------------------------------------
# One-compartment oral PK model (warfarin)
[parameters]
  theta TVCL(0.134, 0.001, 10.0)
  ...
============================================================
DATA SUMMARY
------------------------------------------------------------
  Subjects:       10
  Observations:   110
  Time range:     0.5 to 120
============================================================
PARAMETER ESTIMATES
------------------------------------------------------------
  PARAMETER                            INITIAL           FINAL        SE     %RSE
  TVCL                                   0.134        0.132779   0.01433    10.8%
  TVV                                      8.1         7.69273    0.2899     3.8%
  TVKA                                       1        0.756318   0.03491     4.6%
  ETA_CL                                  0.07       0.0284934  0.006363    22.3%
  ETA_V                                   0.02      0.00959794  0.002161    22.5%
  ETA_KA                                   0.4        0.340694   0.07627    22.4%
  PROP_ERR                                0.01       0.0106401  0.0007838     7.4%
============================================================
ESTIMATION SETTINGS
------------------------------------------------------------
  Optimizer:      bobyqa
  Max iterations: 300
  BLOQ method:    drop
============================================================
OBJECTIVE FUNCTION
------------------------------------------------------------
  OFV   = -280.183
  AIC   = -266.183
  BIC   = -247.28
  Convergence: YES
============================================================
COVARIANCE STEP
------------------------------------------------------------
  Condition number: 12.08
  Eigenvalues:      2.224  1.023  0.9999  0.9921  0.9246  0.6524  0.1841
============================================================
DIAGNOSTICS
------------------------------------------------------------
  ETA shrinkage:  ETA_CL=-0.1%  ETA_V=0.0%  ETA_KA=0.1%
  EPS shrinkage:  15.3%
  DW statistic:   2.609 (negative autocorrelation)
  ETA normality (Shapiro-Wilk):
    ETA_CL: W=0.900 p=0.217
    ETA_V: W=0.849 p=0.057
    ETA_KA: W=0.937 p=0.521
============================================================
FINAL GRADIENT
------------------------------------------------------------
  Gradient not available for this optimizer (BOBYQA / BFGS / GN / SAEM).
============================================================
RUNTIME
------------------------------------------------------------
  Wall time:   0.0 seconds
  Gradient:    finite differences (inner) / N/A (outer)
============================================================

Capture to a file:

log_text <- ferx_runlog(fit, verbose = FALSE)
writeLines(log_text, "warfarin_run1.log")

Gradient convergence check (for gradient-based optimizers such as L-BFGS):

ferx_runlog(fit, gradient_tol = 0.001)  # flag any |grad| > 0.001

The FINAL GRADIENT section shows the gradient components at the best-OFV parameter point and warns if any exceed gradient_tol (default 0.01). BOBYQA, GN, and SAEM do not produce a final gradient.

Other `ferx_fit` fields

fit$ofv                  # objective function value (-2 log-likelihood)
fit$aic                  # AIC = OFV + 2p
fit$bic                  # BIC = OFV + p*ln(N)
fit$converged            # TRUE if the outer optimizer converged
fit$method               # estimation method used
fit$gradient_used        # "dual" (analytic) or "finite_differences"
fit$warnings             # non-fatal warnings from the run

fit$ebe_etas             # empirical Bayes ETAs per subject (data frame)
fit$individual_estimates # individual PK parameters (CL, V, KA, ...) per subject
fit$eta_normality        # Shapiro-Wilk normality test p-value per ETA

fit$shrinkage_eta        # ETA shrinkage per random effect
fit$shrinkage_eps        # EPS shrinkage (scalar or vector for multi-sigma models)
fit$dw_statistic         # pooled Durbin-Watson statistic for IWRES
fit$iwres_lag1_r         # pooled lag-1 Pearson r for IWRES

Shrinkage — values near 1 indicate that the data are not informative about a random effect:

\[\text{ETA shrinkage}_k = 1 - \frac{\text{SD}(\hat{\eta}_k)}{\sqrt{\omega_{kk}}} \qquad \text{EPS shrinkage} = 1 - \text{SD}(\text{IWRES})\]

Convergence diagnostics — `ferx_check_init()`

Before committing to a long fit, run a 5-iteration pilot to confirm the OFV is decreasing from the initial values:

chk <- ferx_check_init(ex$model, ex$data, method = "focei")
chk$summary
#   n_iter  ofv_start  ofv_end  ofv_drop  converged
#        5     -180.2   -186.3      6.1      FALSE

Optimizer trace — `ferx_trace()` / `ferx_plot_trace()`

When optimizer_trace = TRUE is passed to ferx_fit(), ferx writes a per-iteration CSV. Read and plot it:

fit <- ferx_fit(ex$model, ex$data, optimizer_trace = TRUE)

tr <- ferx_trace(fit)          # data frame with one row per outer iteration
head(tr[, c("iter", "ofv", "grad_norm")])

ferx_plot_trace(fit)           # two-panel plot: OFV + method-specific metric
ferx_plot_trace(fit, log_ofv = TRUE)  # OFV on log scale (OFV - min(OFV))

ferx_trace() with no argument recalls the trace from the last ferx_fit() call in the current session:

ferx_trace()    # same as ferx_trace(fit) after the above

See Output Files — Optimizer trace CSV (ferx-core book) for the full column reference.

Parameter estimates — ferx_estimates()

Parameter correlation matrix — ferx_cor_matrix()

ETA-covariate correlations — ferx_eta_cov()

IWRES autocorrelation — check_diagnostics()

Goodness-of-fit plots from fit$sdtab

Structured warnings — ferx_warnings()

Run log — ferx_runlog()

Other ferx_fit fields

Convergence diagnostics — ferx_check_init()

Optimizer trace — ferx_trace() / ferx_plot_trace()

Parameter estimates — `ferx_estimates()`

Parameter correlation matrix — `ferx_cor_matrix()`

ETA-covariate correlations — `ferx_eta_cov()`

IWRES autocorrelation — `check_diagnostics()`

Goodness-of-fit plots from `fit$sdtab`

Structured warnings — `ferx_warnings()`

Run log — `ferx_runlog()`

Other `ferx_fit` fields

Convergence diagnostics — `ferx_check_init()`

Optimizer trace — `ferx_trace()` / `ferx_plot_trace()`