Analysis Toolkit

Analysis Toolkit

Author: Ciprian Ștefan Pleșca

The src/openlongevity/analysis/ package contains four dependency-light statistical utilities. None of them import the API, the database, or the provider layer — each is a pure function or a small stateful estimator operating only on the data structures its caller supplies. This page documents each one, together with the specific methodological limitation the project attaches to it.

flowchart LR
    subgraph analysis["analysis/ package"]
        BA["biological_age.py\nBiologicalAgeModel"]
        SU["survival.py\nkaplan_meier()"]
        PW["pathway.py\nenrich_gene_set()"]
        OM["omics.py\nintegrate_samples()"]
    end
    Caller[Calling code / notebook / future API route] --> BA
    Caller --> SU
    Caller --> PW
    Caller --> OM

Loading

Biological-age regression baseline

BiologicalAgeModel is an ordinary-least-squares regression implemented from first principles — it standardizes features, solves the normal equations via Gauss-Jordan elimination with partial pivoting, and reports mae, rmse, and r_squared.

flowchart TD
    Fit["fit(features, targets)"] --> Standardize["standardize each feature column:\n(x - mean) / std, std=1.0 if zero"]
    Standardize --> Design["prepend intercept column of 1s"]
    Design --> Normal["solve normal equations:\nXᵀX beta = Xᵀy"]
    Normal --> Pivot["Gauss-Jordan elimination\nwith partial pivoting"]
    Pivot --> Singular{"pivot magnitude < 1e-12?"}
    Singular -- yes --> Err["raise ValueError: feature matrix is singular"]
    Singular -- no --> Coeffs["store intercept + coefficients"]
    Coeffs --> Predict["predict(features): de-standardize\nand apply coefficients"]
    Predict --> Evaluate["evaluate(features, targets):\nmae, rmse, r_squared"]

Loading

The docstring is explicit that this is "a small, transparent regression baseline for synthetic biological-age data" and that "prediction is not clinical validation." Structurally, this class has no notion of biomarkers, aging biology, or clinical thresholds — it is a general-purpose linear regression tool that happens to be named for a biological-age use case; its scientific validity depends entirely on the features and targets a caller supplies, none of which are provided by the package itself.

Kaplan–Meier survival estimation

kaplan_meier(times, events) produces the classic step-function survival estimate without any external statistics dependency:

sequenceDiagram
    autonumber
    participant Caller
    participant KM as kaplan_meier()
    Caller->>KM: times[], events[] (events[i]=False means censored)
    KM->>KM: validate: same length, all times >= 0
    KM->>KM: survival = 1.0
    loop each distinct time, ascending
        KM->>KM: at_risk = count(time_j >= t)
        KM->>KM: deaths = count(time_j == t and event_j)
        alt at_risk > 0 and deaths > 0
            KM->>KM: survival *= (1 - deaths/at_risk)
        end
        KM->>KM: append KaplanMeierPoint(t, at_risk, deaths, round(survival, 6))
    end
    KM-->>Caller: list[KaplanMeierPoint]

Loading

The module docstring calls this a "dependency-light Kaplan–Meier estimator for demonstrations" — it returns point estimates only. It does not compute confidence intervals (e.g., via Greenwood's formula), does not perform log-rank tests between groups, and does not handle left-truncation or competing risks. Any of those would need to be layered on top by a caller.

Pathway enrichment

enrich_gene_set(genes, pathways, universe) computes a hypergeometric-test p-value for each pathway's overlap with a query gene set, then applies a Benjamini-Hochberg-style adjustment:

flowchart TD
    In["genes, pathways (name to gene-set mapping), universe"] --> Valid{"genes non-empty,\npathways non-empty,\nuniverse non-empty,\ngenes subset of universe?"}
    Valid -- no --> E[raise ValueError]
    Valid -- yes --> Loop[for each pathway]
    Loop --> Overlap["overlap = |genes ∩ members|\nsize = |members ∩ universe|"]
    Overlap --> Skip{"size == 0 or overlap == 0?"}
    Skip -- yes --> Next[skip this pathway]
    Skip -- no --> Tail["hypergeometric upper-tail sum\nusing math.comb"]
    Tail --> Denom["divide by comb(population, sample)\nclamp to 1.0"]
    Denom --> Collect[collect name, overlap, size, p_value]
    Next --> Loop
    Collect --> Sort["sort tests by ascending p_value"]
    Sort --> BH["adjusted_p = min(1.0, p * n_tests / rank)"]
    BH --> Out["list[EnrichmentResult]"]

Loading

Two implementation details matter for correct interpretation:

• The multiple-testing correction is computed only over pathways that passed the size > 0 and overlap > 0 filter, not over every pathway originally supplied. This means the effective number of tests (n_tests in the BH formula) can be smaller than len(pathways), which the README flags directly: "the pathway routine has an identified multiple-testing adjustment limitation."
• The hypergeometric tail sum guards against invalid comb() arguments with the condition sample_size - k <= population_size - size, preventing a ValueError from math.comb when k would require more "failures" than exist in the remaining population.

Multi-omics integration

integrate_samples(samples) groups MultiOmicsSample records by sample_id, nesting each sample's OmicsLayer value and feature dictionary underneath:

flowchart LR
    S1["MultiOmicsSample(sample_id=S1, layer=GENOMICS, features=...)"] --> R["result: sample S1 maps to layer genomics with its features"]
    S2["MultiOmicsSample(sample_id=S1, layer=PROTEOMICS, features=...)"] --> R2["result: sample S1 maps to layers genomics and proteomics"]

Loading

Two properties are worth calling out explicitly, since they define the boundary of what this function does and does not guarantee:

1. Missing values remain explicit as None rather than being silently dropped or imputed — the docstring states this directly. A caller consuming the integrated structure must handle None values themselves; the function does not average, interpolate, or otherwise fabricate a value.
2. No batch correction or participant-consistency check is performed. The function groups by sample_id string equality only. If two samples share a sample_id but originate from different batch_id values or actually belong to different participant_ids due to an upstream labeling error, integrate_samples will silently merge them — grouping is purely syntactic, and the batch_id field is carried through unused by the grouping logic itself, available only for a caller to inspect afterward.

Why these utilities are kept outside the API layer

None of the four analysis utilities are currently wired into any FastAPI route. This is consistent with the project's separation-of-concerns stance: exposing a regression fit, a survival curve, or a pathway enrichment result over a public API implies a level of curation, validation, and interpretive guidance (confidence intervals, assumption checks, sample provenance) that these first-principles implementations do not yet provide. Keeping them as standalone, importable utilities lets them be used in notebooks or scripts under direct human supervision without implying they are production-grade scientific endpoints.

Next

See Scientific-Methodology-and-Limitations for how each of these utilities' outputs should be framed when reported.