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>)[source]#

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>, force=<typer.models.OptionInfo object>)[source]#

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>)[source]#

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>, force=<typer.models.OptionInfo object>)[source]#

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>, force=<typer.models.OptionInfo object>)[source]#

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>, force=<typer.models.OptionInfo object>)[source]#

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>)[source]#

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>)[source]#

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()[source]#

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

analysis.cli.stratify()[source]#

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

analysis.cli.drift()[source]#

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

analysis.cli.sensitivity()[source]#

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

analysis.cli.report()[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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

analysis.paths.stage_dir(root, stage)[source]#

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

analysis.paths.run_dir(root, stage, run_hash)[source]#

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

analysis.paths.manifest_path(run_dir)[source]#

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)[source]#

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)[source]#

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)[source]#

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

analysis.cache.file_digest(path)[source]#

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)[source]#

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()[source]#

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)[source]#

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)[source]#

Write a dataframe to Parquet, preserving its index.

analysis.cache.load_frame(path)[source]#

Read a dataframe from Parquet.

analysis.cache.save_model(obj, path)[source]#

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

analysis.cache.load_model(path)[source]#

Load a joblib-persisted object.

analysis.cache.save_json(obj, path)[source]#

Write obj to a human-readable JSON file.

analysis.cache.load_json(path)[source]#

Read a JSON file.

analysis.cache.read_manifest(run_dir)[source]#

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)[source]#

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>)[source]#

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)[source]#

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

analysis.run.run_context(stage, params, *, root=None, force=False)[source]#

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)[source]#

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()[source]#

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)[source]#

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()[source]#

Delete the log file if it exists.

analysis.checkpoint.clear_checkpoints(directory)[source]#

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)[source]#

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.

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)[source]#

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)[source]#

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()[source]#

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

supports_timing()[source]#

Return whether the backend can provide diagnosis-timing fields.

analysis.cohort.open_catalogue(root)[source]#

Open the dscat catalogue at the repository root.

analysis.cohort.source_csv(cat, root, dataset, version, table, role='')[source]#

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)[source]#

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

analysis.cohort.read_columns(path, columns)[source]#

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'))[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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()[source]#

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

integrate()[source]#

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)[source]#

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()[source]#

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

integrate()[source]#

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)[source]#

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()[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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.load_category_map(path)[source]#

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)[source]#

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)[source]#

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, progress_bar=1, verbose=1)[source]#

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.
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)[source]#

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)[source]#

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).

analysis.enrich.cohens_d(group, reference)[source]#

Return Cohen’s $d$ of a group against a reference, pooling their variances.

analysis.enrich.feature_enrichment(data, labels, n_classes=4)[source]#

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)[source]#

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)[source]#

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:

ci_low and ci_high (the percentile interval at level), median, the level, and n_valid (the resamples that yielded a defined correlation; a resample that empties or flattens a class is dropped, as in the permutation null).

Return type:

dict

analysis.enrich.contributory_features(enrichment, n_classes=4)[source]#

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)[source]#

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)[source]#

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')[source]#

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)[source]#

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)[source]#

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()[source]#

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

analysis.reference.align_to_named(signature, proportions)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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

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 continues from the first iteration that is missing. Because iteration i always uses seed base_seed + i, a resumed run reproduces what an uninterrupted run would have computed.

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).

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)[source]#

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>)[source]#

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)[source]#

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)[source]#

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'))[source]#

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)[source]#

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

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 continues from the first missing fit and the first missing comparison. 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.

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.

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)[source]#

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

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.

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)[source]#

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)[source]#

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).

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.

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)[source]#

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, 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)[source]#

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'))[source]#

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)[source]#

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)[source]#

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)[source]#

Partition values into ordered, named strata.

spec()[source]#

Return the serialisable policy specification.

class analysis.strata.FixedBands(edges, labels=None, name='fixed')[source]#

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)[source]#

Partition values at the fixed edges.

spec()[source]#

Return the serialisable specification, including the edges.

class analysis.strata.QuantileBins(q, name='quantile')[source]#

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)[source]#

Partition values at its interior quantiles.

spec()[source]#

Return the serialisable specification (the intended q).

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

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)[source]#

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

spec()[source]#

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

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)[source]#

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)[source]#

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)[source]#

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

analysis.strata_data.build_strata_data(root, version, index, age_at_eval, sex)[source]#

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)[source]#

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)[source]#

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)[source]#

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()[source]#

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))[source]#

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