Skip to contents

Changelog

Source: NEWS.md

epicurves 0.0.0.9000

Parametric bootstrap for the general engine

ga_bootstrap() simulates data sets from a ga_fit() result with ga_simulate(), refits them from the estimates and returns the replicates of any statistic of the fit, with percentile and basic intervals (confint.ga_boot()). Missing values of the original data stay missing, and replicate seeds are drawn from the current stream, so results do not depend on the number of cores. To support fits with time-varying rates, ga_simulate() gained a varying argument, and ga_fit() now returns the optim convergence code.

Simpson quadrature of Q_k fixed at the start of each interval

sird_ga_transition() stored, as the first node of every interval after the first, the resolvent of the previous interval instead of the identity: deSolve records the pre-event value at an event time, so the row at t_{k-1} held Phi(t_{k-1}, t_{k-2}). That node enters Simpson’s rule with weight h/3, which made the error of sird_ga_Qk() decrease only in proportion to the step (2.7% with 10 substeps, 0.34% with 80). With the node reset, the error is about 1e-5 against the joint RK4 reference, limited by the lsoda tolerance. Only the prevalence path (sird_ga_loglik(), sird_ga_fit()) used this function; the incidence path and ga_*() never did, so no result of the articles changes.

Generic compartmental engine rewritten (ga_model() and friends)

ga_model(), ga_event(), ga_incidence(), ga_prevalence(), ga_loglik(), ga_simulate(), ga_fit() and ga_profile() replace the earlier closure-based engine (new_jump_model(), general_ga_*(), seird_model(), sird_model()), which is removed. The new engine (src/ga_generic.cpp) handles any mass-action model with all compartments explicit, integrates the mean, the resolvent and the Lyapunov equation of the state-noise covariance jointly by RK4 (the old engine obtained Q_k by Simpson quadrature on the stored ODE grid; the relative errors of 0.3% to 3% that check_Qk.R first reported for it came from the endpoint defect described above, not from the quadrature), observes event incidence through cumulative event counters in an augmented state, and allows piecewise-constant rate parameters. On SIRD with daily cases and deaths it reproduces sird_ga_loglik_inc() to within 2e-10 in the log-likelihood, for constant and time-varying transmission (test-ga-generic.R). ga_fit() uses a Nelder-Mead search with relative tolerance 1e-8 followed by BFGS; on a SEIRD example this reached the same maximum as two long Nelder-Mead runs with about one fifth of the likelihood evaluations.

Ill-conditioned innovation covariance no longer aborts incidence fits

The C++ Kalman log-likelihood (src/sird_ga_fast.cpp) called arma::inv_sympd() on the innovation covariance after only a determinant check. When that matrix passed the check but was ill-conditioned, the call threw and aborted the surrounding optim(), killing whole fits (first seen in the reporting-probability profile of the time-varying Paraíba fit at p = 0.03). The non-throwing overload is now used and the evaluation returns the -Inf log-likelihood sentinel, so optimizers penalize that point and continue. Fits that never hit the exception are unchanged.

Independent verification of Article 1’s inherited simulation tables

The Article 1 Scenario 1-5 SIRD-GA/SIR tables were transcribed from the thesis rather than regenerated by this package. A reduced, disclosed spot check (J=10 replicates, not the manuscript’s J=100, owing to real cost: ~70-130s per fit at N=20,000 since Nelder-Mead needs the full maxit=1000+ budget to converge reliably) was run for Scenario 1 at N=20,000 using sird_ga_fit() directly, with early-extinction trajectories discarded exactly as the manuscript’s own stated protocol requires (an unfiltered first attempt produced a spurious, wildly diverged fit, which traced to an extinct simulated series, not a package bug - see below). Result: mean estimates centred near the true parameters (lambda 0.637 vs true 0.6, gamma 0.294 vs 0.3, alpha 0.131 vs 0.1, p 0.545 vs 0.5, tau2 0.444 vs 0.4), with MSEs roughly one order of magnitude above the manuscript’s J=100 table, consistent with the much higher Monte Carlo noise expected at J=10 (3 of 10 replicates did not fully converge within the iteration cap and were included as-is). This confirms the inherited table is in the right neighbourhood and free of gross implementation bugs, without claiming a full independent replication.

Initial implementation, built to validate and reproduce the results of two companion manuscripts derived from the thesis.

Findings made while validating against the thesis

  • Richards growth-curve sign correction. The thesis’s Richards-curve equation (Chapter 3), transcribed literally, uses exp(+beta*gamma*t); this drives the curve to 0 as t -> Inf instead of converging to the carrying capacity K, inconsistent with the standard Richards (1959) form and with the thesis’s own (correct) Gompertz/Logistic equations. richards_curve() implements the corrected form, exp(-beta*gamma*t). This needs to be corrected in the Article 2 manuscript’s methods section.
  • Partition-point over-segmentation. The flat “50% MSE reduction” split-confirmation rule described in the thesis can over-segment long, structurally stable phases, both in a re-simulation of the thesis’s own Scenario 1 and in the real Paraíba death series. fit_piecewise_ensemble() exposes a confirm_criterion = "aic" option (a BIC-style, parameter-count-penalised alternative) that is markedly more conservative while still recovering the Scenario 1 change points to the same accuracy as the thesis’s own reported values. Recommended for Article 2 as a robustness improvement, and as a direct response to the 2023 referees’ novelty concerns about the (non-piecewise) precursor paper.
  • Generalized compartmental engine + SEIRD-GA. new_jump_model() / general_ga_transition() / general_ga_loglik() / general_ga_fit() generalize the hand-derived SIRD-GA machinery to arbitrary compartmental structure (numerical rather than analytic Jacobian). Cross-validated against the original SIRD-GA engine on identical data (log-likelihoods agree to within 0.03%). seird_model() demonstrates a four-compartment extension (adds an Exposed compartment); MLE on simulated data recovers R0 exactly (2.500 vs true 2.500) and lambda/gamma/alpha within 18%.
  • Permutation-test partition criterion, and its limits. fit_piecewise_ensemble(confirm_criterion = "permutation") implements a formal, multiplicity-corrected hypothesis test (permutation_test_partition()), following Matteson & James (2013)’s e.divisive logic but specialized to parametric growth curves. Found to be underpowered at a computationally feasible replicate count (B=29): p=0.23 for a genuine, obvious change point that both heuristic criteria detect easily. Documented as a genuine cost/power trade-off in Article 2, not hidden.
  • External nonparametric detectors over-segment this data class. Benchmarked ecp::ks.cp3o and ecp::e.divisive against the piecewise ensemble: both are 2-4 orders of magnitude faster but substantially over-segment smoothly-varying growth curves (e.divisive: 8 spurious interior points on Scenario 1’s 3-true-changepoint series), because they assume approximately piecewise-stationary data and growth phases are never actually stationary. A real, reportable finding, not a configuration error.
  • SIRD-GA real-data identifiability: two basins. Fitting SIRD-GA to raw daily case counts is sensitive to day-of-week reporting artefacts; a 7-day centered rolling average (smooth_series()) is recommended before estimation. On the smoothed Paraíba series, a 5-start sweep (sird_ga_fit_multistart()) reveals two distinct optima: a “fast-dynamics” basin (4/5 starts, loglik ~2947, R0~1.03, mean infectious period ~1.4 days, epidemiologically implausible) and a “slow-dynamics” basin (1/5 starts, loglik ~2924, R0~1.03, mean infectious period ~4.5 days, matching the published estimates and known COVID-19 natural history). Within the fast basin, gamma and alpha trade off freely while gamma+alpha stays fixed, mirroring the thesis’s own appendix finding that gamma/alpha are estimated far less precisely than lambda or R0. Recommended for Article 1: report the multi-start summary and prefer the slow-dynamics basin on epidemiological grounds, rather than reporting an unconstrained MLE point estimate. See the SIRD-GA vignette for the full table and discussion.
  • Article 1’s Scenario 2/3 had silently substituted the thesis’s appendix design. The thesis’s main-text Scenario 2/3 (Chapter 2) vary only the reporting rate p in {0.7, 0.9} at fixed (lambda,gamma,alpha,tau)= (0.6,0.3,0.1,0.4), isolating under-reporting severity on its own. The manuscript draft’s “Scenario 2/3” instead reused the thesis appendix’s distinct full-parameter regimes (different R0 values entirely), losing the thesis’s actual finding that reporting-rate severity, on its own, is not the dominant driver of estimation precision. Fixed: the appendix regimes are now correctly labelled Scenario 4/5, and the real p-sensitivity Scenario 2/3 (SIR baseline and SIRD-GA) were inserted with exact numbers read from capitulo2.tex. A genuine table/figure-ordering bug (Scenario 1’s table and confidence-band figure left physically stranded inside the “Scenarios 4-5” subsection) was caught and fixed in the same pass.