Tutorial propensity support - #2
Merged
Merged
Conversation
The SCARF mouse-brain Perturb-seq pilot (58 perturbations, 68-227 cells each, 15,135 controls) exposed two inference defects in LFC: 83% of its discoveries were genes with zero counts in the perturbed arm, and real effects were estimated but not called because the default variance formula was not the variance of the AIPW estimator. Estimator - Influence-function variance var(eta)/n is the only variance estimator; the Welch-by-arm 'unequal' formula (2x the correct SE for equal arms, far more for rare treatments with class-balanced scores) is removed and kept as a warning alias of 'pooled' for one release. - Calibrated propensity scores (class_weight=None) replace 'balanced' as the default in LFC, cross_fitting, estimate_propensity_scores and refit_propensity_scores; ps_clip defaults to a prevalence-aware bound. - Model-based variance floor 1/(n1*tau1) + 1/(n0*tau0) on the log scale; arms with zero observed counts stay estimable at the floor so complete knockouts are reported; var_floored and std_raw columns. - Expression threshold thres_min='auto' (about min_counts=5 expected counts in the smaller arm) applied to observed as well as counterfactual arm means. - Small-sample correction for in-sample nuisance fits: n/(n-d) variance scaling and a t reference with n-d degrees of freedom. - Support columns n_treated, n_control, count_treated, count_control; warnings for floored small arms and for backend='fast' without crispyx. Validation (plan/validation_0.0.10): oracle simulations, fake-perturbation and label-permutation nulls, and old-vs-new comparisons on Perturb-seq, SEA-AD, Adamson, Replogle and SCARF; new tests/test_small_arm_inference.py; Welch tests re-derived. Tutorials - SCARF tutorial added with the full investigation and the 0.0.10 comparison (1,920 -> 12,346 discoveries, Wilcoxon Jaccard 0.028 -> 0.341). - Replogle: prep_tutorial_data.py now downloads the scPerturb raw-count release and subsets with crispyx; the previous subset held log1p-normalised values. Batch results regenerated on counts (r = 30; JIC re-check deferred). - Adamson refit at the r = 30 its JIC table selects, LFC in batches; SEA-AD, Perturb-seq (Python and R) re-rendered; LFC docs and changelog updated. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
- New block-structured batched IRLS for [covariates | one-hot treatments] (causarray/glm_onehot.py): Schur-complement Newton steps, exact vs statsmodels, active-set iteration; routed from fit_glm_auto and estimate_disp_auto. 5 s at d=7 to 132 s at d=231 on Replogle 3,000 x 8,563 versus 902 s (crispyx) and 1,090 s (statsmodels pool). - Fused parallel likelihood kernels nll_mat / grad_genes / grad_cells with fixed reduction order (bitwise identical across thread counts); one GCATE update drops from 2.34 s to 0.50 s on a 3,000 x 3,000 problem. - estimate_r fits the initial GLM once, sorts singular vectors by singular value, passes A_init per r (the previous start was dropped by estimate), and reports wall time per r. - Router constants unchanged; benchmark measurements recorded in docstrings and the changelog. Tests for the solver, the kernels and estimate_r. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…lback - _USE_ONEHOT_FOR_IMPUTE defaults to True: fit_glm_auto(..., A=A, impute=...) fits the joint [W | A] model once with the block-structured IRLS. On the Perturb-seq tutorial: 10.7 s vs 46.8 s (statsmodels pool), tau corr 0.9987, 7,460 of 7,468 / 7,482 discoveries shared. The crispyx per-perturbation fit produced 1,049 coefficients above _FAST_MAX_COEF on the same data and had always fallen back to statsmodels. - estimate_disp_auto returned None when neither the structured solver nor crispyx applied; it now falls back to estimate_disp. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
estimate_r on the 6,000-cell subsample with 200 treatment columns now finishes in 15 minutes and selects r = 10. The notebook keeps its r = 30 results and states that the refit is deferred; run_batch_0.0.10.py reads r from the table. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
crispyx 0.1.5 ships the block-structured `[covariates | one-hot treatments]` IRLS that causarray had been carrying itself, a BLAS Gram matrix in place of the three-operand einsum, internal design preconditioning and `min_mu` defaulting to 0. causarray now calls it instead of duplicating it. Removed: `causarray/glm_onehot.py`; `_fit_glm_fast_single`, `_fit_glm_fast_per_perturbation`, `_scale_design_columns` and the two deviance-residual helpers from `nb_glm_fast.py`; `_FAST_MAX_D`, the `n p / d_eff^2` throughput heuristic, `_USE_ONEHOT_SOLVER`, `_USE_ONEHOT_FOR_IMPUTE`, `_ONEHOT_MIN_GROUPS`, `_moments_dispersion` and both one-hot branches of `fit_glm_auto`. About 1,000 lines of package code; the two GLM modules go from 1,397 to 938 lines. Kept: statsmodels as `backend='original'`, as the fallback when crispyx is missing, and as the last resort when a fit diverges. Benchmarked on 0.1.5 rather than inherited: `_FAST_MIN_P` 50 -> 10 (the batched path is faster at every gene count measured, 533x at p=5 to 2.9x at p=500, agreeing with statsmodels to 1e-5 with the dispersion supplied); `_FAST_MAX_COEF` kept at 1e4 and documented as a non-finite guard, since the largest coefficient on any realistic design measured 17.6. Fixes: - `estimate_disp_auto` returns `None` again when no batched estimate is available, instead of a pooled moments estimate. The pooled estimate cost 0.15 of correlation with the truth on the deconfounding benchmark (0.6179 -> 0.4640; naive 0.6188). - `estimate_disp` and `fit_glm` no longer raise `IndexError` on `offset=False`. - `fit_glm` with more than one treatment column and `impute=False` returned all-zero coefficients: the fitted means were reshaped to `(-1, a)` instead of broadcast across the treatment axis, which raised inside the per-gene `except`. - `mem_limit_gb` is honoured on every imputation route again. - `fit_glm_ondisk` returned all-NaN coefficients when a cell had no counts among the genes read; such cells are dropped with a warning, and `offset=True` elsewhere now names them instead of propagating NaN. Tests: `test_glm_onehot.py` becomes `test_structured_glm.py` and tests causarray's routing and conventions rather than a solver causarray no longer owns. The on-disk tests read the in-repo Adamson subset (or `CAUSARRAY_TEST_H5AD`) instead of an absolute path. Two single-seed knife-edge assertions were rewritten around the claim that holds. 205 passed, 1 skipped. Version stays 0.0.10. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`cache_propensity_batch.py` read `replogle-r.csv`, the pre-0.0.10 log-normalised JIC table, while the notebook reads `replogle-r-0.0.10.csv`, so the focal batch was fitted at a different `r` from the analysis it is meant to diagnose. It also still passed `usevar="unequal"`, which 0.0.10 removed and accepts only as a deprecated alias. `replogle-r-0.0.10.csv` is regenerated on the current GLM path: 13.6 min against a run that had not finished in 14.75 h, selecting r = 10 as before. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
tune_penalty_factor selects the per-treatment L2 penalty for one propensity
covariate and returns the mapping refit_propensity_scores consumes, so the two
compose. It triggers on treatments failing a support check and returns the
smallest factor meeting a target.
Dropping the covariate is the infinite-penalty limit, so it bounds what any
finite factor can achieve. The search evaluates that endpoint first: when the
dropped fit already misses the target, the treatment is reported infeasible
after a single extra fit rather than an exhausted grid. on_infeasible then
decides between the largest factor ('best', the default) and leaving the arm
unpenalized, since the trigger has already said its support is inadequate.
Target a monotone metric. auc and overlap_ratio move monotonically with the
penalty; ess_treated_fraction does not, because a completely separated arm has
near-uniform weights and a deceptively high ESS that falls as the penalty
restores genuine overlap.
gcate_lfc_batch(save_nuisances=True) writes each batch's outcome-model
predictions beside the result cache. The outcome model does not depend on the
propensity design, so they can be fed back to LFC as Y_hat to re-estimate under
a different propensity without refitting the expensive stage.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Every tutorial now follows one layout: N_*.py scripts numbered in the order they run, the notebook unnumbered, data/ for the raw download and prepared inputs, results/ for everything generated, and a README mapping each file to the step that produces and consumes it. Version numbers are gone from filenames. Where stripping created a collision the current artifact takes the plain name and the superseded one gains a -legacy suffix, so replogle-r-0.0.10.csv becomes results/replogle-r.csv and the log-normalised-era table becomes results/replogle-r-legacy.csv. The replogle ignore rules collapse from twelve stale per-file entries to data/ and results/**, with negations for the two small JIC tables that are committed so the notebook runs without a refit. The results/** form matters: git cannot re-include a file whose parent directory is excluded by a directory pattern. 5_refit_propensity.py re-estimates LFC under a tuned propensity, reusing the cached outcome predictions rather than refitting them. sea_ad_lfc.csv is removed: nothing reads it, and it predated the current fit. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The tutorials estimated log-fold changes and then diagnosed the propensity model. That order lets an arm whose inverse-probability weights concentrate on a few cells report a large discovery count before the reader has any reason to distrust it, so diagnostics and the refit now come first. Each notebook runs tune_penalty_factor itself rather than reading a table a script produced, and shows the overlap histograms before and after the tuned penalty on the same arms. The estimate then passes those scores to LFC, so the diagnostics and the results describe the same propensity instead of diagnosing one model and estimating with another. The target has to match the failure mode and it differs by dataset. Replogle and perturb-seq lose histogram overlap while keeping their treated sample, so they target overlap. Adamson's SEC61A1 keeps an overlap of 0.45 while its treated ESS collapses to 0.001, so it targets both; an overlap-only target reported it as already met and left the weights concentrated. Two claims are withdrawn. Library size was described as a consequence of the treatment and as a post-treatment variable; measured library size mixes capture depth with total RNA content and these data do not identify which drives a given arm's shift, so the text now states only what is observed. "Sensitivity" framing is dropped throughout, and the per-treatment covariate filtering is presented as a different operation from penalisation rather than a rival option: a penalty shrinks a coefficient the data cannot support, dropping asserts the covariate does not belong in the model. SEA-AD needed neither reordering nor tuning: its arms are balanced, nothing trips the support check, and its diagnostics already preceded the estimate. Its stale GCATE cache and June-era JIC table are rebuilt against the current engine. Hard-coded discovery counts that no longer matched the fit are removed in favour of the numbers the cells print. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
- gcate_lfc_batch: reject save_nuisances without cache_path up front, write nuisances before the result cache, refit cached batches whose nuisances are missing on resume - Fall back to per-gene NB dispersion instead of a fixed size of 1 - Reject unknown GLM families in fit_gcate and the likelihood kernels - Accept per-treatment ps_clip arrays - tune_penalty_factor: honour mask, score the reported factor, add target_met, fix the docstring trigger key - estimate_r honours backend - Trim docstrings and the 0.0.10 changelog; fold 0.0.11 into 0.0.10 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
![claude[bot]](https://avatars.githubusercontent.com/in/1236702?s=60&v=4){kind=small}
There was a problem hiding this comment.
If your organization's extra usage balance is empty, an organization admin can add extra usage credits at claude.ai/admin-settings/usage. If its monthly spend limit was reached, an admin can raise it on the same page. If neither applies, contact Anthropic support.
Once extra usage is available, comment @claude review on this pull request, or push a new commit, to trigger a review.
The package builds from pyproject.toml, which did not list crispyx (only setup.cfg did). Add it, plus a test extra used by the workflow. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Python 3.8/3.9 are dropped, crispyx becomes required, LFC defaults change and glm_onehot is removed, so this is a minor, not patch, bump. Recorded notebook outputs keep the version they were run with. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The nuisance store is written batch by batch, so an unexpected size was discovered only when the disk filled. Estimate it once before the first write: verbose=True reports the total, and a store above NUISANCE_WARN_BYTES raises a ResourceWarning naming it and the free space on that volume. The estimate is 8 * n_cells * n_genes * n_treatments per batch, which grows with the number of treatments as well as the data. That factor is easy to miss: the SCARF tutorial would write about 112 GB across six batches of ten perturbations, against 58 GB for Replogle's fourteen batches on a third of the cells. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Both deconfounding tests used one seed on which GCATE's stage 1 stops on a plateau before recovering the factors (also at v0.0.9), so they passed or failed on numerical noise. They now assert median MSE, a majority of seeds and a per-seed bound over five seeds. gcate_lfc_batch's cache uses pd.HDFStore, which needs pytables. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.