API reference

The importable analysis API. The package is built stage by stage, so this reference grows as each pipeline stage is added. It currently covers the command-line interface, the run and caching infrastructure, the cohort abstraction, feature typing, the mixture-model wrapper, the enrichment and alignment used to recover and name the reference classes, the model selection, stability, and cross-cohort replication that test the recovered solution, and the stratification axes, binning policies, and acceptance requirements that the stratified analysis is judged against before its bins are frozen.

Command-line interface

analysis command-line interface (Typer).

One subcommand per pipeline stage. Each implemented stage reads named inputs and writes named outputs plus a manifest under artefacts/<stage>/<run-hash>/, so a later run recomputes only what changed (plan section 11). Stages not yet written are grouped under a “planned” panel in analysis --help and exit non-zero; a command leaves that panel when its stage is implemented.

analysis.cli.cohort(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, as_of=<typer.models.OptionInfo object>, sample_n=<typer.models.OptionInfo object>, sample_seed=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Build the harmonised proband-by-feature matrix and its typing manifest.

analysis.cli.fit(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_components=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, as_of=<typer.models.OptionInfo object>, sample_n=<typer.models.OptionInfo object>, sample_seed=<typer.models.OptionInfo object>, no_covariates=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Fit the reference general finite mixture model and predict class labels.

analysis.cli.align(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_components=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, as_of=<typer.models.OptionInfo object>, sample_n=<typer.models.OptionInfo object>, sample_seed=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Compute the seven-category signature and align our classes to Litman’s named classes.

analysis.cli.select(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, k_min=<typer.models.OptionInfo object>, k_max=<typer.models.OptionInfo object>, n_iterations=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, cv=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, as_of=<typer.models.OptionInfo object>, sample_n=<typer.models.OptionInfo object>, sample_seed=<typer.models.OptionInfo object>, workers=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Grid over the number of components and report the information criteria.

analysis.cli.stability(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, mode=<typer.models.OptionInfo object>, n_fits=<typer.models.OptionInfo object>, top_k=<typer.models.OptionInfo object>, n_reps=<typer.models.OptionInfo object>, frac=<typer.models.OptionInfo object>, sub_n_init=<typer.models.OptionInfo object>, n_components=<typer.models.OptionInfo object>, ref_n_init=<typer.models.OptionInfo object>, ref_seed=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, as_of=<typer.models.OptionInfo object>, sample_n=<typer.models.OptionInfo object>, sample_seed=<typer.models.OptionInfo object>, workers=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Summarise multi-initialisation or subsampling stability of the reference fit.

analysis.cli.nmin(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, sizes=<typer.models.OptionInfo object>, n_reps=<typer.models.OptionInfo object>, benchmark=<typer.models.OptionInfo object>, sweep_n_init=<typer.models.OptionInfo object>, n_components=<typer.models.OptionInfo object>, ref_n_init=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, as_of=<typer.models.OptionInfo object>, sample_n=<typer.models.OptionInfo object>, sample_seed=<typer.models.OptionInfo object>, workers=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Find the minimum viable stratum size by refitting at descending sample sizes.

analysis.cli.replicate(version=<typer.models.OptionInfo object>, ssc_version=<typer.models.OptionInfo object>, n_components=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, n_permutations=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, as_of=<typer.models.OptionInfo object>, sample_n=<typer.models.OptionInfo object>, sample_seed=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Project the SPARK model onto the SSC and correlate the category profiles.

analysis.cli.strata_describe(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, quantile_bins=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Characterise the stratification axes and test the candidate binning policies.

A phase-3 (pre-registration) stage. It builds age at diagnosis, the derived diagnostic era, and the measurement-to-diagnosis lag for the modelling cohort, then evaluates each binning policy on both axes against the acceptance requirements (analysis.requirements): the substantive fixed bands, an equal-frequency quantile split, and the max-equal split that is the chosen primary scheme. No model is fitted; the output feeds the frozen bin choice (plan sections 7 and 12).

analysis.cli.strata(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Assign each proband to an age-at-diagnosis and a diagnostic-era stratum.

A phase-4 stage. It builds the two stratifying axes for the modelling cohort and assigns every proband to a stratum on each axis with the frozen primary policy, MaxEqualBins(min_bin_size) (plan section 12a): the finest equal-frequency split that keeps every bin at or above the floor. The per-proband assignments and the realised bin edges are cached for the stratify stage to consume.

analysis.cli.stratify(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, limit=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, workers=<typer.models.OptionInfo object>, shared_with=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Re-estimate the GFMM independently within each stratum of an axis.

A phase-4 stage. For the chosen axis it assigns the modelling cohort to the frozen MaxEqualBins(min_bin_size) strata, then fits the four-class covariate GFMM within each stratum on the same features, typing, and hyperparameters as the reference. Each stratum’s fitted model, hard labels, and class-by-feature centroids are stored so the drift analysis is a pure consumer that never refits. The strata are independent fits, so they run concurrently over a ProcessPoolExecutor (measured: a StepMix fit is single-core, so throughput comes from running many at once, not from one fit using many cores), each pinned to a single BLAS thread (analysis.profiling.single_threaded_blas()) so concurrent workers do not oversubscribe the machine. Every fit is measured (analysis.profiling) and streamed to a resumable checkpoint, so an interrupt continues from the first unfitted stratum. --limit fits only the first few strata, to debug the instrumented pipeline or pilot it before the full run. --shared-with fits on the cross-cohort shared feature set (plan section 8), so the SSC, which provides only a subset of the 238 features, can be re-estimated within its cognitive strata and compared against the reference cohort.

analysis.cli.drift(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, n_permutations=<typer.models.OptionInfo object>, alignment=<typer.models.OptionInfo object>, distance=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, workers=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, reference_scheme=<typer.models.OptionInfo object>, pairwise_mode=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Align stratum classes to the reference and measure drift against the permutation null.

A phase-4 stage (the confirmatory core, frozen in plan section 12a), split so the analysis is decoupled from the fitting. The heavy half (the drift-null sub-stage) refits within size-shuffled pseudo-strata and stores each fit’s centroids and reference contingency, keyed only by the fitting parameters. This cheap half then aligns and measures: it reads the observed per-stratum fits (the stratify stage) and the stored null fits, applies the chosen --alignment and --distance over them, and reads each aligned class’s drift against its size-matched null (beyond the 95th percentile, Benjamini-Hochberg controlled). Because the fits are stored, changing the alignment or distance re-measures without re-fitting. The alignment confidence (per-class Jaccard, overall adjusted Rand index) is reported alongside, so a large shift with low overlap reads as reorganisation, not drift.

analysis.cli.order(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, k_anchor=<typer.models.OptionInfo object>, k_cap=<typer.models.OptionInfo object>, b_screen=<typer.models.OptionInfo object>, b_escalate=<typer.models.OptionInfo object>, escalate_threshold=<typer.models.OptionInfo object>, cv_n_init=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, workers=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Test whether the supported number of latent classes is stable across strata of an axis.

The H0C hypothesis (plan section 7, H0C). Per stratum, the supported number of classes is found by a warm-started parametric bootstrap likelihood-ratio search anchored at four (analysis.order), and read relative to the pooled cohort put through the identical procedure, so the shared over-extraction cancels: an order change is a stratum whose supported order differs from the pooled cohort’s. The strata are the four equal-frequency quantile bins of the axis (a fifth class in a thousand-proband bin is not estimable, so the fine MaxEqualBins scheme is deliberately not used here), plus, for the era axis, the DSM-IV against DSM-5 (2013) split as one targeted secondary pair. The pooled search is cached as its own stage so both axes reuse it. Bootstrap draws run concurrently and stream to a resumable checkpoint. A positive order-change claim needs agreement: the BLRT order differs from the pooled order, the cross-validated elbow knee moves off it, and the adjusted Lo-Mendell-Rubin test agrees.

analysis.cli.sweep(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, scheme=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, n_permutations=<typer.models.OptionInfo object>, alignment=<typer.models.OptionInfo object>, distance=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, workers=<typer.models.OptionInfo object>, no_covariates=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Run every localisation scheme end to end and read them against one shared null.

The phase-4 conductor. For each --scheme (hard bins, the frozen primary, and the kernel LSEM trajectory, by default) it fits the observed sweep and the axis-permutation null, aligns each local fit to the pooled reference, measures its drift with the chosen --distance, and reads it against the null, writing one combined decision table with a scheme and position column so the arms plot as trajectories side by side. Nothing here is frozen in code: the primary scheme is a choice of which row of the table to privilege, so the pipeline stays flexible while the whole battery is computed on every run. This is the heavy stage: the null is n_permutations full sweeps, so the confirmatory both-scheme run is a cluster job; pass --n-permutations 1 to smoke the chain or 100 for an overnight pilot.

analysis.cli.bandwidth(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_points=<typer.models.OptionInfo object>, targets=<typer.models.OptionInfo object>, reduce=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Choose a kernel bandwidth by the effective sample size of its focal fits.

A kernel focal fit is worth the sum of its weights in whole probands, its effective sample size, which sets its power the way a bin’s count sets a hard-bin fit’s. This stage inverts that: for each target effective size it reports the bandwidth whose focal fits reach it, so a kernel window can be set to carry the same power as a hard bin at the recovery floor. With --reduce min (the default) every focal fit clears the target, the same guarantee the hard-bin floor gives every bin; --reduce median fixes the typical focal point instead and lets the edges thin. The bandwidth is in the axis units (years for age at diagnosis, years for the diagnosis year). No model is fitted, so the stage is fast; the reported bandwidths feed a kernel sweep, for example --scheme kernel:<bandwidth>:<n_points>.

analysis.cli.trajectory(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, n_shuffle=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, from_sweep=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Embed the classes and measure how each one’s centroid moves across the strata.

A phase-4 presentation stage that consumes the stratified fits (analysis stratify) and writes only aggregate, class-level outputs. It aligns each stratum’s classes to the named reference by membership (reusing analysis.drift), fits a linear-discriminant embedding of the four pooled classes, projects the reference anchors and the aligned stratum centroids into it, and quantifies each class’s trajectory: a directional test (net young-to-old displacement against an ordering-shuffle null) and a roughness measure (step size against the sampling-noise expectation). The directional test is a pilot on the observed centroids; the confirmatory test is the continuous-trend regression against the refit permutation null (plan section 12a). Nothing per-proband is written, so the outputs may be promoted to the manuscript and the docs.

analysis.cli.attribute(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, alignment=<typer.models.OptionInfo object>, decomposition=<typer.models.OptionInfo object>, contrast=<typer.models.OptionInfo object>, feature_space=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Attribute each stratum class’s movement to features and probands (archived, refit-era).

Archived. This is the refit-era attribution stage. The category attribution ($H_0^F$) is now read from the single cached fit by the block engine (analysis.blocks), so this stage is no longer $H_0^F$’s evidence; it is kept to render the membership-churn and mover figures on the refit-pilot archive page.

A phase-4 interpretation stage and a pure consumer of the reference fit and the stratify fits (no re-fitting). For every stratum it aligns the fit to the reference, splits each aligned class’s centroid shift into per-feature contributions (analysis.attribution), and contrasts the probands that left the class against those that stayed. It writes the per-feature decomposition, the per-category totals, the mover-versus-stayer contrast, and a per-class headline table, so a movement reads as which features and which people carry it. The decomposition and contrast are descriptive readouts of an already-measured drift, so they sit outside the section-12a confirmatory freeze.

analysis.cli.invariance(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_simulations=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, max_grid=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, q=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Test the reference class profiles for stability along an axis, from the single cached fit.

The score-based measurement-invariance test (plan section 7e): for the measurement-only reference, each proband’s casewise score on every class-conditional location parameter is cumulated in axis order into an empirical fluctuation process, standardised to a Brownian bridge under stability. The maxLM and Cramer-von Mises functionals are read against an analytic (simulated-bridge) null, per focal block (each class, and each class crossed with a feature category), with Benjamini-Hochberg control across blocks. No mixture is refitted.

Unlike the drift stage, this consumes the fitted model itself (its parameters and responsibilities), not only the labels, and it reads the marginal (measurement-only) reference, so its estimand matches the kernel and pairwise arms rather than the covariate fit.

analysis.cli.invariance_trajectory(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_points=<typer.models.OptionInfo object>, n_boot=<typer.models.OptionInfo object>, controls=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, q=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Read how each class profile drifts along an axis as a null-free, separation-scaled effect.

The recast of the score-based invariance test (plan section 7e). Freezing the measurement-only reference’s responsibilities, it reads each class’s local centroid as a smooth function of the axis, forms the per-feature displacement from the pooled centroid, and scales magnitudes by the between-class separation. Uncertainty is a family-clustered bootstrap (the tube and per-feature intervals), replacing the saturated bridge null. It reports the in-plane capture fraction so the 2D figure cannot hide out-of-plane drift, a covariance-aware Mahalanobis corroboration, and a control-panel specificity comparison (era and age against household income, area deprivation, and a random ordering). No mixture is refitted; every artefact is class or feature level.

analysis.cli.displacement_atlas(dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_points=<typer.models.OptionInfo object>, min_coverage=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Screen every continuous or ordered non-modelling axis for per-class endpoint displacement.

The generalisation of the specificity panel (plan section 12b). The timing axes are two orderings of the cohort; this stage reads each class’s separation-scaled endpoint displacement along every axis in the catalogue (analysis.axes: the timing axes and the covariate pool), so the atlas figure can sort them from the largest mover to the smallest against the random-order floor. No mixture is refitted and no orthogonality assumption is made about any covariate; the random ordering is the only reference. The 238 clustered features, totals over them, and held-out phenotype instruments are excluded as circular or as a non-null phenotype ceiling. Every artefact is class or axis level.

analysis.cli.demographic_conditioning(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_points=<typer.models.OptionInfo object>, min_coverage=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, include_timing=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Report how much of each class’s timing drift each demographic covariate accounts for.

The conditioning half of the demographic screen (plan section 7g). For every covariate in analysis.demographics, the reference features are residualised on the covariate and the class endpoint displacement along the timing axis is re-read; the shrinkage is the fraction of that class’s drift the covariate linearly accounts for (analysis.blocks.conditioning_shrinkage()). Each covariate’s linear span of the axis is reported beside its shrinkage, because a covariate orthogonal to the axis cannot account for an axis-ordered drift. This is a descriptive partial association, not a causal claim, and it refits nothing; the confirmatory hypotheses are unchanged.

With --include-timing the screen appends age at evaluation and the measurement-to-diagnosis lag (analysis.demographics.TIMING_COVARIATES). These correlate with age at diagnosis by construction, so their shrinkage is not ceilinged near zero as the demographics’ is, and they adjudicate whether the age-at-diagnosis drift reduces to the general age effect the reference model already carries in age at evaluation. The flag changes the run hash, so the demographic-only artefacts stay valid.

analysis.cli.invariance_responsibilities(model_obj, x_values)

Return the frozen posterior responsibilities of the measurement-only reference.

analysis.cli.prevalence(axis=<typer.models.OptionInfo object>, dataset=<typer.models.OptionInfo object>, version=<typer.models.OptionInfo object>, n_points=<typer.models.OptionInfo object>, n_boot=<typer.models.OptionInfo object>, adjusted=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, n_init=<typer.models.OptionInfo object>, q=<typer.models.OptionInfo object>, min_bin_size=<typer.models.OptionInfo object>, force=<typer.models.OptionInfo object>)

Test whether the frozen class proportions trend along an axis (H0B).

The prevalence-drift test. The four classes are held fixed at the measurement-only reference fit (no mixture is refitted); their mixing proportions are regressed on the axis. The rigorous read is a maximum-likelihood three-step correction that removes the classify-analyse bias of a hard-label regression by fixing the classification-error matrix of the frozen posteriors; a naive hard-label multinomial logit is reported beside it as an uncorrected cross-check. It reports, per class, the one-versus-rest axis log-odds slope and odds ratio with a family-clustered bootstrap interval and $p$, the naive Wald and likelihood-ratio $p$-values, the joint likelihood-ratio test of class ~ axis against class ~ 1, the predicted proportion curve with its band, an adjusted axis slope net of sex, the lag, and age at evaluation, and (for era) the DSM-5 pre/post-2013 contrast. Every artefact is class or coefficient level; no per-proband quantity is written.

analysis.cli.sensitivity()

Re-fit under alternative feature sets and within cognitive-level strata.

analysis.cli.report()

Assemble the non-disclosive tables and figures for the manuscript.

Configuration and paths

Run configuration: constants, the reference pipeline’s fixed choices, and input paths.

This module gathers the values that pin a run down: the cohort the reference fit uses, the covariates and age window from Litman et al., the default mixture-model hyperparameters, and the locations of the author-provided reference inputs (the final feature list and the feature-to-category map) and the released typing pickles.

The reference inputs live under the gitignored .literature tree, because they come from the authors’ released code and our correspondence with them and are not ours to redistribute. Each location can be overridden with an environment variable so a run can point at a copy elsewhere.

analysis.config.author_feature_list(root)

Return the path to the authors’ final feature list (mixture_model_columns.csv).

This 238-feature list is the authoritative set the general finite mixture model was fit on. It resolves the ambiguity the released preprocessing left open, because that code drops columns rather than naming the kept set (plan section 5).

Parameters:

root (pathlib.Path) – The monorepo root.

Returns:

Location of the feature-list CSV. Overridable with ANALYSIS_FEATURE_LIST.

Return type:

pathlib.Path

analysis.config.author_category_map(root)

Return the path to the feature-to-category map (feature_to_category_mapping.csv).

The map assigns each feature to one of the seven literature-defined categories and is used only to summarise results, never to fit (plan section 5).

Parameters:

root (pathlib.Path) – The monorepo root.

Returns:

Location of the category-map CSV. Overridable with ANALYSIS_CATEGORY_MAP.

Return type:

pathlib.Path

analysis.config.litman_typing_dir(root)

Return the directory holding the released typing pickles.

The directory holds binary_columns.pkl, categorical_columns.pkl, and continuous_columns.pkl: the feature-type assignments StepMix’s densities depend on. We reconcile a dictionary-derived typing against these (plan section 6, step 2).

Parameters:

root (pathlib.Path) – The monorepo root.

Returns:

Location of the typing-pickle directory. Overridable with ANALYSIS_TYPING_DIR.

Return type:

pathlib.Path

Locating the monorepo root and the cached-artefact paths.

Each pipeline stage writes its outputs and a manifest under <root>/artefacts/<stage>/<run-hash>/, where the run hash is a content hash over the inputs that determine the output (dataset version, feature list, hyperparameters, seed, stratum definition, and the package commit). The directory is gitignored, because it holds participant-derived intermediates that the SFARI consent does not allow into the committed history.

analysis.paths.find_repo_root(start=None)

Walk up from start (or the working directory) to the monorepo root.

Honours ANALYSIS_ROOT if set. Otherwise the root is the nearest ancestor that holds a data/ directory or a pyproject.toml declaring a uv workspace.

Parameters:

start (pathlib.Path, optional) – Directory to start the search from. Defaults to the working directory.

Returns:

The resolved monorepo root.

Return type:

pathlib.Path

analysis.paths.artefacts_dir(root)

Return the cached-artefacts directory, <root>/artefacts.

analysis.paths.stage_dir(root, stage)

Return the directory holding every run of one stage, <root>/artefacts/<stage>.

analysis.paths.run_dir(root, stage, run_hash)

Return one run’s directory, <root>/artefacts/<stage>/<run-hash>.

analysis.paths.manifest_path(run_dir)

Return the manifest path for a run, <run-dir>/manifest.json.

Runs and caching

Content-addressed cache: run hashes, manifests, and artefact serialization.

A stage’s output is identified by a hash over the inputs that determine it (dataset and version, the feature-list digest, model hyperparameters, the covariate set, the seed, and the stratum definition). A run is a cache hit when a manifest with the same hash already exists and finished cleanly, so a later session recomputes only what changed (plan section 11).

The git commit and the resolved package versions are recorded in each manifest for provenance, but neither enters the hash. The repository keeps the manuscript and the docs site beside the analysis code, so its HEAD moves (and the working tree reads -dirty) on edits that touch no analysis input; folding that commit into every hash would discard the expensive fit and stability caches on each unrelated commit. A change the hash should react to is caught at a finer grain instead: a stage that builds its data inline digests the result with frame_digest(), so a harmonisation edit that alters the data invalidates the cache while one that does not, a comment or a refactor, leaves it valid. The replicate stage does this for the integrated SSC frame.

This module holds the hashing, the manifest read and write, the environment capture, and the serialization helpers. The analysis.run module composes them into the run_context lifecycle.

analysis.cache.canonical_json(obj)

Serialise obj to a stable JSON string.

Keys are sorted and whitespace is removed, so the same logical parameters always produce the same string and therefore the same hash.

Parameters:

obj (object) – Any JSON-serialisable object, plus Path (rendered as its string).

Returns:

The canonical JSON encoding.

Return type:

str

analysis.cache.compute_hash(params)

Return the SHA-256 hex digest of a run’s parameters.

Parameters:

params (Mapping) – The inputs that determine the output (plan section 11).

Returns:

The 64-character hex digest.

Return type:

str

analysis.cache.short_hash(full_hash, length=16)

Return the leading length characters of a hash, for directory names.

analysis.cache.file_digest(path)

Return the SHA-256 hex digest of a file’s bytes.

Used to fold the author feature list and the typing manifest into a run hash, so that editing either input invalidates the cache.

Parameters:

path (pathlib.Path) – File to digest.

Returns:

The 64-character hex digest.

Return type:

str

analysis.cache.frame_digest(df)

Return a stable SHA-256 hex digest of a dataframe’s content.

Hashes the column names and the per-row values (the index included), so a different parse, a renamed column, or a dropped proband yields a different digest. The result does not depend on row order: the row hashes are sorted before they are combined, so two frames holding the same probands and values agree.

Use this to fold a matrix built inline into a run hash. The compute_hash() parameters of the cohort stage cover the input files and settings but not the harmonisation code, so a code-only change (the milestone parser, an SSC rename map) leaves them unchanged; digesting the built frame captures that change instead.

Parameters:

df (pandas.DataFrame) – The frame to digest.

Returns:

The 64-character hex digest.

Return type:

str

analysis.cache.environment_versions()

Return the resolved versions of the packages that fits depend on.

Returns:

Package name mapped to its installed version, or "unknown" when the package is not installed.

Return type:

dict of str to str

analysis.cache.git_commit(root)

Return the short git commit of the repository, or "unknown".

Parameters:

root (pathlib.Path) – Repository root.

Returns:

The short commit hash, suffixed with "-dirty" when the working tree has uncommitted changes, or "unknown" when git is unavailable.

Return type:

str

analysis.cache.save_frame(df, path)

Write a dataframe to Parquet, preserving its index.

analysis.cache.load_frame(path)

Read a dataframe from Parquet.

analysis.cache.save_model(obj, path)

Persist a fitted model (or any picklable object) with joblib.

analysis.cache.load_model(path)

Load a joblib-persisted object.

analysis.cache.save_json(obj, path)

Write obj to a human-readable JSON file.

analysis.cache.load_json(path)

Read a JSON file.

analysis.cache.read_manifest(run_dir)

Return a run’s manifest, or None when it does not exist.

Parameters:

run_dir (pathlib.Path) – The run directory artefacts/<stage>/<hash>.

Returns:

The parsed manifest, or None when absent.

Return type:

dict or None

analysis.cache.write_manifest(run_dir, data)

Write a run’s manifest to run_dir/manifest.json.

The run lifecycle: a content-addressed directory, a captured log, and a manifest.

run_context is the single entry point every expensive stage uses. It hashes the run’s parameters, opens (or reuses) the directory artefacts/<stage>/<hash>, captures the run’s standard output into run.log while still showing it on the console, and writes a manifest.json recording the inputs, status, timing, resolved package versions, the repository commit, and the caller’s metrics.

A run whose manifest already exists and finished cleanly is a cache hit: the caller is told so through RunContext.cache_hit and loads the cached artefacts instead of recomputing (plan section 11). Pass force=True to recompute regardless.

class analysis.run.RunContext(stage, run_hash, run_dir, cache_hit, params, metrics=<factory>, log=<factory>)

Handle to one run: where its artefacts live and what to record about it.

stage

The pipeline stage, used as the artefact subdirectory.

Type:

str

run_hash

The full hash over the run’s parameters.

Type:

str

run_dir

The directory holding this run’s artefacts and manifest.

Type:

pathlib.Path

cache_hit

True when a clean manifest already existed, so the caller should load cached artefacts rather than recompute.

Type:

bool

params

The parameters that determined the hash.

Type:

dict

metrics

Scalar results the caller wants recorded in the manifest (final log-likelihood, class proportions, and so on).

Type:

dict

log

Logger writing to both the console and run.log.

Type:

logging.Logger

path(name)

Return the path to a named artefact inside this run’s directory.

analysis.run.run_context(stage, params, *, root=None, force=False)

Open a run: resolve its directory, capture its log, and manage its manifest.

Parameters:

• stage (str) – The pipeline stage (the artefact subdirectory and logger name).

• params (Mapping) – The inputs that determine the output. Their hash names the run directory and is recorded in the manifest, so editing any input invalidates the cache.

• root (pathlib.Path, optional) – Repository root. Defaults to the discovered root.

• force (bool, default False) – Recompute even when a clean manifest already exists.

Yields:

RunContext – The run handle. When cache_hit is True the body should load cached artefacts; otherwise it computes, writes artefacts under run_dir, and may set metrics for the manifest.

Append-only checkpoint logs that let a long iterative stage resume after an interrupt.

The content-addressed cache (analysis.cache, analysis.run) works at the granularity of a whole stage: a run is reused only once it has finished cleanly, and an interrupted run leaves no reusable output, so re-running recomputes the stage from the start. That is fine for the short stages, but the multi-seed loops (model selection, the multi-initialisation and subsampling stability runs, and the minimum-stratum-size sweep) can take tens of minutes to hours, and losing all of it to one interrupt is wasteful.

A CheckpointLog records each unit of work as it completes, one JSON line per unit, appended and flushed to disk. When the stage runs again over the same parameters (hence the same run directory), it reads the completed units back and continues from the first one that is missing. The seeds are derived deterministically from the unit index, so a resumed run reproduces exactly what an uninterrupted run would have computed; the checkpoint changes only how much is recomputed, never the result.

The unit of resumption is one line. Each line holds the whole payload for one unit (for the selection grid, every criterion row for one seeded iteration; for stability, one fit and its comparison). A process killed mid-write can leave a torn final line; CheckpointLog.load() parses up to the first line that does not decode and stops there, so the last, incomplete unit is dropped and recomputed rather than read back half-written. The logs use Python’s json non-finite extension (NaN is written and read back unchanged), since they are read only by this module.

class analysis.checkpoint.CheckpointLog(path)

An append-only log of completed units of work for one resumable loop.

Each call to append() writes one unit’s payload as a JSON line and flushes it to disk; load() reads the completed units back in order. The payload is any JSON-serialisable value: a caller that produces several records per unit stores the list of them, so one line maps to one resumable unit.

Parameters:

path (pathlib.Path) – The log file. It lives inside the stage’s content-addressed run directory, so it is specific to the run’s parameters and is never shared between different parameter sets.

load()

Return the completed unit payloads in the order they were written.

Parsing stops at the first line that does not decode as JSON, which drops a torn final line left by a process killed mid-write. Because units are appended and flushed one at a time, only the last line can be incomplete.

Returns:

One entry per completed unit, in append order. Empty when the log does not exist.

Return type:

list

append(payload)

Append one unit’s payload as a JSON line and flush it to disk.

The write is flushed and fsync-ed so a completed unit survives an interrupt that kills the process before the stage finishes.

Parameters:

payload (object) – Any JSON-serialisable value describing one completed unit.

clear()

Delete the log file if it exists.

analysis.checkpoint.clear_checkpoints(directory)

Remove every checkpoint log in a run directory.

Called when a stage is forced to recompute (so a stale partial run is not resumed) and once it has finished cleanly (so the now-redundant checkpoints do not linger beside the final artefacts).

Parameters:

directory (pathlib.Path) – A stage’s run directory.

A single progress bar per command, with live state in its postfix.

Each pipeline command opens one bar whose total is the sum of all units of work across every loop it runs (for example the sum of initialisations over a grid of component counts, or the features in an enrichment pass). The bar is updated once per unit and its postfix carries the live state: which stage or stratum is running, the best log-likelihood so far, the smallest class proportion, and so on.

The bar is written to the real standard error so it survives the standard-output redirection that analysis.run uses to capture a run’s log, and so its control characters do not pollute run.log.

analysis.progress.task_bar(total, desc, **kwargs)

Open the command’s single progress bar as a context manager.

Parameters:

• total (int) – Total units of work across every loop the command runs.

• desc (str) – Short label shown to the left of the bar.

• **kwargs (Any) – Forwarded to tqdm.tqdm.

Yields:

tqdm.tqdm – The bar. Call update once per unit of work and set_postfix to publish the live state.

Hardware capture and per-unit resource measurement for the fitting stages.

The stratified fits and the permutation null are the heaviest compute in the pipeline, and how feasible they are depends on the machine that runs them: a laptop overnight, or a Harvard O2 (SLURM) job array. This module measures that cost rather than guessing it. It captures the hardware once per run and the resource use of each fitted unit, so a short calibration run yields the per-fit cost that the full run, and any later SLURM resource request, are projected from.

Two pieces, neither of which knows anything about the model, so the same instrumentation serves the stratified fits, the drift null, and any later stage:

• capture_hardware() records the CPU, core counts, memory, BLAS thread pools, and platform. It is written into the run manifest, so every artefact carries the machine that produced it and a laptop run and an O2 run are directly comparable.

• measure() is a context manager around one unit of work (one fit_gfmm call). It times the unit, reads its process CPU time, samples the resident set size on a background thread to catch the peak (including the native BLAS allocations that tracemalloc cannot see), and records the bytes the unit writes, returning a UnitMetrics.

The four measured quantities map onto the four numbers a SLURM job needs: wall time sets --time, peak resident memory sets --mem, CPU utilisation (CPU time over wall time) shows whether a fit uses more than one core and so sets --cpus-per-task, and the output bytes set the scratch storage the run consumes.

Scope. The CPU-time and memory readings are for the calling process and its threads, which is what a StepMix fit is: BLAS-threaded work in one process, no child processes. A stage that forks worker processes would need the children measured too; that is noted where it would matter rather than handled here, since the fitting stages do not fork.

analysis.profiling.single_threaded_blas()

Force every BLAS and OpenMP thread pool to one thread for the enclosed block.

A lone fit is already single-core (measure() shows cpu_utilisation about one), but scikit-learn bundles its own OpenMP runtime (libomp) sized to every logical core by default. Running several fits at once through a ProcessPoolExecutor gives each worker its own copy of that pool, so N concurrent workers compete for N times the machine’s threads the moment any of them calls a parallelised scikit-learn routine. A spawned worker inherits the parent’s environment, and each numerical library reads its thread count once, the first time it is used, so setting these variables before the pool starts keeps every worker to the one core it was given. Not needed around a solitary fit, so it is applied only around a concurrent pool, not a whole CLI invocation.

Yields:

None

analysis.profiling.capture_hardware()

Capture the hardware and threading configuration of the current machine.

The result is JSON-serialisable and recorded once per run, so every artefact carries the machine that produced it. Each probe is guarded: a field that cannot be read is set to None rather than raising, since the capture is diagnostic and partial information is still useful.

Returns:

Keys: cpu_model, architecture, physical_cores, logical_cores, cpu_freq_mhz, total_memory_bytes, available_memory_bytes, hostname, platform, system, python_version, blas (the loaded BLAS thread pools), and thread_env (the threading environment variables).

Return type:

dict

class analysis.profiling.UnitMetrics(wall_s, cpu_s, peak_rss_bytes, start_rss_bytes, n_samples, output_bytes=None)

The measured cost of one unit of work (one fit).

wall_s

Elapsed wall-clock seconds.

Type:

float

cpu_s

Process CPU seconds (user plus system) spent during the unit. Summed across threads, so a BLAS-threaded fit reports more CPU than wall time.

Type:

float

peak_rss_bytes

Highest resident set size seen while the unit ran, the figure a SLURM --mem request must cover.

Type:

int

start_rss_bytes

Resident set size when the unit started, so the unit’s own allocation can be read as a delta against the process baseline.

Type:

int

n_samples

How many memory samples the background thread took, recorded so a peak built from too few samples is visible rather than silently trusted.

Type:

int

output_bytes

Bytes the unit wrote to disk, set by the caller after the artefacts are saved. Feeds the scratch-storage projection. None when not recorded.

Type:

int or None

property cpu_utilisation: float

CPU seconds per wall second, with 1 one fully-used core and above 1 multi-core.

property peak_rss_delta_bytes: int

Peak resident memory above the process baseline at the unit’s start.

to_dict()

Return a JSON-serialisable record, including the derived ratios.

class analysis.profiling.MeasureHandle(output_bytes=None, metrics=None)

The handle measure() yields, so the caller can attach the unit’s output size.

The metrics are filled in when the context exits; output_bytes is set by the caller inside the block (for example to the size of the artefact it just wrote) and is copied into the final UnitMetrics.

analysis.profiling.measure(*, sample_interval_s=0.05)

Measure the wall time, CPU time, and peak resident memory of the enclosed work.

A background thread samples the resident set size every sample_interval_s seconds and keeps the maximum, so the peak captures the native BLAS allocations a Python-level memory tracer would miss. On exit the handle’s metrics holds the UnitMetrics; set handle.output_bytes inside the block to record what the unit wrote.

Parameters:

sample_interval_s (float, optional) – Seconds between resident-memory samples. The default suits fits that run for seconds.

Yields:

MeasureHandle – The handle whose metrics are populated when the block exits.

Examples

>>> with measure() as unit:
...     result = fit_gfmm(inputs)
...     unit.output_bytes = save(result)
>>> unit.metrics.wall_s
2.41

analysis.profiling.path_bytes(path)

Return the size in bytes of a file, or the recursive size of a directory.

Used to record what a unit wrote, so the per-unit output size projects to the storage a full run consumes. A path that does not exist contributes nothing.

Parameters:

path (pathlib.Path) – A file or directory.

Returns:

Total size in bytes.

Return type:

int

analysis.profiling.summarise(metrics)

Aggregate per-unit metrics into the distribution a projection is built from.

The full-run and SLURM projections multiply these measured per-unit figures by the known unit count, so the spread (median to 90th percentile to maximum) is what turns a calibration sample into a bounded estimate rather than a single guess.

Parameters:

metrics (list of UnitMetrics) – The measured units.

Returns:

n_units; total_wall_s and total_cpu_s; the wall-second median / p90 / max; the peak_rss_bytes median and max; the mean cpu_utilisation; and total_output_bytes (None when no unit recorded a size). Empty input returns zeroes.

Return type:

dict

Cohort abstraction

The cohort abstraction: one interface, a SPARK and an SSC backend behind it.

Every analysis is written once against Cohort and the harmonised CohortMatrix it yields, so the reference fit, the stability checks, and the SSC replication run on either cohort without change (plan section 10). Each backend maps its raw tables onto the shared schema; the SPARK-only timing fields (age at diagnosis, era) are exposed as an optional capability that the SSC backend need not provide.

The shared helpers here resolve a table’s source CSV through the dscat catalogue and read only the columns a stage needs, so a backend never loads a whole file into memory.

class analysis.cohort.CohortMatrix(features, covariates, dataset, version)

A harmonised proband-by-feature matrix with its covariates and provenance.

features

Proband-by-feature matrix, indexed by proband id, holding only the requested clustered features (covariates excluded).

Type:

pandas.DataFrame

covariates

Proband-by-covariate matrix on the same index, holding the structural-model covariates.

Type:

pandas.DataFrame

dataset

Cohort the matrix was built from.

Type:

str

version

Dataset version the matrix was built from.

Type:

str

property feature_names: list[str]

Return the feature column names.

property n_probands: int

Return the number of probands (rows).

class analysis.cohort.Cohort(*args, **kwargs)

The contract every cohort backend satisfies.

A backend integrates its instruments into one harmonised, complete-case proband-by-column frame (covariates and the shared feature schema), and declares whether it can provide the SPARK-only diagnosis-timing fields used for the stratification axes (plan sections 5 and 7).

integrate()

Return the harmonised, complete-case proband-by-column frame.

supports_timing()

Return whether the backend can provide diagnosis-timing fields.

axis(name, index, covariates, min_bin_size=1000)

Return a stratification variable and its default binning policy.

The stratified analysis (plan section 7) re-estimates the mixture within strata of an axis. A backend resolves a named axis to the per-proband variable (on index) and the policy that bins it, so a stage names an axis and never reads a cohort-specific column itself. This generalises supports_timing(): the diagnosis-timing axes (age_at_diagnosis, era) are SPARK-only, while the cognitive axes (cognitive_impairment, iq) are shared, so the SSC backend provides the second pair and none of the first (plan sections 5 and 8).

Parameters:

• name (str) – The axis to resolve.

• index (pandas.Index) – The modelling-cohort proband index the variable is built on.

• covariates (pandas.DataFrame) – The cohort covariates, on index; the timing axes read age at evaluation and sex from here.

• min_bin_size (int, default 1000) – The floor passed to a size-based policy (MaxEqualBins); ignored by the fixed dichotomy.

Returns:

The variable (a pandas.Series) and its BinningPolicy, or None when the backend does not provide name.

Return type:

tuple or None

family_ids(index)

Return the per-proband family identifier, or None when the backend has none.

The clustered bootstrap in analysis.trajectory_local resamples families rather than probands, so a backend that groups probands into families exposes the grouping key here. Like axis(), this is an optional capability: a backend that cannot provide a family key returns None.

analysis.cohort.open_catalogue(root)

Open the dscat catalogue at the repository root.

analysis.cohort.source_csv(cat, root, dataset, version, table, role='')

Return the absolute path to a table’s backing CSV.

Parameters:

• cat (dscat.index.Catalogue) – Open catalogue.

• root (pathlib.Path) – Repository root, used to resolve the catalogue’s repo-relative path.

• dataset (str) – The table to locate.

• version (str) – The table to locate.

• table (str) – The table to locate.

• role (str, default "") – Family role for cohorts that split a measure across role folders (SSC). SPARK tables use the empty role.

Returns:

Absolute path to the CSV.

Return type:

pathlib.Path

Raises:

FileNotFoundError – When no source row matches table and role.

analysis.cohort.csv_columns(path)

Return a CSV’s column names without reading its rows.

analysis.cohort.read_columns(path, columns)

Read only the requested columns of a CSV, skipping any that are absent.

Parameters:

• path (pathlib.Path) – CSV to read.

• columns (collections.abc.Sequence of str) – Columns to keep; columns not present in the file are ignored.

Returns:

The requested columns that exist in the file.

Return type:

pandas.DataFrame

analysis.cohort.build_matrix(integrated, feature_names, dataset, version, covariates=('sex', 'age_at_eval_years'))

Split an integrated frame into a feature matrix and a covariate matrix.

Parameters:

• integrated (pandas.DataFrame) – The backend’s harmonised, complete-case frame.

• feature_names (collections.abc.Sequence of str) – The clustered features to keep.

• dataset (str) – Provenance recorded on the matrix.

• version (str) – Provenance recorded on the matrix.

• covariates (collections.abc.Sequence of str, optional) – Covariate columns to split out. Defaults to the structural-model covariates.

Returns:

The feature and covariate matrices on a shared index.

Return type:

CohortMatrix

Raises:

KeyError – When a requested feature or covariate is absent from integrated.

analysis.cohort.get_cohort(dataset, version, root=None, *, as_of=None)

Return the backend for a dataset.

Parameters:

• dataset (str) – "spark" or "ssc".

• version (str) – Dataset version.

• root (pathlib.Path, optional) – Repository root. Defaults to the discovered root.

• as_of (str, optional) – A records cutoff passed to the SPARK backend (for example "2022-12-12"); it restricts the cohort to the probands present at that freeze. Ignored by cohorts that do not carry the timing fields the cutoff needs (the SSC).

Returns:

The backend.

Return type:

Cohort

Raises:

ValueError – When dataset is not a known cohort.

The shared feature schema and the cross-cohort harmonisation maps.

The feature set is the authors’ final 238-feature list. The CBCL competence items arrive as strings in SPARK and are recoded to the ordinal integers the released preprocessing used. The SSC rename maps carry the second cohort’s column names onto the SPARK names, so both cohorts present the same schema to the rest of the pipeline (plan section 10).

analysis.cohort.schema.parse_age_months(value, *, disambiguate=None)

Parse a free-text developmental-milestone age into months.

The SSC records milestone ages as free text, whereas the SPARK features are ages in months. The recognised forms are mapped onto months: a bare number or a number with a month unit ("13 months", "13 mos", "12 mon", "13m") is taken as months; a number with a year unit is multiplied by twelve, with an optional trailing months part ("1 yr 6 mo"); a number with a week unit ("6 weeks") is converted from weeks; a compound ("2 yrs 10 mos") sums its distinct unit parts; a half- or quarter-year fraction ("3 1/2 yrs"), a "y.o." suffix ("3 y/o"), an "age N" phrase, and a trailing "old" are handled; "at birth" is zero; and a range or an “N or M” ("12-14", "18 months to 2 years", "7 or 8 months") is read as its midpoint. A bound ("<3 mos", "before 1 year") has the bound dropped and the stated age taken.

A statement that the milestone was never reached ("never", "not yet", "hasn't") with no age given returns the SPARK “not yet” code of 888, so those severe-delay probands are kept and coded as SPARK codes them, rather than dropped. Other entries are left missing and drop at the complete-case step: text with no numeric age ("normal", "on time", "unknown"); a calendar date entered in the age field ("03/2003"); a regression or loss narrative ("12 mos (lost at 15 mos)"); and a bare number left after a bound ("under 2"), whose scale (years or months) is ambiguous.

A parsed age above the SPARK “over 7 years” code is capped at 85 months, matching the dropdown and discarding mis-parsed outliers. The parsing rules and the forms left missing are set out in the package’s milestone-parsing guide.

A bare number carries no unit, so its scale is ambiguous: “4” on bowel training is four years, but “13” on walking is thirteen months. When disambiguate is given, a unit-less number (and a bare number left after a bound, such as “under 2”) is passed to it to resolve that scale; build_milestone_disambiguator builds such a resolver from the SPARK reference distribution. Without disambiguate a bare number is read as months, as the released code did, and a bare number left after a bound stays missing.

Parameters:

• value (object) – A raw milestone cell: a string, a number, or a missing value.

• disambiguate (callable, optional) – A resolver mapping a unit-less age in x to its value in months, choosing between x (months) and 12 * x (years). When omitted, a bare number is taken as months.

Returns:

The age in months, or None when the entry carries no recognisable age.

Return type:

float or None

analysis.cohort.schema.build_milestone_disambiguator(spark_months)

Build a scale resolver for one milestone from its SPARK distribution.

SPARK records each milestone as an age in months on a fixed dropdown grid, so its distribution is a clean, large-sample reference. The SSC records the same milestones as free text without a consistent unit, so a bare number is ambiguous between months and years. The returned resolver decides, for a unit-less number $x$, between $x$ months and $12x$ months by which reading has the higher likelihood under the SPARK distribution in log-age space. This reads a small number as months for early milestones (a child walks at “13” months) and as years for late ones (a child is bowel trained at “4” years), matching how a human reads the field, rather than applying a blanket per-milestone unit rule.

The years reading is only considered when $12x$ stays within the milestone cap; a value whose months reading already exceeds plausible ages keeps the months reading. A degenerate reference (fewer than two distinct ages) yields an identity resolver.

Parameters:

spark_months (numpy.ndarray) – The SPARK ages in months for one milestone. Non-finite values, the “not yet” code, and ages at or above the 85-month cap are dropped before the distribution is estimated.

Returns:

A function mapping a unit-less age x to its resolved value in months.

Return type:

callable

analysis.cohort.schema.load_feature_list(path)

Read the authors’ feature list from its one-column CSV.

Parameters:

path (pathlib.Path) – Location of mixture_model_columns.csv, which has a feature header.

Returns:

The feature names, in file order.

Return type:

list of str

The SPARK cohort backend.

Reproduces the integration in the released process_integrate_phenotype_data.py: the SCQ, background-history (child and sibling), RBS-R, and CBCL 6-18 instruments are read, screened on age and the per-instrument missingness counter, joined on the proband id, and reduced to complete cases. Two things differ from the released code, both deliberate. The kept columns are pinned to the authors’ final feature list rather than rederived by dropping columns (plan section 5), and only the needed columns of each CSV are read, so a whole instrument file is never loaded (the dscat guardrail).

The released datadf.round() is a fit-time transform, applied in analysis.model rather than here, so the cached cohort matrix matches the unrounded intermediate the authors saved.

The backend takes an optional records cutoff (as_of) that restricts the cohort to the probands present at an earlier SPARK freeze, so a later release can be cut back to (an approximation of) the data Litman et al. fit on. See Subsetting the cohort to the V9 freeze for the method and its limits.

class analysis.cohort.spark.SparkCohort(root, version, *, as_of=None)

Build the harmonised SPARK proband-by-feature frame.

Parameters:

• root (pathlib.Path) – Repository root.

• version (str) – SPARK release version (for example "2026-03-23").

• as_of (str, optional) – A records cutoff. When set (for example "2022-12-12", Litman’s V9 freeze), the cohort is restricted to probands registered by the cutoff year whose every cohort instrument was also completed by then. The cutoff is resolved to its calendar year. When None (the default), the full release is built and the behaviour matches the released preprocessing.

supports_timing()

Return True: SPARK carries the diagnosis-timing fields (plan section 5).

axis(name, index, covariates, min_bin_size=1000)

Resolve a stratification axis to its variable and binning policy.

SPARK provides the two SPARK-only timing axes and the two shared cognitive axes:

• age_at_diagnosis and era reuse analysis.strata_data.build_strata_data() and the frozen MaxEqualBins policy (plan section 12a);

• cognitive_impairment is the binary ml_predicted_cog_impair flag (trained against measured IQ below 80), split by the two-band intellectual-disability policy;

• iq is the medical-record full-scale IQ (iq.fsiq_score), equal-frequency binned like the timing axes.

family_ids(index)

Return the family identifier of each proband, for the clustered bootstrap.

SPARK groups probands into families by core_descriptive_variables.family_sf_id (a de-identified family key). The clustered bootstrap (analysis.trajectory_local) resamples families rather than probands, so siblings move together and the tube respects the within-family correlation. Read as a string so the seven-digit zero-padded key keeps its identity; a proband with no family key is its own singleton family (its proband id), so it is never silently pooled with another.

Parameters:

index (pandas.Index) – The modelling-cohort proband index the identifier is aligned to.

Returns:

The per-proband family identifier, on index.

Return type:

pandas.Series

integrate()

Integrate the instruments into the harmonised, complete-case cohort frame.

Returns:

Proband-by-column frame indexed by proband id, holding the covariates and the pinned feature set, coerced to numeric and reduced to complete cases. When a records cutoff is set, the frame is further restricted to the probands present at the cutoff.

Return type:

pandas.DataFrame

The SSC cohort backend.

Harmonises the SSC proband instruments (CBCL 6-18, RBS-R, SCQ-Lifetime, core descriptive, and background history) onto the shared SPARK feature schema, following the renames in the released generate_ssc_data (plan section 10). Rather than reproduce the authors’ column drop-lists, which targeted their SSC version, features are selected positively: a SPARK feature is provided when its SSC column (after renaming) exists. The SSC backend therefore exposes the subset of the schema the SSC instruments cover.

Two caveats are recorded honestly. The authors read the background-history milestones from a hand-cleaned file that was not released, so both the SSC-to-SPARK mapping (SSC_BH_RENAME) and the parse of the raw free-text ages into months (parse_age_months) are ours, and two SPARK milestone features have no SSC equivalent. The SSC milestone ages are free text without a consistent unit, so a bare number is ambiguous between months and years; the scale is resolved per milestone against the SPARK reference distribution (_milestone_priors and build_milestone_disambiguator), so that a small number reads as months on an early milestone and as years on a late one, as a human cleaner reads it. The fidelity of this backend to the authors’ SSC pipeline, and the exact shared-feature contract, are confirmed in the SSC replication stage (phase 2). The backend does not provide diagnosis-timing fields.

class analysis.cohort.ssc.SscCohort(root, version)

Build the harmonised SSC proband-by-feature frame.

Parameters:

• root (pathlib.Path) – Repository root.

• version (str) – SSC dataset version (for example "15.3").

supports_timing()

Return False: SSC does not expose a clean diagnosis timestamp (plan section 5).

family_ids(index)

Return None: the family-clustered bootstrap is wired for the SPARK timing axes only.

The local-trajectory recast (analysis.trajectory_local) runs on the SPARK-only age and era axes, so the SSC backend does not resolve a family key here. The one-family-per proband it would otherwise need is left unimplemented rather than guessed.

axis(name, index, covariates, min_bin_size=1000)

Resolve a stratification axis to its variable and binning policy.

SSC provides the two shared cognitive axes and neither timing axis (it has no clean diagnosis timestamp, plan section 5). Both cognitive axes read the harmonised full-scale deviation IQ (ssc_diagnosis.fs_deviation_score): cognitive_impairment dichotomises it at the intellectual-disability threshold (config.ID_IQ_THRESHOLD, the same construct as SPARK’s flag), and iq bins it by equal frequency, using the continuous score SPARK lacks. age_at_diagnosis and era return None.

integrate()

Integrate the SSC proband instruments into a harmonised, complete-case frame.

Returns:

Proband-by-column frame indexed by individual id, holding the covariates and the subset of the shared feature schema the SSC instruments provide, coerced to numeric and reduced to complete cases.

Return type:

pandas.DataFrame

Feature typing and the model

Feature typing, reverse-coded items, and the category map.

StepMix fits a Gaussian density to each continuous feature, a Bernoulli to each binary feature, and a multinomial to each categorical feature, so the typing has to be right. We derive it three ways and reconcile them: from the data dictionary (via dscat), from the authors’ released typing pickles, and from the observed cardinality in the cohort. The dictionary rebuild is the primary signal and is required to agree with the pickles; where they disagree, the run defers to the pickle typing (the reproduction target) and records the disagreement (plan section 6, step 2).

This module also carries the 24 reverse-coded SCQ social items, which are flipped before the seven-category summary, and loads the feature-to-category map used only for summaries.

class analysis.features.Typing(continuous, binary, categorical)

A partition of the feature set into the three StepMix density types.

continuous, binary, categorical

The features modelled with a Gaussian, a Bernoulli, and a multinomial density respectively.

Type:

list of str

as_dict()

Return a mapping from each feature to its type.

property counts: dict[str, int]

Return the number of features of each type.

analysis.features.n_value_levels(value_coding)

Count the discrete coded levels in a dictionary value coding.

Parameters:

value_coding (str or None) – The value-coding text, where each coded level is a line containing =.

Returns:

The number of coded levels.

Return type:

int

analysis.features.infer_from_dictionary(field_type, value_coding)

Infer a feature type from its dictionary field type and value coding.

Calculated scores and dropdown age codings are continuous. A radio item with exactly two coded levels is binary, otherwise categorical.

Parameters:

• field_type (str or None) – The dictionary field type (for example "radio" or "calculated").

• value_coding (str or None) – The value-coding text.

Returns:

"continuous", "binary", or "categorical".

Return type:

str

analysis.features.load_pickle_typing(typing_dir, features)

Load the released typing for a set of features.

Parameters:

• typing_dir (pathlib.Path) – Directory holding the three typing pickles.

• features (list of str) – Features to type.

Returns:

Each feature mapped to its pickle type (str), or None when it is absent from all pickles or appears in more than one (ambiguous).

Return type:

dict

analysis.features.observed_cardinality(frame, features)

Return the number of distinct non-null values of each feature in the cohort.

Parameters:

• frame (pandas.DataFrame) – The cohort frame.

• features (list of str) – Features to count.

Returns:

Each feature mapped to its distinct-value count.

Return type:

dict of str to int

analysis.features.dictionary_typing(root, dataset, version, features)

Infer each feature’s type from the data dictionary.

Parameters:

• root (pathlib.Path) – Repository root.

• dataset (str) – Dataset and version whose dictionary to read.

• version (str) – Dataset and version whose dictionary to read.

• features (list of str) – Features to type.

Returns:

Each feature mapped to its dictionary-inferred type.

Return type:

dict of str to str

analysis.features.reconcile(features, dict_typing, pickle_typing, observed)

Reconcile the three typing signals into one typing and a report.

The chosen type defers to the pickle typing where it exists (the reproduction target), otherwise to the dictionary. The report records all three signals, whether the dictionary and pickle agree, and whether the observed cardinality is consistent with the chosen type.

Parameters:

• features (list of str) – Features to reconcile.

• dict_typing (dict of str to str) – Dictionary-inferred types.

• pickle_typing (dict) – Released pickle types, each str or None when absent or ambiguous.

• observed (dict of str to int) – Observed distinct-value counts.

Returns:

The reconciled typing and the per-feature reconciliation report.

Return type:

tuple

analysis.features.build_typing(root, dataset, version, features, frame=None)

Build the reconciled feature typing and its report.

Parameters:

• root (pathlib.Path) – Repository root.

• dataset (str) – Dataset and version whose dictionary to read.

• version (str) – Dataset and version whose dictionary to read.

• features (list of str) – The feature set to type.

• frame (pandas.DataFrame, optional) – The cohort frame, used for the observed-cardinality signal. When omitted, that signal is skipped.

Returns:

The reconciled typing and the reconciliation report.

Return type:

tuple

analysis.features.instrument_map(root, dataset, version, features)

Return each feature’s source instrument, read from the data dictionary.

Iterates the cohort instruments in analysis.config.COHORT_INSTRUMENTS order and reads each one’s feature names from the dscat catalogue, the same describe loop dictionary_typing() uses. A feature is assigned to the first instrument that carries it, so the developmental-milestone columns shared by the child and sibling background-history tables resolve to background_history_child (the proband’s own history, listed first), which the sibling table only duplicates by column name. Resolution fails loudly, mirroring reconcile()’s no-typing-signal guard: a feature carried by none of the instruments raises rather than being dropped, so a gap in the mapping cannot pass silently.

Parameters:

• root (pathlib.Path) – Repository root.

• dataset (str) – Dataset and version whose dictionary to read.

• version (str) – Dataset and version whose dictionary to read.

• features (list of str) – The feature set to map.

Returns:

Each feature mapped to its source instrument.

Return type:

dict of str to str

Raises:

ValueError – When a feature is carried by none of the cohort instruments.

analysis.features.load_category_map(path)

Load the feature-to-category map.

Parameters:

path (pathlib.Path) – Location of feature_to_category_mapping.csv, with category and feature columns.

Returns:

Each feature mapped to its category.

Return type:

dict of str to str

The StepMix general finite mixture model wrapper.

A thin layer over StepMix that builds the mixed-data descriptor from a reconciled typing and fits the one-step covariate parametrisation Litman et al. use: a Gaussian, Bernoulli, or multinomial measurement density per feature, with sex and age at evaluation as structural covariates (plan section 6, step 3). The random restarts are delegated to StepMix’s own n_init, as in the released code; StepMix shows the restart progress and its iteration log is captured into the run log.

The released datadf.round() is applied here, immediately before fitting, so the cached cohort matrix stays unrounded while the model sees the rounded values the authors fit on.

class analysis.model.FitResult(model, labels, measurement_data, metrics)

A fitted GFMM with its labels and selection statistics.

model

The fitted estimator.

Type:

StepMix

labels

The hard class label per proband, indexed by proband id.

Type:

pandas.Series

measurement_data

The descriptor-aligned measurement matrix the model was fit and predicted on.

Type:

pandas.DataFrame

metrics

Selection statistics: average log-likelihood, AIC, BIC, sample-size-adjusted BIC, the number of probands, and the class proportions.

Type:

dict

analysis.model.prepare_inputs(matrix, typing, round_values=True)

Build the StepMix measurement descriptor and the covariate matrix.

Parameters:

• matrix (analysis.cohort.CohortMatrix) – The cohort feature and covariate matrices.

• typing (analysis.features.Typing) – The reconciled feature typing.

• round_values (bool, default True) – Round features and covariates to the nearest integer before fitting, as in the released GFMM.py.

Returns:

The measurement data, the mixed-data descriptor, and the covariate matrix.

Return type:

tuple

analysis.model.fit_gfmm(matrix, typing, *, n_components=4, n_init=200, n_steps=1, random_state=None, sample_weight=None, structural='covariate', progress_bar=1, verbose=1)

Fit the one-step covariate GFMM and predict a hard label per proband.

Parameters:

• matrix (analysis.cohort.CohortMatrix) – The cohort feature and covariate matrices.

• typing (analysis.features.Typing) – The reconciled feature typing that sets each feature’s density.

• n_components (int, optional) – Number of latent classes.

• n_init (int, optional) – Number of random restarts, delegated to StepMix.

• n_steps (int, optional) – StepMix estimation steps (one-step joint estimation by default).

• random_state (int or None, optional) – Seed for reproducible restarts.

• sample_weight (numpy.ndarray or pandas.Series or None, optional) – Per-proband weight on the fit, in the row order of the matrix. StepMix’s expectation maximisation weights each proband’s log-likelihood and its class responsibilities by this value, so a weight of zero excludes a proband and a weight of one includes it in full. This is what a LocalisationScheme supplies: an indicator weight reproduces a hard-bin subset fit, a kernel weight gives a local (LSEM) fit. None leaves every proband at weight one, the pooled fit.

• structural (str or None, optional) – The structural (concomitant) model. "covariate" (the default) is Litman’s one-step covariate parametrisation: the class prior depends on sex and age at evaluation, fitted jointly. None fits the measurement model alone (marginal classes, no covariates). The measurement-only fit exists for the kernel sweep, whose fractional weights make the covariate general linear model diverge; it changes the estimand, so it is compared only against a measurement-only reference, never the covariate one.

• progress_bar (int, optional) – StepMix progress-bar verbosity for the restart loop.

• verbose (int, optional) – StepMix log verbosity; its output is captured into the run log.

Returns:

The fitted model, the predicted labels, the measurement data, and the selection statistics.

Return type:

FitResult

analysis.model.selection_metrics(model, measurement_data, covariates, labels)

Compute the information criteria and class proportions for a fit.

Parameters:

• model (StepMix) – The fitted estimator.

• measurement_data (pandas.DataFrame) – The measurement matrix the model was fit on.

• covariates (pandas.DataFrame) – The covariate matrix.

• labels (pandas.Series) – The predicted hard labels.

Returns:

Average log-likelihood, AIC, BIC, sample-size-adjusted BIC, the proband count, and the class proportions sorted by class id.

Return type:

dict

analysis.model.class_centroids(measurement_data, labels)

Return the per-class feature means (class-by-feature centroid matrix).

Parameters:

• measurement_data (pandas.DataFrame) – The measurement matrix.

• labels (pandas.Series) – The predicted hard labels on the same index.

Returns:

Class-by-feature mean matrix, indexed by class id.

Return type:

pandas.DataFrame

Enrichment and alignment

Per-class feature enrichment and the seven-category signature.

Reproduces the released enrichment pipeline (plan section 6, step 7). Each feature is tested one class against the rest, in both directions: a binomial test for binary features (those with two observed values) and a Welch t-test for the rest. The \(p\)-values are Benjamini-Hochberg corrected within each class and direction; a corrected \(p\) below \(0.05\) marks a feature as enriched or depleted in that class. The 24 reverse-coded SCQ items have their direction flipped, and the features are summarised into the seven literature-defined categories as the signed proportion enriched minus depleted, which is the class signature used to align to the published classes (plan section 6a).

class analysis.enrich.CorrelationInterval

A percentile bootstrap interval on a profile correlation.

class analysis.enrich.BootstrapInterval

The overall interval, plus the same interval for each of the seven categories.

analysis.enrich.cohens_d(group, reference)

Return Cohen’s \(d\) of a group against a reference, pooling their variances.

analysis.enrich.feature_enrichment(data, labels, n_classes=4)

Test each feature for one-versus-rest enrichment in every class.

Parameters:

• data (pandas.DataFrame) – The proband-by-feature measurement matrix used for the fit.

• labels (pandas.Series) – The hard class label per proband, on the same index.

• n_classes (int, default 4) – Number of classes.

Returns:

One row per feature with, per class c: class{c}_dir (+1 enriched, -1 depleted, 0 neither, after correction) and class{c}_effect (fold enrichment for binary features, Cohen’s \(d\) otherwise), plus an is_binary flag.

Return type:

pandas.DataFrame

analysis.enrich.profile_correlation(signature_a, signature_b)

Correlate two seven-category class signatures, overall and per category.

This is the comparison Litman et al. use to declare reproduction and replication: the profile is the signed net-proportion-enriched vector per category, and the correlation is taken over the class-by-category matrix (plan section 6a). The two signatures must already have their classes aligned to a common ordering (for example by analysis.align.greedy_overlap_align() for same-sample comparison or analysis.align.hungarian_align() across cohorts).

Parameters:

• signature_a (pandas.DataFrame) – Class-by-category signed signatures on the same class index and the seven categories in SEVEN_CATEGORIES order.

• signature_b (pandas.DataFrame) – Class-by-category signed signatures on the same class index and the seven categories in SEVEN_CATEGORIES order.

Returns:

The overall Pearson correlation over the flattened class-by-category matrix, and a mapping from each category to its Pearson correlation across classes. A near-constant profile (overall or within a category) has an undefined correlation, returned as None.

Return type:

tuple

Raises:

ValueError – When the two signatures do not share their class index.

analysis.enrich.bootstrap_overall_correlation(measurement, labels, target, category_map, *, n_boot, seed, n_classes=4, level=0.95, reverse_coded=('q02_conversation', 'q09_expressions_appropriate', 'q19_best_friend', 'q20_talk_friendly', 'q21_copy_you', 'q22_point_things', 'q23_gestures_wanted', 'q24_nod_head', 'q25_shake_head', 'q26_look_directly', 'q27_smile_back', 'q28_things_interested', 'q29_share', 'q30_join_enjoyment', 'q31_comfort', 'q32_help_attention', 'q33_range_expressions', 'q34_copy_actions', 'q35_make_believe', 'q36_same_age', 'q37_respond_positively', 'q38_pay_attention', 'q39_imaginative_games', 'q40_cooperatively_games'), keep=None)

Bootstrap the overall profile correlation by resampling probands.

The class labels are held fixed and the probands are resampled with replacement; for each resample the seven-category signature is recomputed and correlated against a fixed target on the same class index. The spread of those correlations is the sampling uncertainty in the reproduction or replication statistic from the finite cohort, with the model fit held fixed. It does not capture uncertainty from refitting the model (the stability stage does that) nor, for the reproduction, the resolution of the figure-read target.

Parameters:

• measurement (pandas.DataFrame) – The proband-by-feature matrix.

• labels (pandas.Series) – The hard class label per proband, positionally aligned with measurement.

• target (pandas.DataFrame) – The fixed signature to correlate against, on the same class index the recomputed signature carries (class ids 0 to n_classes - 1).

• category_map (dict of str to str) – Feature-to-category map for the signatures.

• n_boot (int) – Number of bootstrap resamples.

• seed (int) – Seed for the resampling, so the interval is reproducible.

• n_classes (int, default 4) – Number of classes.

• level (float, default 0.95) – Central probability of the reported percentile interval.

• reverse_coded (tuple of str, optional) – SCQ items whose enrichment direction is flipped before the signature.

• keep (set of str, optional) – The contributory feature set, applied to every resample.

Returns:

The overall interval as ci_low and ci_high (the percentile interval at level), median, the level, and n_valid (the resamples that yielded a defined overall correlation; a resample that empties or flattens a class is dropped, as in the permutation null). A category entry then holds the same interval for each of the seven categories, keyed by category name. A per-category correlation is taken over the four classes alone, so its interval is markedly wider than the overall one and rests on fewer valid resamples.

Return type:

dict

analysis.enrich.contributory_features(enrichment, n_classes=4)

Return the features that contribute to the class signatures.

Reproduces the released non-contributory feature exclusion (plan section 6, step 7). A feature is dropped when it is significantly enriched or depleted in no class, or when its effect size stays below the magnitude threshold in every class: a fold enrichment below \(1.5\) for binary features, an absolute Cohen’s \(d\) below \(0.2\) for the rest. The surviving features are the universe over which the seven-category proportions are computed, so the same set (taken from the reference solution) is applied to both sides of a comparison.

Parameters:

• enrichment (pandas.DataFrame) – The per-feature enrichment from feature_enrichment().

• n_classes (int, default 4) – Number of classes.

Returns:

The contributory feature names, in the enrichment frame’s order.

Return type:

list of str

analysis.enrich.category_signature(enrichment, category_map, n_classes=4, reverse_coded=('q02_conversation', 'q09_expressions_appropriate', 'q19_best_friend', 'q20_talk_friendly', 'q21_copy_you', 'q22_point_things', 'q23_gestures_wanted', 'q24_nod_head', 'q25_shake_head', 'q26_look_directly', 'q27_smile_back', 'q28_things_interested', 'q29_share', 'q30_join_enjoyment', 'q31_comfort', 'q32_help_attention', 'q33_range_expressions', 'q34_copy_actions', 'q35_make_believe', 'q36_same_age', 'q37_respond_positively', 'q38_pay_attention', 'q39_imaginative_games', 'q40_cooperatively_games'), keep=None)

Summarise feature enrichment into the signed seven-category class signature.

For each class and category the signature is the proportion of the category’s features enriched minus the proportion depleted, after flipping the reverse-coded SCQ items.

Parameters:

• enrichment (pandas.DataFrame) – The per-feature enrichment from feature_enrichment().

• category_map (dict of str to str) – Feature-to-category map.

• n_classes (int, default 4) – Number of classes.

• reverse_coded (tuple of str, optional) – SCQ items whose enriched and depleted directions are swapped.

• keep (set of str, optional) – When given, only these features contribute to the proportions (the contributory set from contributory_features()). When omitted, every feature contributes, which is the right behaviour for the reference signature compared to the published figure.

Returns:

Class-by-category signed signature, indexed by class id, columns the seven categories in reporting order.

Return type:

pandas.DataFrame

Hungarian alignment of one set of classes to another by profile similarity.

StepMix assigns class ids arbitrarily on every fit, so a recovered class has no fixed meaning. To compare two solutions (our fit against the published classes, or a stratum against the reference) we align them by matching their profiles with the Hungarian algorithm on a cost matrix of profile distances (Kuhn 1955), as the plan specifies for cross-stratum and cross-cohort comparison (plan section 6, deviations).

class analysis.align.Alignment(mapping, cost, correlations, total_cost)

The result of aligning a source set of classes to a target set.

mapping

Each source class id mapped to the target label it aligns to.

Type:

dict

cost

The source-by-target cost matrix the assignment minimised.

Type:

pandas.DataFrame

correlations

Each source class id mapped to the Pearson correlation of its profile with its assigned target profile.

Type:

dict

total_cost

The minimised total assignment cost.

Type:

float

analysis.align.hungarian_align(source, target, metric='correlation')

Align source classes to target classes by minimising profile distance.

Parameters:

• source (pandas.DataFrame) – Source class-by-profile matrix (for example our four classes by seven categories).

• target (pandas.DataFrame) – Target class-by-profile matrix with the same columns (for example the published named classes).

• metric (str, default "correlation") – "correlation" (cost is one minus Pearson r) or "euclidean".

Returns:

The mapping, the cost matrix, the per-pair correlations, and the total cost.

Return type:

Alignment

Raises:

ValueError – When the source and target columns differ.

analysis.align.greedy_overlap_align(source, target)

Align source class labels to target labels by greedy proband overlap.

This reproduces the released match_class_labels rule, which compares two clusterings of the same probands (across seeds or subsamples). Each source class claims the target class it overlaps most, where overlap is the proportion of the source class that falls in the target class; a collision is resolved in favour of the larger overlap, and any classes left unclaimed are paired in order. The rule needs a shared index, so it applies to same-sample comparison only; for disjoint strata and across cohorts the plan aligns by profile similarity with hungarian_align() instead (plan section 6, deviations).

Parameters:

• source (pandas.Series) – Class label per proband for the solution being aligned.

• target (pandas.Series) – Class label per proband for the reference solution, on an overlapping index.

Returns:

Each source class id mapped to the target class id it aligns to.

Return type:

dict of int to int

Aligning our recovered classes to Litman’s four named classes.

Two of the three alignment routes the authors use are closed to us: we lack SPARK v9 and their per-proband labels, and no reference model was released. We therefore align on the named-class anchors: the substantive, data-driven characteristics the authors use to define each class (the most-developmental class is Mixed ASD with DD; the highest-difficulty and smallest class is Broadly affected; the largest, high-core, no developmental delay class is Social/behavioral; the uniformly lowest is Moderate challenges). The assignment is cross-validated for mutual consistency.

The published seven-category signatures (read from figure 1b) give a profile correlation against each named class and an overall correlation, the analogue of the authors’ own replication measure (their SSC replication reported \(r = 0.927\)). Those values are read to the figure’s resolution, not from a supplementary table (plan section 6a, step 1); the values themselves are tabulated in the reproduction guide.

class analysis.reference.NamedAlignment(mapping, correlations, overall_correlation, anchors, anchors_hold, cost)

The alignment of our classes to the named classes and its validation.

mapping

Each of our class ids mapped to a named class.

Type:

dict

correlations

Per-class Pearson correlation of our signature with the assigned named published signature. A class is None when its published profile is saturated (constant), which makes the correlation undefined.

Type:

dict

overall_correlation

Pearson correlation over the full class-by-category matrix (our classes aligned to the named published classes), the analogue of the authors’ replication measure.

Type:

float

anchors

Each named anchor mapped to whether it held for the assigned class.

Type:

dict

anchors_hold

Whether every anchor held.

Type:

bool

cost

The Hungarian cost matrix (our classes by named classes).

Type:

pandas.DataFrame

analysis.reference.published_signature()

Return the published seven-category signature (figure 1b), classes by category.

analysis.reference.align_to_named(signature, proportions)

Align our recovered classes to the named classes and validate the mapping.

Parameters:

• signature (pandas.DataFrame) – Our class-by-category signature, indexed by class id.

• proportions (dict of int to float) – Our class proportions by class id.

Returns:

The anchor mapping, the per-class and overall profile correlations against the published signature, the anchor consistency checks, and the cost matrix from the (secondary) Hungarian route.

Return type:

NamedAlignment

Selection, stability, and replication

Model selection: the number-of-components grid and its information criteria.

Reproduces the released GFMM_model_validation procedure (plan section 6, step 4). For each candidate number of classes, the one-step covariate GFMM is fitted and scored on a panel of criteria, repeated over several seeds and summarised as a mean and a standard deviation. The criteria are the validation log-likelihood from 3-fold cross-validation, the Akaike, Bayesian, sample-size-adjusted Bayesian, and consistent Akaike information criteria, the approximate weight of evidence, the scaled relative entropy, the average latent-class posterior probability (ALCPP), and the smallest class proportion.

Two faithfulness notes. StepMix 3.0.0 exposes aic, bic, sabic, and caic as methods whose formulas match the authors’ hand-rolled helpers, so those are used directly; the approximate weight of evidence is not a StepMix method and is computed here. The “Lo-Mendell-Rubin likelihood-ratio test” in the released code is a naive chi-square on the cross-validated log-likelihood differences with the degrees of freedom fixed at one, not the analytically correct adjusted test; lmr_lrt_proxy() reproduces that approximation and is documented as such. Litman et al. do not seed the per-iteration fits; we seed them for reproducibility (plan section 11), which is the only deliberate divergence here.

The released decision for four components is asserted visually (a reference line at four on every panel) rather than by an automatic rule, so this module reports the full criteria table and leaves the choice to the methods write-up.

class analysis.selection.SelectionResult(per_iteration, summary, k_values)

The model-selection grid and its summary.

per_iteration

One row per (seed, number of components) with every criterion.

Type:

pandas.DataFrame

summary

Mean and standard deviation of each criterion per number of components, plus the mean Lo-Mendell-Rubin proxy \(p\)-value.

Type:

pandas.DataFrame

k_values

The component counts gridded.

Type:

list of int

analysis.selection.awe(model, measurement, covariates)

Return the approximate weight of evidence for a fitted model.

The penalty sits between the consistent Akaike criterion and a heavier term: \(-2 \ell + k(\ln n + 1.5)\), where \(\ell\) is the total log-likelihood, \(k\) the number of free parameters, and \(n\) the sample size. StepMix does not expose this criterion, so it is computed from score (the average log-likelihood) and n_parameters, as in the authors’ utils.awe.

Parameters:

• model (StepMix) – A fitted estimator.

• measurement (pandas.DataFrame) – The measurement matrix the model was scored on.

• covariates (pandas.DataFrame) – The structural covariate matrix.

Returns:

The approximate weight of evidence (lower is better).

Return type:

float

analysis.selection.per_fit_criteria(model, measurement, covariates)

Compute the per-fit selection criteria for one fitted model.

The information criteria and the average log-likelihood are scored with the covariate channel (as in the released validation code); the hard labels and posteriors used for the ALCPP and the smallest-class proportion come from the measurement posterior without covariates, matching how the reference labels are assigned (plan section 6, step 3).

Parameters:

• model (StepMix) – A fitted estimator.

• measurement (pandas.DataFrame) – The measurement matrix.

• covariates (pandas.DataFrame) – The structural covariate matrix.

Returns:

One value per entry in PER_FIT_CRITERIA.

Return type:

dict of str to float

analysis.selection.validation_log_likelihood(measurement, covariates, descriptor, k_values, *, seed, n_init, cv=3)

Cross-validated mean log-likelihood per number of components.

Uses scikit-learn’s grid search over n_components with the StepMix average log-likelihood as the score, as in the released compute_LL (plan section 6, step 4).

Parameters:

• measurement (pandas.DataFrame) – The measurement matrix.

• covariates (pandas.DataFrame) – The structural covariate matrix, passed as the supervised target so it enters the cross-validated score.

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• k_values (collections.abc.Sequence of int) – The component counts to grid over.

• seed (int) – Random seed for the estimator.

• n_init (int) – Random restarts per fold fit.

• cv (int, default 3) – Number of cross-validation folds.

Returns:

Each number of components mapped to its mean validation log-likelihood.

Return type:

dict of int to float

analysis.selection.lmr_lrt_proxy(val_log_likelihood)

Naive likelihood-ratio-test \(p\)-values between adjacent component counts.

For consecutive counts \(k\) and \(k+1\) the statistic is \(-2(\ell_k - \ell_{k+1})\) referred to a chi-square with one degree of freedom, keyed at \(k+1\). The log-likelihoods are the cross-validated per-sample means. This is the approximation the released code computes; it is not the analytically correct Lo-Mendell-Rubin adjusted test (the degrees of freedom are fixed at one rather than the difference in free parameters), and it is reported as a proxy only.

Parameters:

val_log_likelihood (dict of int to float) – Validation log-likelihood per number of components.

Returns:

Each number of components (from the second upward) mapped to its proxy \(p\)-value.

Return type:

dict of int to float

analysis.selection.run_selection(matrix, typing, *, k_values, n_iterations, n_init, base_seed=0, cv=3, checkpoint_dir=None, workers=0)

Run the component-count grid over several seeds and summarise the criteria.

Iterations are independent (each is its own cross-validation pass plus a fit per K), so they run concurrently over a ProcessPoolExecutor, each pinned to a single BLAS thread (analysis.profiling.single_threaded_blas()) so concurrent workers do not oversubscribe the machine. Each seeded iteration is one resumable unit: when checkpoint_dir is given, an iteration’s criterion rows are appended to a checkpoint as it completes, and a re-run over the same directory recomputes only the iterations missing from it (tracked by iteration index, not by checkpoint length, since concurrent iterations do not finish in submission order). Because iteration i always uses seed base_seed + i, a resumed run reproduces what an uninterrupted run would have computed; the per-iteration rows are reassembled in iteration order regardless of completion order, so a resumed run’s table matches an uninterrupted one exactly.

A single-initialisation fit can fail to converge when the structural covariate M-step hits a near-singular design, most often at higher class counts. Such a fit has its criteria recorded as nan rather than raising, so one bad fit does not lose the whole grid; the per-component summary skips those entries. This matches the cross-validation pass, where scikit-learn already scores a failed fold as nan.

Parameters:

• matrix (analysis.cohort.CohortMatrix) – The cohort feature and covariate matrices.

• typing (analysis.features.Typing) – The reconciled feature typing.

• k_values (collections.abc.Sequence of int) – The component counts to grid over.

• n_iterations (int) – Number of seeded repetitions (the released code uses 200).

• n_init (int) – Random restarts per fit (the released validation fits use one).

• base_seed (int, default 0) – Seeds are base_seed to base_seed + n_iterations - 1.

• cv (int, default 3) – Cross-validation folds for the validation log-likelihood.

• checkpoint_dir (pathlib.Path, optional) – Directory for the resumable checkpoint. When None the run is held in memory and nothing is written, so an interrupt loses the work. The directory must be specific to these parameters (the caller passes the stage’s content-addressed run directory).

• workers (int, default 0) – Concurrent worker processes. 0 uses the logical core count minus one. 1 runs in-process instead of through a pool, which a test that monkeypatches a fitting dependency relies on (a spawned worker process would not see the patch).

Returns:

The per-iteration criteria, the per-component summary, and the gridded counts.

Return type:

SelectionResult

Multi-initialisation and subsampling stability, and the minimum viable stratum size.

Reproduces the released stability and subsampling analyses (plan section 6, step 5) and reuses the same machinery to fix the minimum stratum size for the stratified work (plan section 7b).

Multi-initialisation stability runs many single-initialisation fits from different random starts, ranks them by log-likelihood, and compares the best of them to the reference solution. Subsampling stability refits on random halves of the cohort and compares each back to the reference. The comparison is threefold: the seven-category profile correlation (the authors’ own measure), the class-overlap matrix, and the adjusted Rand index, which the plan adds to the released overlap because it is label-invariant and chance-corrected (plan section 6, deviations). Same-sample comparisons align class labels with the released greedy overlap rule (analysis.align.greedy_overlap_align()); the Rand index needs no alignment.

The released code runs 2,000 single-init fits and reports the best 100, and refits on 100 halves; both counts are configurable here. Litman et al. do not seed these fits; we seed them for reproducibility (plan section 11), the only deliberate divergence.

The minimum viable stratum size is found by refitting at descending sample sizes and recording where four-class recovery degrades: the smallest class proportion, the scaled relative entropy, the average latent-class posterior probability, and the profile correlation to the full-sample reference. The size below which the profile correlation falls past the reproduction benchmark is the floor for the stratification bins.

A fit can fail to converge when the structural covariate M-step meets a near-singular design, most often at higher class counts on a small subsample. Rather than abort the stage, a failed fit is recorded as missing (_try_fit): the multi-initialisation run drops it before ranking, and the subsampling and minimum-size sweeps mark its replicate degenerate, so it falls out of the aggregate means.

class analysis.stability.Comparison(overall_correlation, category_correlation, adjusted_rand_index, smallest_class_proportion, overlap, degenerate)

One fit’s comparison to the reference solution.

overall_correlation

Pearson correlation over the flattened class-by-category profiles, or None when a profile is near-constant.

Type:

float or None

category_correlation

Per-category Pearson correlation across classes.

Type:

dict of str to float or None

adjusted_rand_index

Chance-corrected agreement between the fit and reference labellings.

Type:

float

smallest_class_proportion

The smallest class proportion in the aligned fit, counting a collapsed class as zero (so a four-to-three collapse reads as a proportion near zero, not the smallest surviving class).

Type:

float

overlap

The source-by-reference class-overlap matrix.

Type:

pandas.DataFrame

degenerate

Whether the fit collapsed a class (recovered fewer than n_components classes), in which case the profile correlation is undefined and the fit is dropped from the aggregate means, as in the released code.

Type:

bool

class analysis.stability.StabilitySummary(fits, comparisons, overlap_mean, aggregate=<factory>)

The result of a multi-initialisation or subsampling stability run.

fits

One row per fit with its seed, average log-likelihood, and convergence flag.

Type:

pandas.DataFrame

comparisons

One row per compared fit with the overall and per-category profile correlations, the adjusted Rand index, and the smallest class proportion.

Type:

pandas.DataFrame

overlap_mean

The mean class-overlap matrix over the compared fits (source class on rows, reference class on columns).

Type:

pandas.DataFrame

aggregate

Mean and standard deviation of the overall correlation and the adjusted Rand index, and the mean per-category correlation.

Type:

dict

class analysis.stability.NminResult(per_fit, summary, n_min, benchmark, floor, floor_ci)

The minimum-viable-stratum-size sweep.

per_fit

One row per (target size, replicate) with the recovery metrics.

Type:

pandas.DataFrame

summary

Mean recovery metrics per target size.

Type:

pandas.DataFrame

n_min

The smallest target size whose mean profile correlation holds at or above the benchmark, or None when no swept size clears it. Kept for continuity; it reads one clearing size off a possibly non-monotone curve, so prefer floor below.

Type:

int or None

benchmark

The profile-correlation threshold used.

Type:

float

floor

The recovery floor from a monotone (isotonic) fit of correlation against log-size: the smallest size at which the fitted recovery reaches the benchmark. None when the fitted curve does not reach the benchmark anywhere in the swept range. Pools every fit, so it is robust to the scatter a small replicate count produces.

Type:

int or None

floor_ci

The 90 per cent bootstrap confidence interval (lower, upper) for floor, or None when too few resamples cross to form one. The upper bound is the conservative bin floor.

Type:

tuple of int or None

analysis.stability.class_overlap_matrix(source, target, n_components)

Return the class-overlap matrix between two labellings of shared probands.

Cell (k, j) is the proportion of reference class j whose probands fall in source class k. After the source labels are aligned to the reference, the diagonal is the retention of each class.

Parameters:

• source (pandas.Series) – Source class label per proband (aligned to the reference class ids).

• target (pandas.Series) – Reference class label per proband, on the same index.

• n_components (int) – Number of classes.

Returns:

Source-by-reference overlap proportions.

Return type:

pandas.DataFrame

analysis.stability.compare_to_reference(measurement, fit_labels, reference_labels, reference_enrichment, category_map, *, n_components, reverse_coded=('q02_conversation', 'q09_expressions_appropriate', 'q19_best_friend', 'q20_talk_friendly', 'q21_copy_you', 'q22_point_things', 'q23_gestures_wanted', 'q24_nod_head', 'q25_shake_head', 'q26_look_directly', 'q27_smile_back', 'q28_things_interested', 'q29_share', 'q30_join_enjoyment', 'q31_comfort', 'q32_help_attention', 'q33_range_expressions', 'q34_copy_actions', 'q35_make_believe', 'q36_same_age', 'q37_respond_positively', 'q38_pay_attention', 'q39_imaginative_games', 'q40_cooperatively_games'))

Compare one fit’s labelling to the reference on shared probands.

The fit labels are aligned to the reference by greedy overlap, the seven-category signature is recomputed on the aligned labels, and three statistics are returned: the overall and per-category profile correlations against the reference signature, the adjusted Rand index (on the raw labels, which is label-invariant), and the class-overlap matrix. The contributory feature set is taken from the reference enrichment and applied to both signatures, so the correlation is computed over the same feature universe the authors used (plan section 6, step 7).

Parameters:

• measurement (pandas.DataFrame) – The measurement matrix the fit was made on (full cohort or subsample).

• fit_labels (pandas.Series) – The fit’s hard labels.

• reference_labels (pandas.Series) – The reference solution’s hard labels.

• reference_enrichment (pandas.DataFrame) – The reference solution’s per-feature enrichment, from which the reference signature and the contributory feature set are derived.

• category_map (dict of str to str) – Feature-to-category map for the signature.

• n_components (int) – Number of classes.

• reverse_coded (tuple of str, optional) – SCQ items whose enrichment direction is flipped before the signature.

Returns:

The overall and per-category profile correlations, the adjusted Rand index, the smallest class proportion, the class-overlap matrix, and a degenerate-fit flag.

Return type:

Comparison

analysis.stability.run_multi_init_stability(matrix, typing, reference_labels, reference_enrichment, category_map, *, n_fits, top_k, n_components=4, base_seed=0, checkpoint_dir=None, workers=0)

Run many single-init fits, rank by log-likelihood, and compare the best to the reference.

Both phases parallelise over independent single-init fits, each pinned to a single BLAS thread (analysis.profiling.single_threaded_blas()) so concurrent workers do not oversubscribe the machine. The run has two resumable phases when checkpoint_dir is given. Each fit appends its seed and log-likelihood to a checkpoint as it completes, and each top-top_k comparison appends its result to a second checkpoint; a re-run over the same directory recomputes only the fits and comparisons missing from them (tracked by seed, not by checkpoint length, since concurrent units do not finish in submission order). The per-fit labels are not stored: a comparison whose fit was restored from a prior run refits that seed on demand (the same seed and single initialisation reproduce it), which keeps the checkpoint to scalars while still resuming exactly. Both phases reassemble their rows in a deterministic order regardless of completion order (fits by seed before the final rank, comparisons in best_seeds order), so a resumed run’s tables match an uninterrupted run’s exactly.

A single-initialisation fit that fails to converge (see _try_fit) is kept in the ranked table with a nan log-likelihood and dropped before the top-top_k selection, so one bad fit never crashes the run and only well-formed fits are compared.

Parameters:

• matrix (analysis.cohort.CohortMatrix) – The cohort feature and covariate matrices.

• typing (analysis.features.Typing) – The reconciled feature typing.

• reference_labels (pandas.Series) – The reference solution’s labels.

• reference_enrichment (pandas.DataFrame) – The reference solution’s per-feature enrichment (for the signature and the contributory feature set).

• category_map (dict of str to str) – Feature-to-category map.

• n_fits (int) – Number of single-initialisation fits (the released code uses 2,000).

• top_k (int) – Number of best fits (by log-likelihood) to compare to the reference (released: 100).

• n_components (int, optional) – Number of classes.

• base_seed (int, optional) – Seeds are base_seed to base_seed + n_fits - 1.

• checkpoint_dir (pathlib.Path, optional) – Directory for the resumable checkpoints. When None the run is held in memory and an interrupt loses the work. The directory must be specific to these parameters.

• workers (int, default 0) – Concurrent worker processes per phase. 0 uses the logical core count minus one. 1 runs in-process instead of through a pool, which a test that monkeypatches a fitting dependency relies on (a spawned worker process would not see the patch).

Returns:

The ranked fits, the per-fit comparisons of the top top_k, the mean overlap, and the aggregate statistics.

Return type:

StabilitySummary

analysis.stability.run_subsampling_stability(matrix, typing, reference_labels, reference_enrichment, category_map, *, n_reps, frac=0.5, n_init=20, n_components=4, base_seed=0, checkpoint_dir=None, workers=0)

Refit on random subsamples and compare each back to the reference.

Replicates are independent, so they run concurrently over a ProcessPoolExecutor, each pinned to a single BLAS thread (analysis.profiling.single_threaded_blas()) so concurrent workers do not oversubscribe the machine. Each replicate is one resumable unit: a re-run over the same checkpoint_dir recomputes only the replicates missing from it (tracked by replicate index, not by checkpoint length, since concurrent replicates do not finish in submission order), and the results are reassembled in replicate order regardless of completion order, so a resumed run’s tables match an uninterrupted run’s exactly.

Parameters:

• matrix (analysis.cohort.CohortMatrix) – The cohort feature and covariate matrices.

• typing (analysis.features.Typing) – The reconciled feature typing.

• reference_labels (pandas.Series) – The reference solution’s labels.

• reference_enrichment (pandas.DataFrame) – The reference solution’s per-feature enrichment (for the signature and the contributory feature set).

• category_map (dict of str to str) – Feature-to-category map.

• n_reps (int) – Number of subsample replicates (the released code uses 100).

• frac (float, default 0.5) – Subsample fraction without replacement.

• n_init (int, default 20) – Random restarts per subsample fit (the released code uses 20).

• n_components (int, optional) – Number of classes.

• base_seed (int, optional) – Seeds are base_seed to base_seed + n_reps - 1.

• checkpoint_dir (pathlib.Path, optional) – Directory for the resumable checkpoint. When None the run is held in memory and an interrupt loses the work. The directory must be specific to these parameters.

• workers (int, default 0) – Concurrent worker processes. 0 uses the logical core count minus one. 1 runs in-process instead of through a pool.

Returns:

The per-replicate fits, comparisons, mean overlap, and aggregate statistics.

Return type:

StabilitySummary

analysis.stability.estimate_floor(per_fit, benchmark, *, n_bootstrap=1000, seed=0)

Estimate the recovery floor by isotonic regression with a bootstrap interval.

Fits a monotone (non-decreasing) regression of the per-fit profile correlation on \(\log_{10}\) size, then reads the smallest size at which the fitted recovery reaches benchmark. Because recovery improves with sample size in expectation, the monotone fit irons out the scatter a small replicate count produces, so the estimate is stable where the smallest-clearing-size rule is not. A fit-level bootstrap gives a percentile interval; its upper bound is the conservative bin floor.

Parameters:

• per_fit (pandas.DataFrame) – One row per fit, with size and overall_correlation. Rows whose correlation is missing (a fit that collapsed a class) are dropped.

• benchmark (float) – The profile-correlation threshold that defines recovery.

• n_bootstrap (int, default 1000) – Fit-level bootstrap resamples for the interval.

• seed (int, default 0) – Seed for the bootstrap resampling.

Returns:

(floor, (lower, upper)). floor is the crossing size, or None when the fitted curve does not reach the benchmark in the swept range. The interval is None when fewer than half the resamples cross.

Return type:

tuple

analysis.stability.run_nmin_sweep(matrix, typing, reference_enrichment, reference_labels, category_map, *, sizes, n_reps, benchmark, n_init=20, n_components=4, base_seed=0, checkpoint_dir=None, workers=0)

Refit at descending sample sizes to fix the minimum viable stratum size.

Each target size is fitted n_reps times on a random subsample of that size; the recovery metrics recorded are the smallest class proportion, the scaled relative entropy, the average latent-class posterior probability, and the profile correlation to the full-sample reference. The minimum viable size is the smallest swept size whose mean profile correlation holds at or above benchmark (plan section 7b).

Every (size, replicate) cell is an independent fit, so the sweep runs concurrently over a ProcessPoolExecutor, each worker pinned to a single BLAS thread (analysis.profiling.single_threaded_blas()) so concurrent workers do not oversubscribe the machine. Each cell is one resumable unit: a re-run over the same checkpoint_dir recomputes only the cells missing from it (tracked by (size, replicate), not by checkpoint length, since concurrent cells do not finish in submission order), and the results are reassembled in the sweep’s own grid order regardless of completion order, so a resumed run’s tables match an uninterrupted run’s exactly.

Parameters:

• matrix (analysis.cohort.CohortMatrix) – The cohort feature and covariate matrices.

• typing (analysis.features.Typing) – The reconciled feature typing.

• reference_enrichment (pandas.DataFrame) – The full-sample reference per-feature enrichment (for the signature and the contributory feature set).

• reference_labels (pandas.Series) – The reference solution’s labels (for the profile-correlation alignment).

• category_map (dict of str to str) – Feature-to-category map.

• sizes (collections.abc.Sequence of int) – Target subsample sizes to sweep, largest first.

• n_reps (int) – Replicates per size.

• benchmark (float) – Profile-correlation threshold that defines recovery.

• n_init (int, default 20) – Random restarts per fit.

• n_components (int, optional) – Number of classes.

• base_seed (int, optional) – Base seed; each (size, replicate) gets a distinct derived seed.

• checkpoint_dir (pathlib.Path, optional) – Directory for the resumable checkpoint. When None the run is held in memory and an interrupt loses the work. The directory must be specific to these parameters.

• workers (int, default 0) – Concurrent worker processes. 0 uses the logical core count minus one. 1 runs in-process instead of through a pool, which a test that monkeypatches a fitting dependency relies on (a spawned worker process would not see the patch).

Returns:

The per-fit metrics, the per-size summary, the smallest-clearing-size n_min, and the isotonic recovery floor with its bootstrap interval.

Return type:

NminResult

Cross-cohort replication: train on SPARK, project onto the SSC.

Reproduces the released replication_on_SSC procedure (plan section 6, step 6). The shared feature set is the intersection of the harmonised SPARK and SSC matrices. A fresh GFMM is fitted on SPARK restricted to those shared features, then the fitted model predicts class labels on the SSC. Because both cohorts pass through the one model, the class ids already correspond, so no cross-cohort label alignment is needed (unlike the disjoint-strata case in phase 4). The replication measure is the correlation of the seven-category enrichment profiles between the two cohorts, the same currency Litman et al. use to declare replication.

The plan adds a calibration the released code lacks: a permutation null that breaks the SSC class-to-profile association, so the observed correlation is read against chance rather than asserted (plan section 6, deviations; section 12).

Two faithfulness points. StepMix validates prediction inputs by feature count, not by name, so the SSC measurement matrix is reindexed to the exact SPARK column order before prediction. The SSC harmonisation relies on our own milestone mapping (the authors used a hand-cleaned background-history file that was not released), and the locally held SSC release is small, so the replication is reported with its sample size and these caveats rather than as a clean reproduction of the published value.

class analysis.replicate.ReplicationResult(shared_features, spark_signature, ssc_signature, overall_correlation, category_correlation, null_overall, p_value, correlation_ci, metrics)

The cross-cohort replication and its calibration.

shared_features

The features shared by the two cohorts and used for the fit.

Type:

list of str

spark_signature, ssc_signature

The seven-category signatures of the SPARK fit and the SSC projection.

Type:

pandas.DataFrame

overall_correlation

Pearson correlation over the flattened class-by-category profiles.

Type:

float or None

category_correlation

Per-category Pearson correlation across classes.

Type:

dict

null_overall

The overall correlation under each label permutation.

Type:

list of float

p_value

The proportion of null correlations at or above the observed one, or None when the observed correlation is undefined or no permutations were run.

Type:

float or None

correlation_ci

The proband-bootstrap percentile interval on the overall correlation, with a category block holding the same interval per category, or None when the bootstrap was skipped.

Type:

dict or None

metrics

Sample sizes, class proportions, and the convergence flag.

Type:

dict

analysis.replicate.shared_feature_set(spark, ssc)

Return the features shared by both cohorts, in SPARK feature order.

A stable order (SPARK order) is used rather than a set’s arbitrary order, so the fit and projection see identical column positions (StepMix checks feature count, not name).

Parameters:

• spark (analysis.cohort.CohortMatrix) – The harmonised cohort matrices.

• ssc (analysis.cohort.CohortMatrix) – The harmonised cohort matrices.

Returns:

The shared feature names.

Return type:

list of str

analysis.replicate.run_replication(spark, ssc, typing, category_map, *, n_components=4, n_init=200, n_permutations=200, n_bootstrap=500, seed=0, reverse_coded=('q02_conversation', 'q09_expressions_appropriate', 'q19_best_friend', 'q20_talk_friendly', 'q21_copy_you', 'q22_point_things', 'q23_gestures_wanted', 'q24_nod_head', 'q25_shake_head', 'q26_look_directly', 'q27_smile_back', 'q28_things_interested', 'q29_share', 'q30_join_enjoyment', 'q31_comfort', 'q32_help_attention', 'q33_range_expressions', 'q34_copy_actions', 'q35_make_believe', 'q36_same_age', 'q37_respond_positively', 'q38_pay_attention', 'q39_imaginative_games', 'q40_cooperatively_games'))

Fit a GFMM on the SPARK shared features, project onto the SSC, and correlate profiles.

Parameters:

• spark (analysis.cohort.CohortMatrix) – The harmonised SPARK and SSC matrices.

• ssc (analysis.cohort.CohortMatrix) – The harmonised SPARK and SSC matrices.

• typing (analysis.features.Typing) – The reconciled feature typing (restricted internally to the shared features).

• category_map (dict of str to str) – Feature-to-category map for the signatures.

• n_components (int, optional) – Number of classes.

• n_init (int, optional) – Random restarts for the SPARK fit.

• n_permutations (int, optional) – Number of SSC label permutations for the null. Zero skips the null.

• seed (int, optional) – Random seed for the fit and the permutations.

• reverse_coded (tuple of str, optional) – SCQ items whose enrichment direction is flipped before the signature.

Returns:

The shared features, both signatures, the observed and null correlations, the permutation \(p\)-value, and the sample-size metrics.

Return type:

ReplicationResult

Stratification and pre-registration

Binning policies for the stratification axes.

The stratified analysis (plan section 7) re-estimates the mixture model within strata of a continuous variable: age at diagnosis in years, or the derived calendar year of diagnosis. A binning policy maps that variable onto an ordered set of named strata. The downstream fit consumes a StratumAssignment and never sees which policy produced it, so the analysis is independent of the binning choice and a policy can be swapped without touching the fitting code.

Three policies are defined, a substantive scheme plus two distribution-free ones:

• FixedBands cuts the variable at explicit, substantively motivated edges: clinical age bands, or calendar-era bands anchored on the DSM-5 boundary. The edges are a parameter, frozen at pre-registration (plan section 12) and provisional until then.

• QuantileBins cuts the variable at its own empirical quantiles, giving equal-frequency strata whose edges follow the cohort rather than a fixed rule.

• MaxEqualBins is the same equal-frequency idea but with the bin count chosen from the cohort size and the minimum stratum size, the finest split that keeps every bin above the floor, rather than a fixed number of bins.

Both use left-closed, right-open intervals \([lo, hi)\) with open outer bins, so a value that lands on an interior edge falls in the upper band. For the era axis this places a diagnosis recorded in a boundary year on the later side, matching DSM-5 taking effect in 2013.

class analysis.strata.StratumAssignment(labels, edges, codes, counts, n_missing, spec)

An ordered partition of a continuous variable into named strata.

labels

The stratum names, in ascending order of the variable.

Type:

list of str

edges

The finite interior cut points, of length len(labels) - 1. The outer bins are open, so the full partition is (-inf, edges[0]), [edges[0], edges[1]), …, [edges[-1], +inf).

Type:

list of float

codes

The per-row stratum label, an ordered categorical indexed as the input. Rows whose value is missing are unassigned (categorical NaN).

Type:

pandas.Series

counts

Assigned rows per stratum, in label order.

Type:

dict of str to int

n_missing

Rows dropped because the variable was missing.

Type:

int

spec

The serialisable policy specification, recorded in the run manifest and the frozen pre-registration so a stratification is reproducible from it alone.

Type:

dict

class analysis.strata.BinningPolicy(*args, **kwargs)

The interface every stratification axis depends on.

A policy turns a continuous variable into a StratumAssignment. Downstream code is typed against this protocol, never against a concrete policy, which is what keeps the analysis independent of the binning choice.

assign(values)

Partition values into ordered, named strata.

spec()

Return the serialisable policy specification.

class analysis.strata.FixedBands(edges, labels=None, name='fixed')

Cut at explicit, substantively motivated edges (plan section 7).

The edges are interior cut points in the units of the axis (age at diagnosis in years, or era as a calendar year). They are provisional until frozen at pre-registration.

edges

Strictly ascending interior cut points; k edges give k + 1 bands.

Type:

tuple of float

labels

Band names in ascending order. Defaults to half-open range labels.

Type:

tuple of str, optional

name

Policy name recorded in the spec.

Type:

str

assign(values)

Partition values at the fixed edges.

spec()

Return the serialisable specification, including the edges.

class analysis.strata.QuantileBins(q, name='quantile')

Cut at the variable’s own empirical quantiles, for equal-frequency strata.

The realised edges depend on the cohort, so the spec records the intended number of bins q and the assignment records the edges the cohort produced. Ties can collapse bins, so the realised count of strata can be below q.

q

The number of equal-frequency bins requested (at least 2).

Type:

int

name

Policy name recorded in the spec.

Type:

str

assign(values)

Partition values at its interior quantiles.

spec()

Return the serialisable specification (the intended q).

class analysis.strata.MaxEqualBins(min_bin_size=1000, name='max-equal')

Equal-frequency bins, as many as keep every bin above a minimum size.

The bin count is not fixed: it is the largest q whose equal-frequency split still leaves every bin at or above min_bin_size, starting from floor(n / min_bin_size) and stepping down if ties or skew leave a bin short. The result is the finest equal-frequency partition that clears the floor, so the resolution follows the cohort size and the floor rather than a hand-picked q. Like QuantileBins, the edges are the variable’s own quantiles, so the choice stays on the design side of the pre-registration firewall.

min_bin_size

The size every bin must clear; sets both the starting bin count and the floor the step-down enforces. Defaults to the phase-2 recovery floor.

Type:

int

name

Policy name recorded in the spec.

Type:

str

assign(values)

Partition values into the finest equal-frequency split above the floor.

spec()

Return the serialisable specification (the floor that sets the bin count).

analysis.strata.id_dichotomy(edge, *, low_is_impaired, name='id-dichotomy')

Split a cognitive variable into two intellectual-disability strata.

The cross-cohort cognitive-impairment axis (plan section 8) is a two-band split of a variable that indexes cognitive level: SSC full-scale deviation IQ (impaired band low, so the edge is the IQ threshold), or SPARK’s binary machine-learned impairment flag (impaired band high, so the edge is 0.5). Both cohorts use the same band names, "id" and "no_id", so a stratum matches its counterpart across cohorts regardless of which side of the edge the impaired band sits on.

Parameters:

• edge (float) – The single interior cut point: the IQ threshold for a score, or 0.5 for a 0/1 flag.

• low_is_impaired (bool) – Whether the impaired band is below the edge. True for an IQ score (a low score is impaired), False for an impairment flag (the value 1 is impaired).

• name (str, default "id-dichotomy") – Policy name recorded in the spec.

Returns:

A two-band policy labelled ("id", "no_id") in ascending order of the variable.

Return type:

FixedBands

Construct the stratification variables and the demographics for the cohort.

The phase-3 feasibility stage needs the two stratifying axes, the measurement-to-diagnosis lag, and a covariate frame, all built for the probands in the modelling cohort and on its index. This module reads only the timing and demographic columns each needs through the dscat catalogue (never a whole instrument file) and derives:

• age at diagnosis in years, from core_descriptive_variables.diagnosis_age (months);

• the calendar year of diagnosis, the registration-anchor reconstruction registration_year - age_at_registration_years + age_at_diagnosis (plan section 5);

• the lag in years, age_at_eval_years - age_at_diagnosis, the section 7a confound;

• a numeric demographics frame (sex, age at evaluation, a cognitive-impairment proxy, Hispanic ethnicity, and the race indicators) for the per-bin table and the balance check.

Implausible values are set to missing before any policy is evaluated, and counted in the diagnostics so the cleaning is visible. The age and era axes are gated only on their own validity (a childhood age at diagnosis, a reconstructed year in the cohort window), not on evaluation timing, since age at diagnosis is recalled. Only the lag is gated on the evaluation comparison. Nothing here fits a model: this is the design-side characterisation the binning policy is judged on (plan section 12 firewall).

class analysis.strata_data.StrataData(axes, lag, demographics, instrument_years, diagnostics)

The stratification variables, lag, and demographics for the cohort.

axes

Columns age_at_diagnosis_years and diagnosis_year, cleaned of implausible values, indexed by proband.

Type:

pandas.DataFrame

lag

Measurement-to-diagnosis lag in years, on the same index.

Type:

pandas.Series

demographics

Numeric covariates (sex, age at evaluation, cognitive-impairment proxy, Hispanic ethnicity, race indicators) for the demographic table and balance check.

Type:

pandas.DataFrame

instrument_years

Per-instrument completion year for the dated instruments, for the contemporaneity summary.

Type:

pandas.DataFrame

diagnostics

Per-variable counts (missing, implausible) and quantiles, plus the contemporaneity summary, all recorded in the manifest.

Type:

dict

analysis.strata_data.derive_axes(diagnosis_age_months, registration_year, age_at_registration_years, age_at_eval_years)

Derive the two axes and the lag from the raw timing fields.

Age at diagnosis is a parent-recalled value, valid in the childhood range and independent of when the instruments were completed, so the two axes are gated only on their own validity. The lag alone depends on the evaluation timing. Implausible values set to missing: an age at diagnosis below zero or above the childhood ceiling; a reconstructed year outside the cohort window; and a lag more negative than the one-year rounding tolerance (a diagnosis genuinely postdating evaluation).

analysis.strata_data.summarise(axes, lag, n_total)

Summarise the missingness, implausible drops, and spread of each variable.

analysis.strata_data.build_strata_data(root, version, index, age_at_eval, sex)

Build the stratification variables and demographics for a cohort.

Parameters:

• root (pathlib.Path) – Repository root.

• version (str) – SPARK release version.

• index (pandas.Index) – The modelling-cohort proband index the variables are built for.

• age_at_eval (pandas.Series) – Age at evaluation per proband (from the cohort covariates), on index.

• sex (pandas.Series) – Encoded sex per proband (from the cohort covariates), on index.

Returns:

The axes, lag, demographics, per-instrument years, and diagnostics.

Return type:

StrataData

Acceptance requirements a binning policy must meet before it is frozen.

A binning policy (analysis.strata) is only eligible for the confirmatory stratified fit (plan section 7) if its partition can actually be fitted and is not silently confounded. This module evaluates a concrete policy against a tiered requirement set, computed from the cohort, the stratifying variable, and (optionally) the measurement-to-diagnosis lag and a covariate frame. Nothing here fits a per-stratum mixture model or reads class drift, so the checks stay on the design side of the pre-registration firewall: they decide whether a partition is eligible to be fitted, not whether the fit gives a wanted answer.

Three tiers:

• Tier 1, eligibility gates (hard pass or fail). Every bin clears the empirical \(N_\text{min}\) floor and the projected smallest-class floor, coverage of the modelling cohort is high, and the partition is valid. A policy that fails any Tier 1 gate is ineligible.

• Tier 2, confound and balance (reported with a flag, not a hard fail). Size balance, lag entanglement and small-lag retention for the era axis, covariate balance, and edge robustness. These inform the covariate-versus-subsample decision rather than reject a policy outright.

• Tier 3, demographics. A per-bin summary with standardised differences across the extreme bins, both the manuscript’s per-stratum table and a check that any drift is not trivially a composition artefact.

Thresholds are settable (RequirementThresholds) and recorded in the report, so the frozen pre-registration carries the exact values a policy was judged against. The defaults follow the phase-2 findings: a per-bin floor of 1000 (the practical recovery floor from the \(N_\text{min}\) sweep) and a smallest class near 15 per cent of a bin.

class analysis.requirements.RequirementThresholds(min_bin_size=1000, smallest_class_fraction=0.15, min_projected_smallest_class=150, min_coverage=0.9, min_n_bins=2, max_bin_share=0.65, small_lag_years=2.0, max_lag_correlation=0.3, max_reassigned_fraction=0.05, edge_perturbation=1.0, smd_flag=0.2)

The numeric criteria a policy is judged against.

min_bin_size

Smallest allowed assigned count in any bin (the phase-2 recovery floor).

Type:

int

smallest_class_fraction

Reference smallest-class proportion, used to project a bin’s smallest class.

Type:

float

min_projected_smallest_class

Smallest allowed projected smallest-class count (bin size times the fraction).

Type:

int

min_coverage

Smallest allowed fraction of the modelling cohort assigned (non-missing variable).

Type:

float

min_n_bins

Fewest non-empty bins (a contrast needs at least two).

Type:

int

max_bin_share

Largest share of the assigned cohort a single bin may hold before it is flagged.

Type:

float

small_lag_years

Lag cut for the small-lag subsample used to probe the era axis.

Type:

float

max_lag_correlation

Largest tolerated Spearman correlation between the variable and the lag.

Type:

float

max_reassigned_fraction

Largest share of probands that may move bin under an edge perturbation.

Type:

float

edge_perturbation

Size of the edge perturbation, in the variable’s units.

Type:

float

smd_flag

Standardised difference across the extreme bins above which a covariate is flagged.

Type:

float

class analysis.requirements.RequirementResult(key, tier, description, status, observed, threshold, detail)

The outcome of one requirement.

key

Short identifier.

Type:

str

tier

1, 2, or 3.

Type:

int

description

What the requirement checks.

Type:

str

status

"pass" or "fail" for Tier 1, "ok" or "flag" for Tier 2, "report" for Tier 3, and "skipped" when an input was not supplied.

Type:

str

observed

The measured quantity.

Type:

float or None

threshold

The criterion it was compared against.

Type:

float or None

detail

A human-readable summary.

Type:

str

class analysis.requirements.PolicyReport(spec, n_total, n_assigned, counts, results, demographics=None)

The full evaluation of one policy against the requirement set.

spec

The policy specification (from spec()).

Type:

dict

n_total

Rows in the modelling cohort.

Type:

int

n_assigned

Rows assigned to a bin (the rest are missing on the variable).

Type:

int

counts

Assigned rows per bin, in label order.

Type:

dict of str to int

results

One entry per requirement, across all three tiers.

Type:

list of RequirementResult

demographics

The per-bin covariate summary with the extreme-bin standardised differences, or None when no covariates were supplied.

Type:

pandas.DataFrame or None

property eligible: bool

Whether every Tier 1 gate passed.

property flags: list[str]

Keys of the Tier 2 checks that were flagged.

to_frame()

Return the requirement results as a table, one row per requirement.

analysis.requirements.evaluate_policy(policy, variable, *, lag=None, covariates=None, thresholds=RequirementThresholds(min_bin_size=1000, smallest_class_fraction=0.15, min_projected_smallest_class=150, min_coverage=0.9, min_n_bins=2, max_bin_share=0.65, small_lag_years=2.0, max_lag_correlation=0.3, max_reassigned_fraction=0.05, edge_perturbation=1.0, smd_flag=0.2))

Evaluate a binning policy against the tiered requirement set.

Parameters:

• policy (analysis.strata.BinningPolicy) – The concrete policy to test.

• variable (pandas.Series) – The stratifying variable over the modelling cohort (age at diagnosis or era).

• lag (pandas.Series, optional) – The measurement-to-diagnosis lag in years, on the same index. Enables the Tier 2 lag checks (the era-axis defence).

• covariates (pandas.DataFrame, optional) – Numeric (one-hot encoded where categorical) covariates on the same index. Enables the Tier 2 covariate-balance check and the Tier 3 demographic table.

• thresholds (RequirementThresholds, optional) – The criteria to judge against. Defaults to the phase-2 values.

Returns:

The per-requirement results, the eligibility verdict, and the demographic table.

Return type:

PolicyReport

Stratified analysis

Class drift between a stratum fit and the pooled reference, and its permutation null.

The stratified analysis (plan section 7, frozen in section 12a) asks whether the four reference classes move when the mixture model is re-estimated within a stratum of age at diagnosis or diagnostic era. This module measures that movement and calibrates it, with the expensive part (the fits) separated from the cheap, method-dependent part (alignment and distance), so a different alignment or distance can be tried without re-fitting.

The unit that is fitted and stored is a StratumSummary: per-class feature means (centroids) and dispersions (standard deviations), plus the contingency of the fit’s labels against the pooled reference labels on the same probands. These are the method-independent sufficient statistics. From them:

• an AlignmentMethod maps the fit’s arbitrary class ids to the reference classes. MembershipJaccard (the default) aligns on who is in each class, since a stratum is a subset of the pooled cohort and so carries both labellings on the same probands; this distinguishes a class that moved (same members, shifted centroid) from one that reorganised (different members), which a centroid-only alignment cannot. CentroidHungarian aligns on centroid distance instead.

• a DistanceMethod measures how far each aligned class moved. Mahalanobis (the default) is the covariance-aware distance between centroids, so correlated features count once rather than many times; StandardisedEuclidean and MeanAbsolute are the diagonal (covariance-blind) distances between centroids; JensenShannon compares the class-conditional distributions, treating each feature as Gaussian with the per-class mean and dispersion, so it sees a change in spread that the centroid distances miss.

The drift is read against the between-class separation (the same distance between distinct reference classes) so a shift is on the scale of the partition, and against a permutation null: pseudo-strata of the same sizes from random partitions of the cohort, so the observed shift is read against same-size random partitions (beyond the 95th percentile, then FDR controlled). The alignment also reports its confidence (per-class Jaccard, overall adjusted Rand index), so a large shift with low overlap is flagged as reorganisation, not drift.

summarise_pseudo_stratum() is a top-level, picklable unit of work (fit one subset, return its summary), so the null can be spread across a process pool.

analysis.drift.is_degenerate_fit(fit)

Return whether a fit diverged to non-finite parameters without raising.

Under fractional (kernel) weights the covariate GLM can blow up: it emits overflow and invalid-value warnings and leaves a non-finite log-likelihood rather than raising an exception the workers could catch. Such a fit’s labels are meaningless, so it is treated as degenerate and dropped, the same as a fit that raised.

class analysis.drift.ReferenceModel(centroids, dispersions, pooled_sd, precision, labels)

The pooled reference solution, the fixed target every stratum is compared against.

centroids

Reference class-by-feature centroids (means).

Type:

pandas.DataFrame

dispersions

Reference class-by-feature standard deviations, for the distributional distance.

Type:

pandas.DataFrame

pooled_sd

Per-feature standard deviation across the cohort, the diagonal-distance normaliser.

Type:

pandas.Series

precision

Inverse of the Ledoit-Wolf-shrunk pooled within-class covariance, in the column order of centroids. Shrinkage keeps it well-conditioned at 238 features.

Type:

numpy.ndarray

labels

The reference (pooled) class per proband, used to build a stratum’s contingency.

Type:

pandas.Series

as_stratum()

Return the reference as a stratum summary, for the between-class separation.

analysis.drift.build_reference(measurement_data, labels)

Build the reference model from the pooled fit’s measurement data and labels.

Computes the per-class centroids and dispersions, the per-feature pooled spread, and the Ledoit-Wolf-shrunk precision matrix of the pooled within-class covariance (the residuals of each proband from its class mean). Shrinkage is what makes the 238-feature covariance invertible and stable.

class analysis.drift.StratumSummary(centroids, dispersions, contingency, n)

The method-independent summary of one fit: centroids, dispersions, and contingency.

These are the sufficient statistics for any alignment or distance, so they are stored once and re-measured cheaply when the method changes.

centroids

Fit class-by-feature means.

Type:

pandas.DataFrame

dispersions

Fit class-by-feature standard deviations.

Type:

pandas.DataFrame

contingency

Counts of the fit’s classes (rows) against the reference classes (columns).

Type:

pandas.DataFrame

n

Number of probands.

Type:

int

class analysis.drift.ClassAlignment(mapping, quality, overall)

A mapping from fit classes to reference classes, with its confidence.

mapping

Fit class id to reference class id.

Type:

dict of int to int

quality

Per reference class, the match confidence (Jaccard for membership, a normalised closeness for centroid alignment); higher is more confident.

Type:

dict of int to float

overall

The adjusted Rand index between the two labellings (membership), or the mean per-pair quality (centroid). A low value means the partition reorganised rather than shifted.

Type:

float

analysis.drift.common_columns(source, reference)

Return the feature columns shared by two matrices, in reference order.

analysis.drift.contingency_table(fit_labels, reference_labels, weights=None)

Cross-tabulate a fit’s labels against the reference labels over the shared probands.

With weights each proband contributes its weight to its cell rather than a count of one, so a kernel fit’s contingency is the weighted overlap. Unweighted (the default) gives the plain counts the hard-bin analysis uses.

analysis.drift.adjusted_rand_index(table)

Return the adjusted Rand index between two labellings, from their contingency table.

Chance-corrected agreement (0 is chance, 1 is identical partitions), computed from the counts directly so it needs only the stored contingency, not the per-proband labels.

class analysis.drift.AlignmentMethod(*args, **kwargs)

Map a stratum fit’s classes to the reference classes.

align(stratum, reference)

Return the fit-to-reference class mapping and its confidence.

class analysis.drift.MembershipJaccard(name='membership')

Align on shared membership: pair classes by maximal Jaccard overlap of their probands.

The most direct alignment, since a stratum is a subset of the pooled cohort, so each proband carries both labellings. The Jaccard normalises for the very unequal class sizes, so the largest class does not dominate the match. The overall confidence is the adjusted Rand index of the two labellings.

align(stratum, reference)

Align by Hungarian assignment on one minus the Jaccard overlap.

class analysis.drift.CentroidHungarian(name='centroid')

Align on centroid distance: pair classes by the closest standardised centroids.

A fallback that uses only the centroids, so it cannot tell a class that moved from one that reorganised. Kept to cross-check the membership alignment: a disagreement between the two flags an unsafe mapping.

align(stratum, reference)

Align by Hungarian assignment on the standardised centroid distance.

class analysis.drift.DistanceMethod(*args, **kwargs)

Measure the distance one aligned class moved between a stratum and the reference.

class_distance(stratum, fit_class, reference, ref_class)

Distance between a stratum class and its aligned reference class.

class analysis.drift.StandardisedEuclidean(name='euclidean')

Standardised Euclidean distance: the root-mean-square per-feature shift in SD units.

A diagonal (covariance-blind) distance: it treats the features as independent.

class_distance(stratum, fit_class, reference, ref_class)

Root mean square of the standardised per-feature difference.

class analysis.drift.FullStandardisedEuclidean(name='euclidean-full')

Full (unaveraged) standardised Euclidean distance: the L2 norm in SD units.

The sum-norm counterpart of StandardisedEuclidean: the square root of the summed, not averaged, squared per-feature shift, $lVert delta / sigma rVert$. This is the convention the effect-size trajectory uses for a class’s displacement magnitude (analysis.trajectory_local.grain_magnitude(), an unaveraged norm over the grain), so a between-class separation measured this way puts numerator and denominator on the same scale: a displacement then reads as a genuine fraction of the mean inter-class gap. The averaged StandardisedEuclidean divides by an extra $sqrt{n}$ in the feature count $n$ relative to this, so mixing the two inflates a separation-scaled magnitude by that factor.

class_distance(stratum, fit_class, reference, ref_class)

L2 norm of the standardised per-feature difference (summed, not averaged).

class analysis.drift.MeanAbsolute(name='mean-abs')

Mean absolute per-feature shift in SD units, an outlier-robust diagonal distance.

class_distance(stratum, fit_class, reference, ref_class)

Mean absolute standardised per-feature difference.

class analysis.drift.Mahalanobis(name='mahalanobis')

Mahalanobis distance between centroids, using the shrunk within-class precision.

The covariance-aware distance: correlated features contribute once rather than many times, so a coordinated shift across a correlated block of symptoms is not double-counted. The default, as the statistically proper multivariate distance. Centroids are reindexed to the reference feature order; a feature absent from the stratum contributes no shift.

class_distance(stratum, fit_class, reference, ref_class)

Square root of the precision-weighted squared centroid difference.

class analysis.drift.JensenShannon(name='jsd')

Mean per-feature Jensen-Shannon divergence between the class-conditional distributions.

Each feature’s class-conditional is treated as Gaussian with the per-class mean and dispersion, so the divergence sees a change in spread, not only in location, which the centroid distances cannot. Bounded in [0, 1] per feature and averaged over the shared features. The Gaussian treatment is an approximation for the binary and categorical-coded features.

class_distance(stratum, fit_class, reference, ref_class)

Mean per-feature Jensen-Shannon divergence over the shared features.

analysis.drift.class_distances(stratum, reference, mapping, distance)

Per reference class, the distance its aligned stratum class sits from it.

analysis.drift.class_separation(reference, distance)

Mean distance between distinct reference classes, the drift baseline.

The same distance the drift uses, measured between the reference classes themselves and averaged over pairs, so drift can be read as a fraction of the gap between distinct classes.

class analysis.drift.DriftResult(distances, alignment)

One stratum’s drift: per-class distance plus the alignment that produced it.

analysis.drift.compute_drift(stratum, reference, alignment, distance)

Align a stratum to the reference and measure each aligned class’s drift.

Pure and cheap (no fitting): the method-dependent step run over stored summaries, so a different alignment or distance re-measures without re-fitting.

analysis.drift.null_partition(index, sizes, seed)

Partition index into consecutive random chunks of the given sizes.

Shuffles the proband index with a seeded generator, then splits it into blocks of sizes, so the pseudo-strata have the same sizes as the real strata but no relation to the stratifying axis. The seed is the permutation index, so a resumed null reproduces the same partitions.

analysis.drift.summarise(measurement_data, labels, reference_labels, weights=None)

Build a StratumSummary from a fit’s measurement data and labels.

Computes per-class means and standard deviations and the contingency against the reference labels, the method-independent statistics every distance and alignment is derived from. With weights (a kernel fit) the centroids, dispersions, and contingency are weighted; without them (a hard-bin fit) they are the plain per-class statistics.

analysis.drift.summarise_pseudo_stratum(features, covariates, typing, reference_labels, n_init, seed)

Fit the GFMM on one subset and return its method-independent summary.

A top-level function so it pickles for a process pool. Stores the centroids, dispersions, and reference contingency, not a drift value, so the alignment and distance can be chosen (and changed) afterwards without re-fitting. Returns None if the fit is degenerate (a singular covariate GLM), so the caller drops that pseudo-stratum from the null rather than letting one bad refit abort the whole run.

analysis.drift.serialise_summary(summary, perm, s_idx)

Serialise a stratum summary to a JSON-able record for the null store.

The null fits are stored as their summaries, not their drift, so the alignment and distance can be chosen afterwards. One record per pseudo-stratum, keyed by its permutation and stratum index.

analysis.drift.deserialise_summary(record)

Rebuild a StratumSummary from a serialised null-store record.

analysis.drift.benjamini_hochberg(p_values, q=0.05)

Return a boolean mask of the hypotheses that pass Benjamini-Hochberg FDR control.

Controls the false-discovery rate at q across the strata-by-class drift tests (plan section 12a). A hypothesis is rejected if its p-value is at or below the largest threshold q * rank / m it satisfies, where m is the number of finite p-values; NaN p-values (a degenerate stratum) never pass.

analysis.drift.read_against_null(observed, null_draws)

Read an observed drift against its size-matched null distribution.

Returns the null 95th percentile, whether the observed shift exceeds it, and the permutation p-value with the Phipson-Smyth add-one correction (so the smallest p is 1 / (n + 1) rather than zero). The decision threshold and the FDR step across classes are applied by the caller over these per-class reads.

Attribute a class’s movement between two fits to features and to probands (archived).

Archived. This is the refit-era attribution: it interprets the drift of a mixture re-estimated within a stratum. The category attribution ($H_0^F$) is now read from the single cached fit by the block-attribution engine (analysis.blocks, the additive category decomposition), so this module is no longer $H_0^F$’s evidence. It is kept because it renders the membership-churn and mover-versus-stayer figures on the refit pilot page, which the single-fit engine cannot reproduce (a frozen fit relabels no proband).

The drift stage (analysis.drift) measures how far a reference class moves when the mixture is re-estimated within a stratum, as one distance per class. This module opens that distance up. It asks which features carry the shift, and which probands changed class, so a movement reads as “these features, these people” rather than a single number.

Two families, both cheap readouts over the stored fits (no re-fitting):

• Centroid-shift decomposition splits a class’s distance into signed per-feature contributions that sum back to it. For the Mahalanobis distance the split is the term-by-term expansion of the quadratic form $Deltamu^top P Deltamu = sum_i Deltamu_i (PDeltamu)_i$, so a coordinated shift across correlated features is charged together, matching the distance the drift stage reports. The diagonal split $(Deltamu_i / sigma_i)^2$ is the covariance-blind cross-check.

• Mover and stayer attribution labels each shared proband a class member that stayed or left between the two fits, then contrasts the two groups over a feature frame (the clustered features, the held-out SPARK variables, or, later, the genetic scores). The movement is then explained by what marks the probands the class shed.

The unit both families read is a Comparison: two labellings on the same proband index plus the fit summaries and their alignment. It is agnostic to how the second labelling was produced (a hard stratum bin, a kernel-weighted focal point, a partition-tree node), so the same attribution runs on any of them.

class analysis.attribution.Comparison(reference, stratum, ref_labels, fit_labels, alignment)

Two labellings on a shared proband index, with their summaries and alignment.

reference

The pooled reference: centroids, dispersions, pooled spread, and the within-class precision.

Type:

analysis.drift.ReferenceModel

stratum

The second fit’s method-independent summary (centroids, dispersions, contingency).

Type:

analysis.drift.StratumSummary

ref_labels

Reference class per proband, over the shared index.

Type:

pandas.Series

fit_labels

Second-fit class per proband, over the shared index, in the fit’s own class ids.

Type:

pandas.Series

alignment

The fit-to-reference class mapping and its confidence.

Type:

analysis.drift.ClassAlignment

property shared_index: Index

Probands carrying both labellings.

movements()

One Movement per reference class the alignment mapped, in class order.

class analysis.attribution.Movement(comparison, ref_class, fit_class)

One aligned class pair: how reference ref_class re-expressed as fit_class.

comparison

The two-fit pairing this movement is drawn from.

Type:

Comparison

ref_class

The fixed reference class being tracked.

Type:

int

fit_class

Its aligned partner in the second fit.

Type:

int

analysis.attribution.signed_shift(movement)

Per-feature centroid shift in pooled-SD units, signed (stratum minus reference).

The direction of movement per feature, so a decomposition magnitude can be read together with whether the stratum class sits above or below the reference class on that feature.

class analysis.attribution.DecompositionMethod(*args, **kwargs)

Split a class’s squared distance into signed per-feature contributions.

contributions(movement)

Signed per-feature contributions that sum to the class’s squared distance.

class analysis.attribution.MahalanobisContribution(name='mahalanobis')

Term-by-term split of the Mahalanobis distance: c_i = delta_i (P delta)_i.

The contributions sum to the squared Mahalanobis distance the drift stage reports, so a class’s movement has an additive, covariance-aware feature breakdown. A contribution can be negative when a feature’s shift offsets a correlated block, which the diagonal split cannot show. The default, matching the default drift distance.

contributions(movement)

Return the precision-weighted per-feature contributions to the squared distance.

class analysis.attribution.StandardisedContribution(name='standardised')

Diagonal split c_i = (delta_i / sigma_i)**2: non-negative and covariance-blind.

Sums to the total squared standardised shift, the cross-check on the Mahalanobis split: a feature that ranks high here but low under Mahalanobis is one whose shift is shared with a correlated block rather than its own.

contributions(movement)

Return the squared standardised per-feature shift.

analysis.attribution.category_of(feature, category_map)

Return a feature’s category, or "unmapped" for an absent or blank entry.

The author map leaves a few CBCL composites (for example total_problems_t_score) with a blank category, which loads as NaN, and does not list every clustered feature. Both resolve to "unmapped" here so the category is always a clean string, and so a blank category never silently drops a feature’s contribution from the totals.

analysis.attribution.category_totals(contributions, category_map)

Aggregate per-feature contributions into the literature categories.

Sums the signed contributions within each category (plan section 6, step 7), so a class’s movement reads at the level of the categories, not only the individual features. Every feature is kept (blanks and absentees under "unmapped"), so the totals sum to the same squared distance the per-feature contributions do.

analysis.attribution.membership_counts(movement)

Count the stayers, leavers, and joiners for one class between the two fits.

Returns:

n_stayers (in the class under both fits), n_leavers (a reference member the second fit dropped), and n_joiners (a second-fit member the reference did not assign to the class). Leavers plus joiners over their union is the class churn, one minus the Jaccard overlap the alignment reports.

Return type:

dict

analysis.attribution.movers(movement, kind='either')

Mark the probands that changed class membership between the two fits.

A class can move in two ways: it can shed members (leavers, reference members the second fit drops) and it can absorb members (joiners, second-fit members the reference did not assign there). A class that keeps every member but pulls in new ones still drifts, so a leaver-only view misses it; kind="either" (the default) captures both.

Parameters:

• movement (Movement) – The aligned class pair to score.

• kind (str, default "either") – "either" marks leavers and joiners against the stayers over the union of the two memberships (the class churn). "leavers" restricts to the reference members and marks those the second fit dropped. "joiners" restricts to the second-fit members and marks those the reference did not assign to the class.

Returns:

Boolean moved over the relevant probands: the complement are the stayers the contrast is run against.

Return type:

pandas.Series

class analysis.attribution.AttributionResult(importances, n_movers, n_stayers, method)

Ranked feature attributions distinguishing movers from stayers.

importances

One row per feature: the signed effect (a standardised mean difference or a model coefficient), its magnitude for ranking, and, where the method provides them, a p_value and fdr_significant flag. Sorted by magnitude, most distinguishing first. Empty when the contrast is undegenerate (only movers or only stayers).

Type:

pandas.DataFrame

n_movers

Number of probands that left the class.

Type:

int

n_stayers

Number that stayed.

Type:

int

method

The contrast method’s name.

Type:

str

class analysis.attribution.ContrastMethod(*args, **kwargs)

Rank the features that distinguish movers from stayers.

contrast(moved, features)

Return the ranked feature attributions for one class’s movers against its stayers.

class analysis.attribution.UnivariateContrast(name='univariate')

Per-feature standardised mean difference between movers and stayers.

Cohen’s $d$ per feature (positive when movers score higher), with a Welch $t$-test, Benjamini-Hochberg controlled across features. Model-free, fast, and always defined; the first read on what marks the movers, feature by feature, before any multivariate model. Missing values are handled per feature, so it runs on the held-out variables as it does on the clustered features.

contrast(moved, features)

Contrast movers and stayers with a per-feature effect size and Welch test.

class analysis.attribution.LogisticContrast(name='logistic', c=1.0)

L1-regularised logistic regression of mover status on standardised features.

Signed coefficients rank the features that jointly distinguish movers from stayers, so correlated features share credit rather than each scoring the marginal difference. Features are standardised and missing values median-imputed within the contrasted probands; a constant or fully missing feature is dropped. The L1 penalty keeps the ranking sparse, so a short list of features carries the movement.

contrast(moved, features)

Fit the penalised model and return the signed coefficients, ranked by magnitude.

Reference schemes: what each local fit’s drift is measured against.

The drift stage (analysis.drift) measures how far the reference classes move when the mixture is re-estimated within a stratum, and it compares every stratum to one fixed target, the pooled reference (section 6a). That target is itself a design choice with consequences for what the drift means, so this module turns it into a pluggable axis, orthogonal to the localisation scheme (analysis.localise) that forms the fits. A run picks a reference scheme the same way it picks an alignment or a distance.

A stratum is a proper subset of the pooled cohort, so comparing it to the pooled reference has the stratum contribute to the target it is judged against. The size-matched permutation null (plan section 12a) calibrates that overlap away rather than removing it: a null pseudo-stratum is a subset of the same size, so it carries the same overlap, and the residual pull is conservative, because a stratum’s own members draw the pooled centroids toward it and so understate its drift. Two alternative targets remove the overlap by construction instead:

• PooledReference (the frozen primary) compares each stratum to the pooled reference. The stratum and the reference share probands, so membership alignment applies and a class that moved can be told from one that reorganised.

• PairwiseReference compares each stratum to a neighbouring stratum. The two are disjoint, so nothing is shared and the comparison is independent, at the cost of centroid alignment (the move-versus-reorganise distinction needs shared members). It reads change along the axis rather than distance from a single pooled partition, the honest test when the existence of one reference structure is itself in question.

The work splits in two. Topology (which fit pairs with which) is a Pairing, a pure function of the fit labels and positions, so it is tested without any fit. Resolution turns a pairing into a DriftComparison by building the concrete ReferenceModel for it through a ReferenceResolver, which the caller backs with the cached fits. The measurement itself is unchanged: compute_drift() runs over the query summary and the resolved reference, so no distance or alignment code is duplicated here.

class analysis.reference_scheme.QueryFit(label, position, summary)

The query side of a comparison: one local fit and where it sits on the axis.

label

The fit’s name, unique within a run (a stratum name, or focal=6.5).

Type:

str

position

The fit’s location on the axis, in the axis units, used to order neighbours for the pairwise topology.

Type:

float

summary

The fit’s method-independent summary (centroids, dispersions, contingency).

Type:

analysis.drift.StratumSummary

class analysis.reference_scheme.Pairing(query_label, reference_kind, reference_label=None)

A pure-topology comparison: which fit is the query and where its reference comes from.

Carries no fit data, so a scheme’s topology is decided and tested without building any reference. resolve_comparisons() turns a pairing into a concrete DriftComparison.

query_label

The fit whose drift is measured.

Type:

str

reference_kind

How the reference is obtained: "pooled" (the fixed pooled reference) or "promote" (build a reference from the reference_label fit).

Type:

str

reference_label

The other fit a "promote" pairing draws its reference from; None otherwise.

Type:

str or None

class analysis.reference_scheme.DriftComparison(query_label, position, query, reference, alignment, reference_label=None)

A resolved comparison: a query summary against a concrete reference, with its alignment.

query_label

The query fit’s name.

Type:

str

position

The query fit’s position on the axis.

Type:

float

query

The query fit’s summary.

Type:

analysis.drift.StratumSummary

reference

The target the query is aligned and measured against.

Type:

analysis.drift.ReferenceModel

alignment

The alignment method’s registry name ("membership" when the two share probands, "centroid" when they are disjoint).

Type:

str

reference_label

The fit the reference was built from, for a pairwise or leave-one-out comparison; None for the pooled reference.

Type:

str or None

class analysis.reference_scheme.ReferenceResolver(*args, **kwargs)

Build the concrete reference model a pairing names, backed by the caller’s fits.

A scheme names its references abstractly (the pooled reference, or a named neighbour); the resolver turns each name into a ReferenceModel. The caller implements it over the cached fits, so a pairwise reference is built from a neighbour stratum with no re-fit.

pooled()

Return the fixed pooled reference.

promote(label)

Return the fit named label promoted to a reference (its centroids as a target).

class analysis.reference_scheme.MappingResolver(pooled_reference, promoted)

Resolve references from a fixed pooled reference and a per-label promoted map.

The concrete resolver the drift stage uses. pooled returns the pooled reference; promote looks a fit’s own promoted reference up in promoted, which the caller has built once from each cached fit, so a pairwise comparison reuses a neighbour fit with no re-fit.

pooled_reference

The fixed pooled reference.

Type:

analysis.drift.ReferenceModel

promoted

Each fit’s own promoted reference, keyed by the fit label.

Type:

Mapping of str to analysis.drift.ReferenceModel

pooled()

Return the fixed pooled reference.

promote(label)

Return the promoted reference for the fit named label.

class analysis.reference_scheme.ReferenceScheme(*args, **kwargs)

Decide which fit each drift is measured against.

pairings(ordered)

Return the comparison topology for fits given as (label, position) pairs.

spec()

Return the serialisable scheme specification for the run manifest.

class analysis.reference_scheme.PooledReference(name='pooled', default_alignment='membership')

Compare every stratum to the pooled reference: the frozen confirmatory primary.

The stratum is a subset of the pooled cohort, so both labellings exist on its probands and membership alignment applies. This reproduces the current drift stage exactly, so it is the regression anchor for the abstraction.

pairings(ordered)

One pooled-reference pairing per fit, in the given order.

spec()

Return the specification (the scheme name and its alignment default).

class analysis.reference_scheme.PairwiseReference(mode='adjacent', name='pairwise', default_alignment='centroid')

Compare each stratum to another stratum, so the comparison carries no self-overlap.

The two strata are disjoint, so alignment is by centroid only. mode sets which pairs are formed: "adjacent" (the default) pairs each fit with its successor along the axis, a sequence of local comparisons that reads change along the axis; "all-pairs" compares every earlier fit to every later one. Fits are ordered by position, and each query’s reference is the later fit, so the direction of every comparison runs along the axis.

mode

"adjacent" or "all-pairs".

Type:

str

pairings(ordered)

Pair each fit with a later fit, adjacent by default, all-pairs otherwise.

spec()

Return the specification (the scheme name, its mode, and its alignment default).

analysis.reference_scheme.resolve_comparisons(scheme, queries, resolver)

Turn a scheme’s topology into resolved comparisons, ready to measure.

Reads the topology from scheme.pairings over the query labels and positions, then builds each pairing’s reference through resolver: the pooled reference for a "pooled" pairing, or the named neighbour promoted to a reference for a "promote" pairing. The alignment is the scheme’s default, which the caller may override per comparison.

analysis.reference_scheme.measure(comparisons, distance='mahalanobis')

Measure each resolved comparison’s drift, keyed by the query label.

A thin pass over compute_drift(): each comparison aligns its query to its reference with its own alignment method and measures every aligned class with the shared distance. No fitting happens here, so re-measuring with a different distance is cheap.

Project the classes into a discriminant space and quantify how their centroids move.

The stratified fits give, for each stratum, a class-by-feature centroid aligned to the pooled reference (analysis.drift). This module turns those centroids into the material a trajectory figure needs, and measures the shape of each class’s path, all in aggregate (class-level) terms so nothing per-proband leaves the stage.

Three pieces:

• an Embedding is a linear-discriminant projection fitted on the pooled reference classes. With four classes it spans three axes ($K - 1$), the coordinates in which the classes are maximally separated, so a class moving towards another is read directly. The projection is linear, so distances in it are honest, unlike a nonlinear embedding. It is an illustration; the drift claim rests on the full-dimensional distances of analysis.drift, not on this picture.

• directional_test() asks whether a class moves with the stratifying axis. The statistic is the net displacement from the first third of the strata to the last third, in standardised units. Permuting the stratum order holds the non-directional between-stratum scatter fixed and destroys only the ordering, so a net displacement beyond the shuffled null is movement tied to the axis, not scatter. This is a pilot measure on the observed aligned centroids; the confirmatory test is the continuous-trend regression against the refit permutation null (plan section 12a).

• roughness_metrics() reports the mean step between adjacent strata against the step that independent sampling of a class of that size would produce, so a jagged path can be read as sampling noise rather than movement.

class analysis.trajectory.Embedding(transformer, mean, sd, columns, explained_variance_ratio)

A linear-discriminant projection of the pooled reference classes.

transformer

The fitted transformer, taking standardised feature vectors to discriminant axes.

Type:

sklearn.discriminant_analysis.LinearDiscriminantAnalysis

mean, sd

Per-feature pooled mean and standard deviation used to standardise before projecting, in columns order.

Type:

numpy.ndarray

columns

Feature order the transformer was fitted on.

Type:

list of str

explained_variance_ratio

Fraction of between-class variance carried by each discriminant axis.

Type:

numpy.ndarray

property n_components: int

Return the number of discriminant axes.

analysis.trajectory.fit_embedding(measurement_data, labels, n_components=3)

Fit a linear-discriminant embedding of the pooled classes.

Parameters:

• measurement_data (pandas.DataFrame) – The pooled proband-by-feature matrix.

• labels (pandas.Series) – The reference class per proband, indexed like measurement_data.

• n_components (int, optional) – Discriminant axes to keep, capped at the number of classes minus one. Defaults to 3, the full space for a four-class solution.

Returns:

The fitted projection, with the standardisation it applies before transforming.

Return type:

Embedding

analysis.trajectory.project(embedding, centroids)

Project class-by-feature centroids into the discriminant axes.

Parameters:

• embedding (Embedding) – A fitted embedding.

• centroids (pandas.DataFrame) – Centroids to project, carrying at least the embedding’s feature columns.

Returns:

One row per input centroid, one column per discriminant axis.

Return type:

numpy.ndarray

analysis.trajectory.directional_test(trajectory, *, seed, n_shuffle)

Test whether one class’s trajectory moves with the stratifying axis.

The statistic is the net displacement between the first and last third of the strata, in standardised units. The null permutes the stratum ordering, which preserves the non-directional scatter and removes only the tie to the axis, so the observed value is read as a percentile of the shuffled distribution.

Parameters:

• trajectory (numpy.ndarray) – The class’s standardised centroids, ordered by stratum, shape (n_strata, n_features).

• seed (int) – Seed for the ordering shuffle, for reproducibility.

• n_shuffle (int) – Number of ordering permutations.

Returns:

net (observed net displacement), null95 (95th percentile of the null), p (one-sided, with the Phipson-Smyth add-one), and significant (p < 0.05).

Return type:

dict

analysis.trajectory.roughness_metrics(trajectory, sizes, within_sd)

Measure a class trajectory’s step size against its sampling-noise expectation.

A stratum’s centroid is a mean over that stratum’s members, so two adjacent strata differ by sampling noise even with no real movement. The expected step under sampling alone is $sqrt{sum_f w_f^2 (1/n_i + 1/n_j)}$, where $w_f$ is the within-class standard deviation of feature $f$ (standardised) and $n_i, n_j$ are the adjacent class sizes. A step near this expectation is noise; a step well above it is movement.

Parameters:

• trajectory (numpy.ndarray) – The class’s standardised centroids, ordered by stratum, shape (n_strata, n_features).

• sizes (numpy.ndarray) – The class size in each stratum, in the same order.

• within_sd (numpy.ndarray) – The class’s per-feature within-class standard deviation, standardised.

Returns:

step (mean step between adjacent strata), sampling_noise (mean expected step under sampling), and snr (mean of their per-step ratio).

Return type:

dict

Score-based measurement invariance from a single cached fit (plan section 7e).

The stratified analysis (section 7) asks whether the four reference classes are stable as the mixture is re-estimated within strata of age at diagnosis or diagnostic era. Every refit-based scheme carries a fit cost and a permutation null. This module answers the same question from a single cached fit with an analytic null and no refitting, following the empirical fluctuation process of Merkle and Zeileis (2013, Psychometrika) and Merkle, Fan and Zeileis (2014).

The idea. At the maximum-likelihood estimate every proband contributes a score, the gradient of its log-likelihood with respect to each class-profile parameter. Fisher’s identity gives the casewise score with respect to the class-$k$ value of feature $j$ as $r_{ik},partial_theta log f_j(x_{ij};theta_{jk})$, where $r_{ik}$ is the posterior responsibility (predict_proba). The focal parameters are the class-conditional locations, the profiles the whole analysis measures:

• a Gaussian mean has score $r_{ik}(x_{ij}-mu_{jk})/sigma^2_{jk}$;

• a Bernoulli probability has score $r_{ik}(x_{ij}-p_{jk})/(p_{jk}(1-p_{jk}))$;

• a categorical outcome $l$ has the multinomial-logit score $r_{ik}(mathbb{1}[x_{ij}=l]-p_{jkl})$.

All three are validated against a central finite difference of the per-sample log-likelihood to machine precision (numerical_score(), the correctness gate). The categorical outcomes of a feature sum to a redundant direction (the probabilities are a simplex), which the whitening step below removes, so no reference outcome is dropped by hand.

Ordered by the axis, the standardised running sum of these scores is the empirical fluctuation process $B(t)$, $t in [0,1]$. Standardisation is by the inverse square root of the outer-product-of-gradients covariance of the focal block, which decorrelates its dimensions; the process is pinned to zero at both ends and, under stability, converges to a Brownian bridge. Two functionals read the process: $text{maxLM} = max_t lVert B(t) rVert^2$ (power against an abrupt break, with $argmax_t$ the estimated break position) and $text{CvM} = int_0^1 lVert B(t) rVert^2,mathrm{d}t$ (power against gradual drift). The null is drawn by simulating many $d$-dimensional Brownian bridges on the axis’s own time grid and reading the same functional off each, so the $p$-value is analytic rather than a refit.

The module is a pure consumer of a cached measurement-only fit and a per-proband axis, method independent, mirroring the cheap half of analysis.drift. It does not refit the mixture.

class analysis.invariance.FocalParameter(cls, feature, kind, column, local, outcome=None)

One class-conditional location parameter, the unit a casewise score is taken for.

cls

The latent class the parameter belongs to.

Type:

int

feature

The feature (measurement column) the parameter describes.

Type:

str

kind

The emission type: "gaussian_mean", "bernoulli" or "multinomial_logit".

Type:

str

column

The feature’s column index in the measurement matrix (for reading $x_{ij}$).

Type:

int

local

The feature’s index within its emission sub-model (for reading the fitted parameter).

Type:

int

outcome

For a categorical feature, the outcome $l$ this score is taken for; None otherwise.

Type:

int or None

class analysis.invariance.CasewiseScores(values, parameters, probands, responsibilities)

The casewise focal scores for a fit: one column per focal parameter, one row per proband.

values

The score matrix, shape (n_probands, n_focal).

Type:

numpy.ndarray

parameters

The focal parameter each column scores, in column order.

Type:

list of FocalParameter

probands

The proband ids, in row order (the measurement matrix index).

Type:

pandas.Index

responsibilities

The posterior class responsibilities, shape (n_probands, n_classes).

Type:

numpy.ndarray

analysis.invariance.per_sample_log_likelihood(model, x_values)

Return the marginal log-likelihood of each proband under a measurement-only fit.

The per-sample log-likelihood is $log sum_k pi_k prod_j f_j(x_{ij};theta_{jk})$, the log-sum-exp over classes of the class prior plus the measurement emission. This is the quantity the casewise scores are gradients of, and the finite-difference gate differentiates.

Parameters:

• model (StepMix) – A fitted measurement-only estimator.

• x_values (numpy.ndarray) – The measurement matrix, shape (n_probands, n_features), in the fit’s column order.

Returns:

The per-proband log-likelihood, shape (n_probands,).

Return type:

numpy.ndarray

analysis.invariance.responsibilities(model, x_values)

Return the posterior class responsibilities $r_{ik}$ for a measurement-only fit.

Parameters:

• model (StepMix) – A fitted measurement-only estimator.

• x_values (numpy.ndarray) – The measurement matrix, in the fit’s column order.

Returns:

The responsibilities, shape (n_probands, n_classes), each row summing to one.

Return type:

numpy.ndarray

analysis.invariance.casewise_scores(model, measurement_data, typing)

Compute the casewise focal score of every class-conditional location parameter.

For each class and each feature the analytic score of the fitted location parameter is formed by Fisher’s identity: the responsibility times the emission gradient. Gaussian means and Bernoulli probabilities give one score column per feature; a categorical feature gives one column per observed outcome (the never-observed padding columns are skipped).

Parameters:

• model (StepMix) – A fitted measurement-only estimator.

• measurement_data (pandas.DataFrame) – The measurement matrix the fit was estimated on, in its column order.

• typing (analysis.features.Typing) – The feature typing that assigns each feature its emission.

Returns:

The score matrix, the focal parameter per column, the proband index, and the responsibilities.

Return type:

CasewiseScores

analysis.invariance.numerical_score(model, measurement_data, parameter, *, eps=1e-05)

Return the finite-difference casewise score of one focal parameter (the correctness gate).

Perturbs the single parameter by plus and minus eps, recomputes the per-sample log-likelihood each way, restores the fit, and returns the central difference. A Gaussian mean and a Bernoulli probability are perturbed directly; a categorical outcome is perturbed on the logit scale (the outcome probabilities are renormalised by a softmax), which matches the analytic multinomial-logit score. The analytic casewise_scores() must match this to machine precision per proband; a mismatch means a wrong gradient, which would invalidate every $p$-value.

Parameters:

• model (StepMix) – A fitted measurement-only estimator. Restored to its fitted parameters on return.

• measurement_data (pandas.DataFrame) – The measurement matrix, in the fit’s column order.

• parameter (FocalParameter) – The focal parameter to differentiate.

• eps (float, optional) – The central-difference step.

Returns:

The finite-difference score, shape (n_probands,).

Return type:

numpy.ndarray

class analysis.invariance.TimeGrid(t, dt, index, positions)

The ordinal time grid the fluctuation process is read on.

t

The cumulative sample fraction at each evaluation point, ending at one.

Type:

numpy.ndarray

dt

The spacing between consecutive points (t differenced, with a leading t[0]), the integration weight for the Cramer-von Mises functional and the increment variance for the simulated null.

Type:

numpy.ndarray

index

The proband count included at each evaluation point (the cumulative-sum row to read).

Type:

numpy.ndarray

positions

The axis value at each evaluation point, for reading a break off the process.

Type:

numpy.ndarray

analysis.invariance.build_time_grid(sorted_axis, *, max_points=512)

Build the evaluation grid from the axis-sorted probands, collapsing ties.

The process changes value proband by proband, but within a run of tied axis values the proband order is arbitrary, so the process is only read at the end of each tied run (an order-invariant point). When there are more such points than max_points the grid is thinned to roughly equal ordinal spacing; the observed statistic and the simulated null use the same thinned grid, so the reading stays calibrated.

Parameters:

• sorted_axis (numpy.ndarray) – The axis values of the covered probands, in ascending order.

• max_points (int, optional) – The largest grid the process is evaluated on.

Returns:

The evaluation grid.

Return type:

TimeGrid

analysis.invariance.fluctuation_process(block_values, grid)

Return the standardised, bridge-tied fluctuation process on the grid.

The scores are cumulated in axis order, standardised by the whitening, scaled by the root of the sample size, and tied down linearly so the process is zero at both ends (the residual tilt from an approximate optimiser or from axis coverage below one is removed). Under stability the result converges to a Brownian bridge.

Parameters:

• block_values (numpy.ndarray) – The focal scores of one block, in axis order, shape (n_probands, d).

• grid (TimeGrid) – The evaluation grid.

Returns:

The process, shape (n_grid, d_eff).

Return type:

numpy.ndarray

analysis.invariance.bridge_functionals(process, grid)

Return the maxLM and Cramer-von Mises functionals of a process and the break index.

maxLM is the largest squared norm over the grid, powerful against an abrupt break, and its location is the estimated break. CvM is the grid integral of the squared norm, powerful against a gradual drift.

Parameters:

• process (numpy.ndarray) – The fluctuation process, shape (n_grid, d_eff).

• grid (TimeGrid) – The evaluation grid.

Returns:

maxLM, CvM, and the grid index of the maximum.

Return type:

tuple

analysis.invariance.simulate_bridge_null(d_eff, grid, *, n_sim, seed, chunk=256)

Simulate the null distributions of the two functionals for a d_eff-dimensional bridge.

Draws independent d_eff-dimensional Brownian bridges on the grid and reads the same functionals off each, so the null matches the observed statistic’s grid exactly. This is a Gaussian-process simulation, not a model refit. The draws run in chunks to bound memory.

Parameters:

• d_eff (int) – The effective dimension of the focal block.

• grid (TimeGrid) – The evaluation grid.

• n_sim (int) – The number of simulated bridges.

• seed (int) – The seed for the draw.

• chunk (int, optional) – The number of bridges simulated at once.

Returns:

"maxLM" and "cvm", each an array of n_sim null draws.

Return type:

dict

analysis.invariance.simulate_bridge_band(d_eff, grid, *, n_sim, seed, quantiles=(0.5, 0.95), chunk=256)

Return pointwise quantiles of a simulated bridge’s squared norm, the figure’s null band.

Unlike simulate_bridge_null(), which keeps only the maxLM and CvM of each draw, this keeps the squared norm at every grid point, so the observed process can be drawn against a pointwise null envelope.

Parameters:

• d_eff (int) – The effective dimension of the block.

• grid (TimeGrid) – The evaluation grid.

• n_sim (int) – The number of simulated bridges.

• seed (int) – The seed for the draw.

• quantiles (collections.abc.Sequence of float, optional) – The quantiles of the squared norm to return per grid point.

• chunk (int, optional) – The number of bridges simulated at once.

Returns:

A mapping from each quantile to its per-grid-point curve, shape (n_grid,).

Return type:

dict

analysis.invariance.add_one_pvalue(observed, null_draws)

Return the simulation $p$-value with the Phipson-Smyth add-one correction.

The smallest attainable $p$-value is $1/(B+1)$ rather than zero, so a statistic beyond every draw is not reported as impossible.

analysis.invariance.directional_slopes(block_values, axis_values)

Return each focal parameter’s ordinary-least-squares slope of score against the axis.

A systematic non-zero slope is directional drift: the parameter’s score trends with the ordering variable rather than fluctuating around zero. The sign carries the direction.

Parameters:

• block_values (numpy.ndarray) – The focal scores of one block, shape (n_probands, d).

• axis_values (numpy.ndarray) – The axis value of each proband, in the same row order.

Returns:

The per-parameter slope, shape (d,).

Return type:

numpy.ndarray

class analysis.invariance.BlockResult(label, cls, category, n, d, d_eff, max_lm, cvm, p_max_lm, p_cvm, break_position, break_low, break_high, direction, direction_feature, reject_max_lm=False, reject_cvm=False)

The invariance reading for one focal block.

label

The block identifier (a class, or a class crossed with a category).

Type:

str

cls

The class the block belongs to.

Type:

int

category

The feature category the block is restricted to, or None for a whole-class block.

Type:

str or None

n

The number of probands the process ran over.

Type:

int

d

The number of focal parameters in the block.

Type:

int

d_eff

The effective dimension after dropping redundant directions.

Type:

int

max_lm

The maxLM statistic.

Type:

float

cvm

The Cramer-von Mises statistic.

Type:

float

p_max_lm

The bridge $p$-value of maxLM.

Type:

float

p_cvm

The bridge $p$-value of CvM.

Type:

float

break_position

The axis value at the maxLM maximum, the estimated break.

Type:

float

break_low

The lower edge of the break confidence set.

Type:

float

break_high

The upper edge of the break confidence set.

Type:

float

direction

The signed directional slope of the block’s strongest-trending parameter.

Type:

float

direction_feature

The feature of that strongest-trending parameter.

Type:

str

reject_max_lm

Whether maxLM passes the Benjamini-Hochberg step across blocks (filled by the caller).

Type:

bool

reject_cvm

Whether CvM passes the Benjamini-Hochberg step across blocks (filled by the caller).

Type:

bool

analysis.invariance.test_block(block_values, parameters, sorted_axis, order, grid, *, label, cls, category, n_sim, seed)

Run the fluctuation-process test on one focal block.

Cumulates the block’s scores in axis order, reads the two functionals, draws the analytic null, and estimates the break with a confidence set from the sup statistic. Pure and cheap: no fitting, only the stored scores.

Parameters:

• block_values (numpy.ndarray) – The block’s focal scores in the fit’s row order, shape (n_probands, d).

• parameters (list of FocalParameter) – The focal parameter of each column, for the directional read.

• sorted_axis (numpy.ndarray) – The covered probands’ axis values, ascending.

• order (numpy.ndarray) – The permutation that sorts the probands by the axis.

• grid (TimeGrid) – The evaluation grid for this axis.

• label (str, int, str or None) – The block’s identity.

• cls (str, int, str or None) – The block’s identity.

• category (str, int, str or None) – The block’s identity.

• n_sim (int) – The number of simulated bridges for the null.

• seed (int) – The seed for the null draw.

Returns:

The block’s statistics, $p$-values, break, and directional slope.

Return type:

BlockResult

analysis.invariance.focal_blocks(parameters, category_map, *, by_category)

Group focal-parameter column indices into the blocks the test runs over.

The whole-class blocks pool every parameter of a class. The category blocks restrict a class to one feature category, so a drift can be localised to a class and a symptom domain.

Parameters:

• parameters (list of FocalParameter) – The focal parameter of each score column.

• category_map (Mapping or None) – The feature-to-category mapping; only needed when by_category is set.

• by_category (bool) – Build class-by-category blocks in addition to whole-class blocks.

Returns:

A mapping from (class, category) (category None for a whole-class block) to the list of column indices in that block.

Return type:

dict

class analysis.invariance.TopProcess(label, t, positions, observed, null_q50, null_q95)

The empirical fluctuation process and null band of the most significant block.

Kept so the exploratory figure has a real curve to draw without recomputing the scores.

label

The block the process belongs to.

Type:

str

t

The grid’s cumulative sample fraction.

Type:

numpy.ndarray

positions

The axis value at each grid point.

Type:

numpy.ndarray

observed

The observed squared norm $\lVert B(t) \rVert^2$ at each grid point.

Type:

numpy.ndarray

null_q50

The pointwise median of the simulated bridge’s squared norm.

Type:

numpy.ndarray

null_q95

The pointwise 95th percentile of the simulated bridge’s squared norm.

Type:

numpy.ndarray

class analysis.invariance.InvarianceResult(blocks, n_reference, n_covered, coverage, axis, n_sim, top_process=None, q=0.05)

The invariance reading for a fit against one axis.

blocks

The per-block results, Benjamini-Hochberg decisions filled in.

Type:

list of BlockResult

n_reference

The number of probands in the reference fit.

Type:

int

n_covered

The number with a non-missing axis value, which the test ran over.

Type:

int

coverage

n_covered / n_reference.

Type:

float

axis

The ordering variable.

Type:

str

n_sim

The number of simulated bridges per block.

Type:

int

top_process

The fluctuation process of the most significant block, for the figure.

Type:

TopProcess or None

q

The false-discovery-rate level of the Benjamini-Hochberg step.

Type:

float

analysis.invariance.benjamini_hochberg(p_values, q=0.05)

Return a boolean mask of the hypotheses that pass Benjamini-Hochberg FDR control.

A hypothesis is rejected if its $p$-value is at or below the largest threshold $q,text{rank}/m$ it satisfies, where $m$ is the number of finite $p$-values. Not-a-number entries never pass. This mirrors analysis.drift.benjamini_hochberg(), the repo convention for the strata-by-class tests.

analysis.invariance.run_invariance(model, measurement_data, typing, axis_values, *, axis, category_map=None, by_category=True, n_sim=2000, seed=0, max_grid=512, q=0.05)

Run the score-based invariance test of a fit against an axis, over every focal block.

Joins the axis to the reference probands, drops those with a missing axis value (reporting coverage), computes the casewise scores once, and runs the fluctuation-process test per block. Significance is Benjamini-Hochberg controlled across blocks within the axis.

Parameters:

• model (StepMix) – The fitted measurement-only reference estimator.

• measurement_data (pandas.DataFrame) – The measurement matrix the fit was estimated on.

• typing (analysis.features.Typing) – The feature typing.

• axis_values (pandas.Series) – The per-proband axis value, indexed by proband id.

• axis (str) – The axis name, recorded on the result.

• category_map (Mapping or None, optional) – The feature-to-category mapping for the class-by-category blocks.

• by_category (bool, optional) – Add the class-by-category blocks.

• n_sim (int, optional) – The number of simulated bridges per block.

• seed (int, optional) – The base seed; each block’s null uses a distinct derived seed.

• max_grid (int, optional) – The largest evaluation grid (ties collapsed, then thinned to this many points).

• q (float, optional) – The false-discovery-rate level.

Returns:

The per-block results with FDR decisions, and the axis coverage.

Return type:

InvarianceResult

Local class-profile displacement along an axis, from the single cached fit (plan section 7e).

The score-based invariance test (analysis.invariance) is saturated at the cohort’s sample size: exact measurement invariance is always rejected once the sample runs to the thousands, so the bridge $p$-value cannot discriminate. This module recasts the same question around a null-free effect size. It freezes the pooled responsibilities of the measurement-only reference and reads how each class centroid moves as a smooth function of the axis, with the uncertainty coming from a clustered bootstrap rather than a saturated analytic null. It refits nothing.

The quantity. For the frozen responsibilities $r_{ik}$ (the pooled predict_proba), the local centroid of class $k$ at focal point $f$ is the kernel-and-responsibility-weighted mean

\[\mu_k(f) = \frac{\sum_i w_i(f)\,r_{ik}\,x_i}{\sum_i w_i(f)\,r_{ik}},\]

with $w_i(f)$ the Gaussian kernel weight of proband $i$’s axis value about $f$ (the analysis.localise.gaussian_weights() window at the axis’s chosen bandwidth). With $w equiv 1$ this is the pooled centroid $mu_k$, which equals the fit’s responsibility-weighted class means. The primitive is the per-feature displacement $d_k(f) = mu_k(f) - mu_k$, kept full-dimensional. Magnitudes are divided by the between-class separation (the mean pairwise distance between distinct pooled centroids under the same full standardised-Euclidean norm, separation()), so one separation unit is the mean inter-class gap and a displacement reads as a genuine fraction of that gap, comparable across axes.

Uncertainty is a clustered bootstrap: families are resampled with replacement and the local centroids are recomputed on the resample (re-weighting only, the responsibilities stay frozen), giving a per-focal-point envelope. Resampling families rather than probands respects the within-family correlation, so a block of correlated features carries an honest, wider tube. Everything the module returns is conditional on the pooled fit.

The 2D discriminant plane (analysis.trajectory) is a view of the full-dimensional displacement, not the authority: capture_fraction() reports how much of a class’s displacement lies in that plane, so a drift that is mostly out of plane cannot be hidden by the picture.

analysis.trajectory_local.local_centroids(x_values, responsibilities, weights)

Return the kernel-and-responsibility-weighted class centroids at one focal point.

The centroid of class $k$ is $sum_i w_i r_{ik} x_i / sum_i w_i r_{ik}$, the local weighted mean of the feature matrix under the frozen responsibilities. With weights a Gaussian kernel window this is the local centroid $mu_k(f)$; with unit weights it is the pooled centroid $mu_k$ (pooled_centroids()).

Parameters:

• x_values (numpy.ndarray) – The measurement matrix, shape (n_probands, n_features).

• responsibilities (numpy.ndarray) – The frozen posterior responsibilities $r_{ik}$, shape (n_probands, n_classes).

• weights (numpy.ndarray) – The per-proband kernel weight, shape (n_probands,).

Returns:

The class centroids, shape (n_classes, n_features); a class with no local weight is all not-a-number.

Return type:

numpy.ndarray

analysis.trajectory_local.pooled_centroids(x_values, responsibilities)

Return the pooled (whole-cohort) responsibility-weighted class centroids.

The local centroids with a unit weight on every proband, so this is the frozen-responsibility class mean each local centroid is measured against.

analysis.trajectory_local.separation(reference)

Return the between-class separation, the drift baseline (delegated to analysis.drift).

The full (unaveraged) standardised-Euclidean distance (analysis.drift.FullStandardisedEuclidean), the same sum-norm convention grain_magnitude() uses for a class’s displacement. Numerator and denominator then share a scale, so a separation-scaled magnitude is a genuine fraction of the mean inter-class gap: one separation unit is the mean pairwise distance between distinct reference centroids. The refit-based drift stage keeps its own averaged convention, self-consistent within that stage.

analysis.trajectory_local.grain_magnitude(displacement, pooled_sd, columns, separation_scale)

Return the separation-scaled displacement magnitude of a feature grain, per class.

A grain is a set of feature columns (the 4 classes are read whole, and each class is also read within each of the 7 author categories). Its magnitude is the Euclidean norm of the per-feature displacement in pooled-standard-deviation units over the grain’s features, $lVert d_k / sigma rVert$, divided by the between-class separation. The raw (unaveraged) norm is deliberate: a larger grain carries a larger norm and a wider bootstrap tube, and the tube, not the bare magnitude, calls significance.

Parameters:

• displacement (numpy.ndarray) – The per-feature displacement, shape (..., n_features); the leading axes are kept.

• pooled_sd (numpy.ndarray) – The per-feature pooled standard deviation, shape (n_features,).

• columns (numpy.ndarray) – The integer column indices of the grain’s features.

• separation_scale (float) – The between-class separation the magnitude is divided by.

Returns:

The separation-scaled magnitude, shape displacement.shape[:-1].

Return type:

numpy.ndarray

analysis.trajectory_local.discriminant_plane(embedding)

Return an orthonormal basis of the first two discriminant directions, in standardised space.

The embedding maps a standardised feature vector to the discriminant axes by the linear scalings_; the plane the trajectory figure draws is the span of the first two of those direction vectors. Orthonormalising that span (a thin QR) gives a projector onto the plane, so the in-plane part of a displacement can be measured honestly even though the raw scaling vectors are not orthogonal.

Returns:

The orthonormal basis, shape (n_features, 2).

Return type:

numpy.ndarray

analysis.trajectory_local.capture_fraction(displacement_row, pooled_sd, plane)

Return the fraction of one class’s displacement that lies in the discriminant plane.

The honesty guard on the figure: $lVert P,d_k rVert / lVert d_k rVert$, where $d_k$ is the class’s standardised displacement and $P$ orthogonally projects onto the plane. A value near one means the 2D picture shows essentially all of the movement; a value near zero means the movement is mostly out of plane and the picture understates it.

Parameters:

• displacement_row (numpy.ndarray) – One class’s per-feature displacement, shape (n_features,).

• pooled_sd (numpy.ndarray) – The per-feature pooled standard deviation, shape (n_features,).

• plane (numpy.ndarray) – An orthonormal basis of the plane, shape (n_features, 2) (discriminant_plane()).

Returns:

The in-plane capture fraction, or not-a-number when the displacement is zero.

Return type:

float

analysis.trajectory_local.mahalanobis_magnitude(displacement_row, precision)

Return the Mahalanobis magnitude of one class’s raw displacement.

The covariance-aware corroborating magnitude, $sqrt{d_k^top Sigma^{-1} d_k}$ with the Ledoit-Wolf-shrunk pooled within-class precision, so a coordinated shift across correlated features counts once. This is the analysis.drift Mahalanobis distance evaluated on the local displacement; bootstrap-calibrating it (its clustered-bootstrap band) makes it dimension-fair across grains of different size.

Parameters:

• displacement_row (numpy.ndarray) – One class’s raw per-feature displacement, shape (n_features,).

• precision (numpy.ndarray) – The pooled within-class precision, shape (n_features, n_features).

class analysis.trajectory_local.ObservedTrajectory(focal_points, pooled, displacement, ld, grain_magnitude, mahalanobis, capture, focal_ref, peak_focal)

The observed local-displacement trajectory of a fit against one axis.

focal_points

The axis positions the local centroids were read at, shape (n_focal,).

Type:

numpy.ndarray

pooled

The pooled class centroids, shape (n_classes, n_features).

Type:

numpy.ndarray

displacement

The per-feature displacement $d_k(f)$, shape (n_classes, n_focal, n_features).

Type:

numpy.ndarray

ld

The local centroids’ first two discriminant coordinates, shape (n_classes, n_focal, 2), the trajectory the plane figure draws.

Type:

numpy.ndarray

grain_magnitude

Per grain, the separation-scaled magnitude, shape (n_classes, n_focal).

Type:

dict of str to numpy.ndarray

mahalanobis

The whole-class Mahalanobis magnitude, shape (n_classes, n_focal).

Type:

numpy.ndarray

capture

The per-class in-plane capture fraction of the endpoint displacement, shape (n_classes,).

Type:

numpy.ndarray

focal_ref

The focal index the capture fraction and per-feature inference are anchored at: the endpoint (the last focal point), pre-specified so the per-feature test is not selected on the observed magnitude. The endpoint carries the accumulated drift for a monotone or single-break axis.

Type:

int

peak_focal

Per class, the focal index of the largest whole-class magnitude, reported as where the drift is strongest (informational, not a test anchor), shape (n_classes,).

Type:

numpy.ndarray

analysis.trajectory_local.observed_trajectory(x_values, responsibilities, axis_values, focal_points, bandwidth, *, pooled_sd, separation_scale, grains, embedding, precision, plane)

Compute the observed local-displacement trajectory of a fit against an axis.

Reads the local centroids at each focal point under the frozen responsibilities, forms the per-feature displacement from the pooled centroid, and derives the discriminant-plane coordinates, the separation-scaled grain magnitudes, the Mahalanobis magnitude, and the endpoint capture fraction. Pure and cheap: no fitting, only re-weighting.

Parameters:

• x_values (numpy.ndarray) – The measurement matrix, the frozen responsibilities, and the per-proband axis value.

• responsibilities (numpy.ndarray) – The measurement matrix, the frozen responsibilities, and the per-proband axis value.

• axis_values (numpy.ndarray) – The measurement matrix, the frozen responsibilities, and the per-proband axis value.

• focal_points (numpy.ndarray) – The axis positions to read local centroids at.

• bandwidth (float) – The Gaussian kernel bandwidth, in axis units.

• pooled_sd (numpy.ndarray) – The per-feature pooled standard deviation.

• separation_scale (float) – The between-class separation the magnitudes are divided by.

• grains (dict of str to numpy.ndarray) – The feature-column indices of each grain (whole-class plus per author category).

• embedding (analysis.trajectory.Embedding) – The fixed discriminant embedding of the pooled classes.

• precision (numpy.ndarray) – The pooled within-class precision, for the Mahalanobis magnitude.

• plane (numpy.ndarray) – The orthonormal discriminant plane, for the capture fraction.

Returns:

The observed trajectory and its derived magnitudes.

Return type:

ObservedTrajectory

class analysis.trajectory_local.BootstrapTube(quantiles, ld, grain_bands, mahalanobis_bands, feature_displacement, n_boot, clustered, signed_slope=None, net_trend=None, signed_trajectory=None, break_position=None)

The clustered-bootstrap envelope of a displacement trajectory.

quantiles

The bootstrap quantiles held in each band, in order (typically low, median, high).

Type:

tuple of float

ld

The local centroids’ discriminant coordinates over the replicates, shape (n_boot, n_classes, n_focal, 2), the centroid tube the plane figure draws.

Type:

numpy.ndarray

grain_bands

Per grain, the quantile bands of the separation-scaled magnitude, shape (n_quantiles, n_classes, n_focal).

Type:

dict of str to numpy.ndarray

mahalanobis_bands

The quantile bands of the whole-class Mahalanobis magnitude, shape (n_quantiles, n_classes, n_focal).

Type:

numpy.ndarray

feature_displacement

The standardised per-feature displacement at the reference (endpoint) focal point, over the replicates, shape (n_boot, n_classes, n_features).

Type:

numpy.ndarray

n_boot

The number of bootstrap replicates.

Type:

int

clustered

Whether families (True) or probands (False) were resampled.

Type:

bool

signed_slope

The directional draws: per replicate, each class’s signed net-projected slope (the directional statistic of directional_statistic()), shape (n_boot, n_classes). None when the tube was built without frozen net directions.

Type:

numpy.ndarray or None

net_trend

Per replicate, each class’s separation-scaled net-trend displacement over the focal span, shape (n_boot, n_classes); None as above.

Type:

numpy.ndarray or None

signed_trajectory

Per replicate, each class’s one-dimensional signed trajectory projected onto its frozen net direction, shape (n_boot, n_classes, n_focal); None as above. The band of this is what the directional figure draws.

Type:

numpy.ndarray or None

break_position

Per replicate, each class’s single-break location on the signed trajectory, shape (n_boot, n_classes); None as above.

Type:

numpy.ndarray or None

analysis.trajectory_local.clustered_bootstrap_tube(x_values, responsibilities, axis_values, families, focal_points, bandwidth, *, pooled_sd, separation_scale, grains, embedding, precision, focal_ref, n_boot, seed, clustered=True, quantiles=(2.5, 50.0, 97.5), net_directions=None)

Bootstrap the displacement trajectory by resampling families (or probands).

Each replicate resamples whole families with replacement (so a proband appearing twice contributes twice), recomputes the pooled and local centroids on the resample under the frozen responsibilities, and records the discriminant coordinates, the separation-scaled grain magnitudes, the Mahalanobis magnitude, and the per-feature displacement at each class’s reported focal point. The per-focal-point quantiles of those are the tube. Setting clustered=False resamples individual probands instead, the independent-bootstrap comparison that shows the family clustering is real rather than cosmetic.

When net_directions is given (each class’s frozen unit direction from directional_statistic()), each replicate also records the directional draws: the signed net-projected slope, the separation-scaled net-trend displacement, the one-dimensional signed trajectory, and the single-break location. Their per-class spread is the clustered-bootstrap null the H0E directional test reads against; freezing the direction at the observed value keeps the projected slope a fixed linear functional, so it is signed and its interval can honestly cover zero.

Parameters:

• x_values (numpy.ndarray) – The measurement matrix, frozen responsibilities, and per-proband axis value.

• responsibilities (numpy.ndarray) – The measurement matrix, frozen responsibilities, and per-proband axis value.

• axis_values (numpy.ndarray) – The measurement matrix, frozen responsibilities, and per-proband axis value.

• families (numpy.ndarray) – The per-proband family identifier, shape (n_probands,); the clustering unit.

• focal_points (numpy.ndarray) – The axis positions to read local centroids at.

• bandwidth (float) – The Gaussian kernel bandwidth.

• pooled_sd (ndarray) – As described for observed_trajectory().

• separation_scale (float) – As described for observed_trajectory().

• grains (dict[str, ndarray]) – As described for observed_trajectory().

• embedding (Embedding) – As described for observed_trajectory().

• precision (ndarray) – As described for observed_trajectory().

• focal_ref (int) – The focal index the per-feature displacement is recorded at (the endpoint).

• n_boot (int) – The number of bootstrap replicates.

• seed (int) – The base seed for the resampling.

• clustered (bool, optional) – Resample families (default) or individual probands.

• quantiles (tuple of float, optional) – The bootstrap quantiles kept in each band.

• net_directions (numpy.ndarray, optional) – Each class’s frozen unit net direction, shape (n_classes, n_features). When given, the directional draws are recorded; when None they are left off the tube.

Returns:

The bootstrap replicates and their per-focal-point quantile bands.

Return type:

BootstrapTube

class analysis.trajectory_local.FeatureInference(displacement, ci_low, ci_high, p_value, reject, covers_zero)

The per-feature displacement, its clustered-bootstrap interval, and the FDR decision.

displacement

The observed standardised per-feature displacement at each class’s reported focal point, shape (n_classes, n_features).

Type:

numpy.ndarray

ci_low, ci_high

The bootstrap interval per feature, shape (n_classes, n_features).

Type:

numpy.ndarray

p_value

The two-sided bootstrap $p$-value that the displacement differs from zero, shape (n_classes, n_features).

Type:

numpy.ndarray

reject

The Benjamini-Hochberg decision across the n_classes * n_features tests, shape (n_classes, n_features).

Type:

numpy.ndarray

covers_zero

Whether the bootstrap interval covers zero, shape (n_classes, n_features); most being true is the readable “many features invariant”.

Type:

numpy.ndarray

analysis.trajectory_local.per_feature_inference(observed_displacement, feature_draws, *, q=0.05)

Test each per-feature displacement against zero with a clustered-bootstrap interval and FDR.

Reads the 95 per cent bootstrap interval of each (class, feature) displacement, forms a two-sided bootstrap $p$-value from the fraction of replicates on the far side of zero, and applies Benjamini-Hochberg control across the 4 * n_features tests (the analysis.invariance.benjamini_hochberg() implementation, the repo convention). A displacement whose interval covers zero is invariant at this level; most covering zero is the “many features invariant” reading.

Parameters:

• observed_displacement (numpy.ndarray) – The observed standardised per-feature displacement, shape (n_classes, n_features).

• feature_draws (numpy.ndarray) – The bootstrap replicates, shape (n_boot, n_classes, n_features) (BootstrapTube.feature_displacement).

• q (float, optional) – The false-discovery-rate level.

Returns:

The per-feature displacement, interval, $p$-value, and FDR decision.

Return type:

FeatureInference

class analysis.trajectory_local.ControlComparison(axis_magnitude, control_magnitude, difference, diff_draws, p_value, p_value_greater, n_boot)

A paired-bootstrap comparison of a timing axis against one control variable.

axis_magnitude, control_magnitude

The observed, class-averaged, separation-scaled endpoint magnitude of the timing axis and of the control variable (household income, area deprivation, or a random ordering).

Type:

float

difference

axis_magnitude - control_magnitude.

Type:

float

diff_draws

The paired-bootstrap replicate differences, shape (n_boot,).

Type:

numpy.ndarray

p_value

The two-sided bootstrap $p$-value that the difference is zero, floored at $1/(n_{text{boot}} + 1)$.

Type:

float

p_value_greater

The one-sided bootstrap $p$-value that the axis magnitude exceeds the control’s, floored the same way.

Type:

float

n_boot

The number of paired bootstrap replicates.

Type:

int

analysis.trajectory_local.control_specificity_bootstrap(x_values, responsibilities, families, axis_values, axis_bandwidth, axis_focal, control_values, control_bandwidth, control_focal, *, pooled_sd, separation_scale, n_boot, seed)

Paired family-bootstrap test that a timing axis’s drift exceeds a control’s.

The specificity panel (the invariance-as-an-effect-size guide) reads the timing axis’s endpoint magnitude as larger than a control’s, but as a magnitude comparison only, because the axis and the control were each read from one point estimate. This adds a $p$-value: every bootstrap replicate resamples one set of families and recomputes both the axis and the control magnitude on that same resample, so the two quantities share their sampling variation and the difference is a genuine paired statistic, not the comparison of two separately noisy numbers. The observed difference then acts as its own bootstrap-inverted test: \(p\) is the fraction of replicate differences on the far side of zero (doubled for the two-sided form), the same construction per_feature_inference() and directional_inference() use.

x_values, responsibilities, and families must already be restricted to the probands finite on both axis_values and control_values, so that a family resampled for one quantity is resampled for the other. axis_focal and control_focal are each variable’s own endpoint focal position (its own bandwidth and grid), matching how the specificity panel reads each variable.

Parameters:

• x_values (numpy.ndarray) – The measurement matrix and frozen responsibilities, restricted to the shared rows.

• responsibilities (numpy.ndarray) – The measurement matrix and frozen responsibilities, restricted to the shared rows.

• families (numpy.ndarray) – The per-proband family identifier over the same rows, the clustering unit.

• axis_values (numpy.ndarray) – The timing axis and the control variable, over the same rows.

• control_values (numpy.ndarray) – The timing axis and the control variable, over the same rows.

• axis_bandwidth (float) – Each variable’s own Gaussian kernel bandwidth.

• control_bandwidth (float) – Each variable’s own Gaussian kernel bandwidth.

• axis_focal (float) – Each variable’s own endpoint focal position.

• control_focal (float) – Each variable’s own endpoint focal position.

• pooled_sd (numpy.ndarray) – The per-feature pooled standard deviation.

• separation_scale (float) – The between-class separation the magnitude is divided by.

• n_boot (int) – The number of paired bootstrap replicates.

• seed (int) – The bootstrap seed.

Returns:

The observed magnitudes, their difference, and its bootstrap $p$-values.

Return type:

ControlComparison

analysis.trajectory_local.category_grains(columns, category_map)

Return the column indices of each presentation grain: whole-class and per author category.

The whole-class grain ("class") is every feature; each author category grain ("category:<name>") is the features mapped to that category. A category with no present feature is omitted. The grains are the pre-specified aggregation levels, fixed before the data are seen.

Parameters:

• columns (list of str) – The measurement-matrix feature columns, in order.

• category_map (dict of str to str) – The feature-to-category mapping.

Returns:

The grain name mapped to its integer column indices.

Return type:

dict of str to numpy.ndarray

analysis.trajectory_local.referent_grains(columns, instrument_map, referent_map)

Return the column indices of each referent grain: per instrument and per temporal referent.

Each feature carries the instrument it comes from (instrument_map, derived from the data dictionary) and each instrument carries a pre-registered temporal referent (referent_map, analysis.features.INSTRUMENT_REFERENT). The grains are the per-instrument column sets ("instrument:<name>", the transparent underlay) and the per-referent column sets ("referent:<name>", the two-way headline). Resolution fails loudly, mirroring analysis.features.reconcile()’s no-typing-signal guard: a feature with no instrument, or an instrument with no referent, raises rather than being dropped, so a mapping gap cannot pass silently as an empty grain.

Parameters:

• columns (list of str) – The measurement-matrix feature columns, in order.

• instrument_map (dict of str to str) – The feature-to-instrument mapping.

• referent_map (dict of str to str) – The instrument-to-referent mapping.

Returns:

The grain name mapped to its integer column indices.

Return type:

dict of str to numpy.ndarray

Raises:

ValueError – When a feature resolves to no instrument, or its instrument to no referent.

analysis.trajectory_local.slope_vectors(traj_std, focal_points)

Return each class’s per-feature ordinary-least-squares slope against the axis.

For the standardised displacement trajectory $D_k(f) = d_k(f) / sigma$, the slope of feature $m$ is $sum_j (f_j - bar f),D_k(f_j)[m] / sum_j (f_j - bar f)^2$, the univariate least-squares slope of that feature’s displacement on the axis position. The focal grid is evenly spaced in axis units (analysis.localise.focal_grid()), so an equal weight per focal point is an honest per-axis-unit trend on an irregularly sampled axis. A class is regressed only over the focal points where its local centroid is defined; a class defined at fewer than two focal points has an all-not-a-number slope.

Parameters:

• traj_std (numpy.ndarray) – The standardised displacement trajectory, shape (n_classes, n_focal, n_features).

• focal_points (numpy.ndarray) – The axis positions the trajectory was read at, shape (n_focal,).

Returns:

The per-class slope vector $b_k$, shape (n_classes, n_features).

Return type:

numpy.ndarray

analysis.trajectory_local.net_directions(traj_std)

Return each class’s unit net direction, the direction of its mean displacement.

The net direction $hat u_k$ is the unit vector of the mean standardised displacement across the focal grid, $overline{D_k} / lVert overline{D_k} rVert$. It is the axis the signed directional statistic projects onto. A class whose mean displacement is negligible (a symmetric excursion that cancels, or no drift) has an ill-defined direction and is given the zero vector, so its projected slope is zero rather than a projection onto noise.

Parameters:

traj_std (numpy.ndarray) – The standardised displacement trajectory, shape (n_classes, n_focal, n_features).

Returns:

The per-class unit net direction, shape (n_classes, n_features).

Return type:

numpy.ndarray

analysis.trajectory_local.project_onto(traj_std, directions)

Return each class’s one-dimensional signed trajectory along its net direction.

The projection $s_k(f) = langle D_k(f), hat u_k rangle$ of the standardised displacement onto the class’s frozen net direction, a signed scalar per focal point. This is the one-dimensional signed trajectory the directional figure draws and the changepoint read localises a break on; positive values sit on the net-drift side of the pooled centroid.

Parameters:

• traj_std (numpy.ndarray) – The standardised displacement trajectory, shape (n_classes, n_focal, n_features).

• directions (numpy.ndarray) – The per-class unit net direction, shape (n_classes, n_features).

Returns:

The signed trajectory $s_k(f)$, shape (n_classes, n_focal).

Return type:

numpy.ndarray

analysis.trajectory_local.single_break(positions, series, *, min_segment=3)

Return the single-break location of a one-dimensional signed trajectory.

A descriptive changepoint read: the axis position that best splits the signed trajectory into two independent least-squares segments, minimising the combined residual sum of squares. The break is reported at the midpoint of the two focal points it falls between. It is deliberately two independent lines (a discontinuity is allowed), so a level shift such as a DSM-5 (2013) boundary on the era axis is localised, not smoothed over. It is labelled descriptive: the bridge supremum-LM confidence set saturates at the full sample size, so the break location is read with its bootstrap spread rather than a resolved confidence set.

Parameters:

• positions (numpy.ndarray) – The focal positions, shape (n_focal,).

• series (numpy.ndarray) – The one-dimensional signed trajectory, shape (n_focal,).

• min_segment (int, optional) – The fewest focal points each segment must hold.

Returns:

The break location in axis units, or not-a-number when the series is too short or flat.

Return type:

float

class analysis.trajectory_local.DirectionalResult(slope, net_direction, signed_slope, net_trend, signed_trajectory, slope_norm, break_position, span, focal_points)

The observed per-class directional statistic of a displacement trajectory.

slope

The per-feature slope vector $b_k$, shape (n_classes, n_features).

Type:

numpy.ndarray

net_direction

The per-class unit net direction $\hat u_k$, shape (n_classes, n_features).

Type:

numpy.ndarray

signed_slope

The signed net-projected slope $\langle b_k, \hat u_k\rangle$, in standardised displacement per axis unit, shape (n_classes,). The directional statistic.

Type:

numpy.ndarray

net_trend

The separation-scaled net-trend displacement, the signed slope times the focal span over the between-class separation, shape (n_classes,); the interpretable effect size (how far, in separation units, the linear trend carries the class across the axis).

Type:

numpy.ndarray

signed_trajectory

The one-dimensional signed trajectory $s_k(f)$, shape (n_classes, n_focal).

Type:

numpy.ndarray

slope_norm

The Euclidean norm of the slope vector, shape (n_classes,); reported for context and known to be positively biased, so not the test statistic.

Type:

numpy.ndarray

break_position

The single-break location on each signed trajectory, shape (n_classes,).

Type:

numpy.ndarray

span

The focal span (maximum minus minimum focal position), in axis units.

Type:

float

focal_points

The focal positions, shape (n_focal,).

Type:

numpy.ndarray

analysis.trajectory_local.directional_statistic(displacement, pooled_sd, focal_points, separation_scale)

Compute the observed per-class directional statistic of a displacement trajectory.

Standardises the displacement, fits the per-feature slope against the axis, projects it onto the class’s net direction to get the signed directional statistic, and scales it to a separation-unit net-trend effect size. Also returns the one-dimensional signed trajectory, the (biased) slope norm, and the single-break changepoint location. Pure and cheap: it consumes the trajectory observed_trajectory() already computed and refits nothing.

Parameters:

• displacement (numpy.ndarray) – The per-feature displacement $d_k(f)$, shape (n_classes, n_focal, n_features) (ObservedTrajectory.displacement).

• pooled_sd (numpy.ndarray) – The per-feature pooled standard deviation, shape (n_features,).

• focal_points (numpy.ndarray) – The focal positions, shape (n_focal,).

• separation_scale (float) – The between-class separation the net trend is divided by.

Returns:

The observed directional statistic and its parts.

Return type:

DirectionalResult

class analysis.trajectory_local.DirectionalInference(net_trend, net_trend_lo, net_trend_hi, signed_slope, signed_slope_lo, signed_slope_hi, p_value, reject, break_position, break_lo, break_hi)

The per-class directional test: effect size, clustered-bootstrap interval, and FDR.

net_trend

The observed separation-scaled net-trend displacement per class, shape (n_classes,).

Type:

numpy.ndarray

net_trend_lo, net_trend_hi

The clustered-bootstrap interval of the net trend, shape (n_classes,).

Type:

numpy.ndarray

signed_slope

The observed signed net-projected slope per class, shape (n_classes,).

Type:

numpy.ndarray

signed_slope_lo, signed_slope_hi

The clustered-bootstrap interval of the signed slope, shape (n_classes,).

Type:

numpy.ndarray

p_value

The two-sided bootstrap $p$-value that the signed slope differs from zero, shape (n_classes,).

Type:

numpy.ndarray

reject

The Benjamini-Hochberg decision across the classes at level q, shape (n_classes,); a rejected class is directional along this axis.

Type:

numpy.ndarray

break_position

The observed single-break location per class, shape (n_classes,).

Type:

numpy.ndarray

break_lo, break_hi

The bootstrap spread of the break location, shape (n_classes,); descriptive.

Type:

numpy.ndarray

analysis.trajectory_local.directional_inference(observed, tube, *, q=0.05)

Test each class’s directional statistic against the clustered-bootstrap null.

The signed net-projected slope is a signed scalar, so its clustered-bootstrap distribution (frozen net direction, families resampled) gives a two-sided add-one $p$-value that it differs from zero, floored at one over the replicate count plus one. Benjamini-Hochberg control is applied across the classes at level q; a rejected class is directional along the axis. The net-trend and signed-slope intervals are the bootstrap percentiles, and the break location carries its own bootstrap spread. The clustered bootstrap, not the bare slope, calls significance, because the slope norm is positively biased.

Parameters:

• observed (DirectionalResult) – The observed directional statistic (directional_statistic()).

• tube (BootstrapTube) – The clustered-bootstrap tube built with the frozen net directions, carrying the directional draws.

• q (float, optional) – The false-discovery-rate level across the classes.

Returns:

The per-class effect size, interval, $p$-value, FDR decision, and break spread.

Return type:

DirectionalInference

The block-attribution engine: what carries or co-moves with the class drift (plan section 7f).

The drift stage reports each class’s movement as one separation-scaled number, and the section-7e recast (analysis.trajectory_local) resolves it to a per-class, per-feature displacement trajectory read from the single cached fit. This module asks what that movement is about, and it does so for two kinds of candidate.

A block is a named, proband-indexed set of columns. An internal block is a subset of the reference’s own features, so its sub-displacement is a slice of the drift vector the recast already computed and the per-feature squared magnitudes sum back to the whole-class distance. The seven author categories (H0F) and the two instrument referents (H0G) are the internal blocks the earlier stages already read; this module is the generalisation of their dict[str, ndarray] grain seam. An external block is data the model was never fit on (a held-out phenotype, a polygenic score, a microbiome summary). Its class profile is still defined, because the frozen responsibilities $r_{ik}$ weight any proband-level quantity, but that profile is a separate quantity in its own space that cannot be summed into the phenotype drift, only compared with it.

Two modes follow.

Co-drift is defined for any block. The trap it must avoid is the axis: the phenotype drift is a function of era or age, so any block that varies with the axis correlates with the drift trivially, and beating a random control ordering is no evidence when the block genuinely depends on the real axis. The signal that is not a shared-axis artefact is which classes move. The phenotype drifts hardest in the developmental class along age; a uniform axis effect on a block would move all four classes alike. Co-drift is therefore the correlation between the block’s per-class displacement-magnitude profile and the phenotype’s per-class profile, read at the axis endpoint, with a family-clustered bootstrap tube. A positive, zero-excluding correlation is co-drift (the block moves in the same classes); a correlation covering zero is dissociation. The four-class profile is coarse, so the read stays qualitative, aligned against dissociated, with the tube carrying the uncertainty. This is the general form of the H0J question.

Conditioning asks whether a block accounts for the drift: remove the block’s linear contribution to the features and report whether the class displacement shrinks. It is well posed only for a low-dimensional covariate, so a high-dimensional block earns a conditioning score only through an explicit, named reduction (a composite score, a chosen taxon, a single polygenic score). The read is a descriptive partial association, not a causal mediation claim: a covariate that is itself a consequence of the same developmental process would be over-adjusted.

Everything here is a pure consumer of the cached fit’s frozen responsibilities and the strata axis; it refits nothing. The functions take arrays and return results, leaving the cohort joins and artefact writing to the command layer.

analysis.blocks.per_class_profile(x_values, responsibilities, weights, pooled_sd, separation_scale)

Return each class’s separation-scaled displacement magnitude at one focal point.

The local centroid of every class under the frozen responsibilities and the kernel window, minus the pooled centroid, read as a per-class magnitude over the block’s columns. With the block’s own pooled_sd and separation_scale this is the block’s per-class drift profile; with the reference feature matrix it is the phenotype’s.

Parameters:

• x_values (numpy.ndarray) – The block’s measurement matrix, shape (n_probands, n_columns).

• responsibilities (numpy.ndarray) – The frozen posterior responsibilities $r_{ik}$, shape (n_probands, n_classes).

• weights (numpy.ndarray) – The per-proband Gaussian kernel weight at the focal point, shape (n_probands,).

• pooled_sd (numpy.ndarray) – The block’s per-column pooled standard deviation, shape (n_columns,).

• separation_scale (float) – The between-class separation the magnitude is divided by. A single scalar, so it cancels in the co-drift correlation; it is kept for the magnitude to stay interpretable.

Returns:

The per-class magnitude, shape (n_classes,); a class with no local weight is not-a-number.

Return type:

numpy.ndarray

analysis.blocks.profile_alignment(profile_a, profile_b)

Return the class-resolved co-drift: the correlation of two per-class magnitude profiles.

A Pearson correlation of the two four-class profiles, so a block that peaks in the same class as the phenotype scores near $+1$ (co-drift), a block that peaks in a different class scores negative (dissociation), and a block that moves every class alike, a flat profile, scores near zero. Centring makes it insensitive to the shared positive level of two magnitude vectors, which a raw cosine is not. Classes not-a-number in either profile (no local weight) are dropped; with fewer than two shared finite classes, or with either profile flat, the alignment is zero.

Parameters:

• profile_a (numpy.ndarray) – The two per-class magnitude profiles, each shape (n_classes,).

• profile_b (numpy.ndarray) – The two per-class magnitude profiles, each shape (n_classes,).

Returns:

The correlation in [-1, 1], or 0.0 when it is undefined.

Return type:

float

class analysis.blocks.CoDriftResult(phenotype_profile, block_profile, alignment, ci_low, ci_high, p_value, aligned, n_joint)

The class-resolved co-drift of one block with the phenotype drift.

phenotype_profile, block_profile

The observed per-class magnitude profiles, each shape (n_classes,).

Type:

numpy.ndarray

alignment

The observed profile correlation (profile_alignment()).

Type:

float

ci_low, ci_high

The family-bootstrap 2.5 and 97.5 percentiles of the alignment.

Type:

float

p_value

The two-sided add-one bootstrap $p$ against no alignment, floored at one over the number of bootstrap replicates plus one.

Type:

float

aligned

Whether the interval excludes zero on the positive side (co-drift rather than dissociation).

Type:

bool

n_joint

The number of probands finite on both the axis and the block, the paired sample.

Type:

int

analysis.blocks.co_drift(x_phenotype, x_block, responsibilities, families, axis_values, bandwidth, focal, *, pooled_sd_phenotype, pooled_sd_block, separation_phenotype, separation_block, n_boot, seed)

Test whether a block co-drifts with the phenotype in the same classes.

Both profiles are read along the same axis at the same endpoint focal point, since the block co-drifts along era or age, not along a variable of its own. Every bootstrap replicate resamples one set of families and recomputes both the phenotype and the block profile and their correlation on that same resample, so the two share their sampling variation and the alignment is a genuine paired statistic. The observed alignment then acts as its own bootstrap-inverted test, the construction analysis.trajectory_local.control_specificity_bootstrap() uses.

The inputs must already be restricted to the probands finite on both the axis and the block, so that a family resampled for one profile is resampled for the other.

Parameters:

• x_phenotype (numpy.ndarray) – The reference feature matrix, shape (n_joint, n_features).

• x_block (numpy.ndarray) – The block’s matrix over the same rows, shape (n_joint, n_columns).

• responsibilities (numpy.ndarray) – The frozen responsibilities over the same rows, shape (n_joint, n_classes).

• families (numpy.ndarray) – The per-proband family identifier over the same rows, the clustering unit.

• axis_values (numpy.ndarray) – The timing axis over the same rows.

• bandwidth (float) – The axis Gaussian kernel bandwidth.

• focal (float) – The axis endpoint focal position both profiles are read at.

• pooled_sd_phenotype (numpy.ndarray) – Each matrix’s per-column pooled standard deviation.

• pooled_sd_block (numpy.ndarray) – Each matrix’s per-column pooled standard deviation.

• separation_phenotype (float) – Each matrix’s between-class separation (cancels in the correlation, kept for the magnitude).

• separation_block (float) – Each matrix’s between-class separation (cancels in the correlation, kept for the magnitude).

• n_boot (int) – The number of paired bootstrap replicates.

• seed (int) – The bootstrap seed.

Returns:

The observed profiles, their alignment, and its bootstrap interval and $p$-value.

Return type:

CoDriftResult

class analysis.blocks.Decomposition(partitions, squared_magnitude, whole, share)

An internal block’s per-partition share of a class’s drift.

partitions

The partition names (the grain keys), in a fixed order.

Type:

list of str

squared_magnitude

The separation-scaled squared magnitude of each partition, shape (n_classes, n_partitions).

Type:

numpy.ndarray

whole

The whole-block squared magnitude per class, shape (n_classes,); the partition squared magnitudes sum to this (the partitions tile the columns).

Type:

numpy.ndarray

share

Each partition’s fraction of the whole per class, shape (n_classes, n_partitions).

Type:

numpy.ndarray

analysis.blocks.decompose(displacement, pooled_sd, partitions, separation_scale)

Split an internal block’s class displacement into per-partition squared magnitudes.

For a partition of the block’s columns (the seven author categories, the two referents, or any tiling), the separation-scaled squared magnitude is additive: because the magnitude is a Euclidean norm over standardised per-feature displacements, the squared magnitudes of a set of disjoint column groups that cover the block sum to the whole-block squared magnitude. This is the additive sum-of-squares share the H0G referent decomposition already reports, lifted to an arbitrary partition.

Parameters:

• displacement (numpy.ndarray) – The per-feature class displacement at the reported focal point, shape (n_classes, n_columns).

• pooled_sd (numpy.ndarray) – The per-column pooled standard deviation, shape (n_columns,).

• partitions (dict of str to numpy.ndarray) – The partition name mapped to its integer column indices. The groups must be disjoint and cover every column for the shares to sum to one.

• separation_scale (float) – The between-class separation the magnitude is divided by.

Returns:

The per-partition squared magnitudes, the whole, and the shares.

Return type:

Decomposition

class analysis.blocks.GrainContrast(contrast, ci_low, ci_high, p_value, reject, rms_a, rms_b, share_a, share_b)

The size-fair contrast between two internal grains and its bootstrap test.

The general form of the H0G referent contrast, over any two disjoint column groups of an internal block (current against retrospective instruments, one category against the rest, one instrument against another). The statistic is size-fair: each grain’s intensity is the root-mean-square standardised displacement over its own features, so a grain does not win by holding more features. The additive sum-of-squares shares are returned alongside as the descriptive split, but the contrast, not the share, is the test.

contrast

The observed group-A-minus-group-B root-mean-square contrast per class, shape (n_classes,).

Type:

numpy.ndarray

ci_low, ci_high

The clustered-bootstrap interval of the contrast per class, shape (n_classes,).

Type:

numpy.ndarray

p_value

The two-sided add-one bootstrap $p$-value per class, shape (n_classes,).

Type:

numpy.ndarray

reject

The Benjamini-Hochberg decision across the classes at level q, shape (n_classes,).

Type:

numpy.ndarray

rms_a, rms_b

Each grain’s size-fair root-mean-square intensity per class, shape (n_classes,).

Type:

numpy.ndarray

share_a, share_b

Each grain’s additive sum-of-squares share per class (summing to one over the two grains), shape (n_classes,).

Type:

numpy.ndarray

analysis.blocks.grain_contrast(feature_draws, observed_displacement, group_a_cols, group_b_cols, *, q=0.05)

Test the per-class size-fair contrast between two internal grains from the tube draws.

Reads the size-fair root-mean-square intensity of two disjoint column groups from the standardised endpoint displacement, forms the group-A-minus-group-B contrast per class, and calls significance from the clustered-bootstrap replicates the tube already holds (the standardised per-feature displacement at the endpoint over the same family resamples). No new bootstrap loop runs: the draws are paired (both grains re-read on the same resample) and family-clustered. The two-sided add-one $p$-value is the fraction of replicate contrasts on the far side of zero, doubled and floored, and Benjamini-Hochberg control is applied across the classes. The H0G referent split of the era drift is the instance where group A is the current-state grain and group B the retrospective grain (the invariance-trajectory era stage calls this with those two grains).

Parameters:

• feature_draws (numpy.ndarray) – The standardised per-feature displacement replicates at the endpoint, shape (n_boot, n_classes, n_features).

• observed_displacement (numpy.ndarray) – The observed standardised per-feature displacement at the endpoint, shape (n_classes, n_features).

• group_a_cols (numpy.ndarray) – The column indices of the two grains being contrasted.

• group_b_cols (numpy.ndarray) – The column indices of the two grains being contrasted.

• q (float, optional) – The false-discovery-rate level across the classes.

Returns:

The per-class contrast, its interval, $p$-value, and FDR decision, and the per-grain root-mean-square intensities and additive shares.

Return type:

GrainContrast

analysis.blocks.residualise(x_values, covariate)

Remove a covariate’s linear contribution from every column of a matrix.

The ordinary-least-squares residual of each column of x_values on [1, z], so the returned matrix is the part of the features that a linear reading of the covariate does not explain. Conditioning reads the class displacement on this residual: if the drift was carried by variation the covariate tracks, the residual displacement shrinks.

Parameters:

• x_values (numpy.ndarray) – The feature matrix, shape (n_probands, n_features).

• covariate (numpy.ndarray) – The reduced covariate, shape (n_probands,) or (n_probands, n_reduced).

Returns:

The residualised matrix, shape (n_probands, n_features).

Return type:

numpy.ndarray

class analysis.blocks.ConditioningResult(raw_magnitude, conditioned_magnitude, shrinkage)

How much a class’s drift survives removing a covariate’s linear contribution.

raw_magnitude, conditioned_magnitude

The per-class endpoint magnitude before and after residualising the features on the covariate, each in a fixed pooled-standard-deviation metric, shape (n_classes,).

Type:

numpy.ndarray

shrinkage

1 - conditioned / raw per class, the fraction of the drift the covariate accounts for; near one when the covariate carries the drift, near zero when it is irrelevant.

Type:

numpy.ndarray

analysis.blocks.conditioning_shrinkage(x_values, responsibilities, weights, covariate, pooled_sd, separation_scale)

Report how much of each class’s drift a reduced covariate accounts for.

The class displacement is read twice in the same fixed metric (the raw features’ pooled standard deviation), once on the features and once on the features residualised on the covariate (residualise()). A covariate that carries the along-axis movement leaves a small residual displacement, so the magnitude shrinks; an irrelevant covariate leaves the displacement almost unchanged. The shared metric and the shared separation_scale cancel in the ratio, so the shrinkage reads how much of the movement the covariate linearly explains, not a change of units.

This is a descriptive partial association, H0H’s adjustment generalised to a declared reduction, not a causal mediation claim: a covariate that is itself a downstream consequence of the drift would be over-adjusted.

Parameters:

• x_values (numpy.ndarray) – The reference feature matrix, shape (n_probands, n_features).

• responsibilities (numpy.ndarray) – The frozen responsibilities, shape (n_probands, n_classes).

• weights (numpy.ndarray) – The per-proband kernel weight at the endpoint focal point, shape (n_probands,).

• covariate (numpy.ndarray) – The low-dimensional reduction, shape (n_probands,) or (n_probands, n_reduced).

• pooled_sd (numpy.ndarray) – The fixed per-feature metric both readings use, shape (n_features,).

• separation_scale (float) – The between-class separation both readings divide by (cancels in the shrinkage).

Returns:

The raw and conditioned per-class magnitudes and their shrinkage.

Return type:

ConditioningResult

The ordering-axis catalogue for the displacement atlas (plan section 12b).

The specificity check reads a class’s separation-scaled endpoint displacement along an ordering variable. The timing axes (diagnostic era, age at diagnosis) are two such orderings; this module generalises the ordering to any continuous or ordered variable that is not one of the 238 clustered features, so the same displacement can be screened across many axes and sorted from the largest to the smallest mover. The atlas stage (analysis displacement-atlas) consumes this catalogue.

The axes are the timing pair (diagnostic era and age at diagnosis, the mechanism under test), a covariate pool the cohort already carries or that timing derives (the measurement-to-diagnosis lag, age at evaluation, household income, and the area deprivation index), and a seeded random ordering. The random ordering is the floor: the one control guaranteed to carry no real structure, so an axis whose displacement clears it is above sampling noise. There is no covariate-orthogonality assumption here; the atlas reports every axis’s displacement and leaves the random floor as the only reference.

Two roles are deliberately absent. The 238 clustered features and any total taken over them (an SCQ or RBS-R sum) would order probands by a feature the classes were built from, moving the class centroids by construction, a circular self-drift rather than an external ordering. Held-out phenotype instruments (adaptive behaviour, motor coordination, IQ) are still phenotype, correlated with the clustered features by construction, so they are a non-null ceiling rather than a clean ordering and are left out for the same reason a symptom total is.

class analysis.axes.AxisContext(root, dataset, version, index, covariates, strata, seed)

The inputs an axis loader may draw on, built once per atlas run.

root

Repository root, for resolving the catalogue and source CSVs.

Type:

pathlib.Path

dataset, version

The cohort release the atlas runs on.

Type:

str

index

The modelling-cohort proband index every axis is reindexed onto.

Type:

pandas.Index

covariates

The cohort covariate frame (sex, age at evaluation, and the rest).

Type:

pandas.DataFrame

strata

The derived timing axes (era, age at diagnosis) and the lag.

Type:

analysis.strata_data.StrataData

seed

Base seed for the random ordering.

Type:

int

class analysis.axes.AxisSpec(name, label, kind, load)

One ordering axis for the displacement atlas.

name

The axis key: its column in the atlas artefact and its displacement-atlas handle.

Type:

str

label

A human-readable name for the figure.

Type:

str

kind

"timing", "covariate", "phenotype", or "random". The random axis is the floor; the timing axes are the mechanism under test; the rest are external orderings.

Type:

str

load

Maps an AxisContext to a proband-indexed ordering on the cohort index, with not-a-number where the variable is missing. A higher rank means more of the named quantity.

Type:

callable

The demographic covariate catalogue for the drift-attribution screen (plan section 7f).

The block-attribution engine (analysis.blocks) asks two questions of any proband-level quantity: does the phenotype partition drift along it (co-drift, generalised by the displacement atlas of analysis.axes), and does the timing drift shrink when the phenotype is residualised on it (conditioning). This module supplies the demographic covariates both reads consume: household socioeconomic position, parental education and occupation, family structure, inferred parental age at the child’s birth, perinatal complications, and the individual’s sex and race.

Each covariate is one DemographicSpec. An ordered or scalar covariate (income band, education level, a parental-age year, a complication count) also serves as an ordering axis, so the atlas can read the class displacement along it; a nominal covariate (family type, marital status, race) enters the conditioning read alone, as a low-dimensional one-hot block whose linear contribution is partialled out of the features. The coverage floor and the drop of a near-degenerate or thinly joined covariate are the command layer’s job (plan section 7f); this module only declares the covariates and reads them, one set of columns at a time, under the dscat guardrail.

Two facts about the reads govern how the covariates are built. First, a covariate orthogonal to the timing axis cannot account for a timing-ordered drift, so the conditioning shrinkage of such a covariate is near zero by construction; the screen reports each covariate’s association with the axis alongside its shrinkage so this ceiling is visible. Second, the shrinkage is a descriptive partial association, not a causal claim: a covariate that is itself downstream of the diagnosed phenotype (family living arrangement, say) would be over-adjusted, and the prose says so.

The orderings imposed on the ordinal covariates (the education ladder, the collapse of the occupation and living-arrangement categories) are modelling choices, recorded here so a reader can see and contest them rather than find them buried in a preprocessing step.

class analysis.demographics.DemographicSpec(name, label, kind, coding, load)

One demographic covariate for the drift-attribution screen.

name

The covariate key: its handle on the command line and its stem in the artefacts.

Type:

str

label

A human-readable name for the figure.

Type:

str

kind

The covariate family: "ses", "family", "parental", or "individual".

Type:

str

coding

"ordinal", "scalar", "count", "binary", or "onehot". The first four give a single ordered column that also serves as an atlas ordering axis; "onehot" gives a low-dimensional block for the conditioning read alone.

Type:

str

load

Maps an analysis.axes.AxisContext to a proband-indexed frame on the cohort index, with not-a-number where the covariate is missing. Ordered codings return one column; a one-hot coding returns its indicator columns.

Type:

callable

property ordered: bool

Whether the covariate is a single ordered column the atlas can order probands by.

Prevalence drift: how the frozen class proportions vary along an axis (H0B).

The H0B hypothesis asks whether the four Litman class mixing proportions are constant across diagnostic era and across age at diagnosis, or whether at least one class’s proportion trends along an axis. It is distinct from profile drift (the H0A family): the classes are held fixed at the measurement-only reference fit and only their sizes are read as a function of the axis, so nothing is re-estimated. The estimand is the mixing proportions as a function of the axis.

The rigorous read is a three-step estimator (Vermunt 2010; Bakk, Tekle and Vermunt 2013). The frozen reference gives each proband a posterior over the classes; a naive regression of the hard label on the axis inherits the classification error of that assignment (the classify-analyse bias), which attenuates a real slope and can manufacture a spurious one. The maximum-likelihood (ML) correction removes it: the modal assignment is treated as a single categorical indicator whose class-conditional error probabilities are fixed at the confusion matrix of the frozen posteriors, and the structural model, a multinomial logit of the true latent class on the axis, is fitted by expectation-maximisation with that measurement model held fixed. No mixture is refitted; the correction reuses the same confusion matrix StepMix builds for its own three-step estimator (stepmix.stepmix.compute_bch_matrix).

Because the structural model is a small logit rather than a weighted covariate emission, this implements the correction directly rather than through StepMix’s fit path, which always re-estimates the measurement model in its first step and, on this cohort, is numerically unstable under fractional weights (progress log, 2026-07-05). The naive hard-label regression is reported beside the corrected one as a transparent, uncorrected cross-check.

Uncertainty is a family-clustered bootstrap resampling SPARK families, matching the invariance-trajectory convention: the corrected slope, its odds ratio, and the predicted proportion curve each carry a percentile interval, and the corrected slope carries a two-sided add-one bootstrap $p$. The naive model additionally carries the closed-form Wald and likelihood-ratio $p$-values. Significance is Benjamini-Hochberg controlled across the per-class contrasts within an axis. Everything here is class or coefficient level; no per-proband quantity is returned.

analysis.prevalence.modal_assignment(responsibilities)

Return the hard (modal) class assignment from a posterior responsibility matrix.

Parameters:

responsibilities (numpy.ndarray) – Per-proband posterior over the classes, shape (n, K).

Returns:

The index of the most probable class per proband, shape (n,).

Return type:

numpy.ndarray

analysis.prevalence.bch_confusion(responsibilities, assignment='modal')

Return the classification-error matrix of a frozen posterior.

Wraps stepmix.stepmix.compute_bch_matrix, the confusion matrix at the heart of the three-step corrections. The entry D[c, s] is the probability of assigning a proband to class s given that its true latent class is c, estimated from the posterior itself.

Parameters:

• responsibilities (numpy.ndarray) – Per-proband posterior over the classes, shape (n, K).

• assignment ({"modal", "soft"}, optional) – Whether the predicted class is the modal assignment or the soft posterior, matching StepMix’s assignment option.

Returns:

The confusion matrix D of shape (K, K), rows indexed by true class.

Return type:

numpy.ndarray

class analysis.prevalence.MultinomialFit(coef, loglik)

A fitted multinomial logit with baseline-category coefficients.

coef

Coefficients of shape (p, K) with the baseline class column held at zero.

Type:

numpy.ndarray

loglik

The model log-likelihood. For a naive fit this is the multinomial log-likelihood of the hard labels; for a corrected fit it is the observed-data log-likelihood of the ML correction.

Type:

float

proportions(design)

Return the predicted class proportions for a design matrix.

Parameters:

design (numpy.ndarray) – Design matrix of shape (m, p) sharing the fit’s columns.

Returns:

Predicted proportions of shape (m, K), each row summing to one.

Return type:

numpy.ndarray

analysis.prevalence.fit_soft_multinomial(design, targets, *, ridge=1e-06)

Fit a multinomial logit to soft class targets.

Minimises the cross-entropy between the softmax of design @ coef and the target rows, with the first class held as the baseline (its coefficient column fixed at zero) and a small ridge penalty for identifiability under quasi-separation. The targets are non-negative rows; they are the posterior responsibilities in the correction’s M-step and one-hot label indicators for a naive fit.

Parameters:

• design (numpy.ndarray) – Design matrix of shape (n, p), including an intercept column.

• targets (numpy.ndarray) – Non-negative class weights of shape (n, K).

• ridge (float, optional) – The L2 penalty on the free coefficients.

Returns:

Coefficients of shape (p, K) with the baseline column at zero.

Return type:

numpy.ndarray

analysis.prevalence.fit_naive_multinomial(design, labels, *, n_classes)

Fit the naive multinomial logit of hard labels on a design (the uncorrected cross-check).

Parameters:

• design (numpy.ndarray) – Design matrix of shape (n, p).

• labels (numpy.ndarray) – Hard class labels in 0 .. n_classes - 1, shape (n,).

• n_classes (int) – The number of classes K.

Returns:

The fitted coefficients and multinomial log-likelihood.

Return type:

MultinomialFit

analysis.prevalence.fit_corrected_multinomial(design, responsibilities, *, assignment='modal', max_iter=200, tol=1e-08)

Fit the ML-corrected multinomial logit of latent class on a design.

The modal assignment is treated as a single categorical measurement with class-conditional error probabilities fixed at the confusion matrix of the frozen posteriors. The structural multinomial logit is fitted by expectation-maximisation with that measurement model held fixed, so the classify-analyse bias of the naive fit is removed without refitting the mixture.

Parameters:

• design (numpy.ndarray) – Design matrix of shape (n, p), including an intercept column.

• responsibilities (numpy.ndarray) – Frozen per-proband posterior over the classes, shape (n, K).

• assignment ({"modal", "soft"}, optional) – The assignment rule for the confusion matrix.

• max_iter (int, optional) – The maximum number of expectation-maximisation iterations.

• tol (float, optional) – The absolute log-likelihood tolerance for convergence.

Returns:

The corrected coefficients and the observed-data log-likelihood of the correction model.

Return type:

MultinomialFit

class analysis.prevalence.SlopeResult(ref_class, slope, odds_ratio, ci_low=nan, ci_high=nan, wald_p=nan, lrt_p=nan, boot_p=nan, reject=False)

A per-class one-versus-rest axis slope under one estimator.

ref_class

The class index.

Type:

int

slope

The axis log-odds slope (log odds of membership per axis unit).

Type:

float

odds_ratio

exp(slope).

Type:

float

ci_low, ci_high

The slope confidence interval (family-clustered bootstrap percentiles).

Type:

float

wald_p

The closed-form Wald $p$-value on the slope (naive estimator only; not-a-number for the corrected estimator, which uses the bootstrap $p$).

Type:

float

lrt_p

The closed-form likelihood-ratio $p$-value for the axis term (naive estimator only).

Type:

float

boot_p

The two-sided add-one family-clustered-bootstrap $p$-value.

Type:

float

reject

Whether the contrast survives Benjamini-Hochberg control within the axis.

Type:

bool

analysis.prevalence.corrected_onevsrest(design, responsibilities, ref_class, axis_col, *, assignment='modal', max_iter=200, tol=1e-08)

Return the ML-corrected one-versus-rest axis slope for one class.

Collapses the posterior to the indicator “belongs to ref_class”, forms the two-by-two confusion matrix of that binary assignment, and fits a binary logit of the latent indicator on the design by expectation-maximisation with the confusion matrix held fixed. The returned slope is the coefficient on the axis column.

Parameters:

• design (numpy.ndarray) – Design matrix of shape (n, p) including an intercept column.

• responsibilities (numpy.ndarray) – Frozen posterior over the classes, shape (n, K).

• ref_class (int) – The class whose membership is modelled.

• axis_col (int) – The design column whose coefficient is the axis effect.

• assignment ({"modal", "soft"}, optional) – The assignment rule for the confusion matrix.

• max_iter (int, optional) – The maximum number of expectation-maximisation iterations.

• tol (float, optional) – The absolute log-likelihood tolerance for convergence.

Returns:

The corrected axis log-odds slope for the class.

Return type:

float

analysis.prevalence.naive_onevsrest(design, labels, ref_class, axis_col)

Return the naive one-versus-rest axis slope for one class, with Wald and LRT tests.

A binary logit of the hard-label indicator on the design, fitted by statsmodels so the Wald standard error and the likelihood-ratio test against the axis-free model are exact. It is uncorrected for classification error and reported as a cross-check.

Parameters:

• design (numpy.ndarray) – Design matrix of shape (n, p) including an intercept column.

• labels (numpy.ndarray) – Hard class labels, shape (n,).

• ref_class (int) – The class whose membership is modelled.

• axis_col (int) – The design column whose coefficient is the axis effect.

Returns:

The slope, odds ratio, Wald $p$, and likelihood-ratio $p$ (bootstrap fields left unset).

Return type:

SlopeResult

class analysis.prevalence.ProportionCurve(positions, corrected, naive, band_lo, band_hi, pooled)

Predicted class-proportion curves over an axis grid, with bootstrap bands.

positions

The axis grid, shape (G,).

Type:

numpy.ndarray

corrected, naive

Predicted proportions of shape (K, G) under each estimator.

Type:

numpy.ndarray

band_lo, band_hi

The family-clustered bootstrap percentile band on the corrected curve, shape (K, G).

Type:

numpy.ndarray

pooled

The pooled (axis-free) class proportion, the mean responsibility per class over the cohort, shape (K,). This is the original mixing weight the curve trends away from.

Type:

numpy.ndarray

class analysis.prevalence.JointTest(estimator, lr_stat, df, p_value)

A joint likelihood-ratio test of class ~ axis against class ~ 1.

estimator

"corrected" or "naive".

Type:

str

lr_stat

The likelihood-ratio statistic.

Type:

float

df

The degrees of freedom, (K - 1) times the number of axis terms.

Type:

int

p_value

The chi-square tail probability.

Type:

float

class analysis.prevalence.PrevalenceResult(axis, n, axis_mean, corrected_slopes, naive_slopes, adjusted_slopes, joint_tests, curve, dsm_contrasts=<factory>)

The full H0B read for one axis.

axis

The axis name.

Type:

str

n

The number of probands with a finite axis value.

Type:

int

axis_mean

The centring constant subtracted from the axis before fitting (slopes are per raw unit).

Type:

float

corrected_slopes, naive_slopes

The per-class one-versus-rest axis slopes under each estimator (unadjusted).

Type:

list of SlopeResult

adjusted_slopes

The corrected axis slopes net of sex, the diagnosis-to-measurement lag, and age at evaluation.

Type:

list of SlopeResult

joint_tests

The joint likelihood-ratio tests under each estimator.

Type:

list of JointTest

curve

The predicted proportion curves and their band.

Type:

ProportionCurve

dsm_contrasts

For the era axis, the per-class pre/post-2013 (DSM-5) log-odds contrasts; empty otherwise.

Type:

list of SlopeResult

analysis.prevalence.prevalence_analysis(responsibilities, labels, axis_values, families, *, axis, covariates=None, grid, n_boot=500, seed=0, q=0.05, assignment='modal')

Run the H0B prevalence-drift analysis for one axis.

Fits the corrected and naive per-class one-versus-rest slopes, the joint likelihood-ratio tests, the predicted proportion curves, the family-clustered bootstrap uncertainty, and (for era) the DSM-5 pre/post-2013 contrast. The axis is centred before fitting for numerical conditioning; the reported slopes are per raw axis unit. Significance is Benjamini-Hochberg controlled across the per-class contrasts within the axis, separately for each estimator.

Parameters:

• responsibilities (numpy.ndarray) – Frozen posterior over the classes, shape (n, K).

• labels (numpy.ndarray) – Hard class labels, shape (n,).

• axis_values (numpy.ndarray) – The axis value per proband, shape (n,); must be finite (filter upstream).

• families (numpy.ndarray) – The per-proband family identifier, the bootstrap clustering unit, shape (n,).

• axis (str) – The axis name, used to decide whether the DSM-5 contrast is read.

• covariates (dict of str to numpy.ndarray, optional) – The adjustment covariates (sex, lag, age_at_eval) for the sensitivity model.

• grid (numpy.ndarray) – The axis grid the proportion curves are read at, shape (G,).

• n_boot (int, optional) – The number of family-clustered bootstrap replicates.

• seed (int, optional) – The base seed for the bootstrap.

• q (float, optional) – The Benjamini-Hochberg level.

• assignment ({"modal", "soft"}, optional) – The assignment rule for the confusion matrices.

Returns:

The per-class slopes, joint tests, proportion curves, and DSM contrasts for the axis.

Return type:

PrevalenceResult

Support for the number of latent classes, per stratum, by a bootstrap likelihood-ratio search.

The H0C hypothesis (plan section 7, H0C) asks whether the number of components the Litman four-class general finite mixture model supports is stable across strata of age at diagnosis and diagnostic era, or whether a class splits or merges in some stratum. This module answers that per stratum, and relative to the pooled cohort put through the identical procedure, so the shared over-extraction from feature misspecification (which drove the pooled information criteria to nine classes at this sample size) cancels: an order change is a stratum whose supported order differs from the pooled cohort’s, not a stratum whose raw information criteria differ from four.

The measurement model is fitted measurement-only (structural=None, the fit --no-covariates recipe), on hard subsets of the 238-feature cohort matrix. The confirmatory statistic is a warm-started parametric bootstrap likelihood-ratio test (BLRT):

• the observed statistic for a step is $text{LR} = 2(ell_{K+1} - ell_K)$, the twice log-likelihood gain from one extra class, where the $K$-class fit uses a handful of random restarts and the $K+1$-class fit is warm-started by splitting each of the $K$ classes in turn (so $K$ warm starts) plus a couple of random restarts, keeping the best;

• the null is parametric: datasets are simulated from the fitted $K$-component model at the stratum’s own sample size, and each is put through the identical fitting recipe, so the null LR is not biased low by under-fitting the alternative. The $p$-value is the Phipson-Smyth add-one permutation $p$ (as in analysis.drift).

The search is sequential and anchored at four classes (sequential_search()): the primary (splitting) direction tests $4$ against $5$ and steps outward while each step rejects, up to a cap; the secondary (merging) direction tests $4$ against $3$ and steps down while each step fails to reject, down to a floor. Two corroborators, neither of which decides, sit beside the BLRT: the cross-validated log-likelihood elbow, with the knee found by an inline Kneedle (kneedle_knee()), and the proper adjusted Lo-Mendell-Rubin test (vlmr_test(), Formula 15 of Lo, Mendell and Rubin 2001), not the naive one-degree-of-freedom proxy of analysis.selection.lmr_lrt_proxy().

The fits reuse the StepMix internals the rest of the package already depends on (analysis.model.prepare_inputs(), set_parameters/get_parameters, the measurement emission sample and log_likelihood), and the degeneracy handling of analysis.drift (its DEGENERATE_FIT_ERRORS). bootstrap_lr() is a top-level, picklable unit of work (simulate one null dataset, refit, return its LR), so the bootstrap draws spread across a process pool.

analysis.order.measurement_inputs(matrix, typing)

Return the descriptor-aligned measurement matrix and the mixed-data descriptor.

A thin wrapper over analysis.model.prepare_inputs() that drops the covariate channel, since the H0C fits are measurement-only. The returned frame is rounded as the reference release rounds, so the model sees the same values the covariate reference was fitted on.

Parameters:

• matrix (analysis.cohort.CohortMatrix) – The cohort feature and covariate matrices.

• typing (analysis.features.Typing) – The reconciled feature typing that sets each feature’s density.

Returns:

The measurement matrix and the StepMix mixed-data descriptor.

Return type:

tuple

analysis.order.total_log_likelihood(model, x)

Return the total measurement log-likelihood of x under a fitted model.

The per-sample log-likelihood of a mixture is $log sum_k pi_k prod_j f_j(x_{ij}; theta_{jk})$, read off the emission log-likelihood (model._mm.log_likelihood) and the mixing weights (model.weights_), the same primitives analysis.invariance uses. Summed over samples, this is the quantity the likelihood-ratio statistic differences.

Parameters:

• model (StepMix) – A fitted measurement-only model.

• x (numpy.ndarray) – The measurement matrix in the descriptor column order the model was fitted on.

Returns:

The total log-likelihood, summed over samples.

Return type:

float

analysis.order.fit_k(measurement_data, descriptor, k, *, n_init, seed)

Fit a measurement-only k-class GFMM with random restarts.

Parameters:

• measurement_data (pandas.DataFrame) – The descriptor-aligned measurement matrix.

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• k (int) – The number of latent classes.

• n_init (int) – Random restarts, delegated to StepMix; the best restart is kept.

• seed (int) – Random seed for reproducible restarts.

Returns:

The fitted estimator.

Return type:

StepMix

analysis.order.split_warm_params(model, k, class_id, rng, jitter)

Return warm $K+1$ parameters by splitting one class of a fitted $K$-class model.

The split halves class class_id’s mixing weight into two children and perturbs its class-conditional parameters symmetrically by $pm$ a Gaussian jitter, giving a starting point for the extra class near the class most likely to divide. Applied to each of the $K$ classes in turn, this is the $K$-warm-start half of the alternative-model recipe (plan section 7): a good local optimum for $K+1$ that a purely random restart rarely reaches, so the observed and null alternative fits are on an equal footing.

The measurement parameters are a nested dict of per-emission arrays; the class axis is the one of length k. A block with no such axis (a shared parameter) is carried over unchanged.

Parameters:

• model (StepMix) – The fitted $K$-class model to split.

• k (int) – Its number of classes.

• class_id (int) – The class to split into two children.

• rng (numpy.random.Generator) – The generator for the symmetric jitter.

• jitter (float) – Standard deviation of the Gaussian perturbation applied to the split class.

Returns:

A parameter dict in get_parameters format with $K+1$ classes, for warm_em().

Return type:

dict

analysis.order.warm_em(descriptor, k, warm_params, x, *, max_iter=500, tol=1e-06)

Run the EM algorithm from given warm parameters, without reinitialising.

StepMix’s own fit reinitialises the emission parameters at the start of every EM run, so a warm start cannot go through it. This drives the EM loop directly instead: it sets the given parameters, then alternates the model’s own E-step and M-step until the average log-likelihood stops improving. The E-step and M-step are StepMix’s, so the fit is the same optimiser the random restarts use, only seeded from the split rather than from noise.

Parameters:

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• k (int) – The number of latent classes (one more than the model being split).

• warm_params (dict) – The starting parameters, in get_parameters format (from split_warm_params()).

• x (numpy.ndarray) – The measurement matrix in the descriptor column order.

• max_iter (int, optional) – Maximum EM iterations.

• tol (float, optional) – Absolute tolerance on the average log-likelihood for convergence.

Returns:

The warm-started fitted model.

Return type:

StepMix

analysis.order.fit_k_plus_one(model_k, measurement_data, descriptor, k, x, *, n_random, n_init_random, seed, jitter)

Fit the $K+1$-class alternative by class-splitting warm starts plus random restarts.

The alternative-model recipe (plan section 7), used identically for the observed data and every bootstrap sample: split each of the $K$ classes of model_k in turn (so $K$ warm starts, split_warm_params() then warm_em()), add n_random random restarts, and keep whichever reaches the highest total log-likelihood. A single warm start that fails numerically is skipped rather than fatal; if every start fails the caller reads the returned log-likelihood of minus infinity as a degenerate fit.

Parameters:

• model_k (StepMix) – The fitted $K$-class model whose classes seed the warm starts.

• measurement_data (pandas.DataFrame) – The descriptor-aligned measurement matrix (for the random restarts).

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• k (int) – The number of classes of model_k; the alternative has k + 1.

• x (numpy.ndarray) – The measurement matrix as an array (for the log-likelihood and the warm EM).

• n_random (int) – Random restarts to add to the $K$ warm starts.

• n_init_random (int) – Restarts within each random-restart fit.

• seed (int) – Base seed; the warm starts and random restarts derive their seeds from it.

• jitter (float) – Split-perturbation standard deviation.

Returns:

The best $K+1$ fit (None when none succeeded) and its total log-likelihood (minus infinity when none succeeded).

Return type:

tuple

class analysis.order.FitPair(model_k, model_k1, ll_k, ll_k1, k)

A fitted $K$-class model and its warm-started $K+1$-class alternative.

model_k

The $K$-class (null) fit.

Type:

StepMix

model_k1

The best $K+1$-class (alternative) fit, or None when every start was degenerate.

Type:

StepMix or None

ll_k

Total log-likelihood of the $K$-class fit.

Type:

float

ll_k1

Total log-likelihood of the $K+1$-class fit (minus infinity when degenerate).

Type:

float

k

The null number of classes.

Type:

int

property observed_lr: float

Return the observed likelihood-ratio statistic $2(ell_{K+1} - ell_K)$.

analysis.order.fit_pair(measurement_data, descriptor, k, *, n_init, n_random, seed, jitter)

Fit the $K$-class null and the warm-started $K+1$-class alternative on one dataset.

The observed side of one BLRT step. Applied unchanged to the real stratum and (through bootstrap_lr()) to every simulated null dataset, so the recipe is identical on both sides.

Parameters:

• measurement_data (pandas.DataFrame) – The descriptor-aligned measurement matrix.

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• k (int) – The null number of classes.

• n_init (int) – Random restarts for the $K$-class fit.

• n_random (int) – Random restarts added to the $K$ warm starts for the $K+1$-class fit.

• seed (int) – Base seed for the fits.

• jitter (float) – Split-perturbation standard deviation.

Returns:

The two fits and their log-likelihoods.

Return type:

FitPair

analysis.order.simulate_null(model_k, descriptor, columns, n, seed)

Simulate one parametric-bootstrap dataset from a fitted $K$-class model.

Reconstructs a sampleable StepMix from the fitted model’s parameters and draws n samples from the mixture (the measurement emissions and the mixing weights), the parametric null of the BLRT (plan section 7). The seed sets the emission generator, so a draw is reproducible from its index.

Parameters:

• model_k (StepMix or dict) – The fitted $K$-class model, or its get_parameters dict (for a worker process).

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• columns (pandas.Index) – The descriptor column order, so the sampled array is framed like the real matrix.

• n (int) – The stratum’s sample size.

• seed (int) – Seed for the emission sampler.

Returns:

A simulated measurement matrix of n rows in columns order.

Return type:

pandas.DataFrame

analysis.order.bootstrap_lr(null_params, descriptor, columns, n, k, *, n_init, n_random, seed, jitter)

Simulate one null dataset, refit both models, and return the null LR statistic.

A top-level, picklable unit of work (mirroring analysis.drift.summarise_pseudo_stratum()), so the bootstrap draws spread across a process pool. The worker samples its own dataset from null_params (small, passed once per pool), then puts it through the identical fitting recipe as the observed data. Returns None when a fit is numerically degenerate, so the caller drops that draw and counts it rather than letting one bad refit abort the run.

Parameters:

• null_params (dict) – The fitted $K$-class model’s get_parameters dict, the parametric null to sample from.

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• columns (list of str) – The descriptor column order.

• n (int) – The stratum’s sample size to simulate.

• k (int) – The null number of classes.

• n_init (int) – Random restarts for the $K$-class refit.

• n_random (int) – Random restarts added to the warm starts for the $K+1$-class refit.

• seed (int) – Seed for this draw (sampling and fitting).

• jitter (float) – Split-perturbation standard deviation.

Returns:

The null likelihood-ratio statistic, or None if the draw was degenerate.

Return type:

float or None

analysis.order.phipson_smyth_p(observed, null_draws)

Return the Phipson-Smyth add-one bootstrap $p$-value.

The proportion of null draws at or beyond the observed statistic, with the add-one correction so the smallest attainable $p$ is $1 / (B + 1)$ rather than zero (as in analysis.drift.read_against_null()). Non-finite draws (degenerate refits) are dropped.

Parameters:

• observed (float) – The observed likelihood-ratio statistic.

• null_draws (list of float) – The bootstrap null statistics.

Returns:

The add-one $p$-value, or nan when no finite draw survived.

Return type:

float

analysis.order.kneedle_knee(k_values, scores)

Return the knee of a concave, increasing curve by the Kneedle rule.

The cross-validated log-likelihood rises with the number of classes and flattens; the knee is where the diminishing returns set in, the elbow the H0C corroborator reads. This is the Kneedle construction (Satopaa et al. 2011) for a concave increasing curve: normalise both the class count and the score to the unit interval, take the difference between the normalised score and the normalised class count, and return the class count at its maximum. Ties and a flat curve fall back to the first candidate. This is an inline implementation (no new dependency), an objective substitute for reading the elbow off a plot.

Parameters:

• k_values (list of int) – The number of classes, ascending.

• scores (list of float) – The corresponding score (validation log-likelihood), one per k_values entry.

Returns:

The number of classes at the knee.

Return type:

int

analysis.order.vlmr_test(null_ll, null_params, null_classes, alt_ll, alt_params, alt_classes, n)

Return the adjusted Lo-Mendell-Rubin likelihood-ratio test for $K$ against $K+1$.

The proper VLMR of Lo, Mendell and Rubin (2001, Biometrika 88(3):767-778, Formula 15), not the naive one-degree-of-freedom proxy of analysis.selection.lmr_lrt_proxy(). The likelihood-ratio statistic is corrected by an ad-hoc factor and referred to a chi-square:

\[\text{LR} = 2(\ell_1 - \ell_0), \qquad \text{LMR} = \frac{\text{LR}}{1 + \big[((3k_1 - 1) - (3k_0 - 1)) \ln n\big]^{-1}},\]

with degrees of freedom the difference in free parameters $p_1 - p_0$ and $p = Pr(chi^2_{p_1 - p_0} > text{LMR})$. Here $k_0, k_1$ are the null and alternative class counts and $ell_0, ell_1$ their log-likelihoods. Reproduces calc_lrt of the R package tidyLPA.

Parameters:

• null_ll (float) – Log-likelihoods of the null ($k_0$-class) and alternative ($k_1$-class) fits.

• alt_ll (float) – Log-likelihoods of the null ($k_0$-class) and alternative ($k_1$-class) fits.

• null_params (int) – Numbers of free parameters (StepMix.n_parameters).

• alt_params (int) – Numbers of free parameters (StepMix.n_parameters).

• null_classes (int) – The class counts $k_0$ and $k_1$.

• alt_classes (int) – The class counts $k_0$ and $k_1$.

• n (int) – The sample size.

Returns:

lr (the raw statistic), lmr (the corrected statistic), df, and p.

Return type:

dict of str to float

class analysis.order.StepResult(k_null, k_alt, direction, observed_lr, p_value, b_used, escalated, n_dropped, rejected)

One BLRT step: a $k$-against-$k+1$ comparison with its staged bootstrap $p$-value.

k_null

The null number of classes.

Type:

int

k_alt

The alternative number of classes (k_null + 1).

Type:

int

direction

"split" for an upward (add-a-class) step, "merge" for a downward one.

Type:

str

observed_lr

The observed likelihood-ratio statistic.

Type:

float

p_value

The Phipson-Smyth bootstrap $p$-value over all draws used.

Type:

float

b_used

The number of finite bootstrap draws behind p_value.

Type:

int

escalated

Whether the step escalated from the screen count to the full count.

Type:

bool

n_dropped

Degenerate bootstrap draws dropped for this step.

Type:

int

rejected

Whether the step rejected its null at the step level (p_value <= alpha).

Type:

bool

class analysis.order.SearchResult(supported_k, capped, direction, steps=<factory>, elbow_knee=0, cv_log_likelihood=<factory>, vlmr=<factory>, n_dropped=0)

The sequential search outcome for one dataset (a stratum or the pooled cohort).

supported_k

The supported number of classes.

Type:

int

capped

Whether the search hit the upper cap while still rejecting, so supported_k is a lower bound (reported as ">=cap").

Type:

bool

direction

"split" if the order rose above the anchor, "merge" if it fell below, or "stable" if it stayed at the anchor.

Type:

str

steps

Every step taken, in order.

Type:

list of StepResult

elbow_knee

The cross-validated log-likelihood knee (Kneedle), a corroborator.

Type:

int

cv_log_likelihood

The validation log-likelihood per number of classes behind the knee.

Type:

dict of int to float

vlmr

The adjusted Lo-Mendell-Rubin test at the decisive comparison (empty when the search stayed at the anchor with no decisive step).

Type:

dict of str to float

n_dropped

Total degenerate bootstrap draws dropped across the search.

Type:

int

class analysis.order.Recipe(n_init=5, n_random=2, jitter=0.5)

The identical fitting recipe applied to the observed data and every null draw.

n_init

Random restarts for a $K$-class fit.

Type:

int

n_random

Random restarts added to the $K$ class-splitting warm starts for a $K+1$-class fit.

Type:

int

jitter

Standard deviation of the split perturbation.

Type:

float

spec()

Return the serialisable recipe, folded into the run hash.

class analysis.order.Schedule(b_screen=199, b_escalate=999, escalate_threshold=0.1, alpha=0.05)

The staged bootstrap schedule and the step-level rejection level.

b_screen

Draws every step is screened at.

Type:

int

b_escalate

Draws a step escalates to when its screen $p$ falls below escalate_threshold.

Type:

int

escalate_threshold

The screen-$p$ threshold that triggers escalation.

Type:

float

alpha

The step-level rejection level.

Type:

float

spec()

Return the serialisable schedule, folded into the run hash.

analysis.order.blrt_step(measurement_data, descriptor, k_null, *, direction, dispatch, recipe, schedule, seed)

Run one $k$-against-$k+1$ BLRT step with the staged bootstrap.

Fits the observed null and alternative (fit_pair()), then draws the parametric bootstrap null through dispatch: a screen of schedule.b_screen draws, escalated to schedule.b_escalate only when the screen $p$ falls below the threshold. The $p$-value is the Phipson-Smyth add-one $p$ over every finite draw.

Parameters:

• measurement_data (pandas.DataFrame) – The dataset’s descriptor-aligned measurement matrix.

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• k_null (int) – The null number of classes; the alternative has k_null + 1.

• direction (str) – "split" or "merge", recorded on the step.

• dispatch (callable) – The bootstrap-draw dispatcher (pool-backed in the CLI, in-process in a test).

• recipe (Recipe) – The identical fitting recipe for the observed and null fits.

• schedule (Schedule) – The staged bootstrap schedule and rejection level.

• seed (int) – Base seed for the observed fit (draw seeds are the dispatcher’s concern).

Returns:

The step outcome (StepResult) and the observed fit pair (FitPair, for the VLMR corroborator).

Return type:

tuple

analysis.order.cross_validated_log_likelihood(measurement_data, covariates, descriptor, k_values, *, seed, n_init, cv=3)

Return the 3-fold cross-validated log-likelihood per number of classes.

Reuses analysis.selection.validation_log_likelihood() (the released compute_LL), the elbow corroborator’s input. The elbow only corroborates the BLRT decision, so this uses the covariate scoring the selection stage already implements rather than a separate measurement-only cross-validation.

Parameters:

• measurement_data (pandas.DataFrame) – The dataset’s descriptor-aligned measurement matrix.

• covariates (pandas.DataFrame) – The structural covariate matrix on the same index.

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• k_values (list of int) – The class counts to grid over.

• seed (int) – Random seed for the estimator.

• n_init (int) – Random restarts per fold fit.

• cv (int, default 3) – Cross-validation folds.

Returns:

Each number of classes mapped to its mean validation log-likelihood.

Return type:

dict of int to float

analysis.order.sequential_search(measurement_data, descriptor, *, dispatch, recipe, schedule, anchor=4, cap=7, floor=2, seed=0, cv_scores=None)

Search for the supported number of classes, anchored at anchor.

From the anchor the primary (splitting) direction tests anchor against anchor + 1 and steps up while each step rejects, to the cap (">=cap" when it hits the cap still rejecting). If the first split step does not reject, the secondary (merging) direction tests anchor against anchor - 1 and steps down while each step fails to reject, to the floor. Each step is a warm-started parametric-bootstrap likelihood-ratio test (blrt_step()). The supported order is where the outward stepping stops (plan section 7).

The two corroborators are attached but do not decide: the cross-validated log-likelihood knee (from cv_scores via kneedle_knee()) and the adjusted Lo-Mendell-Rubin test at the decisive comparison (vlmr_test()).

Parameters:

• measurement_data (pandas.DataFrame) – The dataset’s descriptor-aligned measurement matrix.

• descriptor (dict) – The StepMix mixed-data measurement descriptor.

• dispatch (callable) – The bootstrap-draw dispatcher.

• recipe (Recipe) – The identical fitting recipe.

• schedule (Schedule) – The staged bootstrap schedule and rejection level.

• anchor (int, optional) – The anchor number of classes (four).

• cap (int, optional) – The upper cap; the search reports ">=cap" there.

• floor (int, optional) – The lower floor.

• seed (int, optional) – Base seed; each step derives its seed from it and its null class count.

• cv_scores (dict of int to float, optional) – The validation log-likelihood per number of classes for the elbow. When absent the knee is reported as zero.

Returns:

The supported order, the steps taken, and the corroborators.

Return type:

SearchResult

analysis.order.agreement_flag(result, pooled_supported_k)

Return whether a positive order-change claim is corroborated for a stratum.

A positive order-change claim requires agreement (plan section 7): the BLRT supported order differs from the pooled cohort’s, the cross-validated elbow knee has moved off the pooled order too, and the adjusted Lo-Mendell-Rubin test agrees at the decisive comparison. A stratum whose order matches the pooled order is not an order change however far the pooled order sits from four, so it never needs corroboration. BLRT non-rejection alone (the order matching the pooled order) is the stability claim.

Parameters:

• result (SearchResult) – The stratum’s search outcome.

• pooled_supported_k (int) – The pooled cohort’s supported order under the identical procedure.

Returns:

True when the order change is corroborated by all three, False otherwise (including when the order matches the pooled order).

Return type:

bool

Localisation and the sweep

Localisation schemes: how a proband is weighted into a set of local fits.

The stratified analysis (plan section 7) re-estimates the mixture model as a function of a continuous axis (age at diagnosis, or the derived calendar year of diagnosis). A local fit is one re-estimation, weighted toward one region of that axis. A LocalisationScheme turns the axis variable into an ordered set of such local fits.

The generalisation this module adds is that a local fit is a weight vector over the whole cohort, not a subset. Two schemes sit under one interface:

• HardBins wraps any BinningPolicy and emits one local fit per stratum with an indicator weight (one inside the bin, zero outside). Fitting on that weight is identical to fitting on the bin’s subset, so the frozen MaxEqualBins primary reproduces the current stratified fits exactly.

• KernelWindows emits one local fit per focal point on a grid, with a smooth kernel weight $w_i = exp(-(a_i - a_0)^2 / 2h^2)$ that falls off with distance from the focal age $a_0$. Every proband contributes to every nearby fit, weighted by relevance, so no proband is wasted at a boundary. This is the local-likelihood (LSEM) estimator: the class profiles come out as smooth trajectories along the axis rather than a handful of disjoint points.

Downstream (analysis.drift) consumes the fits through the method-independent StratumSummary, so alignment, distance, and the permutation null are unchanged by which scheme produced a fit.

class analysis.localise.LocalFit(label, position, weights)

One local re-estimation: where it sits on the axis and how each proband is weighted.

label

The fit’s name, unique within a scheme (a stratum name, or focal=6.5).

Type:

str

position

The fit’s location on the axis, in the axis units. For a hard bin this is the within-bin median; for a kernel window it is the focal point. It is the abscissa the drift trajectory is plotted and the continuous trend is regressed against.

Type:

float

weights

Per-proband weight over the whole cohort index, in [0, 1]. Zero excludes a proband from this fit; one includes it in full.

Type:

pandas.Series

class analysis.localise.LocalisationScheme(*args, **kwargs)

Turn a continuous axis variable into an ordered set of weighted local fits.

Downstream code is typed against this protocol, never a concrete scheme, so the analysis is independent of whether the fits are hard bins or kernel windows.

locales(values)

Return the ordered local fits for values (the per-proband axis variable).

spec()

Return the serialisable scheme specification for the run manifest.

class analysis.localise.HardBins(policy)

Adapt a BinningPolicy to a localisation scheme.

Emits one local fit per stratum with an indicator weight: one for probands in the bin, zero for the rest. Fitting on that weight is the same as fitting on the bin’s subset, so this reproduces the current stratified analysis while presenting it through the unified interface. The frozen primary is HardBins(MaxEqualBins(1000)).

policy

The binning policy that defines the strata.

Type:

analysis.strata.BinningPolicy

name

Scheme name recorded in the spec, derived from the policy’s name.

Type:

str

locales(values)

One indicator-weighted local fit per stratum, positioned at the within-bin median.

spec()

Return the specification (the scheme name and the wrapped policy’s spec).

analysis.localise.gaussian_weights(values, focal, bandwidth)

Return the Gaussian kernel weight of each value about a focal point.

The weight is $\exp(-(a_i - a_0)^2 / 2h^2)$ for value $a_i$, focal point $a_0$, and bandwidth $h$, so a value at the focal point has weight one and the weight falls off with distance. Missing values get weight zero.

analysis.localise.focal_grid(values, n_points, quantile_span)

Return n_points evenly spaced focal points spanning an inner quantile range.

The grid runs between the quantile_span quantiles of the finite values rather than the full range, so the outermost focal points still sit on a populated part of the axis and the edge fits are not estimated from almost nothing.

class analysis.localise.KernelWindows(bandwidth, grid, kernel='gaussian', quantile_span=(0.025, 0.975))

A local-likelihood (LSEM) scheme: one kernel-weighted fit per focal point.

bandwidth

The Gaussian kernel bandwidth $h$, in axis units, the smoothing knob. A larger bandwidth pulls in more probands per fit (smoother, more biased toward the pooled solution); a smaller one localises more sharply (less biased, noisier).

Type:

float

grid

The focal points, in axis units, or an integer count of evenly spaced points to build with focal_grid().

Type:

tuple of float or int

kernel

The kernel name; only gaussian is defined.

Type:

str

quantile_span

The inner quantile range the integer grid spans, ignored when grid is explicit.

Type:

tuple of float, optional

name

Scheme name recorded in the spec.

Type:

str

locales(values)

One kernel-weighted local fit per focal point.

spec()

Return the serialisable specification (bandwidth, kernel, and the focal grid).

analysis.localise.effective_sample_sizes(values, focal_points, bandwidth)

Return the effective sample size of a kernel fit at each focal point.

The effective sample size at a focal point $a_0$ is the sum of the Gaussian weights, $sum_i exp(-(a_i - a_0)^2 / 2h^2)$, the count of whole probands the weighted fit is worth there. It is what sets a local fit’s power, so it is the natural quantity to fix a bandwidth against. Missing values contribute nothing.

analysis.localise.bandwidth_for_effective_n(values, focal_points, target_n, *, reduce='min', tol=0.001, max_iter=60)

Return the smallest bandwidth whose focal fits reach a target effective sample size.

The effective sample size at a focal point rises with the bandwidth, so a target on it maps to a bandwidth by bisection. reduce sets which focal point the target must hold at: "min" fixes the thinnest focal point at the target, so every focal fit clears it (the same guarantee the hard-bin floor gives every bin); "median" fixes the typical focal point, which lets the edges fall below the target. The bandwidth is in the axis units.

Parameters:

• values (pandas.Series) – The axis variable (age at diagnosis in years, or diagnosis year).

• focal_points (list of float) – The focal grid the bandwidth is chosen for, as focal_grid() builds it.

• target_n (float) – The effective sample size each focal fit should reach (for example the recovery floor).

• reduce (str, optional) – "min" (default) or "median": which focal point the target is held at.

• tol (float and int, optional) – Bisection tolerance (in bandwidth units) and iteration cap.

• max_iter (float and int, optional) – Bisection tolerance (in bandwidth units) and iteration cap.

Returns:

The smallest bandwidth meeting the target.

Return type:

float

analysis.localise.permute_axis(values, seed)

Return the axis values shuffled across probands, the permutation null.

Breaking the pairing between a proband and its axis value, then re-running the same scheme on the shuffled values, removes only the association with the axis while preserving the fit sizes and the weight profile. For HardBins this gives same-size random bins (as analysis.drift.null_partition() does); for KernelWindows it gives kernel windows over randomly relabelled probands. One null definition serves both schemes. The seed is the permutation index, so a resumed null reproduces the same shuffles.

analysis.localise.fit_locale(matrix, typing, locale, *, n_init, random_state, weight_floor=0.0001, structural='covariate', progress_bar=0, verbose=0)

Fit the GFMM for one local fit, dropping negligible-weight probands.

Probands whose weight is at or below weight_floor carry no information and are dropped before fitting, so a hard-bin fit is an exact subset fit and a kernel fit skips the far tail. The remaining probands are passed to StepMix with their weights, so the fit is the weighted maximum-likelihood solution local to this fit’s region of the axis.

Parameters:

• matrix (analysis.cohort.CohortMatrix) – The full cohort matrix; the local fit selects and weights its rows.

• typing (analysis.features.Typing) – The reconciled feature typing.

• locale (LocalFit) – The local fit’s per-proband weights and label.

• n_init (int) – StepMix restarts.

• random_state (int) – Seed for reproducible restarts.

• weight_floor (float, optional) – Probands at or below this weight are excluded from the fit.

• progress_bar (int, optional) – StepMix verbosity, off by default so a sweep of many fits stays quiet.

• verbose (int, optional) – StepMix verbosity, off by default so a sweep of many fits stays quiet.

Returns:

The local fit, its labels on the retained probands, and its metrics.

Return type:

analysis.model.FitResult

The sweep: run one or more localisation schemes end to end and read them together.

The stratified analysis re-estimates the mixture model along an axis and measures how far the reference classes move (plan section 7). analysis.strata and analysis.localise give two ways to localise the re-estimation: hard bins (the frozen primary) and kernel windows (the local-likelihood, LSEM, trajectory). This module runs any such LocalisationScheme through the same machinery, so a single run fits the observed sweep, fits the matched permutation null, aligns each local fit to the reference, measures its drift, and reads it against the null. Because every scheme reduces to the method-independent StratumSummary, the hard-bin and kernel arms are read against one reference with one distance and reported side by side.

The heavy work (the fits) is separated from the cheap work (alignment and distance) exactly as analysis.drift does: fit_local_summaries() returns the stored summaries, and sweep_decision() measures over them, so a different alignment or distance re-measures with no refit.

class analysis.sweep.LocaleSummary(label, position, summary)

One local fit’s position on the axis and its method-independent summary.

label

The local fit’s name (a stratum name, or focal=6.5).

Type:

str

position

The local fit’s location on the axis, the abscissa for the trajectory.

Type:

float

summary

The centroids, dispersions, and reference contingency the drift is measured from.

Type:

analysis.drift.StratumSummary

analysis.sweep.summarise_locale(matrix, typing, locale, reference_labels, *, n_init, seed)

Fit one local fit and return its method-independent summary.

A hard bin (indicator weights that are one on its members) is summarised unweighted, so its summary is byte-identical to the current stratified fit’s. A kernel window (fractional weights) is summarised with those weights, so its centroids are the local weighted means.

analysis.sweep.summarise_local_worker(features, covariates, typing, dataset, version, weights, reference_labels, n_init, seed, structural='covariate')

Fit one local fit from its already-subset rows and return its summary (picklable).

A top-level function so it pickles for a process pool (mirrors analysis.drift.summarise_pseudo_stratum()). The caller drops the negligible-weight rows before submitting, so features holds only the retained probands and weights (or None for a hard bin, where every retained weight is one) is aligned to them. Returns None if the fit is degenerate (a singular covariate GLM), so the caller drops that local fit rather than letting one bad refit abort the sweep.

analysis.sweep.fit_local_summaries(matrix, typing, axis_values, scheme, reference_labels, *, n_init, seed)

Fit the observed sweep: one summary per local fit, in axis order.

Serial by design; the CLI spreads the fits across a process pool. Kept pure so the sweep logic is testable without the run machinery.

analysis.sweep.fit_null_summaries(matrix, typing, axis_values, scheme, reference_labels, *, n_init, n_permutations, seed)

Fit the permutation null: for each permutation, shuffle the axis and re-run the scheme.

Returns (locale_index, summary) pairs across all permutations, so the caller reads each observed local fit against the null draws at its own position. Serial; the CLI parallelises.

analysis.sweep.sweep_decision(scheme_name, observed, null, reference, aligner, distancer)

Align, measure, and read each observed local fit against the null, one row per class.

Mirrors the drift stage’s decision table so the hard-bin and kernel arms share a format: per local fit and reference class, the observed drift, its null 95th percentile and permutation p-value, the drift as a fraction of the between-class separation, the alignment confidence, and the Benjamini-Hochberg and reorganisation flags. The position and scheme columns let the arms be plotted as trajectories and compared.

analysis.sweep.parse_scheme(spec)

Build a localisation scheme from a colon-separated command-line specification.

Grammar (case-insensitive scheme name):

• hardbins:max-equal:<floor> gives HardBins(MaxEqualBins(floor)), the frozen primary.

• hardbins:quantile:<q> gives HardBins(QuantileBins(q)), the coarse sensitivity scheme.

• kernel:<bandwidth>:<n_points> gives KernelWindows(bandwidth, grid=n_points), the local-likelihood trajectory.

Parameters:

spec (str) – The scheme specification.

Returns:

The parsed scheme.

Return type:

analysis.localise.LocalisationScheme

Raises:

ValueError – If the scheme name or its arguments are not recognised.