Repository navigation
Releases: JinHongDu-Lab/causarray
Release list
0.1.1
[0.1.1]
Fixed
LFCno longer excludes genes whose arm means differ by less than
thres_diff, which is now only the floor before the logarithm. The filter
selected on the effect estimate: it dropped the pairs with p-values near 1
(about 10% of expressed pairs on a Perturb-seq negative control), so BH over
the remaining pairs was anti-conservative.- Size factors (
comp_size_factor, used byfit_gcate,LFC(offset=True)
andgcate_lfc_batch) take the median of ratios over genes with mean count
of at leastmin_mean=2instead of every gene. In sparse single-cell data
the old factors captured about a third of the depth variation; the rest
shifted all highly expressed genes of an arm together, so their null
statistics were too spread. On a Replogle negative control (30 fake
perturbations cut from the controls) this cuts false discoveries from 114
to 14, and the null SD for genes above 5 counts per cell from 1.32 to
1.08.min_mean=0restores the old factors;method='libsize'uses total
counts. A cell with no counts in those genes gets its total count, rescaled
to the other cells' factors, instead of a factor near the median. comp_size_factor(method='scale')returns one factor per row (total count
overlib_size); it returned one value per column.fit_gcatenow uses anoffsetarray. The fit reads its offset from
size_factor, which was set only foroffset=True, so a supplied array was
replaced by ones everywhere except the dispersion estimate.LFCshows its per-treatment progress bar only withverbose=True.
Added
select_propensity_factorschooses the propensity covariates per
treatment: for a treatment whose weights concentrate (treated ESS < 0.5,
overlap < 0.3 or AUC > 0.9) it drops the covariates most imbalanced between
that treatment and the controls until its support recovers, and returns a
drop_by_treatmentmapping forrefit_propensity_scores. A treatment that
already meets the target keeps every covariate; when the target cannot be
reached,on_infeasiblekeeps the closest fit ('best', default), drops
every candidate ('all') or none ('none').import causarraywarns when it runs in an Intel (x86_64) Python translated
by Rosetta 2 on an Apple Silicon Mac, where its numerical code runs several
times slower (LFCon a 5,000-cell subsample: 159 s against 32-57 s
natively). The README shows how to create a native environment.
Changed
- crispyx >= 0.1.6 is required. With latent factors in the outcome model,
crispyx 0.1.5's structured GLM fitter was about 2.5 times slower than
0.1.6's (same estimates): 159 s against 62 s for LFC on 5,000 cells, 3,000
genes and 10 treatments. LFC(usevar='unequal')is the by-arm Welch variance again (0.1.0 made it an
alias of'pooled'). It is conservative for small case-control designs: on
85 SEA-AD donors, permuted disease labels gave 0-1 false discoveries with
null SD 0.6-0.9, against SD 1.1-1.8 under'pooled'.'pooled'stays the
default;'unequal'over-inflates the standard errors of rare treatments.
v0.1.0
[0.1.0]
Inference fix for small perturbation arms, and the GLM engine moves to crispyx
(>= 0.1.5 now required).
Changed
- crispyx >= 0.1.5 is required.
causarray/glm_onehot.py, the two-stage
per-perturbation fitter and the design preconditioner are removed in favour
of crispyx'sStructuredGLMBatchFitterand dense batch fitter, which agree
with the statsmodels path to 2e-5 on well-conditioned genes. fit_glm_autoroutes on gene count only (_FAST_MIN_P = 10), with the
_FAST_MAX_COEFdivergence check falling back to statsmodels.- Fitted means are no longer floored at
min_mu = 0.5, which biased genes
below about one count per cell. LFCuses the influence-function variance (usevar='pooled') only; the
Welch'unequal'variance overstated standard errors.- Propensity scores default to calibrated probabilities
(class_weight=None);'balanced'remains available. ps_clipdefaults to'auto', a prevalence-aware bound per treatment;
resolved bounds are returned asestimation['ps_clip_bounds'].LFCdefaultthres_min='auto'tests a gene only when aboutmin_counts
(default 5) counts are expected in the smaller arm.
Added
- A Poisson variance floor
1/(n1*tau1) + 1/(n0*tau0)inLFC, so chance
all-zero arms of sparse genes are not called while complete knockouts stay
significant (var_flooredcolumn). - Small-sample correction for
K=1: variance scaled byn/(n-d)and a t
reference withn-ddegrees of freedom. - Per-pair support columns
n_treated,n_control,count_treated,
count_controlandstd_rawin everyLFCresult. RuntimeWarnings for arms under 200 cells hitting the variance floor, and
forbackend='fast'without crispyx.tune_penalty_factorpicks the smallest per-treatment L2 penalty on one
propensity covariate that meets a support target, checking the
dropped-covariate limit first to detect infeasible treatments.gcate_lfc_batch(save_nuisances=True)saves each batch's outcome-model
predictions to<cache_path stem>.nuisances.h5for re-estimation under a
different propensity model; budget about
8 * n_cells * n_genes * n_treatmentsbytes per batch.estimate_racceptsbackendand returns atime_scolumn.
Performance
- Designs with one-hot treatment blocks are solved by a Schur complement,
flat in the number of treatments (n = p = 3,000: 5.1 s with 10 treatments,
7.7 s with 200). - GCATE's objective and gradients use parallel numba kernels with a fixed
reduction order (results independent of thread count); one update on a
3,000 x 3,000 problem drops from ~2.3 s to 0.5 s. estimate_rfits the initial[X | A]GLM once for all candidater.
Fixed
- NB dispersion falls back to the per-gene
estimate_dispwhen no batched
estimate is available, instead of a fixed size of 1. estimate_dispandfit_glmacceptoffset=False.fit_glmwith several treatments andimpute=Falsereturned zero
coefficients.mem_limit_gbis honoured on every imputation route.fit_glm_ondiskdrops cells with no counts (with a warning) instead of
returningNaNcoefficients.
Deprecated
LFC(usevar='unequal')is an alias of'pooled'and warns.LFC(eps_var=...)is ignored; the variance floor supersedes it.
v0.0.9
[0.0.9] - 2026-07-23
Added
align_test_mask()aligns named treatment-by-gene diagnostic masks to
existing causarray result tables without refitting effects or changing
inference.summarize_treatment_associations()and
plot_treatment_associations()diagnose treatment associations with both
observed covariates and estimated latent factors.refit_propensity_scores()supports treatment-specific covariate removal and
feature-specific logistic L2 penalties while preserving untreated score
columns and reusable outcome predictions.- The Replogle tutorial now compares Wilcoxon and causarray effect estimates in
the extreme negative LFC tail, reports expression support, and demonstrates
diagnostic support rules.
Changed
- Perturb-seq propensity diagnostics now include observed library size and a
treatment-specific library-size regularization sensitivity analysis. - Public diagnostics are documented in the README and LFC API guide, including
the distinction between post-hoc result annotation and formal refitting or
multiple-testing changes. summarize_treatment_associations()acceptsbh_scope='per_treatment'to
adjust within each treatment instead of across the whole treatment-by-covariate
grid, and reportsn_tests_in_familyso the correction is self-describing.plot_treatment_associations()accepts an explicit symmetricvmax.plot_propensity_scores()accepts and forwardsclip_bounds, which was
previously fixed at its default in the summary it returns.refit_propensity_scores()reportsdegenerate_designandscore_stdper
refitted treatment.
Fixed
refit_propensity_scores()documented that carried-over score columns are
preserved exactly, whileclipwas in fact applied to the whole returned
matrix. The uniform behaviour is kept, sinceLFCconsumes a single bound,
but it is now stated and covered by a test. Passclip=Noneto leave
carried-over scores untouched.refit_propensity_scores()now warns when covariate filtering leaves a
constant design, which previously collapsed to a covariate-free0.5score
in silence and quietly turned AIPW into an unweighted contrast.- A scalar
clipraisedTypeError: object of type 'float' has no len()
instead of the documentedValueErrorinestimate_propensity_scores()and
refit_propensity_scores(). plot_treatment_associations()scaled the colour bar to the largest observed
magnitude even forspearman_rho, which is bounded in [-1, 1]. Tiny
correlations rendered as fully saturated and panels were not comparable
across subsets; bounded statistics now use a fixed(-1, 1)range.summarize_propensity_scores()reportedclipped_fraction = 0.0for raw
scores by inferring clipping from a defaultclip_boundsthe caller may
never have applied.clip_bounds=Nonenow yieldsNaN.- The Replogle propensity cache script derived log-library size from raw counts
whileprep_causarray_dataderives it from capped counts, so its design
matrix did not match the tutorial's. The script now caps first and asserts
the preprocessing contract before fitting.
v0.0.8
[0.0.8] - 2026-07-20
Added
LFC()andgcate_lfc_batch()results now includelog2fcand
log2fc_se, the base-2 equivalents of the existing natural-logtauand
stdcolumns. Existing statistics, p-values, and discoveries are unchanged.
Changed
- Older compatible
gcate_lfc_batch()caches are upgraded in memory with the
new log2-scale columns when resumed. - The Python environment now requires NumPy >2. A separate
causarray-r
environment supplies R 4.4 and a NumPy-2-compatiblereticulate, avoiding
changes to the Python numerical stack for R interoperability.
v0.0.7
[0.0.7] - 2026-07-18
Added
estimate_propensity_scores()estimates in-sample or out-of-fold treatment
probabilities without rerunning the outcome model.summarize_propensity_scores()andplot_propensity_scores()report and
visualize common support, tail mass, and inverse-weight effective sample size.LFC()now returns both clippedpi_hatand diagnosticpi_hat_rawscores.LFC()results include raw AIPWmean_control,mean_treated, and
estimablediagnostics.- Propensity and latent-factor association diagnostics in the Perturb-seq
tutorial.
Changed
estimate_propensity_scores()now defaults to class-balanced logistic
weighting, matchingLFC()and historical causarray behavior. Pass
class_weight=Noneexplicitly for calibrated treatment probabilities.
Fixed
LFC()now always uses standard, unclipped AIPW pseudo-outcomes while
retaining its established class-balanced propensity nuisance fit to limit
unrelated result drift. This removes directional bias from clipping negative
cell-level influence-function values under treatment-group imbalance.
Calibrated scores remain available withps_class_weight=None; the legacy
cell-level clipping option has been removed.- Valid aggregate counterfactual means are floored only for the final log-ratio
and delta-method calculation, after averaging the unchanged pseudo-outcomes. - Nonfinite or nonpositive aggregate counterfactual means are reported as
non-estimable rather than silently floored before taking a log ratio. cross_est=Truenow enables two-fold nuisance cross-fitting, with an explicit
Ktaking precedence.- Intercept-only propensity models now respect the requested class weighting:
balanced fits return 0.5 and unweighted fits return the case prevalence.
Masks are indexed correctly within cross-fitting folds.
v0.0.6
[0.0.6] - 2026-05-31
Added
-
Batch-wise fitting API (
causarray/gcate.py,causarray/DR_learner.py,causarray/utils.py)fit_gcate_batch(Y, X, A, r, batch_size=10, max_cells=2000, n_ctrl=2000, ...):
Fits GCATE independently on batches of perturbations. A shared control
subsample ofn_ctrlcells is reused across batches; dispersion is
pre-estimated once on the control pool. Supportsskip_batchesto resume
interrupted runs; reports per-batch wall time and ETA whenverbose=True.gcate_lfc_batch(Y, X, A, r, batch_size=10, max_cells=2000, n_ctrl=2000, cache_path=None, ...):
End-to-end batch pipeline — runs GCATE and LFC per batch, freeing large
intermediate arrays after each batch.cache_pathenables HDF5 disk
caching viapandas.HDFStoreso interrupted runs resume from the last
completed batch. Returns a concatenated DataFrame with a'batch'column.LFC_batch(...): deprecated alias forgcate_lfc_batch; emits
DeprecationWarningand will be removed in a future release.n_batchesparameter forfit_gcate_batchandgcate_lfc_batch:
specifies total number of batches instead of per-batch count; overrides
batch_sizewhen set.estimate_r(max_cells=N, random_state=0): new parameter that
automatically subsamples to at mostNcells before running JIC
selection, prioritising control cells.
-
Fast GLM backend via crispyx (
nb_glm_fast.py,gcate_glm.py)fit_glm_fast(): Batch NB-GLM fitting using crispyx'sNBGLMBatchFitter,
replacing per-gene statsmodels IRLS with vectorized batch IRLS.estimate_disp_fast(): Vectorized method-of-moments dispersion estimation.fit_glm_ondisk(): On-disk streaming GLM fitting for large h5ad files.- Per-perturbation fitting (
_fit_glm_fast_per_perturbation): for
multi-treatment data, fits binary (ctrl vs. treatment_k) models
independently, then assembles the full coefficient matrix. fit_glm_auto(): Routes tofit_glm_fast()when crispyx is available and
the effective design dimension is small; falls back to statsmodels otherwise.estimate_disp_auto(): Routes toestimate_disp_fast()for large gene
counts; falls back to statsmodels otherwise.
Fixed
- Numba TBB fork warning: Set
NUMBA_THREADING_LAYER_PRIORITYto prefer
OpenMP over TBB in__init__.py, eliminating fork warnings when Joblib forks
after Numba parallel execution. Addedllvm-openmpto conda dependencies. - Fast-path threshold (
gcate_glm.py): Raised the effective design-dimension
ceiling so largerrvalues and wide batch designs correctly use the crispyx path. - Backend toggle (
gcate_glm.py): Re-added_USE_FAST_BACKENDmodule flag
and_backend_override()context manager for reliable statsmodels fallback. - Weighted dispersion (
nb_glm_fast.py): Dispersion averaging is now
cell-count-weighted; low-coverage perturbations contribute proportionally less. - Control-cell residuals (
nb_glm_fast.py): Fixed last-perturbation overwrite
bug; control-cell deviance residuals and fitted values are now initialised from
the global covariate model. - Module-qualified imports (
gcate_opt.py,gcate.py,DR_estimation.py):
Backend toggles now propagate correctly at call time. estimate_rbare name (gcate.py): FixedNameErrorcaused by a bare
fit_glm_autoreference.- crispyx availability check (
gcate_glm.py): Users without crispyx now get
a transparent fallback to statsmodels instead of a traceback.
Changed
-
⚠️ alter_min()early-stopping defaults (gcate_opt.py):- Default
kwargs_es['max_iters']reduced from 500 → 50. - Default
tolerancereduced from1e-3→0.0; new scale-invariant
rel_tol=2e-4introduced. To reproduce pre-v0.0.6 behavior, pass
kwargs_es_1=dict(max_iters=500)andkwargs_es_2=dict(max_iters=500).
- Default
-
⚠️ BREAKING —LFC()variance and defaultusevar(DR_learner.py):- Default
usevarchanged from'pooled'to'unequal'(Welch). Revert
withLFC(..., usevar='pooled')if reproducing pre-v0.0.6 results. 'unequal'formula corrected: variance is nows₀²/n₀ + s₁²/n₁
(standard Welch); the prior version used(s₀²/n₀ + s₁²/n₁)/2
("half-Welch"), under-estimating the standard error by √2.- p-values now use the t-distribution with Welch-Satterthwaite degrees of
freedom per gene; the prior version used a Normal approximation.
- Default
-
alter_min()initialisation,_check_input(),estimate_r(), and
cross_fitting()now use the auto-dispatch GLM/dispersion paths. -
LFC()now acceptsbackend: str = "auto"("fast"forces crispyx,
"original"forces statsmodels). -
comp_size_factor()vectorized withnp.nanmean/np.nanmedian.
Performance
Benchmarked on Perturb-seq data (n = 2,926 cells, p = 3,221 genes, 29 perturbations):
| Component | Original | Fast | Speedup |
|---|---|---|---|
| GCATE | 331.6 s | 298.5 s | 1.1× |
| LFC | 87.8 s | 65.7 s | 1.3× |
| Total | 419.3 s | 364.2 s | 1.2× |
On synthetic data (n = 500, p = 200): 61.5× GLM fit speedup, 7.1× imputation speedup.
Latent factor recovery: mean canonical correlation 0.998. LFC correlation: 0.856.
Additional LFC throughput improvements on the Replogle tutorial dataset
(79,865 cells × 8,563 genes, 200 perturbations, 14 batches):
| Change | Speedup contribution | Accuracy impact |
|---|---|---|
Stage 1 max_iter 50 → 5 (NB) / 10 (Poisson) |
−10 min | identical (r=1.000) |
| Stage 1 ≤3,000-cell mixed subsample | −55 min | tau r=0.992, Jaccard=0.80 |
| Stage 2 joint fit | −5 min | tau r=0.9994, Jaccard=0.975 |
| Combined | −70.6 min / 1.48× | tau r=0.9994, Jaccard=0.975 |
Full-run: 217.5 min → 146.9 min (1.48× faster); sig pairs −0.2%, perts with ≥1 hit −0.6%.
v0.0.5
- Add more tuning options for random-forest classifiers
- Improve the uninformative tests thresholding
Full Changelog: v0.0.4...v0.0.5