spacr.sp_stats¶

Provide statistical tests and multiple-comparison helpers.

Group-comparison functions delegate test selection to spacr.figures.stats while preserving this module’s established call signatures and result keys. Two groups use Student’s t, Welch’s t, or Mann–Whitney U as supported by the data; larger designs use one-way ANOVA, Welch’s ANOVA, or Kruskal–Wallis. Results identify the selected test and include the assumption checks used to select it.

perform_normality_tests() reports underpowered checks as uninformative, and perform_levene_test() uses the median-centred Brown–Forsythe statistic. Imports of the plotting-backed statistical engine remain local so callers that only need adjustment or contingency-table helpers avoid loading the plotting stack.

Arrayed-screen hit statistics¶

score_arrayed_screen() scores every well of an arrayed screen against its negative control: robust z (median/MAD), SSMD with the method-of-moments, UMVUE and robust estimators for designs with and without replicates (Zhang 2011, J Biomol Screen 16:775), and the B-score (Brideau et al. 2003), whose row and column effects come from median_polish(), a port of R’s stats::medpolish. Scores are computed per plate or against the pooled negative control, hits are called at chosen thresholds, and write_hit_report() writes the ranked hit table as CSV and one plate heatmap per statistic through spacr.plot.save_figure(). Wells are located by spacr.plate_qc and controls named in the spacr.well_spec notation, and the per-plate Z’ comes from the control-chart screen’s own zprime_frame().

The constants: MAD_SCALE is 1 / Phi^-1(0.75) = 1.4826, which makes the MAD of a normal sample estimate its standard deviation (the factor R’s mad() and so cellHTS2 apply). SSMD_ESTIMATORS are mm (method of moments), umvue (uniformly minimal variance unbiased) and robust (median/MAD, SSMD*). DEFAULT_HIT_THRESHOLDS are 3 for every statistic: SSMD 3 is Zhang’s “strong” effect, and 3 robust sigma is the usual cut-off for robust z and the B-score. MEDIAN_POLISH_MAX_ITER and MEDIAN_POLISH_EPS are R’s medpolish defaults (10 and 0.01), kept so the residuals agree with the B-score cellHTS2 computes.

Image-based profiling¶

The profiling helpers turn Measure’s per-object tables into well and treatment profiles and score them, following pycytominer and copairs so the numbers agree with those tools on the same features. Each object table is aggregated to a per-well median, a plate map adds the treatment annotations, features are normalised plate by plate against the negative-control wells (robust MAD by default), uninformative and redundant features are removed (variance, frequency, correlation, missing values and outliers), and replicate wells are collapsed into consensus profiles. Replicate reproducibility is scored as mean average precision (mAP) with copairs’ permutation null and Benjamini-Hochberg correction: phenotypic activity (replicates against controls) and, given a phenotype label, phenotypic consistency (treatments that share the label). Percent replicating is reported beside it. Profiles are written as CSV and Parquet with Metadata_ columns, the consensus also as GCT, with an mAP plot and a consensus-similarity heatmap. Measure runs the whole recipe at the end of a run when profiling is on.

Exceptions¶

HitScoringError

Raised when a screen cannot be scored, with the reason and the way out.

Classes¶

`ArrayedHitResult`	Per-well scores, per-treatment SSMD and per-plate QC of one screen.
`MedianPolish`	Tukey's median polish of one plate: `x = overall + row + column + residual`.

Functions¶

`b_scores`(→ Tuple[numpy.ndarray, MedianPolish, float])	B-score of every well on one plate (Brideau et al. 2003).
`call_hits`(→ numpy.ndarray)	Which scores pass `threshold` in `direction`. NaN is never a hit.
`chi_pairwise`(raw_counts[, verbose])	Run pairwise chi-square (or Fisher's exact) tests across group pairs.
`choose_p_adjust_method`(num_groups, num_data_points)	Recommend a multiple-comparison correction method for the given design.
`hit_heatmap`(result, method, *[, target])	Every plate's `method` scores on one diverging scale, hits outlined.
`hit_table`(→ pandas.DataFrame)	The ranked hit table: sample wells, strongest first.
`mad`(→ float)	Median absolute deviation of the finite entries of `values`.
`median_polish`(→ MedianPolish)	Tukey's two-way median polish, NaN-aware, following R's `medpolish`.
`perform_levene_test`(df, grouping_column, data_column)	Levene's test for equal variance, MEDIAN-centred.
`perform_normality_tests`(df, grouping_column, data_columns)	Report per-group normality, and say when the check had no power.
`perform_posthoc_tests`(df, grouping_column, ...)	Run pairwise post-hoc tests across groups with p-value adjustment.
`perform_statistical_tests`(df, grouping_column, ...[, ...])	Run a supported group comparison for each data column.
`robust_z_scores`(→ numpy.ndarray)	Robust z against a reference: `(x - median) / (1.4826 MAD)`.
`score_arrayed_screen`([positive_levels, ...])	`screen_wells()` then `score_screen()`, in one call.
`score_screen`(→ ArrayedHitResult)	Score every well of an arrayed screen against its negative control.
`screen_wells`([positive_levels, negative_wells, ...])	Collapse a measurement table to one row per well, with each well's role.
`ssmd_replicated`(→ float)	SSMD of one treatment from its replicate paired differences.
`ssmd_unreplicated`(→ numpy.ndarray)	SSMD of single wells against a negative reference, without replicates.
`treatment_ssmd`(→ pandas.DataFrame)	Replicate SSMD per treatment from the sample wells' paired differences.
`write_hit_report`(→ Dict[str, str])	Write the scored screen: CSV tables and one plate heatmap per method.

Module Contents¶

exception spacr.sp_stats.HitScoringError[source]¶

Bases: ValueError

Raised when a screen cannot be scored, with the reason and the way out.

Initialize self. See help(type(self)) for accurate signature.

class spacr.sp_stats.ArrayedHitResult[source]¶

Per-well scores, per-treatment SSMD and per-plate QC of one screen.

Parameters:

wells – one row per well with robust_z, ssmd, b_score, difference, the hit_<method> flags, hit (by options['rank_by']) and rank.
treatments – replicate SSMD per treatment; empty without a treatment column.
plates – per-plate summary including Z’.
options – the settings the scores were computed with.
notes – sentences about what could not be scored and why.

hits() → pandas.DataFrame[source]¶: The called sample wells, strongest first.

report() → str[source]¶: The result in sentences, for a text panel or a log.

class spacr.sp_stats.MedianPolish[source]¶

Tukey’s median polish of one plate: x = overall + row + column + residual.

Parameters:

overall – the fitted grand effect.
row – one effect per plate row; NaN for a row with no fitted well.
column – one effect per plate column; NaN likewise.
residuals – the residual matrix; NaN where the input was NaN.
iterations – sweeps run.
converged – whether the relative change fell under the tolerance.

fitted() → numpy.ndarray[source]¶: The additive fit overall + row + column on the full grid.

spacr.sp_stats.b_scores(matrix, *, fit_mask=None, scale: float | None = None) → Tuple[numpy.ndarray, MedianPolish, float][source]¶

B-score of every well on one plate (Brideau et al. 2003).

Row and column effects are fitted by median_polish() on the wells in fit_mask (normally the sample wells, so a control column does not become a column effect), then removed from every well on the plate. The residuals are divided by the scaled MAD of the fitted wells’ residuals, as cellHTS2 does; Brideau’s unscaled MAD differs by the constant 1.4826 only and ranks the wells identically. A row or column with no fitted well, such as a column holding only controls, has no effect to remove, so its wells are scored against the overall and the other effect alone.

Parameters:

matrix – the plate, NaN for an absent well.
fit_mask – boolean matrix of the wells the effects are fitted on; None fits on every finite well.
scale – divide by this instead of the plate’s own residual MAD, which is how the pooled scope applies the screen-wide MAD.

Returns:

(scores, polish, scale_used); a score is NaN for an absent well, and everywhere when the fitted residuals have no spread.

spacr.sp_stats.call_hits(scores, threshold: float, direction: str = 'both') → numpy.ndarray[source]¶

Which scores pass threshold in direction. NaN is never a hit.

Parameters:

scores – the statistic per well.
threshold – the cut-off; its sign is ignored.
direction – both (|s| >= t), up (s >= t) or down (s <= -t).

Returns:

a boolean array.

Raises:

HitScoringError – for an unknown direction.

spacr.sp_stats.chi_pairwise(raw_counts, verbose=False)[source]¶

Run pairwise chi-square (or Fisher’s exact) tests across group pairs.

Uses Fisher’s exact for 2x2 contingency tables and chi-square otherwise, then applies a multiple-comparison correction selected via choose_p_adjust_method().

Two degenerate inputs used to crash rather than report, and both are routine for a sparse per-well contingency table:

Fewer than two groups. There is no pair to compare, so the p-value correction was handed an empty list and raised ZeroDivisionError. An empty result frame is the correct answer, not an exception.
A category no group observed, or a group with no observations at all. chi2_contingency computes an expected frequency of zero and raises ValueError. A category with zero counts on both sides of a pair carries no information about that pair, so it is dropped before testing – which is the standard handling, not a fudge. If dropping leaves fewer than two categories, or either group is empty, the test is genuinely undefined and the pair is reported with a NaN p-value and a reason instead of being silently omitted.

Parameters:

raw_counts – Contingency-table DataFrame indexed by group.
verbose – When True, print the resulting DataFrame.

Returns:

DataFrame with Group 1, Group 2, Test Name, p-value, p-value_adj, adj and note. Empty (with those columns) when there is no pair to compare.

spacr.sp_stats.choose_p_adjust_method(num_groups, num_data_points)[source]¶

Recommend a multiple-comparison correction method for the given design.

Parameters:

num_groups – Number of unique groups being compared.
num_data_points – Number of data points per group (balanced groups assumed).

Returns:

One of 'holm', 'fdr_bh', 'sidak', or 'bonferroni'.

spacr.sp_stats.hit_heatmap(result: ArrayedHitResult, method: str, *, target: str | None = None)[source]¶

Every plate’s method scores on one diverging scale, hits outlined.

Drawn by spacr.figures.plates.build_plates(), the house plate small multiple, with spacr.figures.plates.score_ramp() and a scale symmetric about zero, so the same colour is the same score on every plate and on either side of the negative control.

Parameters:

result – the scored screen.
method – one of HIT_METHODS.
target – 'screen' or 'print'; default the preference.

Returns:

(figure, panel) as build_plates returns them.

Raises:

HitScoringError – for an unknown method.

spacr.sp_stats.hit_table(scored: pandas.DataFrame, *, hits_only: bool = True) → pandas.DataFrame[source]¶

The ranked hit table: sample wells, strongest first.

Parameters:

scored – the wells of an ArrayedHitResult.
hits_only – keep only the called wells.

Returns:

the ranked rows with the identifying columns first.

spacr.sp_stats.mad(values, *, scale: bool = True) → float[source]¶

Median absolute deviation of the finite entries of values.

Parameters:

values – numbers; NaN and infinities are ignored.
scale – multiply by MAD_SCALE so the result estimates a normal standard deviation, as R’s mad() does.

Returns:

the MAD, or NaN with no finite value.

spacr.sp_stats.median_polish(matrix, *, max_iter: int = MEDIAN_POLISH_MAX_ITER, eps: float = MEDIAN_POLISH_EPS) → MedianPolish[source]¶

Tukey’s two-way median polish, NaN-aware, following R’s medpolish.

The sweep order, the centring of the effects and the stopping rule are R’s, so the fit matches stats::medpolish(x, na.rm = TRUE) and the cellHTS2 B-score built on it. A row or column with no finite well gets a NaN effect instead of a guessed one.

Parameters:

matrix – 2-D array, NaN for a well left out of the fit.
max_iter – most sweeps.
eps – stop when |sum|r| - previous| < eps * sum|r|.

Returns:

the MedianPolish.

Raises:

HitScoringError – for an input that is not two-dimensional.

spacr.sp_stats.perform_levene_test(df, grouping_column, data_column)[source]¶

Levene’s test for equal variance, MEDIAN-centred.

Delegates to spacr.figures.stats.check_equal_variance(). Two things moved when it did, and both change the number a caller writes into a CSV:

The centring is the median (Brown-Forsythe), not SciPy’s default mean. Median centring is less sensitive to non-normal data, and this function is called before the normality verdict is known.
Below spacr.figures.stats.MIN_N_FOR_ASSUMPTIONS observations in the smallest group the result is (nan, nan). On three replicates Levene has almost no power, so “p = 0.7, variances are equal” means “we could not tell”, and printing 0.7 into a results table invites exactly the reading that publishes a difference that is not there.

Parameters:

df – Input DataFrame containing the grouping and value columns.
grouping_column – Column name identifying the group of each row.
data_column – Numeric column to test.

Returns:

Tuple (statistic, p_value), both NaN when the check had no power.

spacr.sp_stats.perform_normality_tests(df, grouping_column, data_columns)[source]¶

Report per-group normality, and say when the check had no power.

The VERDICT and the reported ROWS both come from spacr.figures.stats.check_normality(), so the summary and the detail cannot drift apart. That check is Shapiro-Wilk against a Bonferroni threshold across the groups, and it refuses – reporting NaN and Informative=False – when the smallest group is below spacr.figures.stats.MIN_N_FOR_ASSUMPTIONS. A row whose statistic is NaN is not a failed computation; it is the check saying it could not see.

This module used to run D’Agostino-Pearson or Shapiro per group and read “not rejected” as “normal”, which on three replicates is a decision the data cannot support. The p-values it printed for such groups looked perfectly reasonable, which is why the defect survived.

Groups with fewer than three observations are still reported as 'Skipped': Shapiro-Wilk genuinely cannot run on two points.

Parameters:

df – Input DataFrame containing the grouping and value columns.
grouping_column – Column name identifying the group of each row.
data_columns – Iterable of numeric column names to test.

Returns:

Tuple (is_normal, results). is_normal is True only when every requested column passes – it used to be the verdict for the LAST column examined, so a two-column call answered about the wrong one. results is a list of per-group dicts carrying Comparison, Test Statistic, p-value, Test Name, Column, n, Informative and Verdict.

spacr.sp_stats.perform_posthoc_tests(df, grouping_column, data_column, is_normal)[source]¶

Run pairwise post-hoc tests across groups with p-value adjustment.

Uses Tukey HSD when data is normal, Dunn’s test otherwise with a correction method chosen by choose_p_adjust_method().

is_normal should come from perform_normality_tests(), which is the one engine’s verdict. Passing a hand-computed one puts the omnibus test and the pairwise tests on different footing – Kruskal-Wallis across the groups followed by Tukey between them is two different assumptions about one dataset.

Parameters:

df – Input DataFrame containing the grouping and value columns.
grouping_column – Column name identifying the group of each row.
data_column – Numeric column to compare across groups.
is_normal – Whether the data satisfy the normality assumption.

Returns:

List of dicts with pairwise comparison metadata and p-values.

spacr.sp_stats.perform_statistical_tests(df, grouping_column, data_columns, paired=False)[source]¶

Run a supported group comparison for each data column.

Parameters:

df (pandas.DataFrame) – Data containing the grouping and numeric value columns.
grouping_column (str) – Column identifying each observation’s group.
data_columns (iterable of str) – Numeric columns to test.
paired (bool, default=False) – Request paired analysis. Paired analysis is not implemented; when enabled, no result rows are returned.

Returns:

list of dict – Per-column test name, statistic, p-value, sample counts, effect size, and selection rationale. Refused comparisons use Test Name='not testable' and include the reason.

Notes

spacr.figures.stats.compare() selects Student’s t, Welch’s t, Mann-Whitney U, one-way ANOVA, Welch’s ANOVA, or Kruskal-Wallis from the available groups and informative assumption checks.

spacr.sp_stats.robust_z_scores(values, reference) → numpy.ndarray[source]¶

Robust z against a reference: (x - median) / (1.4826 MAD).

With the negative-control wells as reference this is Zhang’s z* score: how far a well is from the negative control, in robust standard deviations of the negative control.

Parameters:

values – the wells to score.
reference – the reference wells, normally the negative controls.

Returns:

one score per value; NaN everywhere when the reference has fewer than two finite wells or no spread.

spacr.sp_stats.score_arrayed_screen(frame: pandas.DataFrame, value_col: str, *, plate_column: str | None = None, control_column: str | None = None, negative_levels=(), positive_levels=(), negative_wells=None, positive_wells=None, treatment_column: str | None = None, grouping: str = 'mean', min_count: int = 0, **options) → ArrayedHitResult[source]¶

screen_wells() then score_screen(), in one call.

Example

from spacr.sp_stats import score_arrayed_screen, write_hit_report
result = score_arrayed_screen(
    df, 'cell_area', negative_wells='c1', positive_wells='c24',
    treatment_column='gene', rank_by='ssmd', scope='plate')
write_hit_report(result, 'results/hits')

Parameters:

frame – the measurement table.
value_col – the measurement to score.
options – passed to score_screen() (scope, ssmd_estimator, replicate_estimator, thresholds, direction, rank_by).

Returns:

an ArrayedHitResult.

spacr.sp_stats.score_screen(wells: pandas.DataFrame, *, scope: str = 'plate', ssmd_estimator: str = 'mm', replicate_estimator: str = 'umvue', thresholds: Dict[str, float] | None = None, direction: str = 'both', rank_by: str = 'ssmd') → ArrayedHitResult[source]¶

Score every well of an arrayed screen against its negative control.

Robust z and SSMD are taken against the negative-control wells of the well’s own plate (scope='plate') or of every plate together (scope='pooled'). The B-score is always fitted per plate, because row and column effects belong to a plate; the pooled scope divides by the screen-wide residual MAD instead of each plate’s own. Controls are scored too, which is how a positive control shows the assay window, but only sample wells are called as hits.

Parameters:

wells – the table screen_wells() returns.
scope – one of HIT_SCOPES.
ssmd_estimator – per-well SSMD estimator, one of SSMD_ESTIMATORS.
replicate_estimator – SSMD estimator for treatments with replicates.
thresholds – cut-off per method; missing ones use DEFAULT_HIT_THRESHOLDS.
direction – one of HIT_DIRECTIONS.
rank_by – the method that decides hit and rank.

Returns:

an ArrayedHitResult.

Raises:

HitScoringError – for an unknown option or a screen without negative-control wells.

spacr.sp_stats.screen_wells(frame: pandas.DataFrame, value_col: str, *, plate_column: str | None = None, control_column: str | None = None, negative_levels=(), positive_levels=(), negative_wells=None, positive_wells=None, treatment_column: str | None = None, grouping: str = 'mean', min_count: int = 0) → pandas.DataFrame[source]¶

Collapse a measurement table to one row per well, with each well’s role.

Wells are located by the plate-QC reader (spacr.plate_qc: prc, a rowID/columnID pair or a well column), so this and the Plate Viewer agree about where every object sits. Controls are named by a column and its levels, as the control-chart screen names them, or by plate position in the negative_control_wells / positive_control_wells notation of spacr.well_spec (c1, r1, A01), or both.

Parameters:

frame – per-object or per-well table.
value_col – the measurement to score.
plate_column – column naming the plate; default the reader’s choice.
control_column – column holding the control labels.
negative_levels – level(s) of control_column that are the negative control.
positive_levels – level(s) that are the positive control.
negative_wells – well spec of the negative-control wells.
positive_wells – well spec of the positive-control wells.
treatment_column – column naming what is in each well, for the replicate SSMD; each well keeps its most common value.
grouping – mean or median of the objects in a well.
min_count – drop wells with fewer objects than this.

Returns:

one row per well: plateID, well, row_index, column_index, prc, n, value, role and, when asked, treatment.

Raises:

HitScoringError – when the table cannot be scored.

spacr.sp_stats.ssmd_replicated(differences, estimator: str = 'umvue') → float[source]¶

SSMD of one treatment from its replicate paired differences.

differences are D_j = x_j - median(negative reference) for each replicate j, the reference taken on the replicate’s own plate (Zhang 2011, screens with replicates):

mm: mean(D) / sd(D);
umvue: sqrt(2/(n-1)) Gamma((n-1)/2) / Gamma((n-2)/2) mean(D) / sd(D), needing n >= 3;
robust: median(D) / (1.4826 MAD(D)).

Parameters:

differences – the replicate differences of one treatment.
estimator – one of SSMD_ESTIMATORS.

Returns:

the SSMD, or NaN with too few replicates or no spread.

Raises:

HitScoringError – for an unknown estimator.

spacr.sp_stats.ssmd_unreplicated(values, reference, estimator: str = 'mm') → numpy.ndarray[source]¶

SSMD of single wells against a negative reference, without replicates.

Zhang (2011), assuming a tested well has the variability of the negative reference N:

mm: (x - mean_N) / (sqrt(2) s_N);
umvue: (x - mean_N) / (sqrt(2 (n_N - 1) / K) s_N) with K = 2 (Gamma((n_N-1)/2) / Gamma((n_N-2)/2))^2, needing n_N >= 3;
robust (SSMD*): (x - median_N) / (1.4826 sqrt(2) MAD_N).

Parameters:

values – the wells to score.
reference – the negative-reference wells.
estimator – one of SSMD_ESTIMATORS.

Returns:

one SSMD per value; NaN where the reference cannot support the estimator (too few wells, or no spread).

Raises:

HitScoringError – for an unknown estimator.

spacr.sp_stats.treatment_ssmd(scored: pandas.DataFrame, *, estimator: str = 'umvue', threshold: float = DEFAULT_HIT_THRESHOLDS['ssmd'], direction: str = 'both') → pandas.DataFrame[source]¶

Replicate SSMD per treatment from the sample wells’ paired differences.

Parameters:

scored – the wells of an ArrayedHitResult; without a treatment column the result is empty.
estimator – one of SSMD_ESTIMATORS.
threshold – SSMD cut-off for the hit call.
direction – one of HIT_DIRECTIONS.

Returns:

one row per treatment, strongest first: treatment, n_replicates, plates, mean_difference, ssmd_mm, ssmd_umvue, ssmd_robust, ssmd (the chosen estimator), median_robust_z, median_b_score, hit and rank.

spacr.sp_stats.write_hit_report(result: ArrayedHitResult, out_dir, *, methods: Sequence[str] = HIT_METHODS, target: str | None = None) → Dict[str, str][source]¶

Write the scored screen: CSV tables and one plate heatmap per method.

Figures go through spacr.plot.save_figure(), so their format, resolution and print repaint follow the user’s figure preferences. The low-contrast colour warning is off for these: the centre of a diverging score map is near the page colour on purpose, because a score of zero is the negative control and is meant to recede.

Parameters:

result – the scored screen.
out_dir – folder to write into; created if absent.
methods – which heatmaps to draw.
target – figure target passed to hit_heatmap().

Returns:

{name: path} of everything written: hit_table (the ranked hits), hit_scores_wells, hit_plates, hit_treatments with a treatment column, and hit_heatmap_<method> per figure.