"""
CodonAdaptPy Comparative and Multivariate Analysis
This module performs group summaries, Welch tests, effect-size estimation,
Benjamini-Hochberg correction, correlation matrices, PCA, linear discriminant
classification, hierarchical clustering, and pairwise distances. Lightweight
summaries work with the standard library; multivariate operations require the
analysis extra.
Classes:
- ComparativeAnalyzer: Analyze collections of scalar sequence metrics.
:Created: July 20, 2026
:Updated: July 20, 2026
:Author: Naveen Duhan
:Version: 1.0.2
"""
from __future__ import annotations
import math
import random
import statistics
from collections import defaultdict
from collections.abc import Mapping, Sequence
from typing import Any
from .exceptions import OptionalDependencyError
[docs]
class ComparativeAnalyzer:
"""Provide group-level statistical and multivariate analysis helpers."""
[docs]
def summarize(self, values: Sequence[float]) -> dict[str, float]:
"""Return count, mean, median, standard deviation, minimum, and maximum."""
clean = [float(value) for value in values if math.isfinite(float(value))]
if not clean:
return {key: float("nan") for key in ("count", "mean", "median", "sd", "min", "max")}
return {
"count": float(len(clean)),
"mean": statistics.fmean(clean),
"median": statistics.median(clean),
"sd": statistics.stdev(clean) if len(clean) > 1 else 0.0,
"min": min(clean),
"max": max(clean),
}
[docs]
def group_summaries(
self, records: Sequence[Mapping[str, Any]], *, group_key: str, metric_keys: Sequence[str]
) -> dict[str, dict[str, dict[str, float]]]:
"""Summarize selected numeric metrics for each metadata group."""
groups: dict[str, list[Mapping[str, Any]]] = defaultdict(list)
for record in records:
groups[str(record[group_key])].append(record)
return {
group: {metric: self.summarize([row[metric] for row in rows]) for metric in metric_keys}
for group, rows in groups.items()
}
[docs]
def compare_two_groups(
self, left: Sequence[float], right: Sequence[float], *, test: str = "welch"
) -> dict[str, float]:
"""Run Welch or Mann-Whitney testing and calculate Hedges' g."""
if len(left) < 2 or len(right) < 2:
raise ValueError("At least two observations per group are required.")
try:
from scipy import stats
except ImportError as exc:
raise OptionalDependencyError("Statistical tests require CodonAdaptPy[analysis].") from exc
if test == "welch":
test_result = stats.ttest_ind(left, right, equal_var=False, nan_policy="omit")
elif test == "mannwhitney":
test_result = stats.mannwhitneyu(left, right, alternative="two-sided")
else:
raise ValueError("test must be 'welch' or 'mannwhitney'.")
n_left, n_right = len(left), len(right)
pooled = math.sqrt(
((n_left - 1) * statistics.variance(left) + (n_right - 1) * statistics.variance(right))
/ (n_left + n_right - 2)
)
d = (statistics.fmean(left) - statistics.fmean(right)) / pooled if pooled else 0.0
correction = 1 - 3 / (4 * (n_left + n_right) - 9)
return {
"statistic": float(test_result.statistic),
"p_value": float(test_result.pvalue),
"hedges_g": d * correction,
}
[docs]
def automatic_test(
self,
left: Sequence[float],
right: Sequence[float],
*,
paired: bool = False,
alpha: float = 0.05,
estimand: str = "mean",
confidence: float = 0.95,
bootstrap_replicates: int = 2000,
random_seed: int = 1,
) -> dict[str, Any]:
"""Run an estimand-first paired or independent two-sample analysis.
``estimand="mean"`` tests a mean difference using a paired t test or
Welch's unequal-variance t test. ``estimand="distribution"`` tests a
paired shift with Wilcoxon's signed-rank test or an independent
stochastic-ordering difference with the Mann--Whitney U test. This
choice is made from the scientific estimand and sampling design, not
from a preliminary normality-test decision.
Normality diagnostics remain available as model checks. The result
also includes a percentile-bootstrap confidence interval for the
corresponding location contrast and an effect size appropriate to the
selected design. Missing and non-finite observations are removed; in
paired mode, pairs are removed together to preserve correspondence.
"""
if not 0 < alpha < 1:
raise ValueError("alpha must be between zero and one.")
if estimand not in {"mean", "distribution"}:
raise ValueError("estimand must be 'mean' or 'distribution'.")
if not 0 < confidence < 1:
raise ValueError("confidence must be between zero and one.")
if bootstrap_replicates < 100:
raise ValueError("bootstrap_replicates must be at least 100.")
if paired and len(left) != len(right):
raise ValueError("Paired samples must contain the same number of observations.")
try:
from scipy import stats
except ImportError as exc:
raise OptionalDependencyError("Automatic statistical testing requires CodonAdaptPy[analysis].") from exc
if paired:
pairs = [
(float(first), float(second))
for first, second in zip(left, right, strict=True)
if math.isfinite(float(first)) and math.isfinite(float(second))
]
if len(pairs) < 3:
raise ValueError("At least three complete pairs are required for automatic test selection.")
clean_left = [first for first, _ in pairs]
clean_right = [second for _, second in pairs]
differences = [first - second for first, second in pairs]
normality = {"paired_differences": self._normality(differences, alpha=alpha)}
if estimand == "mean":
result = stats.ttest_rel(clean_left, clean_right, nan_policy="omit")
spread = statistics.stdev(differences)
effect = statistics.fmean(differences) / spread if spread else 0.0
test_name = "paired_t"
effect_name = "cohens_dz"
contrast_name = "mean_difference"
contrast = statistics.fmean(differences)
else:
result = stats.wilcoxon(clean_left, clean_right, alternative="two-sided")
nonzero = [value for value in differences if value != 0]
ranked = stats.rankdata([abs(value) for value in nonzero])
positive = sum(rank for rank, value in zip(ranked, nonzero, strict=True) if value > 0)
negative = sum(rank for rank, value in zip(ranked, nonzero, strict=True) if value < 0)
total = positive + negative
effect = (positive - negative) / total if total else 0.0
test_name = "wilcoxon_signed_rank"
effect_name = "matched_rank_biserial"
contrast_name = "median_paired_difference"
contrast = statistics.median(differences)
else:
clean_left = [float(value) for value in left if math.isfinite(float(value))]
clean_right = [float(value) for value in right if math.isfinite(float(value))]
if len(clean_left) < 3 or len(clean_right) < 3:
raise ValueError("At least three observations per group are required for automatic test selection.")
normality = {
"left": self._normality(clean_left, alpha=alpha),
"right": self._normality(clean_right, alpha=alpha),
}
if estimand == "mean":
result = stats.ttest_ind(clean_left, clean_right, equal_var=False, nan_policy="omit")
n_left, n_right = len(clean_left), len(clean_right)
pooled = math.sqrt(
((n_left - 1) * statistics.variance(clean_left) + (n_right - 1) * statistics.variance(clean_right))
/ (n_left + n_right - 2)
)
d_value = (statistics.fmean(clean_left) - statistics.fmean(clean_right)) / pooled if pooled else 0.0
effect = d_value * (1 - 3 / (4 * (n_left + n_right) - 9))
test_name = "welch_t"
effect_name = "hedges_g"
contrast_name = "mean_difference"
contrast = statistics.fmean(clean_left) - statistics.fmean(clean_right)
else:
result = stats.mannwhitneyu(clean_left, clean_right, alternative="two-sided")
effect = 2 * float(result.statistic) / (len(clean_left) * len(clean_right)) - 1
test_name = "mann_whitney_u"
effect_name = "rank_biserial"
contrast_name = "median_difference"
contrast = statistics.median(clean_left) - statistics.median(clean_right)
confidence_interval = self._bootstrap_location_interval(
clean_left,
clean_right,
paired=paired,
location="mean" if estimand == "mean" else "median",
confidence=confidence,
replicates=bootstrap_replicates,
random_seed=random_seed,
)
return {
"test": test_name,
"paired": paired,
"alpha": alpha,
"normality": normality,
"estimand": estimand,
"selection_policy": "test selected from the prespecified estimand and paired or independent design",
"n_left": len(clean_left),
"n_right": len(clean_right),
"statistic": float(result.statistic),
"p_value": float(result.pvalue),
"effect_size": float(effect),
"effect_size_name": effect_name,
"contrast": float(contrast),
"contrast_name": contrast_name,
"confidence_level": confidence,
"confidence_interval": confidence_interval,
"bootstrap_replicates": bootstrap_replicates,
"random_seed": random_seed,
}
@staticmethod
def _bootstrap_location_interval(
left: Sequence[float],
right: Sequence[float],
*,
paired: bool,
location: str,
confidence: float,
replicates: int,
random_seed: int,
) -> list[float]:
"""Return a deterministic percentile-bootstrap interval for a contrast."""
statistic = statistics.fmean if location == "mean" else statistics.median
generator = random.Random(random_seed)
estimates: list[float] = []
if paired:
differences = [first - second for first, second in zip(left, right, strict=True)]
for _ in range(replicates):
sample = generator.choices(differences, k=len(differences))
estimates.append(float(statistic(sample)))
else:
for _ in range(replicates):
first = generator.choices(list(left), k=len(left))
second = generator.choices(list(right), k=len(right))
estimates.append(float(statistic(first) - statistic(second)))
estimates.sort()
tail = (1 - confidence) / 2
lower = estimates[max(0, math.floor(tail * replicates))]
upper = estimates[min(replicates - 1, math.ceil((1 - tail) * replicates) - 1)]
return [lower, upper]
@staticmethod
def _normality(values: Sequence[float], *, alpha: float) -> dict[str, float | str | bool]:
"""Return a named normality test and decision for finite values."""
try:
from scipy import stats
except ImportError as exc:
raise OptionalDependencyError("Normality testing requires CodonAdaptPy[analysis].") from exc
if len(values) < 3:
raise ValueError("Normality testing requires at least three observations.")
if len(values) <= 5000:
result = stats.shapiro(values)
method = "shapiro_wilk"
else:
result = stats.normaltest(values)
method = "dagostino_pearson"
return {
"test": method,
"statistic": float(result.statistic),
"p_value": float(result.pvalue),
"normal": bool(result.pvalue >= alpha),
}
[docs]
@staticmethod
def adjust_bh(p_values: Sequence[float]) -> list[float]:
"""Apply Benjamini-Hochberg false-discovery-rate correction."""
size = len(p_values)
ordered = sorted(enumerate(p_values), key=lambda item: item[1])
adjusted = [1.0] * size
previous = 1.0
for rank, (index, value) in reversed(list(enumerate(ordered, start=1))):
previous = min(previous, float(value) * size / rank)
adjusted[index] = min(1.0, previous)
return adjusted
[docs]
def multivariate(self, matrix: Sequence[Sequence[float]], *, components: int = 2) -> dict[str, Any]:
"""Return standardized PCA scores, loadings, variance, and distances."""
try:
import numpy as np
from sklearn.decomposition import PCA
from sklearn.metrics import pairwise_distances
from sklearn.preprocessing import StandardScaler
except ImportError as exc:
raise OptionalDependencyError("PCA and distances require CodonAdaptPy[analysis].") from exc
array = np.asarray(matrix, dtype=float)
scaled = StandardScaler().fit_transform(array)
model = PCA(n_components=min(components, *scaled.shape)).fit(scaled)
with np.errstate(divide="ignore", invalid="ignore"):
corr_matrix = np.nan_to_num(np.corrcoef(array, rowvar=False), nan=0.0).tolist()
return {
"scores": model.transform(scaled).tolist(),
"loadings": model.components_.tolist(),
"explained_variance_ratio": model.explained_variance_ratio_.tolist(),
"euclidean_distances": pairwise_distances(scaled).tolist(),
"correlation_matrix": corr_matrix,
}
[docs]
def hierarchical_clustering(
self,
matrix: Sequence[Sequence[float]],
*,
method: str = "average",
metric: str = "euclidean",
) -> dict[str, Any]:
"""Return SciPy linkage and leaf order for reproducible clustering."""
try:
import numpy as np
from scipy.cluster.hierarchy import leaves_list, linkage
except ImportError as exc:
raise OptionalDependencyError("Clustering requires CodonAdaptPy[analysis].") from exc
linked = linkage(np.asarray(matrix, dtype=float), method=method, metric=metric)
return {"linkage": linked.tolist(), "leaf_order": leaves_list(linked).tolist()}
[docs]
def linear_discriminant(
self,
matrix: Sequence[Sequence[float]],
labels: Sequence[str],
*,
test_size: float = 0.3,
random_seed: int = 1,
) -> dict[str, Any]:
"""Train and evaluate stratified linear discriminant classification.
The returned projection contains every input observation, while the
confusion matrix and classification report are calculated only from
the held-out test subset. The split indices and random seed are
retained to make reported accuracy exactly reproducible.
"""
if len(matrix) != len(labels) or len(matrix) < 3:
raise ValueError("matrix and labels must contain the same three or more observations.")
if not 0 < test_size < 1:
raise ValueError("test_size must be between zero and one.")
try:
import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.metrics import accuracy_score, classification_report, confusion_matrix
from sklearn.model_selection import train_test_split
except ImportError as exc:
raise OptionalDependencyError("LDA requires CodonAdaptPy[analysis].") from exc
values = np.asarray(matrix, dtype=float)
categories = np.asarray([str(label) for label in labels])
classes, counts = np.unique(categories, return_counts=True)
if len(classes) < 2 or np.any(counts < 2):
raise ValueError("LDA requires at least two classes with two observations each.")
indices = np.arange(len(categories))
train_indices, test_indices = train_test_split(
indices,
test_size=test_size,
random_state=random_seed,
stratify=categories,
)
model = LinearDiscriminantAnalysis().fit(values[train_indices], categories[train_indices])
predicted = model.predict(values[test_indices])
dimensions = min(2, len(classes) - 1)
projection = model.transform(values)[:, :dimensions]
return {
"classes": classes.tolist(),
"accuracy": float(accuracy_score(categories[test_indices], predicted)),
"confusion_matrix": confusion_matrix(categories[test_indices], predicted, labels=classes).tolist(),
"classification_report": classification_report(
categories[test_indices],
predicted,
labels=classes,
output_dict=True,
zero_division=0,
),
"projection": projection.tolist(),
"predictions": predicted.tolist(),
"train_indices": train_indices.tolist(),
"test_indices": test_indices.tolist(),
"test_labels": categories[test_indices].tolist(),
"test_size": test_size,
"random_seed": random_seed,
}
[docs]
def repeated_grouped_discriminant(
self,
matrix: Sequence[Sequence[float]],
labels: Sequence[str],
groups: Sequence[str | int],
*,
folds: int = 5,
repeats: int = 20,
random_seed: int = 1,
shrinkage: bool = True,
) -> dict[str, Any]:
"""Evaluate LDA with repeated stratified, group-exclusive folds.
All observations sharing a group are kept in the same fold, preventing
exact duplicates, isolates, or sequence clusters from appearing in
both training and test data. Metrics are reported per fold and as
mean, standard deviation, and percentile intervals across all held-out
predictions. Shrinkage LDA is used by default for high-dimensional
codon-feature matrices.
"""
if len(matrix) != len(labels) or len(labels) != len(groups) or len(matrix) < 4:
raise ValueError("matrix, labels, and groups must contain the same four or more observations.")
if folds < 2 or repeats < 1:
raise ValueError("folds must be at least two and repeats at least one.")
try:
import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.metrics import accuracy_score, balanced_accuracy_score, confusion_matrix, f1_score
from sklearn.model_selection import StratifiedGroupKFold
except ImportError as exc:
raise OptionalDependencyError("Repeated grouped LDA requires CodonAdaptPy[analysis].") from exc
values = np.asarray(matrix, dtype=float)
categories = np.asarray([str(label) for label in labels])
group_values = np.asarray([str(group) for group in groups])
classes = np.unique(categories)
if len(classes) < 2:
raise ValueError("Grouped LDA requires at least two classes.")
groups_per_class = [len(np.unique(group_values[categories == category])) for category in classes]
if min(groups_per_class) < folds:
raise ValueError("Each class must contain at least as many distinct groups as folds.")
fold_results: list[dict[str, Any]] = []
confusion = np.zeros((len(classes), len(classes)), dtype=int)
for repeat in range(repeats):
seed = random_seed + repeat
splitter = StratifiedGroupKFold(n_splits=folds, shuffle=True, random_state=seed)
for fold, (train_indices, test_indices) in enumerate(
splitter.split(values, categories, groups=group_values), start=1
):
train_groups = set(group_values[train_indices])
test_groups = set(group_values[test_indices])
if train_groups & test_groups:
raise RuntimeError("A group crossed the training/test boundary.")
parameters = {"solver": "lsqr", "shrinkage": "auto"} if shrinkage else {"solver": "svd"}
model = LinearDiscriminantAnalysis(**parameters).fit(values[train_indices], categories[train_indices])
predicted = model.predict(values[test_indices])
observed = categories[test_indices]
fold_confusion = confusion_matrix(observed, predicted, labels=classes)
confusion += fold_confusion
fold_results.append(
{
"repeat": repeat + 1,
"fold": fold,
"random_seed": seed,
"accuracy": float(accuracy_score(observed, predicted)),
"balanced_accuracy": float(balanced_accuracy_score(observed, predicted)),
"macro_f1": float(f1_score(observed, predicted, average="macro", zero_division=0)),
"train_indices": train_indices.tolist(),
"test_indices": test_indices.tolist(),
"train_groups": sorted(train_groups),
"test_groups": sorted(test_groups),
"observed": observed.tolist(),
"predicted": predicted.tolist(),
}
)
summary: dict[str, dict[str, float | list[float]]] = {}
for metric in ("accuracy", "balanced_accuracy", "macro_f1"):
metric_values = np.asarray([row[metric] for row in fold_results], dtype=float)
summary[metric] = {
"mean": float(metric_values.mean()),
"sd": float(metric_values.std(ddof=1)) if len(metric_values) > 1 else 0.0,
"interval_95": np.percentile(metric_values, [2.5, 97.5]).tolist(),
"min": float(metric_values.min()),
"max": float(metric_values.max()),
}
return {
"classes": classes.tolist(),
"folds": folds,
"repeats": repeats,
"random_seed": random_seed,
"shrinkage": shrinkage,
"group_exclusive": True,
"summary": summary,
"aggregate_confusion_matrix": confusion.tolist(),
"fold_results": fold_results,
}
[docs]
def correspondence_analysis(self, matrix: Sequence[Sequence[float]], *, components: int = 2) -> dict[str, Any]:
"""Perform correspondence analysis on a non-negative count/RSCU matrix."""
try:
import numpy as np
except ImportError as exc:
raise OptionalDependencyError("Correspondence analysis requires CodonAdaptPy[analysis].") from exc
values = np.asarray(matrix, dtype=float)
if values.ndim != 2 or values.size == 0 or np.any(values < 0) or values.sum() <= 0:
raise ValueError("Correspondence analysis requires a non-empty non-negative matrix.")
probabilities = values / values.sum()
row_masses = probabilities.sum(axis=1)
column_masses = probabilities.sum(axis=0)
expected = np.outer(row_masses, column_masses)
residuals = np.divide(
probabilities - expected,
np.sqrt(expected),
out=np.zeros_like(probabilities),
where=expected > 0,
)
left, singular, right = np.linalg.svd(residuals, full_matrices=False)
count = min(components, len(singular))
row_coordinates = np.divide(
left[:, :count] * singular[:count],
np.sqrt(row_masses)[:, None],
out=np.zeros((values.shape[0], count)),
where=row_masses[:, None] > 0,
)
column_coordinates = np.divide(
right[:count, :].T * singular[:count],
np.sqrt(column_masses)[:, None],
out=np.zeros((values.shape[1], count)),
where=column_masses[:, None] > 0,
)
inertia = singular**2
return {
"row_coordinates": row_coordinates.tolist(),
"column_coordinates": column_coordinates.tolist(),
"explained_inertia": (inertia[:count] / inertia.sum()).tolist() if inertia.sum() else [0.0] * count,
}