Bootstrap

Maturity: beta — see Feature Maturity for what this means.

The non-parametric case bootstrap estimates parameter uncertainty by resampling subjects with replacement, refitting the model to each replicate dataset, and reading bias, standard errors and confidence intervals off the spread of the estimates.

It lives in ferx-tools, not in the estimation engine: a bootstrap is many fits of one model, which is precisely the workspace boundary rule — if it calls fit() more than once, it is a tool. Nothing about the likelihood, the optimizer or the diagnostics changes; only the data each fit sees.

ferx bootstrap warfarin.ferx --data warfarin.csv --samples 200 --seed 12345 --threads 8

Why bootstrap instead of the covariance step

ferx’s standard errors are the pure R⁻¹ matrix — the same thing NONMEM reports under $COVARIANCE MATRIX=R. They are asymptotic: they assume the log-likelihood is locally quadratic and that the sampling distribution of every parameter is normal. That assumption is symmetric by construction, so an R⁻¹ interval is always centred on the estimate and always the same width on both sides.

The bootstrap assumes none of that. It also keeps working when the covariance step fails or returns a non-positive-definite matrix, which is common exactly where uncertainty matters most.

Here is the shipped examples/warfarin.ferx fit both ways — 200 samples, seed 12345:

parameter	estimate	SE (`R⁻¹`)	SE (bootstrap)	95% percentile CI	95% normal CI
`TVCL`	0.1330	0.0066	0.0081	0.1185 – 0.1506	0.1171 – 0.1488
`TVV`	7.7307	0.2341	0.2444	7.237 – 8.127	7.252 – 8.210
`TVKA`	0.7252	0.1245	0.1405	0.4675 – 1.0641	0.4499 – 1.0006
`OMEGA(ETA_CL,ETA_CL)`	0.02859	0.0128	0.0105	0.0063 – 0.0474	0.0081 – 0.0491
`OMEGA(ETA_V,ETA_V)`	0.00958	0.0043	0.0037	0.0016 – 0.0158	0.0023 – 0.0169
`OMEGA(ETA_KA,ETA_KA)`	0.34896	0.1608	0.2034	0.1092 – 0.8622	−0.0498 – 0.7477
`PROP_ERR`	0.01075	0.0009	0.0010	0.00885 – 0.01259	0.00879 – 0.01271

Most parameters agree closely, which is the expected — and reassuring — result. The last OMEGA is the interesting row, and it is the argument for the whole tool:

its bootstrap standard error is 26% larger than the asymptotic one;
its bootstrap distribution is strongly right-skewed (0.109 – 0.862 around 0.349), which no symmetric interval can express;
and the normal-approximation interval reaches below zero — an impossible value for a variance. The percentile interval cannot, because every endpoint it reports is an estimate that some replicate actually produced.

That is the diagnostic: where the two intervals agree, the asymptotic standard errors are fine; where they disagree, believe the bootstrap.

What is resampled

Whole subjects. For each ID, every record is either in or out of a replicate — never split. This is what keeps the case bootstrap valid for a hierarchical model.

A subject drawn more than once enters the fit as independent individuals: that is the point of sampling with replacement, and it means the copies must carry distinct identifiers. ferx keeps the original ID and appends #2, #3, … so a diagnostic can still be traced back. (PsN does the same thing by renumbering IDs while writing the replicate data file.)

Which parameters are bootstrapped

Every estimated population parameter: each theta, the free lower triangle of the between-subject Omega, each sigma, and — when the model declares kappa — the free lower triangle of the inter-occasion Omega, reported as OMEGA_IOV(KAPPA_CL,KAPPA_CL) and so on. Structural zeros (an off-diagonal between two etas that are not in the same block_omega) are not parameters and are not reported.

The same vector is what --update-inits hands to the replicates and what --dofv evaluates, so a model’s IOV variance is started from the base fit’s estimate and evaluated at the replicate’s, exactly like every other parameter.

With --keep-covariance, each parameter also gets a per-replicate se_ column. A cell there is empty when core reports no standard error for that coordinate — an IOV covariance, for instance, since only the IOV variances have one. Empty means “not reported”, never zero.

Mixture models cannot be bootstrapped — use SIR

A [mixture] model’s classes are identified only up to relabelling: replicate k’s “class 1” need not be the original fit’s class 1. Averaging estimates across replicates without resolving that mixes two different things, inflating the standard error toward the between-class separation — a property of the model, not a sampling uncertainty — while every replicate converges and the table looks entirely reasonable. The error grows with how well separated the classes are, so it is worst exactly when the mixture is best identified and most likely to be believed.

Relabelling rules exist and would work most of the time. That is the wrong standard for the artefact a reader trusts instead of redoing the analysis, so ferx bootstrap refuses a mixture model rather than reporting a plausible-looking table. This is a settled decision, not a missing feature.

For parameter uncertainty on a mixture model, use SIR (sir = true). SIR never re-optimises — it samples around the maximum-likelihood estimates and reweights by the likelihood — so there is no second basin to fall into and no labelling to resolve. The covariance step is likewise well defined at a single fit’s fixed labelling.

A model that reads ID as a covariate cannot be bootstrapped

Because the duplicate copies are necessarily renamed, a model whose predictions depend on the subject identifier would silently see different values than it did in the base fit. ferx refuses such a run with an error naming the fix. PsN documents the same hazard as a known bug and leaves it to the user to notice.

Reproducibility

Each replicate’s draw is derived from the master --seed and its own index, so the resampled datasets depend only on --seed, --samples, --sample-size and the strata. They do not depend on --threads, on the order replicates happen to finish in, or on whether the base model was fitted first.

This is a deliberate departure from PsN

The PsN bootstrap user guide documents the opposite as a known wart: “the results of two runs will be different even if the seed is the same if the lst-file of the base model is present at the start of one run but not the other”, because running the base model advances one shared random-number generator. Deriving per replicate removes that coupling, and makes the draw stream something a test can pin.

Stratified resampling

If the dataset mixes populations — the classic case is 10 subjects with rich profiles and 90 with sparse steady-state samples — resample within each group so that every replicate keeps the same composition:

ferx bootstrap run12.ferx --data d.csv --stratify-on STUDY

The stratification column is read straight from the dataset, so it does not have to be a covariate the model declares, and its values may be non-numeric group labels. It must take exactly one value per subject; if it does not, the run stops and names the offending subject, because an individual that belongs to two strata cannot be resampled from either.

By default each stratum contributes as many subjects as it has. --sample-size changes that in one of two ways:

# One number: strata are allocated in proportion to their size (rounded).
ferx bootstrap run12.ferx --stratify-on STUDY --sample-size 50

# Explicit per stratum — PsN's syntax. Every stratum must be listed.
ferx bootstrap run12.ferx --stratify-on STUDY --sample-size "1001=>12,1002=>24,1003=>10"

With proportional allocation the rounded per-stratum counts need not sum exactly to the requested total; in the default case (sample size = number of subjects) they always do.

Which replicates count

A replicate whose fit did not converge, or whose estimate landed on a bound, is not a draw from the sampling distribution of a converged fit — including it biases the interval. Two filters are therefore on by default, matching PsN:

filter	default	drops replicates where…
`--no-skip-minimization-terminated`	on	the fit did not converge
`--no-skip-estimate-near-boundary`	on	an estimate sits on a bound
`--skip-covariance-step-terminated`	off	the covariance step failed
`--skip-with-covstep-warnings`	off	the covariance step warned

The two covariance filters read a diagnostic that only exists when the covariance step actually ran, and replicate fits skip it by default (it is most of the cost, and the bootstrap standard error comes from the spread of the estimates, not from any per-replicate R⁻¹). Combining them without --keep-covariance is an error rather than a silent no-op.

Filters are applied when the statistics are computed, not when a replicate finishes. That is what makes it possible to change your mind afterwards:

ferx bootstrap --summarize --directory run12-bootstrap --no-skip-minimization-terminated

--summarize re-reads raw_results.csv — which carries each replicate’s diagnostics alongside its estimates — and rewrites the statistics. It refits nothing.

Watching a run

A 200-sample bootstrap is minutes to hours of fitting behind a single command, so the run reports each fit as it completes. On the command line that is a progress bar on stderr:

  [00:01:23] [===============>              ] 97/200 replicates (ETA 00:01:28)

The base model gets a spinner of its own before the bar appears — it is one fit, and a bar sitting at 0/200 through it reads as a run that is stuck. --dofv adds a second pass over every replicate, so it gets its own bar rather than extending the first one past what it promised. A --resume run counts what is left: 10 replicates missing from a 200-sample directory is a bar of 10.

The bar is drawn only when stderr is a terminal, so a run in a log file, a pipe or a scheduler is unaffected. --no-progress suppresses it for an interactive run you want quiet.

ferx bootstrap run12.ferx --data data.csv --samples 200 --no-progress

In R the same events drive a cli progress bar:

bs <- ferx_bootstrap(model, data, samples = 200, threads = 8, progress = TRUE)

The progress stream is presentation only. The draws, the fits and the artefacts are identical whether or not anything is watching, so a run in a terminal and the same run in a cron job produce the same numbers from the same seed.

Recovering an interrupted run

A 200-sample bootstrap is hours of fitting, and the two ways it goes wrong need different answers. PsN’s user guide treats them as one recovery section; ferx splits them by flag.

what happened	what to do
too many samples were filtered out	`--summarize` — nothing needs refitting, only the exclusion criteria change
the process was killed, the machine crashed, the job hit a wall clock	`--resume` — refit only the samples that never finished
some samples failed for a transient reason (out of memory, full disk) that has since gone away	`--resume --retry-failed`

# The run dies at sample 137 of 200.
ferx bootstrap run12.ferx --data data.csv --samples 200 --directory run12-bootstrap

# Pick it up. Samples 1..136 are read back; 137..200 are fitted.
ferx bootstrap run12.ferx --data data.csv --samples 200 --directory run12-bootstrap --resume

The per-replicate files are written as each replicate finishes, not at the end, so a run killed at any point leaves everything it had already completed.

A resumed run is the uninterrupted run

Not “statistically equivalent” — identical. Each replicate’s draw comes from the master seed and its own index, so replicate 137 of the resumed run is fitted on exactly the dataset replicate 137 of the uninterrupted run would have seen, from exactly the same starting estimates: the base model is reused rather than refitted, because --update-inits starts every replicate from its results and a second base fit could land somewhere marginally different. raw_results.csv comes out identical either way, wall clock aside, and that is asserted as a test rather than assumed.

PsN cannot promise this. Its guide documents one shared, order-dependent random-number generator, so a resumed PsN run draws a different set of datasets than an uninterrupted one would have.

Failed replicates are kept unless you say otherwise

A replicate whose fit errored is recorded with the error text and, by default, carried forward on a resume rather than refitted. This is PsN’s answer — its recovery instructions have you delete the run files you want redone — and it is the right default because a fit that failed on this dataset usually fails again, and paying for it a second time reaches the same place.

--retry-failed inverts that. Use it when the failure was a transient resource problem — an out-of-memory kill, a full disk, a lost network mount — rather than the model.

Note that an errored replicate is not the same thing as a filtered one. A replicate that merely failed to converge fitted fine, is in raw_results.csv with minimization_successful = 0, and is dropped by a filter at summary time; --retry-failed does not touch it. Change your mind about those with --summarize.

Resuming the wrong directory is refused

Every run writes bootstrap_run.json beside its CSVs, recording the seed, sample count, sample size, stratification column, --keep-covariance / --dofv, what the replicates were started from (--update-inits / --run-base-model), the parameter names, and SHA-256 hashes of the model file and the dataset. A --resume whose options or inputs disagree is refused, naming the field that differs:

--resume: `run12-bootstrap` was created with --seed 1, but this run asks for 42.
Resuming would mix two different runs' replicates into one raw_results.csv.
Use a fresh --directory, or drop the conflicting option.

Without that check the resumed file would carry one model’s estimates under another model’s column names — a corruption nothing downstream could detect, because the file would still parse and the statistics would still compute.

The subtlest case the check covers is --update-inits. Interrupt a default run and resume it with --no-update-inits, and every reused replicate was started from the base fit while every refitted one starts from the model file’s estimates. Both halves are perfectly good fits; the file holding them is a bootstrap of neither. So the starting point is pinned like everything else, and changing it means a fresh --directory.

A run killed mid-write leaves a partial trailing row. It is dropped and that one sample refitted; only the tail is treated this way, and a malformed row anywhere else is an error, because silently dropping a row in the middle would change the statistics without saying so.

raw_results.csv is written at full precision

Its numbers are printed as the shortest decimal that reads back to the same bits, rather than rounded to a fixed number of places, because the file is the run’s checkpoint: a resumed run starts its replicates from the base fit’s estimates as read back from it. Rounding would perturb that starting point and the resumed run would only match the uninterrupted one to about 1e-9. The statistics files (bootstrap_results.csv, bootstrap_diagnostics.csv) keep their fixed ten decimals — nothing is ever fitted from those.

Confidence intervals

Percentile intervals use PsN’s estimator, the weighted average at x₍ₙ₊₁₎ₚ:

n = number of samples + 1
p = percentile / 100
i = integer part of n·p
f = decimal part of n·p
percentile = (1−f)·xᵢ + f·xᵢ₊₁

A tail can only be reported when there is an order statistic at or below it — i.e. when (n+1)·p ≥ 1. That single condition reproduces every minimum PsN tabulates:

interval	samples needed
5% – 95%	19
2.5% – 97.5%	39
0.5% – 99.5%	199
0.05% – 99.95%	1999

Below the threshold the percentile interval is left empty rather than extrapolated past the smallest observation, and only the normal-approximation interval (estimate ± z·SE_bootstrap) is reported. PsN’s rule of thumb — and the default — is 200 samples for standard errors.

Δofv

--dofv evaluates each replicate’s parameter vector on the original dataset with no estimation (outer_maxiter = 0, exactly NONMEM’s MAXEVAL=0) and reports ofv_bootstrap − ofv_original in delta_ofv.csv.

Every value is necessarily positive: the original fit’s parameters minimise the objective on the original data, so any other vector scores worse. The useful comparison is against a χ² distribution with degrees of freedom equal to the number of estimated parameters (written to bootstrap_diagnostics.csv as chi_square_df, counted over the free parameter coordinates — so a free 2x2 block_omega contributes three, not two). In theory the bootstrap Δofv distribution should sit at or below that reference; if it does not, the bootstrap is not describing the uncertainty well and another method — SIR — is the better choice.

Output files

Written to --directory (default {model}-bootstrap). Names follow PsN so existing scripts and habits transfer.

file	contents
`raw_results.csv`	one row per fit, first row the original dataset; estimates, OFV and termination diagnostics (plus `se_` columns under `--keep-covariance`)
`bootstrap_results.csv`	one row per parameter: original, mean, bias, standard error, median, percentile CI, normal CI
`bootstrap_diagnostics.csv`	run counts, exclusion tallies, and the mean of each per-replicate diagnostic
`included_individuals1.csv`	per replicate, the subject IDs drawn (repeats included), in dataset order
`included_keys1.csv`	per replicate, the internal 1..N ordinals of those subjects
`sample_keys1.csv`	per replicate, one count per original subject: how many times it was drawn
`all_individuals1.csv`	every subject of the original dataset, in original order
`delta_ofv.csv`	`--dofv` only
`bootstrap_run.json`	the seed, options and input hashes the directory was created with — what `--resume` checks itself against

Row j refers to the same replicate in raw_results, included_individuals, included_keys and sample_keys, as in PsN.

The first four are written incrementally, one row per replicate as it finishes, so an interrupted run leaves what it had already done; the rest are written once at the end. A completed run’s files are always in sample order regardless of how its replicates were scheduled, or whether it was resumed. Starting a fresh run into a directory that already holds one overwrites it — add --resume to continue it instead.

Two files depart from PsN’s layout on purpose. bootstrap_results.csv is tidy — one row per parameter — rather than PsN’s stacked blocks of statistics, and the run-level diagnostics live in their own file instead of as a second block inside the results file. A single CSV carrying two different header shapes cannot be read by a data frame, which defeats the point of writing it.

Options

option	default	meaning
`--samples N`	200	number of bootstrap datasets
`--seed N`	1	master seed
`--threads N`	Rayon default	replicates fitted concurrently (each fit uses one thread); an exact ceiling, so peak memory is `N` fits’ worth. When unset, honors the ambient Rayon pool, including `RAYON_NUM_THREADS`.
`--sample-size M`	dataset size	subjects per replicate, or a per-stratum map
`--stratify-on COL`	—	resample within strata from a dataset column
`--directory DIR`	`{model}-bootstrap`	where the CSVs go
`--ci LEVEL`	95	two-sided confidence level, percent
`--no-update-inits`	—	start replicates from the model file’s inits instead of the base fit’s estimates
`--no-run-base-model`	—	skip the fit on the original data (disables bias, the normal CI, `--update-inits` and `--dofv`)
`--keep-covariance`	off	run the covariance step for every replicate
`--dofv`	off	compute Δofv against the original dataset
`--summarize`	off	recompute statistics from an existing run directory
`--resume`	off	continue an interrupted run in `--directory`, refitting only the samples it does not already hold
`--retry-failed`	off	with `--resume`, also refit the samples recorded as failed

Options in the PsN guide with no ferx analogue — -allow_ignore_id, -copy_data, -mceta, -rplots, -clean — are NONMEM run-directory plumbing and are not implemented. -bca (bias-corrected accelerated intervals, which need an additional jackknife pass) is not implemented yet.

Library API

use ferx_core::prepare_run;
use ferx_tools::bootstrap::{run_bootstrap, BootstrapOptions};

let prepared = prepare_run("warfarin.ferx", Some("warfarin.csv"))?;
let result = run_bootstrap(
    &prepared,
    &BootstrapOptions { samples: 200, seed: 12345, ..Default::default() },
)?;

for p in &result.summary.parameters {
    println!("{}: {:.4} (SE {:.4})", p.name, p.mean, p.standard_error);
}

prepare_run is the “load a model and its data, but do not fit” entry point: it resolves the data path against the model’s [data] block, applies [data_selection], reads the population and binds any theta NAME[...] level blocks — everything ferx model.ferx does up to the fit itself.

Cancelling a run

A 200-sample bootstrap is minutes to hours of fitting behind one call, so a library caller gets a cooperative abort: BootstrapOptions::cancel takes a CancelFlag, and setting it from another thread ends the run at the next replicate boundary.

use ferx_core::CancelFlag;
use ferx_tools::bootstrap::{run_bootstrap, BootstrapError, BootstrapOptions};

let flag = CancelFlag::new();
let options = BootstrapOptions { cancel: Some(flag.clone()), ..Default::default() };

// … from another thread, on a Ctrl-C or a user's "stop" button:
flag.cancel();

match run_bootstrap(&prepared, &options) {
    Ok(result) => report(result),
    Err(BootstrapError::Cancelled) => eprintln!("stopped; --resume picks it up"),
    Err(e) => eprintln!("Error: {e}"),
}

Three things follow from where the flag is checked.

A cancelled run is not a failed one. The flag also reaches every replicate’s own fit, so a replicate that was mid-fit unwinds with an error of its own rather than running to completion. Those are dropped, never written to raw_results.csv — a recorded failure is one that a resume carries forward instead of refitting, which would silently cost the run a replicate. What the flag stops is scheduling: the draws that had not started are not fitted at all.

Everything already finished is on disk. The journal flushes each replicate as it lands, so a cancelled run leaves the same directory an interrupted one does, and --resume continues it into exactly the run that was cancelled. The statistics files are the exception: a cancelled run does not write an interval over a truncated set of replicates.

A cancel before the run starts is a no-op on the directory. The flag is read before the journal is opened, and opening it rewrites the per-replicate files — so cancelling a mistyped command cannot destroy the artefacts of the run that is already there.

ferx bootstrap does not install a signal handler yet, so the CLI has nothing to set the flag from; on the R side, ferx_bootstrap() currently runs the bootstrap on the R main thread and so cannot poll for an interrupt either. Both are tracked separately — this is the engine half.