Real-Data Trajectory Analysis

motco.analysis runs the measurement the simulation studies validated on your own data: pooled per-feature standardization, a stage-supervised PLS-DA latent space with the double-CV rank rule, delta / angle / shape between two groups in that space with the net-displacement orientation estimator, and RRPP p-values. analyze_trajectories is the Python entry point and motco analyze the command-line one; both write the same analysis artifact.

Use this entry point rather than assembling motco plsr + motco de by hand. Those commands are a different pipeline: they skip the M-value contract, choose the rank without the parsimony tie-break, measure angle with PC1 by default, and record neither the selected rank, the null summary nor the spectrum.

Input contract

• Blocks. One or more named omic blocks, each a samples × features table with string feature names. Blocks are concatenated in the order you declare them, and every joint feature is named "{block}__{feature}". Rows must equal the metadata's sample ids in the same order; every block needs at least one feature and only finite values.
• Units. Each block declares a unit transform: identity (the default) or beta_to_mvalue (the logit of methylation beta values; values must lie in the open interval (0, 1)). MOTCO does not normalize counts, filter features or impute values. Supply blocks in the continuous units you want measured, for example log-normalized pseudobulk expression. The record states each block's declared units.
• Metadata. A sample-id column, a group column with exactly two labels and a stage column, all named by you. Every group × stage cell must be observed.

• Stage order. Required and explicit. Net displacement runs from the first to the last declared stage, and shape, the spectrum and every table follow this order. Stage order is never inferred by sorting labels, because High < Intermediate < Low alphabetically. Group order is optional and defaults to the labels sorted as strings.

• Covariates (optional). covariate_cols names metadata columns, in model order; see Covariates.

Violations raise TrajectoryAnalysisError (a ValueError) naming the block, labels or cell at fault.

What is measured and recorded

The latent space is fitted once on all samples, with the rank chosen by double cross-validation (the modal rank across outer repeats, the smaller rank winning ties) from the candidate ranks at or above the rank floor, min_components=2 by default (CLI --min-components), so angle and shape are measured in a space of at least two dimensions. min_components=1 restores the unfloored 0.9.0 rule and reproduces every 0.9.0 record and table byte for byte; a floor below 1 is rejected before anything is fitted. The CV defaults are the ones the operating studies ran: cv2_splits=4, cv1_splits=3, n_repeats=5, max_components=20, random_state=1203, each clamped to what the smallest stage can support (max_components to at least the floor plus one; a sample that cannot satisfy the floor is an error naming the counts). forced_components fixes the rank and skips cross-validation for the dimensionality sensitivity analysis, ignoring the floor; the record then says component_selection: "forced".

record (written as result.json) holds:

Key	Contents
statistics, p_values	delta, angle, shape and their RRPP p-values (shape is null with two stages)
null_summary	per statistic: retained draws, mean, sd, 50/90/95/99 % quantiles of the permutation null
config_spectrum	eigenspectrum of the latent stage-mean configuration, pooled and per group, plus the permutation eigengap summary
latent	selected_lv, component_selection, CV settings (min_components beside them when the floor is above 1), cv_mean_auroc, feature counts, feature_order_signature
design	group and ordered stage levels, cell sizes per group × stage, contrast, analysed sample count
covariates	only when covariates are declared: names, kinds, factor levels and reference levels, design columns, the LS-mean reference row, missing policy, dropped samples, coefficients per latent component
blocks	each block's name, feature count, transform and units
measurement	orientation estimator, selected rank, permutations, seed, runtime
params, software	resolved parameters; MOTCO, numpy, scipy, scikit-learn, pandas and Python versions
attribution	the compact attribution record, or a not_requested marker

The selected rank sits beside the angle result because orientation is only interpretable relative to the space it was measured in. The null summary shows how wide the observed statistic's own null was, and the relative eigengap in the spectrum indicates whether a group's orientation is resolvable at all.

The floor guarantees two dimensions, not a second axis with signal. Under the unfloored rule a weak stage signal makes every candidate rank tie near an AUROC of 0.5 and parsimony selects rank 1, where the two groups' orientations are parallel or antiparallel by construction, angle is 0° or 180°, its RRPP null is a two-point distribution, and collinear stage means carry no bend for shape; both SEA-AD cohorts selected rank 1. With the floor, when the stage signal supports one component the second is the best of the remaining candidates, and angle's null widens accordingly. Read the record's config_spectrum relative eigengap and null_summary.angle together: they say how resolvable orientation is in each analysis. The floor was validated on a small-n construction calibrated to select rank 1 by scripts/rank_floor_check.py (rank_floor_decision.json), with the floor-1 arm reported beside it: the 2026-10-23 run returned adopt_floor at both calibrated strengths (report), with power at the null in both arms, so the floor defines angle and shape without making them powerful.

Equivalence to the simulation harness

analyze_trajectories and evaluate_semisynthetic_trajectory(..., integration_method="pls") call the same block-scaling cores and the same fit_pls_latent_space. On identical data, stage order and settings, the rank floor included, they return the same latent matrix, statistics, p-values, selected rank, null summary and spectrum. tests/test_analysis.py checks this for the default estimator, PC1 and a two-block layer selection.

Attribution

With attribution=AttributionDiagnosticSettings(enabled=True, ...) (CLI: --attribution), frozen-model orientation attribution runs on the fitted pooled estimator and the standardized joint matrix it was fitted to, with your stage order and group labels. The latent space is not refitted. Real data carry no generator truth, so the truth-recovery block is reported unavailable rather than invented. The full tables are written under attribution/.

Analysis artifact

write_analysis_outputs(result, out_dir) writes:

• result.json: the record above.
• latent_scores.csv: samples × pls_*, indexed by sample id.
• pls_loadings.csv: block, feature, pls_* (the PLS x_rotations_, in joint order).
• stage_means.csv: group, stage, pls_*, the fitted latent LS-mean vectors in declared stage order. Pass these to motco.viz.plot_trajectories.
• attribution/*.csv, when attribution ran.

It refuses to replace an existing result.json unless overwrite=True (CLI: --overwrite). The fitted scikit-learn estimator is not pickled, because pickles are version-fragile. The scores, loadings and recorded feature order signature are enough to audit and re-project the fit.

Covariates

Donor covariates such as age at death, PMI or RIN enter the trajectory regression in the latent space, so the stage means, the RRPP reduced model and therefore what is permuted all account for them. Declare them on the dataset:

dataset = TrajectoryDataset(
    ...,
    covariate_cols=["age_at_death", "pmi", "education_level"],
    categorical_covariates=["education_level"],   # optional; non-numeric columns are categorical anyway
    missing_covariates="complete_case",           # default "error"
)

CLI: --covariate age_at_death --covariate pmi --categorical-covariate education_level --missing-covariates complete_case.

• Coding. A numeric (non-boolean) column enters as one design column holding its values. Any other column, or one listed in categorical_covariates (for numeric-coded factors), enters as drop-first indicators name[level], with levels ordered as strings (numerically when every label is an integer); the first level is the reference.
• Model. The full model is group × stage + covariates and the reduced model group + stage + covariates. Covariates are additive: no interaction with group or stage. RRPP permutes the reduced model's residuals as before, so every permuted dataset keeps the fitted covariate effects.
• LS-means. Every group × stage mean is evaluated at the pooled mean of each covariate design column (for a factor, its observed level proportions). The covariate term then adds the same vector to every stage mean, and delta, angle and shape are translation-invariant, so the statistics do not depend on the reference. Only stage_means.csv does, and the record states the reference row.
• Missing values. missing_covariates="error" (the default) rejects any missing covariate value, naming the column and up to ten samples. "complete_case" drops those samples from the metadata and every block before standardization, the latent fit and the design, and records the dropped ids. A deletion that empties a group × stage cell is an error. Covariates are never imputed; blocks must still be complete.
• Guards. A covariate that is also the sample, group or stage column, a constant numeric covariate, a factor with one observed level, a non-finite numeric value, and a model with covariates that is not of full column rank (for example a covariate collinear with group × stage) are rejected with an error naming the covariate(s).

By default the latent space is covariate-blind, and the adjustment is still exact. PLS-DA is fitted on the standardized features against stage alone, exactly as without covariates, so the latent scores, loadings and selected rank do not change. OLS commutes with a fixed linear map: regressing the scores T = XW on [design, covariates] gives the coefficients of regressing X on the same design, projected by W. The stage means are therefore adjusted exactly as in feature space; only the basis W was chosen without the covariates. If a covariate effect is large and aligned with stage, PLS may spend a component on it; latent adjustment chooses the basis with the covariates instead.

Attribution is not available with covariates. Frozen-model attribution works from unadjusted cell means, so a request for both raises an error before any fitting.

Keep the covariate set parsimonious. Small cells (the SEA-AD MTG design has 9–28 donors per cell) and many covariate columns can make the model near-singular or the null unstable. The record's design_columns states the degrees of freedom used.

Type I validation. scripts/covariate_type_i_check.py measures size at the Phase 6 design with covariates confounded with group, stage and their interaction (adjusted vs unadjusted arms on a shared latent fit); see the report. Report covariate-adjusted results only under its covariates_validated decision.

Without covariates the record, the artifact and every output are identical to MOTCO 0.7.0 (except the runtime); the covariates key is absent and the schema version is unchanged (tests/test_analysis_identity.py).

Latent adjustment

With covariates declared, the PLS-DA basis is still fitted on the features alone. A large technical axis aligned with stage (library quality, for example) can then take a whole component, and at a selected rank of 1 or 2 that leaves little room for the biology. covariate_latent_adjustment (CLI: --latent-adjustment) fits the basis on features residualized on the declared covariates instead:

params = TrajectoryAnalysisParams(permutations=999, seed=1, covariate_latent_adjustment="stage")

motco analyze ... --covariate age_at_death --covariate pmi --latent-adjustment stage

• Conditional residualization. Every pooled-standardized joint feature is regressed on [1, D, C], where C holds the declared covariate design columns and D the conditioning design, and only the fitted covariate part C γ is subtracted. Under stage, D is the drop-first stage indicators; under group_stage, the group, stage and group × stage indicators. By Frisch–Waugh–Lovell, γ is the covariates' within-design effect, so stage signal that is merely correlated with a covariate stays in the features. Residualizing on the covariates alone would remove it: label-blind residualization upstream of MOTCO dropped the case study's held-out CV AUROC from 0.63 to 0.41. stage keeps the group label out of the basis; under group_stage the group enters only through the pooled within-cell slope. Use stage: it is the conditioning the Type I check admitted (see below).
• Leak-free cross-validation. Inside the double CV, every outer rest set and inner training set estimates its own γ and residualizes its held-out samples as X_held_out − C_held_out γ, so held-out labels never enter the residualization. Both CV levels stratify their folds on the conditioning cells (stages under stage, group × stage cells under group_stage), so every training set can estimate the conditioning design. Under stage the folds equal those of the unadjusted analysis. The final basis uses γ estimated on all samples.
• The pooled scale is kept. The residualized features are not re-standardized: every PLS fit on this path uses scale=False. A feature whose variance was mostly covariate is down-weighted, not inflated back to unit variance, and when the covariate effects are weak the basis approaches the unadjusted one.
• The raw features are measured in the new basis. The latent scores are the pooled-standardized features (not the residuals), centered at the residuals' mean and projected onto the basis. The trajectory model, the LS-means, the RRPP reduced model and the permutations are exactly those of the covariate-adjusted analysis above. With the covariates in both RRPP models, projecting the raw or the residualized features gives identical statistics, p-values and null draws, because the difference (C − C̄) γ W lies in the models' column space. So only the choice of basis (and of CV rank) changes. Measuring the raw features keeps the covariates block's coefficients and stage_means.csv on the same footing as without the setting.
• The same covariates in both places. The residualization covariates are exactly covariate_cols, which also enter the trajectory model; that equality is what makes the projection choice immaterial. To residualize on technical variables, declare them as covariates.
• Guards. The setting without covariates is an error before any fitting. The full-sample residualization design must have full column rank, otherwise the error names the covariates and the conditioning. Every CV training design must too, otherwise the error names the repeat and fold; this happens when a rare factor level drops out of a training set, and depends on random_state. Merge rare levels upstream rather than relying on a minimum-norm fit, whose held-out shift would be arbitrary. Attribution stays unavailable with covariates.
• Record. result.json gains a latent_adjustment block after covariates, and latent.integration_params.latent_adjustment names the conditioning. The block holds conditioning, covariates, design_columns (conditioning then covariate columns, beyond the intercept), n_samples, residual_df (n − 1 − k − q), cv_min_training_df (the smallest training-set residual df seen in CV; null with forced_components), feature_scale (pooled_standardized) and feature_r2. That last field summarizes each feature's partial R² of the covariates given the intercept and conditioning (5/25/50/75/95 % quantiles, mean, fraction above 0.5), the share of variance the residualization removes. params.covariate_latent_adjustment is recorded only when set.

Mediator limitation. Residualizing a variable that mediates the stage effect removes real signal, and no simulation can tell a confounder from a mediator. Choosing and justifying the covariate set is the analysis's job, not MOTCO's: a technical variable that is itself downstream of disease stage (a cell-quality metric that worsens with pathology, say) is a mediator candidate.

Validation. scripts/latent_adjustment_type_i_check.py measures size at the Phase 6 design under a covariate-aligned nuisance axis calibrated so that the unadjusted basis spends its first component on it (median |Spearman| ≥ 0.8). It runs three arms per replicate (unadjusted basis, stage, group_stage) at 8 declared covariates, and its decision file names the admissible conditionings (adopt_stage | adopt_group_stage | diagnose_latent_excess). The 2026-10-20 run returned adopt_stage (report). stage held delta 0.031, angle 0.067 and shape 0.058 against 0.0706. group_stage exceeded the bound on angle (0.082) and is not admissible. About 10 % of replicates lost an arm to a rare factor level. scripts/latent_adjustment_power_comparison.py compares power with and without such an axis against the unadjusted basis. Its 2026-10-20 run (report) returned candidate_fails for both conditionings. stage raised orientation power (angle 0.92–0.93 against 0.64–0.78 at the top effect) but lost shape power when the nuisance axis was present, and small-magnitude power without it. group_stage lost most orientation and shape power. Residualizing on covariates confounded with the design removes signal even when they have no effect. Report latent-adjusted results only with stage conditioning, as a sensitivity analysis beside the covariate-adjusted analysis rather than as a replacement, under the conditions of both reports.

Validation at the case study's configuration. scripts/case_study_type_i_check.py measures size where the SEA-AD case study is actually run: the rank floor of 2, six numeric covariate design columns (one confounded with group × stage and injected into the features, five pure noise), a stage signal weakened to a held-out CV AUROC of 0.69–0.70 and a nuisance the covariate-blind basis follows (median |Spearman| 0.92), with three arms per unit on the same data and permutation draws: covariates in the trajectory model only (adjusted, the primary analysis), none (unadjusted, descriptive) and the stage latent adjustment (latent_stage, the fallback). Its decision names primary_validated | adopt_latent_stage | diagnose_configuration_excess over the anchor and the translation control. The 2026-10-29 run (500 × 999 at δ = 0.1 and 0.05, conjunction) returned diagnose_configuration_excess (report): the adjusted model held every statistic at δ = 0.05 and exceeded the 0.0695 bound on translation delta at δ = 0.1 (0.072; delta 0.060–0.072 in every reading, 0.066 pooled), and the stage adjustment exceeded at the anchor on angle (0.084–0.090) and shape (0.078–0.092) at both strengths. At that configuration, read covariate-adjusted delta p-values against a realized size of about 0.066, and do not interpret stage-adjusted angle or shape rejections, until the diagnosis the report names is done.

Diagnosis of that excess. scripts/case_study_diagnosis.py reran the same units, with the committed adjusted and latent_stage p-values reproduced exactly (0 mismatches over 6,000 per strength), under arms that separate the candidate causes:

• on the covariate-blind basis, the confounded covariate estimated or its injection removed exactly, crossed with 0, 2 or 5 noise covariates;
• each of those rerun with the Freedman–Lane residuals rescaled by √(n / (n − p_red));
• the stage adjustment at a forced rank of 2.

The 2026-10-30 run (report) returned two verdicts:

• residual_df_shrinkage for the adjusted delta. RRPP's permuted residuals shrink with the reduced model's column count, so pooled delta rises from 0.0545 to 0.066 between one and six covariate columns at n = 80. The corrected permutation holds every reading.
• residualized_basis for the stage adjustment. It still exceeds on anchor angle at a fixed rank of 2 (0.084 / 0.082).

With many covariates at small n, read covariate-adjusted delta p-values against that realized size until the corrected permutation is offered and validated (roadmap, Next changes item 8).

Without the setting, every output is identical to MOTCO 0.8.0 except the runtime, record["params"] included (tests/test_analysis_identity.py).

Example: SEA-AD-style design

Two blocks (RNA and ATAC donor pseudobulks), sex as the group, three ADNC stages:

motco analyze \
  --block rna=rna_logcpm.csv --block atac=atac_logcpm.csv \
  --metadata donors.csv --sample-col donor_id --group-col Sex --stage-col ADNC \
  --stage-order "Not AD + Low,Intermediate,High" \
  --permutations 999 --seed 20261015 \
  --out-dir results/motco_mtg_astro

For the dimensionality sensitivity analysis, re-run with --forced-components 2, 3, … into separate output directories.

from motco.analysis import (
    OmicsBlock, TrajectoryAnalysisParams, TrajectoryDataset,
    analyze_trajectories, write_analysis_outputs,
)

dataset = TrajectoryDataset(
    blocks=[OmicsBlock("rna", rna), OmicsBlock("atac", atac)],
    metadata=donors,
    sample_col="donor_id", group_col="Sex", stage_col="ADNC",
    stage_order=["Not AD + Low", "Intermediate", "High"],
)
result = analyze_trajectories(dataset, TrajectoryAnalysisParams(permutations=999, seed=20261015))
write_analysis_outputs(result, "results/motco_mtg_astro")

API

motco.analysis

Real-data trajectory analysis: the validated MOTCO measurement on your own data.

:func:analyze_trajectories runs the production pipeline the simulation harness characterizes — pooled per-feature standardization, a stage-supervised PLS-DA latent space whose rank is chosen by double cross-validation, delta / angle / shape between two groups in that space, and RRPP p-values — on caller-supplied omic blocks and sample metadata. It calls the same cores as :func:motco.simulations.evaluate_semisynthetic_trajectory (block scaling in :mod:motco.simulations.preprocessing, :func:~motco.simulations.evaluation.fit_pls_latent_space, the trajectory design and the record helpers), so the measurement exists once.

Two inputs are load-bearing and therefore explicit:

• Stage order. Net displacement runs from the first to the last declared stage, and shape and the spectrum follow that order. It is never inferred by sorting labels.
• Units. Each block declares its unit transform (identity or beta_to_mvalue). MOTCO does not normalize counts, filter features or impute values; blocks are analysed in the units supplied.

Declared covariates enter the trajectory model. By default the latent space is built without them; TrajectoryAnalysisParams.covariate_latent_adjustment ("stage" or "group_stage") instead fits the basis on features residualized on the covariates conditional on that design, and measures the raw features in it.

:func:write_analysis_outputs writes the result as plain JSON/CSV (record, latent scores, PLS loadings, latent stage means and, when requested, the attribution tables) — the fitted estimator is never pickled.

TrajectoryAnalysisError

Bases: ValueError

Raised when real-data analysis inputs or parameters are invalid.

Source code in src/motco/analysis.py

class TrajectoryAnalysisError(ValueError):
    """Raised when real-data analysis inputs or parameters are invalid."""

OmicsBlock dataclass

One named omic block: a samples × features table and its unit transform.

Source code in src/motco/analysis.py

@dataclass(frozen=True)
class OmicsBlock:
    """One named omic block: a samples × features table and its unit transform."""

    name: str
    values: pd.DataFrame
    transform: UnitTransform = "identity"

TrajectoryDataset dataclass

Aligned omic blocks plus sample metadata, validated on construction.

Every block's rows must equal metadata[sample_col] in order. The group column must hold exactly two labels; stage_order must list every observed stage label exactly once (labels are compared as strings), and every group × stage cell must be observed.

covariate_cols names metadata columns that enter the trajectory model additively, in declared order. A covariate is categorical when listed in categorical_covariates or when its dtype is not numeric (booleans are categorical); otherwise numeric. missing_covariates decides what a missing covariate value does: "error" rejects it, "complete_case" drops the sample from the metadata and every block before anything else is validated or fitted (dropped_sample_ids lists them). Covariates are never imputed.

Source code in src/motco/analysis.py

@dataclass(frozen=True)
class TrajectoryDataset:
    """Aligned omic blocks plus sample metadata, validated on construction.

    Every block's rows must equal ``metadata[sample_col]`` in order. The group
    column must hold exactly two labels; ``stage_order`` must list every
    observed stage label exactly once (labels are compared as strings), and
    every group × stage cell must be observed.

    ``covariate_cols`` names metadata columns that enter the trajectory model
    additively, in declared order. A covariate is categorical when listed in
    ``categorical_covariates`` or when its dtype is not numeric (booleans are
    categorical); otherwise numeric. ``missing_covariates`` decides what a
    missing covariate value does: ``"error"`` rejects it, ``"complete_case"``
    drops the sample from the metadata and every block before anything else is
    validated or fitted (``dropped_sample_ids`` lists them). Covariates are
    never imputed.
    """

    blocks: Sequence[OmicsBlock]
    metadata: pd.DataFrame
    sample_col: str
    group_col: str
    stage_col: str
    stage_order: Sequence[str]
    group_order: Sequence[str] | None = None
    covariate_cols: Sequence[str] = ()
    categorical_covariates: Sequence[str] = ()
    missing_covariates: MissingCovariatePolicy = "error"
    _dropped_sample_ids: tuple[str, ...] = field(default=(), init=False, repr=False, compare=False)
    _covariate_design: CovariateDesign | None = field(default=None, init=False, repr=False, compare=False)

    def __post_init__(self) -> None:
        object.__setattr__(self, "blocks", tuple(self.blocks))
        object.__setattr__(self, "stage_order", tuple(str(stage) for stage in self.stage_order))
        if self.group_order is not None:
            object.__setattr__(self, "group_order", tuple(str(group) for group in self.group_order))
        object.__setattr__(self, "covariate_cols", tuple(self.covariate_cols))
        object.__setattr__(self, "categorical_covariates", tuple(self.categorical_covariates))
        self._validate_metadata()
        self._validate_covariate_declarations()
        cells_before = self._observed_cells()
        self._apply_missing_covariate_policy()
        self._validate_blocks(self.sample_ids)
        self._build_covariate_design()
        self._validate_design(cells_before)

    @property
    def sample_ids(self) -> list[str]:
        return self.metadata[self.sample_col].astype(str).tolist()

    @property
    def groups(self) -> tuple[str, ...]:
        """Group labels in analysis order (declared, else sorted as strings)."""

        if self.group_order is not None:
            return tuple(self.group_order)
        return tuple(sorted(set(self.metadata[self.group_col].astype(str))))

    @property
    def stages(self) -> tuple[str, ...]:
        return tuple(self.stage_order)

    @property
    def dropped_sample_ids(self) -> tuple[str, ...]:
        """Samples removed by complete-case deletion (empty under ``"error"``)."""

        return self._dropped_sample_ids

    @property
    def covariate_design(self) -> CovariateDesign | None:
        """The covariate design columns on the retained samples, or ``None`` without covariates."""

        return self._covariate_design

    def _observed_cells(self) -> set[tuple[str, str]]:
        return set(zip(self.metadata[self.group_col].astype(str), self.metadata[self.stage_col].astype(str)))

    def _validate_covariate_declarations(self) -> None:
        if self.missing_covariates not in MISSING_COVARIATE_POLICIES:
            raise TrajectoryAnalysisError(
                f"Unknown missing_covariates policy {self.missing_covariates!r}; allowed: "
                f"{list(MISSING_COVARIATE_POLICIES)}."
            )
        covariates = list(self.covariate_cols)
        repeated = sorted({name for name in covariates if covariates.count(name) > 1})
        if repeated:
            raise TrajectoryAnalysisError(f"Covariates must be declared once; repeated: {repeated}.")
        reserved = [name for name in covariates if name in {self.sample_col, self.group_col, self.stage_col}]
        if reserved:
            raise TrajectoryAnalysisError(
                f"Covariates must differ from the sample, group and stage columns; got {reserved}."
            )
        absent = [name for name in covariates if name not in self.metadata.columns]
        if absent:
            raise TrajectoryAnalysisError(f"metadata is missing covariate column(s): {absent}.")
        undeclared = [name for name in self.categorical_covariates if name not in covariates]
        if undeclared:
            raise TrajectoryAnalysisError(
                f"categorical_covariates must also be declared in covariate_cols; not declared: {undeclared}."
            )

    def _apply_missing_covariate_policy(self) -> None:
        if not self.covariate_cols:
            return
        missing = self.metadata[list(self.covariate_cols)].isna()
        if not missing.to_numpy().any():
            return
        sample_ids = self.metadata[self.sample_col].astype(str)
        if self.missing_covariates == "error":
            column = next(name for name in self.covariate_cols if missing[name].any())
            affected = sample_ids[missing[column].to_numpy()].tolist()
            raise TrajectoryAnalysisError(
                f"Covariate {column!r} has missing values for {len(affected)} sample(s): {affected[:10]}. "
                "Pass missing_covariates='complete_case' to drop them; covariates are never imputed."
            )
        drop = missing.any(axis=1).to_numpy()
        keep = ~drop
        object.__setattr__(self, "_dropped_sample_ids", tuple(sample_ids[drop].tolist()))
        object.__setattr__(self, "metadata", self.metadata.loc[keep].reset_index(drop=True))
        retained = []
        for block in self.blocks:
            values = block.values
            if isinstance(values, pd.DataFrame) and values.shape[0] == keep.size:
                values = values.loc[keep]
            retained.append(OmicsBlock(block.name, values, block.transform))
        object.__setattr__(self, "blocks", tuple(retained))

    def _build_covariate_design(self) -> None:
        if not self.covariate_cols:
            return
        try:
            design = build_covariate_matrix(self.metadata, self.covariate_cols, self.categorical_covariates)
        except ValueError as exc:
            raise TrajectoryAnalysisError(str(exc)) from exc
        object.__setattr__(self, "_covariate_design", design)

    def _validate_metadata(self) -> None:
        if not isinstance(self.metadata, pd.DataFrame):
            raise TrajectoryAnalysisError("metadata must be a pandas DataFrame.")
        columns = {"sample_col": self.sample_col, "group_col": self.group_col, "stage_col": self.stage_col}
        missing = [f"{role}={name!r}" for role, name in columns.items() if name not in self.metadata.columns]
        if missing:
            raise TrajectoryAnalysisError(f"metadata is missing column(s): {', '.join(missing)}.")
        if len({self.sample_col, self.group_col, self.stage_col}) != 3:
            raise TrajectoryAnalysisError("sample_col, group_col and stage_col must name three different columns.")
        if self.metadata.shape[0] == 0:
            raise TrajectoryAnalysisError("metadata must contain at least one sample.")
        for name in columns.values():
            if self.metadata[name].isna().any():
                raise TrajectoryAnalysisError(f"metadata column {name!r} contains missing values.")
        sample_ids = self.sample_ids
        duplicated = sorted({sid for sid in sample_ids if sample_ids.count(sid) > 1})
        if duplicated:
            raise TrajectoryAnalysisError(f"metadata sample ids are not unique: {duplicated[:10]}.")

    def _validate_blocks(self, sample_ids: list[str]) -> None:
        if not self.blocks:
            raise TrajectoryAnalysisError("At least one omic block is required.")
        names = [block.name for block in self.blocks]
        repeated = sorted({name for name in names if names.count(name) > 1})
        if repeated:
            raise TrajectoryAnalysisError(f"Block names must be unique; repeated: {repeated}.")
        for block in self.blocks:
            if not isinstance(block.name, str) or not block.name:
                raise TrajectoryAnalysisError(f"Block names must be non-empty strings, got {block.name!r}.")
            if block.transform not in UNIT_TRANSFORMS:
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} declares unknown transform {block.transform!r}; "
                    f"allowed: {list(UNIT_TRANSFORMS)}."
                )
            values = block.values
            if not isinstance(values, pd.DataFrame):
                raise TrajectoryAnalysisError(f"Block {block.name!r} must be a pandas DataFrame.")
            if values.shape[0] != len(sample_ids):
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} has {values.shape[0]} rows but metadata has {len(sample_ids)} samples."
                )
            if values.index.astype(str).tolist() != sample_ids:
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} rows are not aligned to metadata column {self.sample_col!r} "
                    "(same sample ids in the same order are required)."
                )
            if values.shape[1] == 0:
                raise TrajectoryAnalysisError(f"Block {block.name!r} must contain at least one feature.")
            features = values.columns.astype(str).tolist()
            if len(set(features)) != len(features):
                raise TrajectoryAnalysisError(f"Block {block.name!r} has duplicate feature names.")
            try:
                matrix = values.to_numpy(dtype=float)
            except (TypeError, ValueError) as exc:
                raise TrajectoryAnalysisError(f"Block {block.name!r} holds non-numeric values.") from exc
            if not np.all(np.isfinite(matrix)):
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} contains missing or infinite values; complete matrices are "
                    "required and no imputation is performed."
                )
            if block.transform == "beta_to_mvalue" and not np.all((matrix > 0.0) & (matrix < 1.0)):
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} declares beta_to_mvalue but holds values outside the open "
                    "interval (0, 1)."
                )

    def _validate_design(self, cells_before: set[tuple[str, str]]) -> None:
        emptied = sorted(cells_before - self._observed_cells())
        if emptied:
            raise TrajectoryAnalysisError(
                f"Complete-case deletion of {len(self._dropped_sample_ids)} sample(s) with missing covariates "
                f"emptied group/stage combination(s): {emptied}."
            )
        observed_stages = sorted(set(self.metadata[self.stage_col].astype(str)))
        declared = list(self.stage_order)
        if len(set(declared)) != len(declared) or set(declared) != set(observed_stages):
            raise TrajectoryAnalysisError(
                f"stage_order must list every observed stage label exactly once; declared {declared}, "
                f"observed {observed_stages}."
            )
        if len(declared) < 2:
            raise TrajectoryAnalysisError(f"At least two stages are required; observed {observed_stages}.")
        observed_groups = sorted(set(self.metadata[self.group_col].astype(str)))
        if len(observed_groups) != 2:
            raise TrajectoryAnalysisError(
                f"Exactly two groups are required in {self.group_col!r}; found {len(observed_groups)}: "
                f"{observed_groups}."
            )
        if self.group_order is not None:
            group_order = list(self.group_order)
            if len(set(group_order)) != len(group_order) or set(group_order) != set(observed_groups):
                raise TrajectoryAnalysisError(
                    f"group_order must list both observed groups exactly once; declared {group_order}, "
                    f"observed {observed_groups}."
                )
        present = self._observed_cells()
        missing = [(group, stage) for group in self.groups for stage in declared if (group, stage) not in present]
        if missing:
            raise TrajectoryAnalysisError(f"metadata has no sample for group/stage combination(s): {missing}.")

groups property

Group labels in analysis order (declared, else sorted as strings).

dropped_sample_ids property

Samples removed by complete-case deletion (empty under "error").

covariate_design property

The covariate design columns on the retained samples, or None without covariates.

TrajectoryAnalysisParams dataclass

Measurement settings; the defaults are the simulation harness's.

forced_components fixes the latent rank (dimensionality sensitivity) and skips cross-validation; None selects it by double CV. seed seeds RRPP exactly as the harness does.

covariate_latent_adjustment (None, "stage" or "group_stage") opts into covariate-aware latent construction: the basis is fitted on features residualized on the declared covariates, conditional on the stage or the group × stage indicators. It requires declared covariates and is omitted from :meth:to_dict while None.

min_components is the rank floor of the cross-validated latent space (default 2, so angle and shape are measured in at least two dimensions). 1 restores the unfloored 0.9.0 rule and normalizes to None, which is omitted from :meth:to_dict and from the record's latent block, so a floor-1 analysis reproduces every 0.9.0 artifact. A floor below 1 is rejected before any fitting. The floor guarantees two dimensions, not a second axis with signal.

Source code in src/motco/analysis.py

@dataclass(frozen=True)
class TrajectoryAnalysisParams:
    """Measurement settings; the defaults are the simulation harness's.

    ``forced_components`` fixes the latent rank (dimensionality sensitivity) and
    skips cross-validation; ``None`` selects it by double CV. ``seed`` seeds
    RRPP exactly as the harness does.

    ``covariate_latent_adjustment`` (``None``, ``"stage"`` or
    ``"group_stage"``) opts into covariate-aware latent construction: the
    basis is fitted on features residualized on the declared covariates,
    conditional on the stage or the group × stage indicators. It requires
    declared covariates and is omitted from :meth:`to_dict` while ``None``.

    ``min_components`` is the rank floor of the cross-validated latent space
    (default 2, so ``angle`` and ``shape`` are measured in at least two
    dimensions). ``1`` restores the unfloored 0.9.0 rule and normalizes to
    ``None``, which is omitted from :meth:`to_dict` and from the record's
    ``latent`` block, so a floor-1 analysis reproduces every 0.9.0 artifact.
    A floor below 1 is rejected before any fitting. The floor guarantees two
    dimensions, not a second axis with signal.
    """

    permutations: int = 0
    seed: int | None = None
    n_jobs: int = 1
    orientation_estimator: str = "net_displacement"
    forced_components: int | None = None
    cv1_splits: int = 3
    cv2_splits: int = 4
    n_repeats: int = 5
    max_components: int = 20
    random_state: int = 1203
    attribution: AttributionDiagnosticSettings = field(default_factory=AttributionDiagnosticSettings)
    progress: bool = False
    covariate_latent_adjustment: str | None = None
    min_components: int | None = 2

    def __post_init__(self) -> None:
        floor = self.min_components
        if floor is not None:
            if isinstance(floor, bool) or not isinstance(floor, (int, np.integer)) or floor < 1:
                raise TrajectoryAnalysisError(f"min_components must be an integer of at least 1, got {floor!r}.")
            # 1 is the unfloored 0.9.0 rule: held as ``None`` so the record stays byte-identical.
            object.__setattr__(self, "min_components", None if int(floor) == 1 else int(floor))
        if self.covariate_latent_adjustment is not None and self.covariate_latent_adjustment not in CONDITIONING_KINDS:
            raise TrajectoryAnalysisError(
                f"Unknown covariate_latent_adjustment {self.covariate_latent_adjustment!r}; allowed: "
                f"{list(CONDITIONING_KINDS)} or None."
            )
        if self.orientation_estimator not in EVALUATION_ORIENTATION_ESTIMATORS:
            raise TrajectoryAnalysisError(
                f"Unknown orientation_estimator {self.orientation_estimator!r}; allowed: "
                f"{list(EVALUATION_ORIENTATION_ESTIMATORS)}."
            )
        if self.permutations < 0:
            raise TrajectoryAnalysisError("permutations must be greater than or equal to 0.")
        if self.n_jobs == 0:
            raise TrajectoryAnalysisError("n_jobs must be -1 or a non-zero integer.")
        settings = self.attribution
        if settings.enabled:
            if settings.bootstrap_replicates < 0:
                raise TrajectoryAnalysisError("attribution.bootstrap_replicates must be greater than or equal to 0.")
            if settings.top_k < 1:
                raise TrajectoryAnalysisError("attribution.top_k must be at least 1.")

    def to_dict(self) -> dict[str, Any]:
        data = asdict(self)
        if self.covariate_latent_adjustment is None:
            # Absent means the 0.8.0 analysis: its resolved parameters stay byte-identical.
            data.pop("covariate_latent_adjustment")
        if self.min_components is None:
            # Floor 1 means the 0.9.0 rank rule: its resolved parameters stay byte-identical.
            data.pop("min_components")
        return data

    @property
    def resolved_min_components(self) -> int:
        """The rank floor the cross-validation runs with (1 when omitted)."""

        return self.min_components or 1

resolved_min_components property

The rank floor the cross-validation runs with (1 when omitted).

TrajectoryAnalysisResult dataclass

The analysis record plus the fitted representation.

record is JSON-safe. latent_scores are the per-sample PLS scores, pls_loadings the joint-feature rotations (block/feature columns first), stage_means the fitted latent LS-mean vector per group × stage in declared stage order, and attribution the full attribution result when it was requested.

Source code in src/motco/analysis.py

@dataclass(frozen=True)
class TrajectoryAnalysisResult:
    """The analysis record plus the fitted representation.

    ``record`` is JSON-safe. ``latent_scores`` are the per-sample PLS scores,
    ``pls_loadings`` the joint-feature rotations (``block``/``feature`` columns
    first), ``stage_means`` the fitted latent LS-mean vector per group × stage
    in declared stage order, and ``attribution`` the full attribution result
    when it was requested.
    """

    record: dict[str, Any]
    latent_scores: pd.DataFrame
    pls_loadings: pd.DataFrame
    stage_means: pd.DataFrame
    attribution: OrientationAttributionResult | None = None

    def to_json(self, **kwargs: Any) -> str:
        return json.dumps(self.record, **kwargs)

analyze_trajectories(dataset, params=None)

Measure two groups' trajectories across ordered stages on real data.

Standardizes every block with pooled per-feature scalers (after its unit transform), concatenates the blocks in declared order as f"{block}__{feature}", builds the stage-supervised PLS-DA latent space (double-CV rank, or params.forced_components), estimates delta, angle (net displacement by default) and shape between the groups in that space and, when params.permutations > 0, their RRPP p-values. Records the selected rank, null summary, stage-mean spectrum, cell sizes, block units and software versions beside the statistics. Optionally runs frozen-model orientation attribution on the fitted estimator.

Declared covariates (dataset.covariate_cols) are appended additively to the full (group × stage) and reduced (group + stage) trajectory models, so RRPP permutes residuals of the reduced model with covariates. Every LS-mean is evaluated at the pooled covariate mean; delta, angle and shape do not depend on that reference. The latent space is fitted without the covariates unless params.covariate_latent_adjustment is set: then the pooled-standardized features are residualized on the covariates conditional on the stage ("stage") or group × stage ("group_stage") indicators before the basis is fitted (per training set inside the cross-validation, at the pooled scale), the raw features are projected onto that basis, and the trajectory model, LS-means and RRPP are unchanged. Attribution is not available with covariates.

Source code in src/motco/analysis.py

def analyze_trajectories(
    dataset: TrajectoryDataset,
    params: TrajectoryAnalysisParams | None = None,
) -> TrajectoryAnalysisResult:
    """Measure two groups' trajectories across ordered stages on real data.

    Standardizes every block with pooled per-feature scalers (after its unit
    transform), concatenates the blocks in declared order as
    ``f"{block}__{feature}"``, builds the stage-supervised PLS-DA latent space
    (double-CV rank, or ``params.forced_components``), estimates ``delta``,
    ``angle`` (net displacement by default) and ``shape`` between the groups
    in that space and, when ``params.permutations > 0``, their RRPP p-values.
    Records the selected rank, null summary, stage-mean spectrum, cell sizes,
    block units and software versions beside the statistics. Optionally runs
    frozen-model orientation attribution on the fitted estimator.

    Declared covariates (``dataset.covariate_cols``) are appended additively to
    the full (group × stage) and reduced (group + stage) trajectory models, so
    RRPP permutes residuals of the reduced model with covariates. Every LS-mean
    is evaluated at the pooled covariate mean; ``delta``, ``angle`` and
    ``shape`` do not depend on that reference. The latent space is fitted
    without the covariates unless ``params.covariate_latent_adjustment`` is
    set: then the pooled-standardized features are residualized on the
    covariates conditional on the stage (``"stage"``) or group × stage
    (``"group_stage"``) indicators before the basis is fitted (per training set
    inside the cross-validation, at the pooled scale), the raw features are
    projected onto that basis, and the trajectory model, LS-means and RRPP are
    unchanged. Attribution is not available with covariates.
    """

    params = params or TrajectoryAnalysisParams()
    covariates = dataset.covariate_design
    if params.covariate_latent_adjustment is not None and covariates is None:
        raise TrajectoryAnalysisError(
            f"Latent adjustment ({params.covariate_latent_adjustment!r}) requires declared covariates; "
            "declare covariate_cols on the dataset or drop the setting."
        )
    if covariates is not None and params.attribution.enabled:
        raise TrajectoryAnalysisError(
            "Attribution is not available with covariates: it works from unadjusted cell means. "
            "Disable attribution or drop the covariates."
        )
    start = perf_counter()
    blocks = {block.name: block.values for block in dataset.blocks}
    order = [block.name for block in dataset.blocks]
    transforms = {block.name: block.transform for block in dataset.blocks}

    scalers = fit_block_scalers(blocks, transforms)
    standardized = transform_blocks(blocks, scalers, transforms)
    X = concatenate_blocks(standardized, order).reset_index(drop=True)

    groups = list(dataset.groups)
    stages = list(dataset.stages)
    group_codes = {label: str(index) for index, label in enumerate(groups)}
    stage_codes = {label: str(index) for index, label in enumerate(stages)}
    working = pd.DataFrame(
        {
            _GROUP: dataset.metadata[dataset.group_col].astype(str).map(group_codes).to_numpy(),
            _STAGE: dataset.metadata[dataset.stage_col].astype(str).map(stage_codes).to_numpy(),
        }
    )

    integration_params: dict[str, Any] = {
        "cv1_splits": params.cv1_splits,
        "cv2_splits": params.cv2_splits,
        "n_repeats": params.n_repeats,
        "max_components": params.max_components,
        "random_state": params.random_state,
        "n_jobs": params.n_jobs,
    }
    if params.forced_components is not None:
        integration_params["forced_components"] = params.forced_components
    nuisance = _latent_nuisance(working, covariates, params.covariate_latent_adjustment)
    sample_ids = dataset.sample_ids
    try:
        latent = fit_pls_latent_space(
            X,
            working[_STAGE],
            integration_params,
            stage_col=dataset.stage_col,
            sample_ids=sample_ids,
            layer_feature_counts={name: int(frame.shape[1]) for name, frame in blocks.items()},
            layers_metadata={},
            original_scales=_joint_original_scales(scalers, order, X.columns),
            methylation_units=_methylation_units(transforms),
            nuisance=nuisance,
            min_components=params.resolved_min_components,
        )
        design = build_simulation_trajectory_design(working, group_col=_GROUP, stage_col=_STAGE)
    except SimulationEvaluationError as exc:
        raise TrajectoryAnalysisError(str(exc)) from exc
    group_levels = [groups[int(code)] for code in design.group_levels]
    stage_levels = [stages[int(code)] for code in design.stage_levels]
    if stage_levels != stages or group_levels != groups:
        raise TrajectoryAnalysisError(
            f"Internal level mapping mismatch: design levels {group_levels}/{stage_levels}, "
            f"declared {groups}/{stages}."
        )
    model_full, model_reduced, ls_means = _covariate_adjusted_design(
        design.model_full, design.model_reduced, design.ls_means, covariates
    )

    observed_deltas, observed_angles, observed_shapes, observed_spectra = estimate_difference(
        Y=latent.matrix,
        model_matrix=model_full,
        LS_means=ls_means,
        contrast=design.contrast,
        return_spectra=True,
        orientation_estimator=params.orientation_estimator,
    )
    config_spectrum = _config_spectrum_block(observed_spectra, group_levels)
    shape_available = len(stage_levels) >= 3
    pair_statistics = _extract_pair_statistics(
        observed_deltas, observed_angles, observed_shapes, shape_available=shape_available
    )

    p_values: dict[str, float] = {}
    null_summary: dict[str, dict[str, float]] = {}
    if params.permutations > 0:
        dist_delta, dist_angle, dist_shape, dist_eigengap = RRPP(
            latent.matrix,
            model_full,
            model_reduced,
            ls_means,
            design.contrast,
            permutations=params.permutations,
            n_jobs=params.n_jobs,
            progress=params.progress,
            seed=params.seed,
            return_eigengaps=True,
            orientation_estimator=params.orientation_estimator,
        )
        null_values = _extract_null_distributions(
            dist_delta, dist_angle, dist_shape, shape_available=shape_available
        )
        p_values = {
            statistic: _empirical_p_value(null_values[statistic], observed)
            for statistic, observed in pair_statistics.items()
            if statistic in null_values and np.isfinite(observed)
        }
        null_summary = _summarize_null_distributions(null_values)
        config_spectrum["permutation_pooled_eigengap"] = _summarize_permutation_eigengaps(dist_eigengap)

    latent_scores = latent.matrix.copy()
    latent_scores.index.name = dataset.sample_col
    artifacts = latent.artifacts
    assert artifacts is not None  # the PLS core always returns its fitted artifacts
    pls_loadings = _loadings_frame(artifacts.model, X.columns, list(latent.matrix.columns))
    betas = np.asarray(estimate_betas(model_full, latent.matrix.to_numpy(dtype=float)), dtype=float)
    stage_means = _stage_means_frame(
        ls_means,
        betas,
        group_levels,
        stage_levels,
        list(latent.matrix.columns),
    )

    attribution_record, attribution_result = _run_attribution(dataset, params, latent, X, transforms)
    cell_sizes = {
        group: {
            stage: int(
                (
                    (dataset.metadata[dataset.group_col].astype(str) == group)
                    & (dataset.metadata[dataset.stage_col].astype(str) == stage)
                ).sum()
            )
            for stage in stage_levels
        }
        for group in group_levels
    }
    record: dict[str, Any] = {
        "schema_version": ANALYSIS_SCHEMA_VERSION,
        "statistics": pair_statistics,
        "p_values": p_values,
        "null_summary": null_summary,
        "config_spectrum": config_spectrum,
        "latent": latent.metadata,
        "design": {
            "group_col": dataset.group_col,
            "stage_col": dataset.stage_col,
            "sample_col": dataset.sample_col,
            "group_levels": group_levels,
            "stage_levels": stage_levels,
            "cell_sizes": cell_sizes,
            "contrast": design.contrast,
            "n_samples": len(sample_ids),
        },
        **(
            {"covariates": _covariate_record(dataset, covariates, betas, list(latent.matrix.columns))}
            if covariates is not None
            else {}
        ),
        **(
            {"latent_adjustment": _latent_adjustment_record(dataset, artifacts.latent_adjustment)}
            if artifacts.latent_adjustment is not None
            else {}
        ),
        "blocks": [
            {
                "name": block.name,
                "n_features": int(block.values.shape[1]),
                "transform": block.transform,
                "units": _TRANSFORM_UNITS[block.transform],
            }
            for block in dataset.blocks
        ],
        "measurement": {
            "orientation_estimator": params.orientation_estimator,
            "selected_lv": latent.metadata["selected_lv"],
            "component_selection": latent.metadata["component_selection"],
            "permutations": params.permutations,
            "seed": params.seed,
            "n_jobs": params.n_jobs,
            "shape_available": shape_available,
            "runtime_seconds": perf_counter() - start,
        },
        "params": params.to_dict(),
        "software": _software_versions(),
        "attribution": attribution_record,
    }
    return TrajectoryAnalysisResult(
        record=_json_safe(record),
        latent_scores=latent_scores,
        pls_loadings=pls_loadings,
        stage_means=stage_means,
        attribution=attribution_result,
    )

write_analysis_outputs(result, out_dir, *, overwrite=False)

Write the analysis artifact to out_dir and return the written paths.

Writes result.json, latent_scores.csv, pls_loadings.csv, stage_means.csv and, when attribution ran, its tables under attribution/. Refuses to touch a directory that already holds result.json unless overwrite is set.

Source code in src/motco/analysis.py

def write_analysis_outputs(
    result: TrajectoryAnalysisResult,
    out_dir: str | Path,
    *,
    overwrite: bool = False,
) -> dict[str, Path]:
    """Write the analysis artifact to ``out_dir`` and return the written paths.

    Writes ``result.json``, ``latent_scores.csv``, ``pls_loadings.csv``,
    ``stage_means.csv`` and, when attribution ran, its tables under
    ``attribution/``. Refuses to touch a directory that already holds
    ``result.json`` unless ``overwrite`` is set.
    """

    directory = Path(out_dir)
    result_path = directory / "result.json"
    if result_path.exists() and not overwrite:
        raise FileExistsError(f"{result_path} already exists; pass overwrite=True to replace the analysis.")
    directory.mkdir(parents=True, exist_ok=True)
    paths: dict[str, Path] = {}
    result_path.write_text(json.dumps(result.record, indent=2, allow_nan=False) + "\n", encoding="utf-8")
    paths["result"] = result_path
    paths["latent_scores"] = directory / "latent_scores.csv"
    result.latent_scores.to_csv(paths["latent_scores"])
    paths["pls_loadings"] = directory / "pls_loadings.csv"
    result.pls_loadings.to_csv(paths["pls_loadings"], index=False)
    paths["stage_means"] = directory / "stage_means.csv"
    result.stage_means.to_csv(paths["stage_means"], index=False)
    if result.attribution is not None:
        for name, path in write_attribution_outputs(result.attribution, directory / "attribution").items():
            paths[f"attribution/{name}"] = path
    return paths