Evaluation
==========

Abstract
--------

**pyglmGamPoi** is a pure Python/C++ reimplementation of
``sctransform::fit_glmGamPoi_offset``, the negative-binomial offset model at the
core of Seurat SCTransform v2. We evaluate numerical parity against R on both
synthetic count matrices (100–2000 genes × 200–2000 cells) and three real
10x Genomics samples from GSE288946 (~2000 × 2000 SCT step-1 subsamples each).

Across all benchmarks, finite-gene ``theta`` estimates agree with R at
machine-precision levels (median |Δ| ≈ 2×10\ :sup:`-7`), intercepts match at
**~10\ :sup:`-9`** (100% within 0.01 with shrinkage enabled), and
Poisson-boundary ``theta = Inf`` assignments match R exactly (**463/463**
pooled). With OpenMP (8 threads) and a native CSR sparse path, Python completes
the real-data step-1 workload in **7.1 s** vs R **46.5 s** (~**6.4×** faster),
while eliminating the R runtime dependency for
`trackcell <https://github.com/seqyuan/trackcell>`_ SCT v2 pipelines.

Motivation
----------

Single-cell RNA-seq normalization via **SCTransform v2** fits, for every gene, a
negative-binomial generalized linear model with a sequencing-depth offset:

.. math::

   y_{gc} \sim \mathrm{NB}(\mu_{gc}, \theta_g), \quad
   \log \mu_{gc} = \beta_g + \log(10^{\mathrm{log10\_umi}_c})

Seurat implements this through the R package **glmGamPoi**
(`Hafemeister & Satija, 2019 <https://github.com/const-ae/glmGamPoi>`_).
Python pipelines such as **trackcell** historically bridged to R via
``Rscript`` subprocesses, adding environment complexity (conda ``st``, CRAN/Bioconductor
versions, ``TRACKCELL_RSCRIPT`` overrides) and blocking fully native workflows.

**pyglmGamPoi** reimplements the glmGamPoi algorithms in C++17 (Newton–Raphson
intercept fit, Cox–Reid overdispersion MLE, optional cross-gene
``overdispersion_shrinkage``) and exposes a drop-in API for trackcell:

.. code-block:: python

   from pyglmGamPoi import fit_offset_model, is_available

   if is_available():
       model_pars = fit_offset_model(umi, regressor_data, gene_index)

This page documents the **replacement effect**: what R call is replaced, how we
measure parity, and what accuracy/speed users should expect.

What is replaced (and what is not)
-----------------------------------

.. list-table::
   :header-rows: 1
   :widths: 35 15 50

   * - Component
     - pyglmGamPoi
     - Notes
   * - ``fit_glmGamPoi_offset``
     - **Replaced**
     - Per-gene NB offset fit; primary scope
   * - Cox–Reid overdispersion MLE
     - **Replaced**
     - ``allow_inf_theta`` → ``theta = Inf`` at Poisson boundary
   * - Cross-gene ``overdispersion_shrinkage``
     - **Replaced**
     - ``loc_median_fit`` trend + intercept refit (R ``glm_gp`` default)
   * - ``vst()`` subsampling (cells/genes step1)
     - trackcell
     - RNG differs between R ``sample()`` and NumPy
   * - ``_reg_model_pars`` ksmooth / outlier logic
     - trackcell
     - Kernel smoothing after initial fit
   * - Pearson residuals, HVG, integration
     - trackcell
     - Downstream of ``model_pars``

Replacing only the NB fit isolates a well-defined numerical target and avoids
conflating RNG subsampling differences with algorithmic ones.

Benchmark design
----------------

Two complementary suites:

**Synthetic scaling** (``scripts/generate_benchmark_report.py``)
  Poisson-distributed UMI matrices at five sizes, same API as production.
  Validates that parity holds from smoke tests (100 genes) up to SCT step-1 scale
  (2000 × 2000).

**Real data — Tier B** (``scripts/run_tcl_benchmark.py``)
  Three GSE288946 samples (GSM8779707–709), filtered 10x counts aligned to R
  ``cells_step1`` / ``genes_step1`` exports from Seurat ``vst(seed = 1448145)``.
  This **Tier-B** protocol (see trackcell ``tcl_test`` README) removes
  subsampling RNG as a confounder so differences reflect the fitter only.
  Count matrices are passed as **CSR sparse** (typical 10x input) without
  densifying.

Reference implementation
  Fresh ``sctransform::fit_glmGamPoi_offset(umi, model_str = "y ~ log_umi", ...)``
  in conda env ``st`` (Seurat 4.4 + glmGamPoi). Cached ``model_pars_step1.csv``
  from ``vst()`` is **not** used as ground truth — R itself differs from
  ``fit_glmGamPoi_offset`` on the intercept scale (~0.5).

Environment
  Benchmarks run with ``OMP_NUM_THREADS=8`` unless noted. Regenerate with
  ``python scripts/run_all_benchmarks.py``.

Metrics
  Per gene: ``theta``, ``Intercept``, ``log_umi`` (= ``log(10)``).
  Reported: median absolute difference, Pearson *r*, fraction within 0.01,
  Inf-``theta`` match count, wall time.

Results
-------

Synthetic scaling
~~~~~~~~~~~~~~~~~

.. include:: _generated/benchmark_table.rst
   :start-line: 0

At SCT step-1 scale (2000×2000 synthetic), Python is **~8.8×** faster than R with
identical ``theta`` / ``Intercept`` (med |Δ| ≈ 10\ :sup:`-11` / 10\ :sup:`-7`).

Real data (GSE288946)
~~~~~~~~~~~~~~~~~~~~~

.. include:: _generated/tcl_benchmark_table.rst
   :start-line: 0

Performance summary
~~~~~~~~~~~~~~~~~~~

.. list-table:: Speed vs R (``fit_glmGamPoi_offset``, 8 OpenMP threads)
   :header-rows: 1
   :widths: 28 18 18 18

   * - Workload
     - Python
     - R (conda ``st``)
     - Speedup
   * - Synthetic 2000×2000
     - 2.0 s
     - 17.5 s
     - 8.8×
   * - Real GSM8779707 (1995×2000)
     - 2.5 s
     - 15.5 s
     - 6.2×
   * - Real 3-sample pooled
     - 7.1 s
     - 46.5 s
     - 6.4×

Implementation milestones (same Tier-B data, GSM8779707):

.. list-table::
   :header-rows: 1
   :widths: 40 20 20

   * - Configuration
     - Wall time
     - vs R
   * - Initial single-thread C++ (dense)
     - ~15 s
     - ~1.0×
   * - + OpenMP (8 threads)
     - ~2.2 s
     - ~7×
   * - + CSR sparse path (no ``toarray()``)
     - ~2.2 s
     - ~7× (parity unchanged)

The CSR path avoids materializing ~90% zeros in 10x matrices; numerical results
match dense after ``sort_indices()`` (required when gene/cell slicing leaves
unsorted CSR rows).

Representative gene-level agreement (GSM8779707, vs R live)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. list-table::
   :header-rows: 1
   :widths: 20 15 15 15

   * - Gene
     - Python θ
     - R θ
     - |Δ Intercept|
   * - RAB7A
     - 4.859
     - 4.859
     - ~10\ :sup:`-9`
   * - RPS7
     - 7.559
     - 7.559
     - ~10\ :sup:`-9`
   * - TTK (multimodal OD)
     - 0.00774
     - 0.00774
     - ~10\ :sup:`-8`
   * - AL392046.1 (Poisson)
     - Inf
     - Inf
     - exact

Downstream effect (trackcell integration)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

When Python uses R-exported ``model_pars_fit`` (ksmooth applied in R) and
computes Pearson residuals, residual correlation vs R reaches
*r* ≈ 1.0 (median |diff| ≈ 10\ :sup:`-13` on GSM8779707). This confirms
that **once model parameters match, downstream SCT quantities are unchanged**
within floating-point noise.

Discussion
----------

**Accuracy.** Finite ``theta`` values are effectively identical to R (median |Δ|
~ 2×10\ :sup:`-7` on real data). With default cross-gene shrinkage, intercepts
match R at machine precision (median |Δ| ~ 2.5×10\ :sup:`-9`; **100%** of genes
within 0.01). Shrinkage ``gene_means`` now follow R ``glm_gp`` first-pass
``rowMeans(Mu)`` (rough-dispersion beta), not the post-OD intercept. ``theta =
Inf`` genes (Poisson boundary under ``allow_inf_theta = TRUE``) match R exactly
(217/217, 188/188, 58/58 across the three samples). Pearson *r* for ``theta``
on real data is not informative when many genes share Inf status; use median
|diff| instead.

**Speed.** OpenMP parallelizes the per-gene C++ loop (``schedule(dynamic, 16)``).
The native CSR path reads 10x count matrices in place (implicit zeros still
contribute to the likelihood). On GSE288946 step1 (2000×2000, 8 threads), Python
is **~6–7× faster** than R while preserving identical ``theta`` / ``Intercept``
parity. Set ``OMP_NUM_THREADS`` before importing ``pyglmGamPoi``.

**Why not compare to ``vst()`` exports directly?** Seurat's ``vst()`` pipeline
applies additional steps before/after ``fit_glmGamPoi_offset``; cached
``model_pars_step1.csv`` shows intercept offset ~0.5 vs a direct
``fit_glmGamPoi_offset`` call even in R. Our comparison target is the function
trackcell actually replaces.

**Integration path.** In trackcell ``fit_offset_model``, priority becomes:
(1) pyglmGamPoi native, (2) R ``fit_glmGamPoi_offset``, (3) statsmodels
``nb_offset`` fallback. Users install with ``pip install pyglmGamPoi`` — no
``conda env st`` required for the NB fit step.

Conclusion
----------

pyglmGamPoi provides a **drop-in, R-free replacement** for
``sctransform::fit_glmGamPoi_offset`` with demonstrated numerical parity on
synthetic and real scRNA-seq data, and **~6–7×** wall-clock speedup on production-scale
SCT step-1 workloads (8 threads). For trackcell SCT v2 users, this removes the
largest R dependency in the normalization path while preserving downstream
compatibility with Seurat-derived workflows.

Reproducing
-----------

.. code-block:: bash

   pip install -e ".[test]"
   export OMP_NUM_THREADS=8
   python scripts/run_all_benchmarks.py
   pytest tests/test_theta_od_parity.py -q
   pytest tests/test_parity_r.py -m parity_r   # optional, needs R

Environment variable ``PYGLMGAMPOI_RSCRIPT`` overrides the Rscript path for
parity tests (default: ``conda run -n st Rscript``).