v0.4.0

[jtimonen]

Summary

This PR prepares bmstate 0.4.0, a correctness-focused minor release addressing statistical, prediction, numerical, pharmacokinetic, and performance issues found during a package audit.

The most consequential defect placed the empirical baseline-hazard prior location on the wrong scale. Other fixes affect joint posterior-draw alignment, spline identifiability, transition-probability prediction, state-visit probabilities, path simulation, and PK calculations.

Rerun requirement

Fits and experiments produced with affected bmstate versions must be rerun before their numerical results are reported or cited.

The main reasons are:

• natural-scale baseline-hazard reference rates were passed directly to a parameter representing a log hazard;
• spline centering changes the identifiable representation of the baseline hazard;
• posterior prediction could combine parameters from different joint posterior draws; and
• prediction and path-simulation defects could change estimated event risks.

For example, a natural-scale reference rate of 0.001 was previously used as the location of log_w0, centering the conditional median hazard near exp(0.001) = 1.001 rather than 0.001. This was a real scale error, not merely a misleading variable name.

Statistical and prediction corrections

• Correct the baseline-hazard prior scale. The public R API stores positive natural-scale reference rates, whereas Stan models the baseline intercept on the log-hazard scale. Stan-data construction now applies log() exactly once at the R–Stan boundary, and the data field is renamed from mu_w0 to mu_log_w0.

• Make empirical hazard calibration explicit. The former Cox/Breslow-based summary could silently omit transitions with no events and then average transition-specific log rates within a transition type. Calibration now pools observed event counts and exact at-risk time within transition type and uses the resulting count/person-time rate as the natural-scale prior reference location. It is documented as a reference rate, not an arithmetic prior mean.

• Reject unsupported transitions before fitting. Every transition retained in a non-PK-only hazard model must have at least one observed occurrence and positive at-risk time. Validation now stops before Stan is called and names the affected transitions. There is no half-event regularization or prior-only estimation of an unobserved transition.

• Keep prediction-data validation separate from fitting. Out-of-sample subjects need not experience or be at risk for every transition represented by the already fitted model. Generic Stan-data construction therefore permits absent transition rows for prediction, while fit_stan() applies the stricter support validation before fitting. Stan fitting retains positive event and risk counts and contains no zero-count likelihood branches.

• Preserve joint posterior-draw alignment. Baseline-hazard arrays were previously repeated in a different draw/subject order from hazard multipliers and prediction metadata. This could combine spline coefficients and intercepts from one posterior draw with covariate effects from another. All prediction inputs now use one explicit draw-major ordering.

• Return explicit zero state-visit probabilities. p_state_visit() previously constructed results only from observed event rows, silently omitting requested subjects or groups with no visit to an event state. It now completes the state-by-group grid, computes denominators separately, and returns explicit zero counts and probabilities.

• Separate spline shape from baseline level. Because the intercept-containing B-spline basis partitions unity, adding a common value to every spline coefficient and subtracting it from log_w0 left the likelihood unchanged. Realized spline coefficients are now centered to sum to zero, leaving the baseline intercept to carry the constant level. The raw hierarchy remains full-dimensional, so this removes likelihood-level confounding rather than every prior-only nuisance direction.

Numerical and simulation corrections

• Keep simulated paths inside their requested horizon. Applying min_t_step could move a genuine event to or beyond t_max, after which the row was labeled as censoring while retaining a nonzero transition index. Genuine pre-horizon events are now preserved, censoring rows always use transition index zero, reversed horizons are rejected, and zero-length horizons return a valid one-row censored path.

• Retain the intended fast thinning approximation. Path simulation continues to use the maximum hazard evaluated on a 100-point grid, inflated by 1.05, and skips transitions whose inflated envelope is at most 1e-9. Bounds are now calculated only for transitions possible from the current state, non-finite bounds are rejected, and the approximation is documented explicitly. This is not described as a mathematically guaranteed envelope.

• Construct exactly the declared integration grid. For noninteger model horizons, the old seq() construction could create more values than the declared N_grid dimension. The grid is now formed directly from exactly N_grid midpoint indices.

• Honor delayed prediction start times. solve_trans_prob_matrix() previously initialized the forward equation at the first requested output time rather than at t_start. It now includes t_start in the integration times, returns only the requested output rows, and validates dimensions even for a zero-duration identity result.

• Return the requested PK interval width. The PK summary helper divided the tail probability by two twice, so an 80% request returned a 90% interval. Quantiles now use (1 - ci_alpha) / 2 and its complement directly.

Pharmacokinetic corrections

• Stabilize the oral one-compartment solution near equal rates. Direct Bateman-function evaluation divided by ka - ke and could produce cancellation, 0 / 0, or 0 * Inf. The implementation now uses sign-aware expm1() factorizations and the continuous equal-rate limit for transient and steady-state calculations.

• Limit parameter-dependent branching to the numerical PK guard. Structural and validation branches introduced elsewhere depend only on data or model dimensions. The equal-rate PK branch is parameter-dependent, but its limiting expression is value- and first-derivative-matched at equality to minimize discontinuity while avoiding the singular direct formula.

• Tighten PK input validation. The package now rejects nonpositive absorption, clearance, volume, dosing-interval, concentration-bound, and exposure-normalization values before they reach numerically invalid operations.

Stan gradient performance

Hazard covariates are time invariant, but the previous Stan program expanded them over every interval and rebuilt the complete hazard-multiplier matrix inside each transition loop.

The revised program differentiates one predictor per subject and transition type and then gathers the required values at event and at-risk indices. This is algebraically equivalent but constructs a substantially smaller reverse-mode graph.

In a representative benchmark with 500 subjects, 10,000 intervals, six transitions, two transition types, and four covariates, repeated gradient time decreased from approximately 1.93 seconds to 1.19 seconds, a reduction of about 38%. Log densities agreed exactly and maximum gradient differences were below 5e-13. No measurable regression was found when every interval belonged to a different subject.

Additional edge-case fixes

• Base data frames without covariates no longer lose dimensions through single-column subsetting.
• Paths without terminal events are returned unchanged instead of failing during truncation.
• Full posterior draws are returned correctly when no draw name is supplied.
• Malformed delayed-start requests and events on nonpositive-duration intervals are rejected explicitly.
• Simulated paths are accumulated in a list and bound once, avoiding repeated output reallocation.

Compatibility

• Package version is bumped from 0.3.3 to 0.4.0.
• The public R API continues to accept natural-scale baseline-hazard reference rates.
• Custom Stan programs selected through options(bmstate_stan_file = ...) must rename mu_w0 to mu_log_w0 and interpret it as a log-scale location.
• Fits produced with the incorrect prior transformation must be rerun.

Validation

• The complete package test suite passes.
• After the final transition-validation adjustment, the targeted correctness regression suite passes 61/61 tests.
• Stan parses successfully with pedantic warnings limited to the documented numerical guards and existing hard parameter bounds.
• git diff --check passes.
• Package vignette and site rendering are left to CI.

Known limitations

• The retained 100-point thinning envelope is a performance approximation and may miss a sufficiently narrow spline peak. The runtime check detects encountered violations but is not a proof that the envelope dominates everywhere.
• The absolute 1e-9 simulation cutoff is deliberately retained for speed and is meaningful only relative to the package’s intended time scale.
• Midpoint hazard integration can be inaccurate on a coarse grid when hazards vary rapidly; substantive analyses should include a finer-grid sensitivity check.
• The half-normal (0, 3) log-intercept scale remains very broad and implies no finite marginal arithmetic mean hazard after integrating over the scale.
• Centering removes the spline/intercept likelihood redundancy, but a reduced-rank contrast parameterization could further remove prior-only common-shift directions.