Blood Cancer Projections

Incidence and prevalence projection methods for Australian blood cancers, using age–period–cohort incidence models and a modified counting-method for prevalence, with parameter-level Monte Carlo uncertainty. The framework is disease-agnostic; the current implementation covers lymphoma and is designed to extend to other blood cancers (each cancer added as its own module).

---

Lymphoma (current implementation)

Code and data accompanying the manuscript:

Irving A, Luo Q, Petrie D, Fanning L, Li JJ, Ghijben P, Chung E, Wellard C, Waters N, McQuilten ZK, Wood EM, Williams J, Watt R, Winton S, Opat S, Barraclough A, Cheah CY, El-Galaly TC, Bishton MJ, Hawkes EA. Lymphoma incidence and prevalence projections in Australia to 2045: a statistical modelling study. [Journal]; submitted [date].

This repository contains the R code, source data, and modelled outputs for projections of incidence (to 2045) and prevalence (2-, 3-, 5-, 10-, and 40-year horizons, 2012 to 2045) for Australian lymphoma, in a two-tier structure: a headline tier of Hodgkin lymphoma (HL) and aggregate non-Hodgkin lymphoma (NHL) — disjoint groups whose sum is the total-lymphoma denominator — and a decomposition tier of three NHL subtypes, diffuse large B-cell lymphoma (DLBCL), follicular lymphoma (FL), and mantle cell lymphoma (MCL), that are a subset of aggregate NHL. All by sex (and persons for the totals). The aggregate NHL/HL models are fitted from 1990 (selected by a start-year holdout sensitivity) and the subtype models from 2003; observed incidence is shown from 1982.

Uncertainty intervals throughout are 95% credible intervals from parameter-level Monte Carlo simulation (B = 1,000 draws). Methodology is summarised below; detailed implementation is in the comments of code/apc_model.R, code/prev_model.R, and code/supplement.R.

Citation

If you use this code or these projections, please cite both the paper and the archive:

[Paper citation, with DOI on publication]

Irving A, Luo Q, Petrie D, Fanning L, Li JJ, Ghijben P, Chung E, Wellard C, Waters N, McQuilten ZK, Wood EM, Williams J, Watt R, Winton S, Opat S, Barraclough A, Cheah CY, El-Galaly TC, Bishton MJ, Hawkes EA. Blood Cancer Projections (lymphoma): incidence and prevalence projections for Australia, 2022–2045. Zenodo. DOI: 10.5281/zenodo.22785391

Repository structure

.
├── code/                              R scripts: the pipeline (run in this order)
│   ├── _setup.R                       Shared constants, colours, figure helpers, save_fig, summarise_draws
│   ├── import_data.R                  AIHW/ABS workbook → CSV transformation (raw workbooks not in repo)
│   ├── apc_model.R                    Fit APC models, project incidence with MC uncertainty; Figure 1, Tables 1 and S2, Table 1 change CrIs
│   ├── prev_model.R                   Compute prevalence (2/3/5/10/40-year horizons); Figures 2–3, Table 2
│   ├── supplement.R                   AIHW validation, SI tables, SI figures
│   ├── _tests.R                       Structural test harness (run after edits)
│   ├── reproduce.R                    One-click reproduction for external reproducers: snapshot the committed files, run, compare
│   ├── supporting/                    Standalone analyses whose results are reported in the paper or SI
│   │   ├── holdout.R                  Temporal incidence holdout, also the scoring window for selection (fit to 2011, predict 2012–2021); SI §S1.1
│   │   ├── knot_selection.R           APC knot-count grid, selected by holdout error + parsimony; SI §S1.1
│   │   ├── start_year_sensitivity.R   Aggregate fit-start-year selection by holdout error; SI §S1.3
│   │   ├── damping_sensitivity.R      Damping-rate robustness for the aggregate tier; SI §S1.3
│   │   ├── covid_sensitivity.R        COVID (2020–21) exclusion robustness check; Methods and SI §S1.3
│   │   └── hl_ageing.R                HL age-at-diagnosis shift; Discussion
│   └── diagnostics/                   One-off checks for co-author questions; nothing in the paper depends on them
│       ├── fan_chart_fl_5yr.R         Full uncertainty fan for 5-year FL prevalence
│       └── male_nhl_decomposition.R   Male aggregate NHL rise decomposed into modelled subtypes and the remainder
├── data/                              Source data (derived from public AIHW tables)
│   ├── incidence_agg.csv              Aggregate NHL/HL incidence, 1982–2021
│   ├── incidence_agg_proj.csv         AIHW aggregate projections, 2026–2035
│   ├── incidence_subtype.csv          Subtype incidence, 2003–2021
│   ├── incidence_subtype_back.csv     Back-estimated subtype incidence, 1982–2002 (regenerated by prev_model.R)
│   ├── pop_hist.csv                   Australian population, 1982–2021
│   ├── pop_proj.csv                   Australian population projections, 2022–2045
│   ├── prevalence_agg.csv             AIHW aggregate prevalence (Book 6): 1/5/10-yr for 2000–2025 (2022–2025 are AIHW projections, input only; validation filters to 2012–2021) and 40-yr for 2021
│   ├── prevalence_subtype.csv         AIHW subtype prevalence (5/10-yr), 2012–2021
│   └── survival.csv                   AIHW observed survival by subtype + aggregate NHL, 2007–2011 / 2012–2016 / 2017–2021
├── output/                            Modelled outputs (generated by the scripts)
│   ├── apc_results.rds                Fitted APC models with β draws and projected case arrays
│   ├── incidence_projections.csv      Annual incidence projections, 2022–2045 (with 95% credible intervals)
│   ├── incidence_historical_asr.csv   Historical observed ASR, 1982/2003–2021
│   ├── apc_effects.csv                APC age/period/cohort effect estimates with credible intervals
│   ├── prevalence_projections.csv     Prevalence projections, all durations × series (incl. aggregate NHL + total lymphoma) × sexes/persons × years (base scenario)
│   ├── prevalence_sensitivity.csv     Prevalence under three survival improvement scenarios
│   ├── prevalence_validation_10yr.csv 10-year prevalence 2012–2021 (used in Table S8)
│   ├── table_3_prevalence_combined.csv Combined-duration prevalence (persons), HL / NHL / DLBCL / FL / MCL / Total lymphoma, with draw-level 95% CrIs
│   ├── table_*.csv                    Manuscript tables (1–3, plus the change-CrI tables) and SI tables (S1–S12); table_s_*.csv and table_holdout/hl_ageing are repository-only diagnostics
│   └── figure_*.{pdf,png}             Manuscript figures (1–3) and SI figures (S1–S13); figure_holdout_validation and figure_hl_ageing are repository-only
├── dashboard/                         Static Quarto dashboard published to GitHub Pages
│   ├── index.qmd                      Dashboard source (Observable JS; runs entirely in the browser)
│   ├── _quarto.yml                    Quarto project and format configuration
│   ├── styles.css                     Dashboard styling
│   ├── export_dashboard_data.R        Builds the public dashboard dataset from committed outputs
│   ├── data/                          Generated public release CSV and its data dictionary — do not edit by hand
│   ├── assets/logos/                  Organisation and funder logos
│   └── README.md                      Dashboard development, preview and review notes
├── .github/workflows/                 GitHub Actions
│   └── publish-dashboard.yml          Renders dashboard/ and deploys it to GitHub Pages on push to main
└── README.md                          This file

How to reproduce

1. Install R and required packages

R ≥ 4.3 (developed on R 4.5.2). Required packages:

install.packages(c("tidyverse", "Epi", "MASS", "cowplot", "gtable", "scales"))

The published results were produced with R 4.5.2 and Epi 2.61 (seed 20260507, B = 1,000 draws). Other versions run the pipeline, but a refit under a different R or Epi version can move the Monte Carlo results slightly: with B = 1,000 the Monte Carlo standard error is roughly ±4 on a median and ±8 on a credible-interval bound for the Table 1 totals, and larger for MCL. code/reproduce.R reports whether a difference is of that size.

2. Set the working directory

All scripts use paths relative to the repository root (e.g. data/incidence_subtype.csv). Set R's working directory to the repository root before running anything.

3. Run the pipeline

Reproducing the published results (recommended): on a fresh download, restart R (Session → Restart R in RStudio), set the working directory to the repository root, open code/reproduce.R and click Source. It moves the committed output/ (and copies data/) into _reproduce/committed/, runs the whole pipeline and the test harness, and prints a comparison against the snapshot: md5 identity of the key outputs, the point rate ratios in apc_effects.csv (which involve no Monte Carlo), a side-by-side table of the headline Table 1 and Table 3 values, and an interpretation. Set rebuild_data <- TRUE at the top to rebuild data/*.csv from the raw workbooks first (they must be in data/raw/; see "Files required to recreate the analysis from raw data").

Two traps this script exists to avoid. First, the committed output/ is the published baseline: a fresh download already contains the published numbers, so checking the manuscript against a download that has not been re-run proves nothing. Second, each pipeline file follows a library + driver pattern: sourcing the file defines functions silently with no side effects, while running it as a script (or calling its run_*() driver) executes the pipeline. In RStudio, clicking Source or Run All on apc_model.R, prev_model.R or supplement.R therefore does nothing visible; call the driver, use Rscript, or use reproduce.R. All three drivers read upstream inputs from disk, so a fresh R session is supported.

Batch:

Rscript code/apc_model.R     # APC fits + MC incidence + Figure 1 + Tables 1, S2
Rscript code/prev_model.R    # Back-estimation + MC prevalence + Figures 2, 3 + Table 2
Rscript code/supplement.R    # SI Tables S1, S3–S12 + SI Figures S1–S13 (Table S2 is written by apc_model.R)

Or in a single command:

Rscript code/apc_model.R && Rscript code/prev_model.R && Rscript code/supplement.R

On Windows Rscript may not be on the path; use its full path, for example "C:\Program Files\R\R-4.5.2\bin\Rscript.exe" code/apc_model.R.

Interactive (RStudio or R REPL): sourcing alone defines the functions; the driver must be called.

source("code/apc_model.R");  res <- run_apc_model()    # B = 1000 by default; also writes the Table 1 change CrIs
source("code/prev_model.R"); pr  <- run_prev_model()
source("code/supplement.R"); sp  <- run_supplement()

The scripts in code/supporting/ are standalone: each sources apc_model.R (and, for the sensitivity grids, supporting/holdout.R) for the model library and is run from the repository root, for example Rscript code/supporting/holdout.R. They are not part of the apc_model → prev_model → supplement chain; their results are fixed inputs to the paper (the holdout validation, the knot, start-year, damping and COVID sensitivities, and the HL ageing analysis) and are re-run only when the model changes. code/diagnostics/ holds one-off checks made for co-author questions; nothing in the manuscript or SI is generated from them. dashboard/export_dashboard_data.R builds the public dashboard dataset from committed outputs.

Pipeline order matters: apc_model.R produces output/apc_results.rds that prev_model.R reads; prev_model.R writes the prevalence CSVs that supplement.R reads. Each driver also accepts a B parameter controlling the number of Monte Carlo draws — B = 1000 is the production setting; B = 200 reduces every working array proportionally and brings the full pipeline runtime to roughly 15 seconds for fast iteration on figures or table styling.

4. Run time

Approximately 90 seconds for the full pipeline at B = 1,000 on a current laptop. The bulk of the time is in prev_model.R (8 chunks × 3 scenarios × 5 durations of per-draw prevalence). Peak memory ~660 MB; output disk usage ~17 MB.

5. Test script

After any change to apc_model.R, prev_model.R, or supplement.R:

Rscript code/_tests.R

The script runs structural assertions (file existence, output dimensions, sign and ordering of CIs, parameter reconstruction, runtime budget, reproducibility under fixed seed, and the shape and interval ordering of table_1_change_cri.csv). The first failing stopifnot() halts execution with a readable message. From a clean cache the script re-runs the relevant pipeline phase as needed; with up-to-date outputs it completes in ~5 seconds.

code/import_data.R is included to document how the CSVs in data/ were derived from the raw AIHW and ABS workbooks. It is not part of the reproduction pipeline — the raw workbooks are excluded from the repository (see .gitignore). Unlike the pipeline files it has no auto-run guard: sourcing it runs it and rewrites data/*.csv, which is what reproduce.R relies on when rebuild_data <- TRUE.

Methodology

Incidence is modelled with age–period–cohort (APC) models using the Epi package (apc.fit() with natural cubic splines at case-weighted quantiles, Poisson log-linear link with log-person-years offset). Knot counts for age, period and cohort are selected by holdout error on 2012–2021 with a parsimony preference (code/supporting/knot_selection.R; the same window also selects the aggregate start year, so the holdout results reported in SI §S1.1 are backtesting used for selection, not an independent test): 5/5/5 for the aggregate tier (NHL, HL; the package default) and 4/5/5 for the subtype tier (DLBCL, FL, MCL; fewer age knots suit the shorter series). Period and cohort effects are damped toward a relative risk of 1.0 at 8% per year beyond the last observed knot, following Sasieni and Luo et al.

Prevalence is computed by the modified counting method:

P_T(Y) = sum_{t=0..T-1} I(Y-t) × S(t)

with S(0) := √S(1) (geometric mean of S(0)=1 and S(1) under constant within-year hazard). Survival is observed for years 1–5 from AIHW and extrapolated to year 40 by constant conditional survival anchored on √(S(5)/S(3)).

Uncertainty is quantified by parameter-level Monte Carlo simulation:

• APC parameters sampled from MVN(coef(fit$Model), vcov(fit$Model)) per fit;
• Survival sampled from cell-specific logit-normal distributions centred at the clipped AIHW point estimate, with scale inferred from the published logit confidence-limit span, so the simulated marginal limits approximate rather than reproduce the published bounds. Suppressed cells are filled by .fill_nearest() from one donor age per (subtype, sex, period) cell, the nearest age group complete across all five horizons: a fully suppressed age group takes the donor's whole curve and shares its draws by reference; a partially suppressed age group keeps its published values and extends the tail by the donor's conditional ratios anchored on its own last published value, with limits scaled by the same donor point ratio (an approximation, not a calibrated interval), and is sampled from its own point and limits. Each follow-up year is sampled independently, so sampled survival paths are not constrained to decline with follow-up; _tests.R reports the share of rising paths (71.55%) and the count of rising deterministic point curves (1,631 of 5,760) as regression guards rather than gates, and the manuscript Limitations disclose both;
• Prevalence assembled from per-draw incidence and survival via the counting method, producing 1,000 prevalence trajectories per (subtype, sex, scenario, duration).

These intervals propagate parameter uncertainty conditional on the fitted Poisson model and the specified survival distributions. They are intervals for expected counts rather than prediction intervals for realised counts, and they are not adjusted for extra-Poisson variation (deviance/df reaches 3.9 in output/table_s2_apc_fit_stats.csv) or for model-selection uncertainty. Reported point estimates are medians across draws and need not equal deterministic plug-in estimates. The mean is also written to output CSVs (as *_mean) for analyses that want the bias-aware expectation; it is biased upward by exp(σ²/2) on small-count cells. 95% credible intervals are 2.5/97.5 percentiles across draws. Fixed seed: 20260507.

Prevalence is reported in a two-tier structure, consistent with incidence. The headline tier is Hodgkin lymphoma and aggregate NHL (using aggregate NHL survival from AIHW Book 11f1); these are disjoint and together are all lymphoma, so the headline total lymphoma = HL + aggregate NHL, formed at the draw level (total_b = HL_b + NHL_b per Monte Carlo draw, then median and 2.5/97.5 percentiles; the two come from independent APC fits). The decomposition tier is DLBCL, FL and MCL, a subset of aggregate NHL (they cover ~53–73% of it, falling over time) that decomposes it and is never summed into the total (that would double-count). Aggregate NHL is therefore a headline result; its AIHW validation is retained but reclassified as validation of a headline result.

The pipeline is validated against AIHW published values in Supporting Information Tables S4, S8 and S9 (Table S4 is a comparison with AIHW's own projections, not a validation; the incidence holdout in SI §S1.1 is temporal backtesting whose assessment window also informed model selection, not an independent test). Incidence agrees within a few per cent of AIHW for both case counts and ASRs; 5-year prevalence agrees within ±2% (HL, FL), ±3% (DLBCL, aggregate NHL) and ±6% (MCL); 10-year within ±2% (HL, FL, aggregate NHL), ±3% (DLBCL) and ±5% (MCL); Hodgkin lymphoma 40-year prevalence is within 1% of the AIHW aggregate value and aggregate NHL 40-year within ~4–9%.

Three asymptotic approximations are worth noting for anyone reading the code. The multivariate normal distribution used for APC parameter draws is an asymptotic approximation to the maximum-likelihood posterior. The logit-normal distribution used for survival cells is centred at the clipped published point estimate and cannot in general reproduce both asymmetric published limits: among the 1,095 published cells, about 6% have an implied limit more than one percentage point from its published counterpart (at most 9.1 points), concentrated in boundary cells such as extreme-low-survival cells in the oldest mantle-cell age band and saturated cells where AIHW reports a 95% upper bound of 1.000. Per-draw clipping in the long-horizon extrapolation produces a Monte Carlo median roughly 10% higher than a single-shot deterministic computation for the smallest subtype × sex × duration cell (40-year mantle-cell female prevalence); reported estimates are Monte Carlo medians and need not equal deterministic plug-in estimates. These limitations are discussed in the manuscript Discussion (Limitations).

Output column conventions

MC summaries in the output CSVs use the suffixes *_mid (median; the headline number), *_p025 and *_p975 (2.5th and 97.5th percentiles of the 95% credible interval), and *_mean (mean across draws; not used as the headline). Applies to cases, asr, rr, and prev. By construction, *_p025 ≤ *_mid ≤ *_p975.

Reproducing tables and figures

Manuscript / SI item	Generated by	Output file
Manuscript Figure 1 (incidence)	apc_model.R	output/figure_1_combined_projections.{pdf,png}
Manuscript Figure 2 (5-year prevalence)	prev_model.R	output/figure_2_prevalence_5yr.{pdf,png}
Manuscript Figure 3 (40-year prevalence)	prev_model.R	output/figure_3_prevalence_40yr.{pdf,png}
Manuscript Table 1 (incidence summary, incl. Total lymphoma = HL + NHL)	apc_model.R	output/table_1_incidence_summary.csv
Manuscript Table 2 (prevalence summary, all durations)	prev_model.R	output/table_2_prevalence_summary.csv
Manuscript Table 3 (combined-duration prevalence, incl. Total lymphoma)	prev_model.R	output/table_3_prevalence_combined.csv
Manuscript Table 1 change CrIs	apc_model.R (step 10 of run_apc_model(); run_table1_change_cri() can also be called after source("code/apc_model.R") to rebuild the table from output/apc_results.rds without refitting)	output/table_1_change_cri.csv
Manuscript Table 2 change CrIs	prev_model.R	output/table_2_change_cri.csv
SI Table S1 (disease definitions, ICD-10 and ICD-O-3)	supplement.R	output/table_s1_disease_definitions.csv
SI Table S2 (APC fit statistics)	apc_model.R	output/table_s2_apc_fit_stats.csv
SI Table S3 (back-estimation validation)	supplement.R	output/table_s3_back_estimation.csv
SI Table S4 (aggregate incidence vs AIHW projections, a forecast comparison)	supplement.R	output/table_s4_incidence_aihw.csv
SI Table S5 (subtype sum as % of aggregate NHL incidence)	supplement.R

Data sources

All source data are derived from publicly available Australian Institute of Health and Welfare (AIHW) tables:

• Cancer Data in Australia (incidence, mortality, AIHW projections): https://www.aihw.gov.au/reports/cancer/cancer-data-in-australia
• Blood Cancer Framework (subtype incidence and survival): https://www.aihw.gov.au/reports/cancer/aihw-cancer-statistics-framework
• Australian Bureau of Statistics (historical and projected populations): https://www.abs.gov.au/statistics/people/population

Specific source tables are referenced in the manuscript Methods and Supporting Information §S1.1–S1.2. The AIHW workbook table identifiers (S1a.1, S1e.1, S11a.1, S11e.1, S11f.1, S11i.1, Book 6) are the AIHW's own and are distinct from the SI's Table S1–S12 numbering.

Files required to recreate the analysis from raw data

code/import_data.R is the only script that reads data/raw/; it transforms the workbooks below into the small CSVs in data/ that the rest of the pipeline uses. data/raw/ is gitignored (large, third-party), so to reproduce the analysis end-to-end, download these files, place them exactly as shown, and run Rscript code/import_data.R before the apc_model.R → prev_model.R → supplement.R pipeline.

AIHW — Cancer Data in Australia 2025 (the "Data downloads" on the Cancer Data in Australia page above), placed in data/raw/AIHW/:

File	Sheet	Provides
CDiA-2025-Book-1a-Cancer-incidence-age-standardised-rates-5-year-age-groups.xlsx	Table S1a.1	Aggregate NHL and HL incidence, 1982–2021
CDiA-2025-Book-1e-Long-term-cancer-incidence-projections.xlsx	Table S1e.1	AIHW published incidence projections (comparison)
CDiA-2025-Book-6-Cancer-prevalence.xlsx	Table S6.2	Aggregate NHL/HL prevalence (validation)
CDiA-2025-Book-11a-Blood-cancer-incidence-by-histology-main-framework.xlsx	Table S11a.1	NHL-subtype incidence, 2003–2021
CDiA-2025-Book-11e-Blood-cancer-survival-and-age-group-proportions-by-histology-main-framework.xlsx	Table S11e.1	NHL-subtype and HL survival
CDiA-2025-Book-11f1-Part1-Blood-cancer-survival-by-histology-ICD10-framework-10-year-age-groups.xlsx	Table S11f.1	Aggregate NHL survival (part 1)
CDiA-2025-Book-11f1-Part2-Blood-cancer-survival-by-histology-ICD10-framework-10-year-age-groups.xlsx	Table S11f.1	Aggregate NHL survival (part 2)
CDiA-2025-Book-11i-Blood-cancer-prevalence-main-framework.xlsx	Table S11i.1	NHL-subtype prevalence (validation)

ABS, placed in data/raw/ABS/:

File (expected local name)	Sheet	Provides
3222_Table_B9.xlsx — Population Projections, Australia (cat. 3222.0), Table B9, Series B (medium growth)	Data1	Projected population, 2022–2050
3101059.xlsx — National, State and Territory Population (cat. 3101.0), Table 59	Data1	Historical population, 1982–2021

The committed CSVs in data/ already contain these processed inputs, so the pipeline runs without the raw workbooks; the list above is only needed to rebuild those CSVs from source.

What is and is not in this repository

Included: R code, processed input CSVs (a curated subset of AIHW tables in a structure suitable for direct use in the analysis), all modelled outputs, figures, and the source of the public Quarto dashboard.

The dashboard: dashboard/ is a static site built from the committed outputs by dashboard/export_dashboard_data.R and published to GitHub Pages by .github/workflows/publish-dashboard.yml on each push to main. It runs entirely in the reader's browser, requires no server, and contains only aggregate modelled results. See dashboard/README.md to preview or refresh it.

Not included: Original AIHW Excel workbooks (the included CSVs are a processed selection of the relevant rows and columns); the manuscript or Supporting Information PDFs (these are obtained from the journal). The R binary file apc_results.rds containing fitted model objects is regenerated by apc_model.R; if size becomes a problem for the public archive, it can be excluded and rebuilt from source.

Licence

• Code (code/ and dashboard/ directories, including the dashboard source and its export script): MIT licence — see LICENSE.
• Data and modelled outputs (data/ and output/ directories): Creative Commons Attribution 4.0 (CC-BY-4.0). Source AIHW data are released by AIHW under their own terms; consult the AIHW website for current licensing.

Contact

Adam Irving, Centre for Health Economics, Monash University — [email protected] (corresponding author).

For methodological detail, refer to the manuscript Methods / Supporting Information. The R scripts in code/ are extensively commented and cross-reference the relevant Supplementary Information sections.