Spike-counting benchmark

Benchmarks methods for estimating the number of significant principal components ("spikes") in high-dimensional data, using a synthetic signal-plus-noise model where the true number of spikes is known.

Methods compared

Method	What it is	Where
scree	Largest-gap elbow rule (naive baseline)	R/baselines.R
parallel	Parallel analysis, noise-edge variant	R/baselines.R
GetRank	RMT4DS sequential rank test (gap-ratio)	R/rmt_methods.R
ScreeNot	RMT4DS denoising threshold (Donoho et al.)	R/rmt_methods.R

RMT4DS is Prof. Xiucai Ding's package (github.com/xcding1212/RMT4DS).

Reproduce everything

R (>= 4.6 used here). From the project root, in R:

renv::restore()                       # rebuild the exact package environment
testthat::test_dir("tests/testthat")  # run the tests
source("scripts/run_benchmark.R")     # produce the table + plot in results/

Outputs land in results/:

• results_table.csv — mean estimate and mean abs error per method per true-k
• results_raw.csv — every (scenario, seed, method) run
• comparison_plot.png — estimated vs. true, with the y=x reference line

Data model

simulate_spiked() (R/data_sim.R) builds Y = S + noise, where the noise entries are N(0, 1/n) and S is a rank-k signal. The k spike strengths are the singular values of S; the spike count k is the ground truth. This is the signal-plus-noise model RMT4DS's GetRank/ScreeNot expect (they SVD Y directly), so all four methods work off the singular values of Y.

Findings

• ScreeNot and the noise-edge parallel analysis recover the true count accurately across all tested k, including weak spikes near the detection edge.
• GetRank under-counts once spikes are graded/closely spaced (k >= 4 here), because its gap-ratio statistic has no clear gap to latch onto. So ScreeNot is the recommended RMT4DS estimator for this regime; GetRank's behaviour is kept in the benchmark as a documented result.
• scree (largest-gap) fails from k=3 on — the top gap dominates.

Notes

• RMT4DS lives in a subdirectory of its repo, so it installs with renv::install("xcding1212/RMT4DS/RMT4DS") (note the trailing /RMT4DS). It does not declare its dependencies, so MASS, pracma, rARPACK, nleqslv, lpSolve and rootSolve are installed alongside it.
• SpikeDetection (github.com/CharbelAbiYounes/SpikeDetection) was not integrated: it is a Julia codebase, not an R package. Bridging to Julia would break the renv-based reproducibility of this project. A future option is to port its BEMA method to R as a fifth comparison.