Statistics

The statistics module performs cluster-based permutation testing on MNI-space NIfTI volumes or fsaverage surface caches produced by the simulation pipeline. It supports two analysis types: group comparison (responders vs non-responders) and correlation (voxelwise correlation with a continuous outcome measure, ACES-style).

graph LR
    NIFTI([MNI NIfTI Volumes]) --> TEST[Voxelwise Test]
    CSV([Subject CSV]) --> TEST
    TEST --> PERM[Permutation Engine]
    PERM --> REPORT([Summary Report])
    PERM --> PLOTS([Null Distribution Plots])
    PERM --> MAPS([NIfTI Output Maps])
    style NIFTI fill:#1a3a5c,stroke:#48a,color:#fff
    style CSV fill:#1a3a5c,stroke:#48a,color:#fff
    style TEST fill:#2d5a27,stroke:#4a8,color:#fff
    style PERM fill:#2d5a27,stroke:#4a8,color:#fff
    style REPORT fill:#1a5c4a,stroke:#4a8,color:#fff
    style PLOTS fill:#1a5c4a,stroke:#4a8,color:#fff
    style MAPS fill:#1a5c4a,stroke:#4a8,color:#fff

Group Comparison

Compare two groups (e.g., responders vs non-responders) using cluster-based permutation testing with voxelwise t-tests. Supports both unpaired (independent samples) and paired designs.

from tit.stats import GroupComparisonConfig, run_group_comparison

# Load subjects from CSV (columns: subject_id, simulation_name, response)
subjects = GroupComparisonConfig.load_subjects("/data/my_project/subjects.csv")

config = GroupComparisonConfig(
    analysis_name="responder_comparison",
    subjects=subjects,
    test_type=GroupComparisonConfig.TestType.UNPAIRED,
    alternative=GroupComparisonConfig.Alternative.TWO_SIDED,
    n_permutations=5000,
    alpha=0.05,
    cluster_threshold=0.05,
    cluster_stat=GroupComparisonConfig.ClusterStat.MASS,
    tissue_type=GroupComparisonConfig.TissueType.GREY,
    group1_name="Responders",
    group2_name="Non-Responders",
)

result = run_group_comparison(config)
print(f"Significant clusters: {result.n_significant_clusters}")
print(f"Significant voxels:   {result.n_significant_voxels}")
print(f"Analysis time:        {result.analysis_time:.1f}s")
print(f"Output directory:     {result.output_dir}")

Correlation

Test for voxelwise correlation between electric field magnitude and a continuous outcome measure (e.g., clinical effect size). Supports Pearson and Spearman correlation with optional subject-level weights.

from tit.stats import CorrelationConfig, run_correlation

# Load subjects from CSV (columns: subject_id, simulation_name, effect_size; optional: weight)
subjects = CorrelationConfig.load_subjects("/data/my_project/correlation_subjects.csv")

config = CorrelationConfig(
    analysis_name="efield_outcome_correlation",
    subjects=subjects,
    correlation_type=CorrelationConfig.CorrelationType.PEARSON,
    n_permutations=5000,
    alpha=0.05,
    cluster_threshold=0.05,
    cluster_stat=CorrelationConfig.ClusterStat.MASS,
    tissue_type=CorrelationConfig.TissueType.GREY,
    use_weights=True,
    effect_metric="Clinical Improvement",
)

result = run_correlation(config)
print(f"Significant clusters: {result.n_significant_clusters}")
print(f"Significant voxels:   {result.n_significant_voxels}")

Configuration

Subject Definitions

Both analysis types use a nested Subject dataclass to define per-subject metadata. Subjects are typically loaded from CSV files via the load_subjects class method.

Group ComparisonCorrelation

CSV columns: subject_id, simulation_name, response (0 or 1).

subjects = GroupComparisonConfig.load_subjects("/data/my_project/subjects.csv")

Field	Type	Description
subject_id	str	Subject identifier (e.g., "001")
simulation_name	str	Name of the simulation to load
response	int	Group assignment: 1 = group 1, 0 = group 2

CSV columns: subject_id, simulation_name, effect_size; optional: weight.

subjects = CorrelationConfig.load_subjects("/data/my_project/subjects.csv")

Field	Type	Default	Description
subject_id	str		Subject identifier
simulation_name	str		Name of the simulation to load
effect_size	float		Continuous outcome measure
weight	float	1.0	Subject-level weight

Enums

Enum	Values	Used By
GroupComparisonConfig.TestType	UNPAIRED, PAIRED	Group comparison only
GroupComparisonConfig.Alternative	TWO_SIDED, GREATER, LESS	Group comparison only
CorrelationConfig.CorrelationType	PEARSON, SPEARMAN	Correlation only
ClusterStat	MASS, SIZE	Both (nested as GroupComparisonConfig.ClusterStat / CorrelationConfig.ClusterStat)
TissueType	GREY, WHITE, ALL	Both (nested as GroupComparisonConfig.TissueType / CorrelationConfig.TissueType)

Statistical Parameters

These parameters are shared by both GroupComparisonConfig and CorrelationConfig:

Parameter	Type	Default	Description
analysis_name	str		Name for this analysis run
subjects	list[Subject]		List of subject definitions
cluster_threshold	float	0.05	Uncorrected p-value threshold for cluster formation
cluster_stat	ClusterStat	MASS	Cluster statistic: "mass" (sum of t-values) or "size" (voxel count)
n_permutations	int	1000	Number of permutations for null distribution
alpha	float	0.05	Family-wise error rate for cluster significance
n_jobs	int	-1	Number of parallel jobs (-1 = all CPUs)
tissue_type	TissueType	GREY	Tissue mask: "grey", "white", or "all"
nifti_file_pattern	str \| None	None	Custom NIfTI filename pattern (auto-resolved from tissue_type if None)
atlas_files	list[str]	[]	Atlas filenames for overlap analysis

NIfTI File Patterns

The tissue_type field determines which NIfTI files are loaded from each subject's simulation output:

TissueType	Resolved Pattern
GREY	grey_{simulation_name}_TI_MNI_MNI_TI_max.nii.gz
WHITE	white_{simulation_name}_TI_MNI_MNI_TI_max.nii.gz
ALL	{simulation_name}_TI_MNI_MNI_TI_max.nii.gz

Set nifti_file_pattern to override the auto-resolved pattern.

Surface Analysis

Both configuration classes accept space="fsaverage" to analyze cortical surface caches instead of MNI volumes. fsaverage_field selects the cached field (default "TI_max") and fsaverage_spacing selects ico spacing 5, 6, or 7 (default 5). The corresponding caches must exist for every input simulation. NIfTI patterns and tissue-mask options apply to the default space="mni" workflow only.

Output

Results are saved to derivatives/ti-toolbox/stats/<analysis_type>/<analysis_name>/ within the project directory. For MNI-volume analysis, both analysis types produce the following:

File	Description
significant_voxels_mask.nii.gz	Binary mask of significant voxels
pvalues_map.nii.gz	Negative log10 p-value map
permutation_null_distribution.pdf	Null distribution with observed clusters
cluster_size_mass_correlation.pdf	Cluster size vs mass scatter plot
analysis_summary.txt	Text summary of results
permutation_details.txt	Per-permutation log
*_analysis_*.log	Timestamped run log

Group comparison additionally produces:

File	Description
average_responders.nii.gz	Mean field map for group 1
average_non_responders.nii.gz	Mean field map for group 2
difference_map.nii.gz	Group 1 minus group 2 difference map

Correlation additionally produces:

File	Description
correlation_map.nii.gz	Full voxelwise correlation map
correlation_map_thresholded.nii.gz	Correlation map masked to significant voxels
t_statistics_map.nii.gz	Voxelwise t-statistic map
average_efield.nii.gz	Mean electric field across all subjects

Result Dataclasses

GroupComparisonResult

Field	Type	Description
success	bool	Whether the analysis completed
output_dir	str	Path to the output directory
n_responders	int	Number of group 1 subjects
n_non_responders	int	Number of group 2 subjects
n_significant_voxels	int	Count of significant voxels
n_significant_clusters	int	Count of significant clusters
cluster_threshold	float	Cluster statistic threshold from null distribution
analysis_time	float	Total runtime in seconds
clusters	list	Cluster details (size, MNI center)
log_file	str	Path to the analysis log

CorrelationResult

Field	Type	Description
success	bool	Whether the analysis completed
output_dir	str	Path to the output directory
n_subjects	int	Number of subjects
n_significant_voxels	int	Count of significant voxels
n_significant_clusters	int	Count of significant clusters
cluster_threshold	float	Cluster statistic threshold from null distribution
analysis_time	float	Total runtime in seconds
clusters	list	Cluster details (size, MNI center, mean/peak r)
log_file	str	Path to the analysis log

API Reference

tit.stats.config.GroupComparisonConfig dataclass

GroupComparisonConfig(analysis_name: str, subjects: list[Subject], test_type: TestType = UNPAIRED, alternative: Alternative = TWO_SIDED, cluster_threshold: float = 0.05, cluster_stat: ClusterStat = MASS, n_permutations: int = 1000, alpha: float = 0.05, n_jobs: int = -1, seed: int | None = None, tissue_type: TissueType = GREY, nifti_file_pattern: str | None = None, space: AnalysisSpace = MNI, fsaverage_field: str = 'TI_max', fsaverage_spacing: int = 5, group1_name: str = 'Responders', group2_name: str = 'Non-Responders', value_metric: str = 'Electric Field Magnitude', atlas_files: list[str] = list())

Configuration for cluster-based permutation testing between two groups.

Compares voxelwise field intensities between responders and non-responders using a t-test with cluster-based permutation correction for multiple comparisons.

Attributes

analysis_name : str Human-readable name for this analysis run. subjects : list of Subject Subject entries, each labelled as responder (1) or non-responder (0). test_type : TestType or str "unpaired" (default) or "paired" t-test. Strings are coerced to :class:TestType. alternative : Alternative or str Sidedness: "two-sided" (default), "greater" or "less". cluster_threshold : float Uncorrected p-value threshold for forming clusters. Default 0.05. cluster_stat : ClusterStat or str Cluster-level statistic used for permutation testing, "mass" (default) or "size". n_permutations : int Number of permutations for the null distribution. Default 1000; the smallest reportable p-value is 1 / (n + 1). alpha : float Family-wise error rate for significance. Default 0.05. n_jobs : int Number of parallel workers (-1, the default, for the global CPU limit -- :func:tit.cpu.cpu_limit; larger values are clamped to it). seed : int or None Seed of the permutation generator; None (default) draws a fresh one. Either way the seed used is written in the run's summary (analysis_summary.txt; surface_stats.json for fsaverage). tissue_type : TissueType or str Which tissue compartment to analyze: "grey" (default), "white" or "all". Ignored when space is "fsaverage". nifti_file_pattern : str or None Filename pattern for subject NIfTI files, with a {simulation_name} placeholder. If None (default), derived from tissue_type, e.g. "grey_{simulation_name}_TI_MNI_MNI_TI_max.nii.gz". space : AnalysisSpace or str Where the statistics run: "mni" (default, voxelwise on the MNI-space NIfTIs) or "fsaverage" (vertexwise on the per-subject fsaverage projections written by the simulator). fsaverage_field : str Surface field for space="fsaverage"; one of :data:tit.constants.FSAVG_FIELD_NAMES ("TI_max", "TI_normal", "hf_peak", "hf_rms"). Default "TI_max". fsaverage_spacing : int fsaverage ico spacing for space="fsaverage": 5 (default), 6 or 7. group1_name : str Display label for the responder group. Default "Responders". group2_name : str Display label for the non-responder group. Default "Non-Responders". value_metric : str Label for the field values in the summary. Default "Electric Field Magnitude". atlas_files : list of str Atlas filenames for overlap analysis (looked up in the bundled atlas directory). Default empty.

Raises

ValueError If subjects lacks at least one responder and one non-responder, if a response is not 0 or 1, if a paired test has unequal group sizes, if alpha or cluster_threshold is not strictly between 0 and 1, if n_permutations is not a whole number or is below floor(1 / alpha) (20 at alpha 0.05: with fewer, no p-value could fall below alpha), if a string value is not a member of its enum, or if fsaverage_field/fsaverage_spacing are invalid for space="fsaverage".

Examples

from tit.stats import GroupComparisonConfig subjects = [ ... GroupComparisonConfig.Subject("ernie", "L_Insula", response=1), ... GroupComparisonConfig.Subject("101", "L_Insula", response=0), ... ] cfg = GroupComparisonConfig( ... analysis_name="active_vs_sham", subjects=subjects, ... test_type="unpaired", alternative="two-sided", ... cluster_stat="mass", n_permutations=1000, tissue_type="grey", ... ) cfg.test_type is GroupComparisonConfig.TestType.UNPAIRED True cfg.nifti_file_pattern 'grey_{simulation_name}_TI_MNI_MNI_TI_max.nii.gz'

Subjects can also come from a CSV with columns subject_id, simulation_name, response:

subjects = GroupComparisonConfig.load_subjects("subjects.csv") # doctest: +SKIP

See Also

GroupComparisonResult : Result container returned by the analysis. run_group_comparison : Orchestration function that consumes this config.

TestType

Bases: StrEnum

Type of statistical test for group comparison.

Attributes

UNPAIRED : str "unpaired" -- independent-samples t-test. PAIRED : str "paired" -- paired-samples t-test (groups must be the same size and ordered pairwise).

Alternative

Bases: StrEnum

Sidedness of the test hypothesis.

Attributes

TWO_SIDED : str "two-sided". GREATER : str "greater" -- responders > non-responders. LESS : str "less" -- responders < non-responders.

Subject dataclass

Subject(subject_id: str, simulation_name: str, response: int)

A single subject in a group comparison analysis.

Attributes

subject_id : str Subject identifier (without sub- prefix). simulation_name : str Name of the simulation to load for this subject. response : int Group label -- 1 for responder, 0 for non-responder.

load_subjects classmethod

load_subjects(csv_path: str) -> list[Subject]

Load group comparison subjects from a CSV file.

Expected columns: subject_id, simulation_name, response (0 or 1). The sub- prefix is stripped from subject IDs automatically.

Parameters

csv_path : str Path to a CSV file with the required columns.

Returns

list of Subject Subject instances parsed from the CSV rows.

Raises

ValueError If required columns are missing, or a row lacks one of their values.

Source code in tit/stats/config.py

@classmethod
def load_subjects(cls, csv_path: str) -> list["GroupComparisonConfig.Subject"]:
    """Load group comparison subjects from a CSV file.

    Expected columns: ``subject_id``, ``simulation_name``, ``response``
    (0 or 1).  The ``sub-`` prefix is stripped from subject IDs
    automatically.

    Parameters
    ----------
    csv_path : str
        Path to a CSV file with the required columns.

    Returns
    -------
    list of Subject
        Subject instances parsed from the CSV rows.

    Raises
    ------
    ValueError
        If required columns are missing, or a row lacks one of their values.
    """
    subjects = []
    for line, row in _read_participant_rows(
        csv_path, {"subject_id", "simulation_name", "response"}
    ):
        if None in (row["subject_id"], row["simulation_name"], row["response"]):
            raise ValueError(
                f"CSV row {line}: subject_id, simulation_name and response are "
                "all required"
            )
        subjects.append(
            cls.Subject(
                subject_id=_participant_id(row["subject_id"]),
                simulation_name=row["simulation_name"],
                response=int(float(row["response"])),
            )
        )
    return subjects

tit.stats.config.GroupComparisonResult dataclass

GroupComparisonResult(success: bool, output_dir: str, n_responders: int, n_non_responders: int, n_significant_voxels: int, n_significant_clusters: int, cluster_threshold: float, analysis_time: float, clusters: list, log_file: str)

Result of a group comparison permutation test.

Attributes

success : bool Whether the analysis completed without error. output_dir : str Absolute path to the directory containing all outputs (NIfTI maps, plots, summary text, log). n_responders : int Number of responder subjects included. n_non_responders : int Number of non-responder subjects included. n_significant_voxels : int Total voxels surviving cluster-corrected threshold. n_significant_clusters : int Number of spatially contiguous clusters that survived permutation correction. cluster_threshold : float Cluster-level statistic threshold derived from the permutation null distribution at the requested alpha. analysis_time : float Wall-clock duration of the full analysis in seconds. clusters : list of dict One entry per significant cluster, containing size, mass, peak coordinates, and atlas overlap info. log_file : str Absolute path to the analysis log file.

See Also

GroupComparisonConfig : Configuration that produced this result. run_group_comparison : Function that returns this result.

tit.stats.config.CorrelationConfig dataclass

CorrelationConfig(analysis_name: str, subjects: list[Subject], correlation_type: CorrelationType = PEARSON, cluster_threshold: float = 0.05, cluster_stat: ClusterStat = MASS, n_permutations: int = 1000, alpha: float = 0.05, n_jobs: int = -1, seed: int | None = None, use_weights: bool = True, tissue_type: TissueType = GREY, nifti_file_pattern: str | None = None, space: AnalysisSpace = MNI, fsaverage_field: str = 'TI_max', fsaverage_spacing: int = 5, effect_metric: str = 'Effect Size', field_metric: str = 'Electric Field Magnitude', atlas_files: list[str] = list())

Configuration for correlation-based cluster permutation testing.

Tests voxelwise correlation between brain field intensities and a continuous behavioral or clinical measure (effect size) across subjects, with cluster-based permutation correction for multiple comparisons.

Attributes

analysis_name : str Human-readable name for this analysis run. subjects : list of Subject Subject entries with associated effect sizes. correlation_type : CorrelationType Pearson or Spearman rank correlation. cluster_threshold : float Uncorrected p-value threshold for forming clusters. cluster_stat : ClusterStat Cluster-level statistic used for permutation testing ("mass" or "size"). n_permutations : int Number of permutations for the null distribution. alpha : float Family-wise error rate for significance. n_jobs : int Number of parallel workers (-1 for the global CPU limit -- :func:tit.cpu.cpu_limit; larger values are clamped to it). seed : int or None Seed of the permutation generator; see :class:GroupComparisonConfig. use_weights : bool Whether to apply per-subject weights during correlation. Default True. tissue_type : TissueType or str Which tissue compartment to analyze: "grey" (default), "white" or "all". nifti_file_pattern : str or None Filename pattern for subject NIfTI files. If None, derived automatically from tissue_type. space : AnalysisSpace or str "mni" (default) or "fsaverage"; see :class:GroupComparisonConfig. fsaverage_field : str Surface field for space="fsaverage". Default "TI_max". fsaverage_spacing : int fsaverage ico spacing (5, 6 or 7). Default 5. effect_metric : str Label for the behavioral/clinical variable in plots. field_metric : str Label for the field intensity axis in plots. atlas_files : list of str Atlas filenames for overlap analysis (looked up in the bundled atlas directory).

Raises

ValueError If fewer than 3 subjects are given, if an effect size is not a finite number or every effect size is equal (the correlation is undefined), if use_weights is on and a weight is not a positive finite number, or for the alpha/cluster_threshold/ n_permutations rules of :class:GroupComparisonConfig.

See Also

CorrelationResult : Result container returned by the analysis. run_correlation : Orchestration function that consumes this config.

CorrelationType

Bases: StrEnum

Type of correlation coefficient to compute.

Subject dataclass

Subject(subject_id: str, simulation_name: str, effect_size: float, weight: float = 1.0)

A single subject in a correlation analysis.

Attributes

subject_id : str Subject identifier (without sub- prefix). simulation_name : str Name of the simulation to load for this subject. effect_size : float Continuous behavioral or clinical measure to correlate with field intensity. weight : float Per-subject weight (default 1.0).

load_subjects classmethod

load_subjects(csv_path: str) -> list[Subject]

Load correlation subjects from a CSV file.

Expected columns: subject_id, simulation_name, effect_size. Optional column: weight. Rows with NaN subject_id or effect_size are silently skipped. The sub- prefix is stripped from subject IDs automatically.

Parameters

csv_path : str Path to a CSV file with the required columns.

Returns

list of Subject Subject instances parsed from valid CSV rows.

Raises

ValueError If required columns are missing, a kept row has no simulation_name, or no valid subjects are found.

Source code in tit/stats/config.py

@classmethod
def load_subjects(cls, csv_path: str) -> list["CorrelationConfig.Subject"]:
    """Load correlation subjects from a CSV file.

    Expected columns: ``subject_id``, ``simulation_name``,
    ``effect_size``.  Optional column: ``weight``.  Rows with NaN
    ``subject_id`` or ``effect_size`` are silently skipped.  The ``sub-``
    prefix is stripped from subject IDs automatically.

    Parameters
    ----------
    csv_path : str
        Path to a CSV file with the required columns.

    Returns
    -------
    list of Subject
        Subject instances parsed from valid CSV rows.

    Raises
    ------
    ValueError
        If required columns are missing, a kept row has no simulation_name, or
        no valid subjects are found.
    """
    subjects = []
    for line, row in _read_participant_rows(
        csv_path, {"subject_id", "simulation_name", "effect_size"}
    ):
        if row["subject_id"] is None or row["effect_size"] is None:
            continue
        if row["simulation_name"] is None:
            raise ValueError(f"CSV row {line}: simulation_name is required")
        weight = row.get("weight")
        subjects.append(
            cls.Subject(
                subject_id=_participant_id(row["subject_id"]),
                simulation_name=row["simulation_name"],
                effect_size=float(row["effect_size"]),
                weight=1.0 if weight is None else float(weight),
            )
        )

    if not subjects:
        raise ValueError("No valid subjects found in CSV")
    return subjects

tit.stats.config.CorrelationResult dataclass

CorrelationResult(success: bool, output_dir: str, n_subjects: int, n_significant_voxels: int, n_significant_clusters: int, cluster_threshold: float, analysis_time: float, clusters: list, log_file: str)

Result of a correlation-based cluster permutation test.

Attributes

success : bool Whether the analysis completed without error. output_dir : str Absolute path to the directory containing all outputs (NIfTI maps, plots, summary text, log). n_subjects : int Number of subjects included in the analysis. n_significant_voxels : int Total voxels surviving cluster-corrected threshold. n_significant_clusters : int Number of spatially contiguous clusters that survived permutation correction. cluster_threshold : float Cluster-level statistic threshold derived from the permutation null distribution at the requested alpha. analysis_time : float Wall-clock duration of the full analysis in seconds. clusters : list of dict One entry per significant cluster, containing size, mass, peak coordinates, mean/peak correlation coefficients, and atlas overlap info. log_file : str Absolute path to the analysis log file.

See Also

CorrelationConfig : Configuration that produced this result. run_correlation : Function that returns this result.

tit.stats.permutation.run_group_comparison

run_group_comparison(config: GroupComparisonConfig, callback_handler: Handler | None = None, stop_callback: Callable[[], bool] | None = None) -> GroupComparisonResult

Run cluster-based permutation testing for a two-group comparison.

Loads responder and non-responder field maps, performs voxelwise (or, for config.space == "fsaverage", vertexwise) t-tests, applies cluster-based permutation correction, generates diagnostic plots, and saves all outputs under <project>/derivatives/ti-toolbox/stats/group_comparison/<analysis_name>/.

Parameters

config : GroupComparisonConfig Fully specified group comparison configuration. callback_handler : logging.Handler or None, optional Extra handler attached to the run-scoped logger (used by the app to stream log lines). stop_callback : callable or None, optional Zero-argument callable returning True to request early termination. Checked between pipeline stages.

Returns

GroupComparisonResult Summary of the outcome: output directory, subject counts, significant voxel/cluster counts, the permutation-derived cluster threshold (NaN when no cluster formed at the cluster-forming threshold, so no permutation ran) and the per-cluster details.

Raises

KeyboardInterrupt If stop_callback returns True during execution. ValueError If the subjects' NIfTIs do not share one grid (the message names the offending subject). FileNotFoundError If a subject's MNI-space NIfTI (per config.nifti_file_pattern) or fsaverage projection is missing.

Notes

Permutation p-values are (b + 1) / (m + 1), so with 1000 permutations the floor is 1/1001; a result is never exactly 0.

Re-using an analysis_name replaces the earlier run: its folder is deleted before the inputs load, together with any file put there by hand, so a re-run that then fails leaves nothing under that name. The app asks before it submits such a run; a script does not, so choose a new name to keep an earlier result.

Examples

from tit.stats import GroupComparisonConfig, run_group_comparison subjects = GroupComparisonConfig.load_subjects("subjects.csv") # doctest: +SKIP cfg = GroupComparisonConfig( ... analysis_name="active_vs_sham", subjects=subjects, ... test_type="unpaired", n_permutations=1000, tissue_type="grey", ... ) # doctest: +SKIP res = run_group_comparison(cfg) # doctest: +SKIP res.success, res.n_significant_clusters, res.output_dir # doctest: +SKIP (True, 2, '.../derivatives/ti-toolbox/stats/group_comparison/active_vs_sham')

See Also

GroupComparisonConfig : The configuration consumed here. GroupComparisonResult : The returned container. run_correlation : Brain-behaviour correlation with the same correction.

Source code in tit/stats/permutation.py

def run_group_comparison(
    config: GroupComparisonConfig,
    callback_handler: logging.Handler | None = None,
    stop_callback: Callable[[], bool] | None = None,
) -> GroupComparisonResult:
    """Run cluster-based permutation testing for a two-group comparison.

    Loads responder and non-responder field maps, performs voxelwise (or,
    for ``config.space == "fsaverage"``, vertexwise) t-tests, applies
    cluster-based permutation correction, generates diagnostic plots, and
    saves all outputs under
    ``<project>/derivatives/ti-toolbox/stats/group_comparison/<analysis_name>/``.

    Parameters
    ----------
    config : GroupComparisonConfig
        Fully specified group comparison configuration.
    callback_handler : logging.Handler or None, optional
        Extra handler attached to the run-scoped logger (used by the app
        to stream log lines).
    stop_callback : callable or None, optional
        Zero-argument callable returning ``True`` to request early
        termination.  Checked between pipeline stages.

    Returns
    -------
    GroupComparisonResult
        Summary of the outcome: output directory, subject counts,
        significant voxel/cluster counts, the permutation-derived cluster
        threshold (NaN when no cluster formed at the cluster-forming
        threshold, so no permutation ran) and the per-cluster details.

    Raises
    ------
    KeyboardInterrupt
        If *stop_callback* returns ``True`` during execution.
    ValueError
        If the subjects' NIfTIs do not share one grid (the message names
        the offending subject).
    FileNotFoundError
        If a subject's MNI-space NIfTI (per ``config.nifti_file_pattern``)
        or fsaverage projection is missing.

    Notes
    -----
    Permutation p-values are ``(b + 1) / (m + 1)``, so with 1000
    permutations the floor is ``1/1001``; a result is never exactly 0.

    Re-using an ``analysis_name`` replaces the earlier run: its folder is
    deleted before the inputs load, together with any file put there by
    hand, so a re-run that then fails leaves nothing under that name. The
    app asks before it submits such a run; a script does not, so choose a
    new name to keep an earlier result.

    Examples
    --------
    >>> from tit.stats import GroupComparisonConfig, run_group_comparison
    >>> subjects = GroupComparisonConfig.load_subjects("subjects.csv")  # doctest: +SKIP
    >>> cfg = GroupComparisonConfig(
    ...     analysis_name="active_vs_sham", subjects=subjects,
    ...     test_type="unpaired", n_permutations=1000, tissue_type="grey",
    ... )  # doctest: +SKIP
    >>> res = run_group_comparison(cfg)  # doctest: +SKIP
    >>> res.success, res.n_significant_clusters, res.output_dir  # doctest: +SKIP
    (True, 2, '.../derivatives/ti-toolbox/stats/group_comparison/active_vs_sham')

    See Also
    --------
    GroupComparisonConfig : The configuration consumed here.
    GroupComparisonResult : The returned container.
    run_correlation : Brain-behaviour correlation with the same correction.
    """
    from tit.telemetry import track_operation
    from tit import constants as _const

    with track_operation(_const.TELEMETRY_OP_STATS):
        try:
            if config.space == GroupComparisonConfig.AnalysisSpace.FSAVERAGE:
                from tit.stats.surface import run_surface_group_comparison

                return run_surface_group_comparison(
                    config, callback_handler, stop_callback
                )
            return _run_group_comparison_inner(config, callback_handler, stop_callback)
        except KeyboardInterrupt:
            raise
        except Exception as exc:
            _log_failure("group_comparison", exc)
            raise

tit.stats.permutation.run_correlation

run_correlation(config: CorrelationConfig, callback_handler: Handler | None = None, stop_callback: Callable[[], bool] | None = None) -> CorrelationResult

Run cluster-based permutation testing for a brain-behaviour correlation.

Correlates each voxel's field value with a continuous per-subject measure (Subject.effect_size), Pearson or Spearman, with cluster-based permutation correction (ACES-style). Outputs go under <project>/derivatives/ti-toolbox/stats/correlation/<analysis_name>/.

Parameters

config : CorrelationConfig Fully specified correlation configuration. callback_handler : logging.Handler or None, optional Extra handler attached to the run-scoped logger. stop_callback : callable or None, optional Zero-argument callable returning True to request early termination.

Returns

CorrelationResult Output directory, subject count, significant voxel/cluster counts, cluster threshold (NaN when no cluster formed at the cluster-forming threshold, so no permutation ran) and per-cluster details.

Raises

KeyboardInterrupt If stop_callback returns True during execution. FileNotFoundError If a subject's NIfTI is missing.

Notes

Re-using an analysis_name replaces the earlier run, as in :func:run_group_comparison: its folder is deleted before the inputs load, so a re-run that then fails leaves nothing under that name.

Examples

from tit.stats import CorrelationConfig, run_correlation subjects = CorrelationConfig.load_subjects("subjects_continuous.csv") # doctest: +SKIP cfg = CorrelationConfig(analysis_name="dose_response", subjects=subjects, ... correlation_type="spearman") # doctest: +SKIP run_correlation(cfg).n_significant_clusters # doctest: +SKIP 1

See Also

CorrelationConfig : The configuration consumed here. run_group_comparison : Two-group comparison with the same correction.

Source code in tit/stats/permutation.py

def run_correlation(
    config: CorrelationConfig,
    callback_handler: logging.Handler | None = None,
    stop_callback: Callable[[], bool] | None = None,
) -> CorrelationResult:
    """Run cluster-based permutation testing for a brain-behaviour correlation.

    Correlates each voxel's field value with a continuous per-subject
    measure (``Subject.effect_size``), Pearson or Spearman, with
    cluster-based permutation correction (ACES-style).  Outputs go under
    ``<project>/derivatives/ti-toolbox/stats/correlation/<analysis_name>/``.

    Parameters
    ----------
    config : CorrelationConfig
        Fully specified correlation configuration.
    callback_handler : logging.Handler or None, optional
        Extra handler attached to the run-scoped logger.
    stop_callback : callable or None, optional
        Zero-argument callable returning ``True`` to request early
        termination.

    Returns
    -------
    CorrelationResult
        Output directory, subject count, significant voxel/cluster counts,
        cluster threshold (NaN when no cluster formed at the cluster-forming
        threshold, so no permutation ran) and per-cluster details.

    Raises
    ------
    KeyboardInterrupt
        If *stop_callback* returns ``True`` during execution.
    FileNotFoundError
        If a subject's NIfTI is missing.

    Notes
    -----
    Re-using an ``analysis_name`` replaces the earlier run, as in
    :func:`run_group_comparison`: its folder is deleted before the inputs
    load, so a re-run that then fails leaves nothing under that name.

    Examples
    --------
    >>> from tit.stats import CorrelationConfig, run_correlation
    >>> subjects = CorrelationConfig.load_subjects("subjects_continuous.csv")  # doctest: +SKIP
    >>> cfg = CorrelationConfig(analysis_name="dose_response", subjects=subjects,
    ...                         correlation_type="spearman")  # doctest: +SKIP
    >>> run_correlation(cfg).n_significant_clusters  # doctest: +SKIP
    1

    See Also
    --------
    CorrelationConfig : The configuration consumed here.
    run_group_comparison : Two-group comparison with the same correction.
    """
    try:
        return _run_correlation_inner(config, callback_handler, stop_callback)
    except KeyboardInterrupt:
        raise
    except Exception as exc:
        _log_failure("correlation", exc)
        raise

tit.stats.engine.PermutationEngine

PermutationEngine(*, cluster_threshold: float = 0.05, n_permutations: int = 1000, alpha: float = 0.05, cluster_stat: str = 'mass', alternative: str = 'two-sided', n_jobs: int = -1, seed: int | None = None, log: Logger | None = None)

Master class for cluster-based permutation testing.

Stores all control parameters as self.* so methods only need data arrays.

Source code in tit/stats/engine.py

def __init__(
    self,
    *,
    cluster_threshold: float = 0.05,
    n_permutations: int = 1000,
    alpha: float = 0.05,
    cluster_stat: str = "mass",
    alternative: str = "two-sided",
    n_jobs: int = -1,
    seed: int | None = None,
    log: logging.Logger | None = None,
):
    self.cluster_threshold = cluster_threshold
    self.n_permutations = n_permutations
    self.alpha = alpha
    self.cluster_stat = cluster_stat
    self.alternative = alternative
    self.n_jobs = n_jobs
    self.seed = seed  # None: fresh OS entropy; an int reproduces the null
    self._log = log or logger

correct_groups

correct_groups(responders, non_responders, *, p_values, t_statistics, valid_mask, test_type: str = 'unpaired', perm_log_file: str | None = None, subject_ids_resp: list | None = None, subject_ids_non_resp: list | None = None) -> tuple

Cluster-based permutation correction for group comparison.

Returns (sig_mask, threshold, sig_clusters, null_dist, observed_clusters, correlation_data).

Source code in tit/stats/engine.py

def correct_groups(
    self,
    responders,
    non_responders,
    *,
    p_values,
    t_statistics,
    valid_mask,
    test_type: str = "unpaired",
    perm_log_file: str | None = None,
    subject_ids_resp: list | None = None,
    subject_ids_non_resp: list | None = None,
) -> tuple:
    """Cluster-based permutation correction for group comparison.

    Returns ``(sig_mask, threshold, sig_clusters, null_dist, observed_clusters,
    correlation_data)``.
    """
    from .io_utils import save_permutation_details

    if self.cluster_stat == "mass" and t_statistics is None:
        raise ValueError("t_statistics required when cluster_stat='mass'")

    self._log.info(
        "Cluster-based permutation correction (%s, %s)",
        self.cluster_stat,
        self.alternative,
    )

    initial_mask = (p_values < self.cluster_threshold) & valid_mask
    labeled_array, n_clusters = label_signed(
        initial_mask, t_statistics, self.alternative
    )
    self._log.info(
        "Clusters at p<%.3f (uncorrected, sign-separated): %d",
        self.cluster_threshold,
        n_clusters,
    )

    empty = {"sizes": np.array([]), "masses": np.array([])}
    if n_clusters == 0:
        self._log.warning(
            "No clusters found. Consider increasing cluster_threshold."
        )
        return (
            np.zeros_like(p_values, dtype=int),
            0,
            [],
            np.array([]),
            [],
            empty,
        )

    # Pre-extract voxel data
    all_data = np.concatenate([responders, non_responders], axis=-1)
    n_resp = responders.shape[-1]
    n_total = n_resp + non_responders.shape[-1]

    test_coords = np.argwhere(valid_mask)
    n_test = len(test_coords)
    self._log.info("Pre-extracting %d voxels, %d subjects", n_test, n_total)

    # float64: a paired relabelling (mean +/- diff) of float32 values is then exact, so a
    # constant group stays constant instead of gaining rounding-noise variance.
    test_data = np.zeros((n_test, n_total), dtype=np.float64)
    for idx, (i, j, k) in enumerate(test_coords):
        test_data[idx, :] = all_data[i, j, k, :]
    del all_data
    gc.collect()

    self._log.info("Test data: %.1f MB", test_data.nbytes / (1024**2))

    actual_jobs = resolve_n_jobs(self.n_jobs)
    self._log.info(
        "Running %d permutations on %d cores",
        self.n_permutations,
        actual_jobs,
    )

    seeds = np.random.default_rng(self.seed).integers(
        0, 2**31, size=self.n_permutations
    )
    track = (
        perm_log_file is not None
        and subject_ids_resp is not None
        and subject_ids_non_resp is not None
    )

    if actual_jobs == 1:
        results = [
            _run_single_permutation(
                test_data,
                test_coords,
                n_resp,
                n_total,
                self.cluster_threshold,
                valid_mask,
                p_values.shape,
                test_type=test_type,
                alternative=self.alternative,
                cluster_stat=self.cluster_stat,
                seed=seeds[i],
                return_indices=track,
            )
            for i in range(self.n_permutations)
        ]
    else:
        results = Parallel(n_jobs=actual_jobs, verbose=0)(
            delayed(_run_single_permutation)(
                test_data,
                test_coords,
                n_resp,
                n_total,
                self.cluster_threshold,
                valid_mask,
                p_values.shape,
                test_type=test_type,
                alternative=self.alternative,
                cluster_stat=self.cluster_stat,
                seed=seeds[i],
                return_indices=track,
            )
            for i in range(self.n_permutations)
        )
        gc.collect()

    del test_data
    gc.collect()

    if track:
        null_stats = np.array([r[0] for r in results])
        perm_indices = [r[1] for r in results]
        null_sizes = [r[2] for r in results]
        null_masses = [r[3] for r in results]
    else:
        null_stats = np.array([r[0] for r in results])
        perm_indices = None
        null_sizes = [r[1] for r in results]
        null_masses = [r[2] for r in results]

    # Determine threshold
    sorted_null = np.sort(null_stats)[::-1]
    ti = max(1, min(int(self.alpha * self.n_permutations), len(sorted_null)))
    cluster_stat_threshold = sorted_null[ti - 1]

    stat_unit = "voxels" if self.cluster_stat == "size" else "mass units"
    self._log.info(
        "Threshold (p<%.3f): %.2f %s  " "(null min=%.2f, mean=%.2f, max=%.2f)",
        self.alpha,
        cluster_stat_threshold,
        stat_unit,
        np.min(null_stats),
        np.mean(null_stats),
        np.max(null_stats),
    )

    if perm_log_file is not None and perm_indices is not None:
        info = [
            {
                "perm_num": i,
                "perm_idx": perm_indices[i],
                "max_cluster_size": null_stats[i],
            }
            for i in range(self.n_permutations)
        ]
        save_permutation_details(
            info,
            perm_log_file,
            subject_ids_resp,
            subject_ids_non_resp,
            paired=test_type == "paired",
        )
        self._log.info("Permutation log saved: %s", perm_log_file)

    # Identify significant clusters (MNE-style per-cluster p-values)
    sig_mask, sig_clusters, all_observed = _identify_significant_clusters(
        labeled_array,
        n_clusters,
        t_statistics,
        null_stats,
        self.cluster_stat,
        self.alpha,
        self.alternative,
        self._log,
    )

    self._log.info(
        "Significant: %d clusters, %d voxels",
        len(sig_clusters),
        np.sum(sig_mask),
    )

    correlation_data = {
        "sizes": np.array(null_sizes),
        "masses": np.array(null_masses),
    }
    return (
        sig_mask,
        cluster_stat_threshold,
        sig_clusters,
        null_stats,
        all_observed,
        correlation_data,
    )

correct_correlation

correct_correlation(subject_data, effect_sizes, *, r_values, t_statistics, p_values, valid_mask, correlation_type: str = 'pearson', weights=None, perm_log_file: str | None = None, subject_ids: list | None = None) -> tuple

Cluster-based permutation correction for correlation analysis.

Returns (sig_mask, threshold, sig_clusters, null_dist, observed_clusters, correlation_data).

Source code in tit/stats/engine.py

def correct_correlation(
    self,
    subject_data,
    effect_sizes,
    *,
    r_values,
    t_statistics,
    p_values,
    valid_mask,
    correlation_type: str = "pearson",
    weights=None,
    perm_log_file: str | None = None,
    subject_ids: list | None = None,
) -> tuple:
    """Cluster-based permutation correction for correlation analysis.

    Returns ``(sig_mask, threshold, sig_clusters, null_dist, observed_clusters,
    correlation_data)``.
    """
    from .io_utils import save_correlation_permutation_details

    effect_sizes = np.asarray(effect_sizes, dtype=np.float64)
    n_subjects = len(effect_sizes)

    self._log.info(
        "Correlation cluster correction (%s, %s, %s)",
        correlation_type,
        self.cluster_stat,
        self.alternative,
    )

    # Form initial clusters based on alternative
    initial_mask = (p_values < self.cluster_threshold) & valid_mask
    labeled_array, n_clusters = label_signed(
        initial_mask, t_statistics, self.alternative
    )
    self._log.info("Clusters at p<%.3f: %d", self.cluster_threshold, n_clusters)

    empty = {"sizes": np.array([]), "masses": np.array([])}
    if n_clusters == 0:
        self._log.warning("No clusters found.")
        return (
            np.zeros_like(p_values, dtype=int),
            0,
            [],
            np.array([]),
            [],
            empty,
        )

    # Pre-extract voxel data
    valid_coords = np.argwhere(valid_mask)
    n_valid = len(valid_coords)
    voxel_data = np.zeros((n_valid, n_subjects), dtype=np.float64)
    for idx, (i, j, k) in enumerate(valid_coords):
        voxel_data[idx, :] = subject_data[i, j, k, :]

    # Pre-rank for Spearman
    preranked = False
    if correlation_type == "spearman":
        self._log.info("Pre-ranking voxel data for Spearman (%d voxels)", n_valid)
        voxel_data = np.apply_along_axis(rankdata, 1, voxel_data)
        preranked = True

    actual_jobs = resolve_n_jobs(self.n_jobs)
    self._log.info(
        "Running %d permutations on %d cores",
        self.n_permutations,
        actual_jobs,
    )

    seeds = np.random.default_rng(self.seed).integers(
        0, 2**31, size=self.n_permutations
    )
    track = perm_log_file is not None and subject_ids is not None

    if actual_jobs == 1:
        results = [
            _run_single_correlation_permutation(
                voxel_data,
                effect_sizes,
                valid_coords,
                self.cluster_threshold,
                valid_mask,
                p_values.shape,
                correlation_type=correlation_type,
                weights=weights,
                cluster_stat=self.cluster_stat,
                alternative=self.alternative,
                seed=seeds[i],
                return_indices=track,
                voxel_data_preranked=preranked,
            )
            for i in range(self.n_permutations)
        ]
    else:
        results = Parallel(n_jobs=actual_jobs, verbose=0)(
            delayed(_run_single_correlation_permutation)(
                voxel_data,
                effect_sizes,
                valid_coords,
                self.cluster_threshold,
                valid_mask,
                p_values.shape,
                correlation_type=correlation_type,
                weights=weights,
                cluster_stat=self.cluster_stat,
                alternative=self.alternative,
                seed=seeds[i],
                return_indices=track,
                voxel_data_preranked=preranked,
            )
            for i in range(self.n_permutations)
        )
        gc.collect()

    del voxel_data
    gc.collect()

    if track:
        null_stats = np.array([r[0] for r in results])
        perm_indices = [r[1] for r in results]
        null_sizes = [r[2] for r in results]
        null_masses = [r[3] for r in results]
    else:
        null_stats = np.array([r[0] for r in results])
        perm_indices = None
        null_sizes = [r[1] for r in results]
        null_masses = [r[2] for r in results]

    sorted_null = np.sort(null_stats)[::-1]
    ti = max(1, min(int(self.alpha * self.n_permutations), len(sorted_null)))
    cluster_stat_threshold = sorted_null[ti - 1]

    stat_unit = "voxels" if self.cluster_stat == "size" else "mass units"
    self._log.info(
        "Threshold (p<%.3f): %.2f %s  " "(null min=%.2f, mean=%.2f, max=%.2f)",
        self.alpha,
        cluster_stat_threshold,
        stat_unit,
        np.min(null_stats),
        np.mean(null_stats),
        np.max(null_stats),
    )

    if perm_indices is not None:
        info = [
            {
                "perm_num": i,
                "perm_idx": perm_indices[i],
                "max_cluster_size": null_stats[i],
            }
            for i in range(self.n_permutations)
        ]
        save_correlation_permutation_details(info, perm_log_file, subject_ids)
        self._log.info("Permutation log saved: %s", perm_log_file)

    # Identify significant clusters
    sig_mask, sig_clusters, all_observed = _identify_significant_clusters(
        labeled_array,
        n_clusters,
        t_statistics,
        null_stats,
        self.cluster_stat,
        self.alpha,
        self.alternative,
        self._log,
        r_values=r_values,
    )

    self._log.info(
        "Significant: %d clusters, %d voxels",
        len(sig_clusters),
        np.sum(sig_mask),
    )

    correlation_data = {
        "sizes": np.array(null_sizes),
        "masses": np.array(null_masses),
    }
    return (
        sig_mask,
        cluster_stat_threshold,
        sig_clusters,
        null_stats,
        all_observed,
        correlation_data,
    )

tit.stats.engine.cluster_analysis

cluster_analysis(sig_mask, affine, log=None)

Connected-component analysis with MNI coordinate mapping.

Returns a list of cluster dicts sorted by size (descending).

Source code in tit/stats/engine.py

def cluster_analysis(sig_mask, affine, log=None):
    """Connected-component analysis with MNI coordinate mapping.

    Returns a list of cluster dicts sorted by size (descending).
    """
    import nibabel as nib

    _log = log or logger
    labeled_array, num_clusters = label(sig_mask)

    if num_clusters == 0:
        _log.info("No clusters found in significance mask")
        return []

    clusters = []
    for cid in range(1, num_clusters + 1):
        coords = np.argwhere(labeled_array == cid)
        com_voxel = np.mean(coords, axis=0)
        com_mni = nib.affines.apply_affine(affine, com_voxel)
        clusters.append(
            {
                "cluster_id": cid,
                "size": len(coords),
                "center_voxel": com_voxel,
                "center_mni": com_mni,
            }
        )

    clusters.sort(key=lambda c: c["size"], reverse=True)

    _log.info("Found %d clusters", num_clusters)
    for c in clusters[:10]:
        _log.info(
            "  Cluster %d: %d voxels, MNI (%.1f, %.1f, %.1f)",
            c["cluster_id"],
            c["size"],
            c["center_mni"][0],
            c["center_mni"][1],
            c["center_mni"][2],
        )

    return clusters

tit.stats.nifti.load_subject_nifti_ti_toolbox

load_subject_nifti_ti_toolbox(subject_id: str, simulation_name: str, nifti_file_pattern: str = 'grey_{simulation_name}_TI_MNI_MNI_TI_max.nii.gz', dtype=float32) -> tuple[ndarray, Nifti1Image, str]

Load a single subject's NIfTI file from TI-Toolbox BIDS structure.

Parameters

subject_id : str Subject identifier (e.g. '070'). simulation_name : str Simulation folder name (e.g. 'ICP_RHIPPO'). nifti_file_pattern : str, optional Filename pattern with {subject_id} / {simulation_name} placeholders, in the 2-pair spelling; an mTI simulation's file is read from mTI/niftis/ (see :func:tit.analyzer.field_selector.simulation_nifti_path). dtype : numpy dtype, optional Data type for the returned array. Default is np.float32.

Returns

data : numpy.ndarray 3-D array of voxel values. img : nibabel.Nifti1Image The loaded NIfTI image (useful for affine / header). filepath : str Absolute path of the loaded file.

Raises

FileNotFoundError If the resolved NIfTI path does not exist and is not one the simulation can produce now; a simulator NIfTI it did not save is written first, from the field mesh (:func:tit.volume_fields.ensure_nifti_file).

Source code in tit/stats/nifti.py

def load_subject_nifti_ti_toolbox(
    subject_id: str,
    simulation_name: str,
    nifti_file_pattern: str = "grey_{simulation_name}_TI_MNI_MNI_TI_max.nii.gz",
    dtype=np.float32,
) -> tuple[np.ndarray, nib.Nifti1Image, str]:
    """Load a single subject's NIfTI file from TI-Toolbox BIDS structure.

    Parameters
    ----------
    subject_id : str
        Subject identifier (e.g. ``'070'``).
    simulation_name : str
        Simulation folder name (e.g. ``'ICP_RHIPPO'``).
    nifti_file_pattern : str, optional
        Filename pattern with ``{subject_id}`` / ``{simulation_name}``
        placeholders, in the 2-pair spelling; an mTI simulation's file is read
        from ``mTI/niftis/`` (see
        :func:`tit.analyzer.field_selector.simulation_nifti_path`).
    dtype : numpy dtype, optional
        Data type for the returned array.  Default is ``np.float32``.

    Returns
    -------
    data : numpy.ndarray
        3-D array of voxel values.
    img : nibabel.Nifti1Image
        The loaded NIfTI image (useful for affine / header).
    filepath : str
        Absolute path of the loaded file.

    Raises
    ------
    FileNotFoundError
        If the resolved NIfTI path does not exist and is not one the simulation
        can produce now; a simulator NIfTI it did not save is written first,
        from the field mesh (:func:`tit.volume_fields.ensure_nifti_file`).
    """
    from tit.volume_fields import ensure_nifti_file

    # TI/niftis/ for a 2-pair run, mTI/niftis/ (mTI file names) for an mTI run.
    path = simulation_nifti_path(subject_id, simulation_name, nifti_file_pattern)
    nifti_dir, filepath = str(path.parent), str(path)
    ensure_nifti_file(subject_id, simulation_name, path, say=logger.info)

    # Load the file (inline basic loading)
    if not os.path.exists(filepath):
        # Provide extra context to make debugging path/layout issues easier
        if os.path.isdir(nifti_dir):
            try:
                existing = sorted(os.listdir(nifti_dir))
            except OSError:
                existing = []
            preview = existing[:20]
            suffix = ""
            if len(existing) > len(preview):
                suffix = f" (showing first {len(preview)} of {len(existing)})"
            raise FileNotFoundError(
                f"NIfTI file not found: {filepath}. "
                f"Directory exists: {nifti_dir}. "
                f"Files in directory: {preview}{suffix}"
            )
        raise FileNotFoundError(f"NIfTI file not found: {filepath}")

    img = nib.load(filepath)
    data = img.get_fdata(dtype=dtype)

    # Ensure 3D data (squeeze out extra dimensions if present)
    while data.ndim > 3:
        data = np.squeeze(data, axis=-1)

    return data, img, filepath

tit.stats.nifti.load_group_data_ti_toolbox

load_group_data_ti_toolbox(subject_configs: list[dict], nifti_file_pattern: str = 'grey_{simulation_name}_TI_MNI_MNI_TI_max.nii.gz', dtype=float32, reference: tuple | None = None) -> tuple[ndarray, Nifti1Image, list[str]]

Load and stack multiple subjects into a 4-D array.

Parameters

subject_configs : list of dict Each dict must contain 'subject_id' and 'simulation_name' keys (e.g. {'subject_id': '070', 'simulation_name': 'ICP_RHIPPO'}). nifti_file_pattern : str, optional Filename pattern forwarded to :func:load_subject_nifti_ti_toolbox. dtype : numpy dtype, optional Data type for the returned arrays. Default is np.float32. reference : tuple or None, optional (shape, affine, where) of a grid every subject must be on -- another group's, when two groups are compared voxel by voxel; where names it in the error. None (default) makes the first subject the reference.

Returns

data_4d : numpy.ndarray Shape (X, Y, Z, n_subjects). template_img : nibabel.Nifti1Image Image from the first subject (affine / header reference). subject_ids : list of str Subject identifiers in the same order as the last axis of data_4d.

Raises

ValueError If no subjects could be loaded.

Source code in tit/stats/nifti.py

def load_group_data_ti_toolbox(
    subject_configs: list[dict],
    nifti_file_pattern: str = "grey_{simulation_name}_TI_MNI_MNI_TI_max.nii.gz",
    dtype=np.float32,
    reference: tuple | None = None,
) -> tuple[np.ndarray, nib.Nifti1Image, list[str]]:
    """Load and stack multiple subjects into a 4-D array.

    Parameters
    ----------
    subject_configs : list of dict
        Each dict must contain ``'subject_id'`` and ``'simulation_name'``
        keys (e.g. ``{'subject_id': '070', 'simulation_name': 'ICP_RHIPPO'}``).
    nifti_file_pattern : str, optional
        Filename pattern forwarded to :func:`load_subject_nifti_ti_toolbox`.
    dtype : numpy dtype, optional
        Data type for the returned arrays.  Default is ``np.float32``.
    reference : tuple or None, optional
        ``(shape, affine, where)`` of a grid every subject must be on -- another
        group's, when two groups are compared voxel by voxel; *where* names it in
        the error. ``None`` (default) makes the first subject the reference.

    Returns
    -------
    data_4d : numpy.ndarray
        Shape ``(X, Y, Z, n_subjects)``.
    template_img : nibabel.Nifti1Image
        Image from the first subject (affine / header reference).
    subject_ids : list of str
        Subject identifiers in the same order as the last axis of
        *data_4d*.

    Raises
    ------
    ValueError
        If no subjects could be loaded.
    """
    if len(subject_configs) == 0:
        raise ValueError("No subjects could be loaded successfully")

    data_4d = None
    subject_ids: list[str] = []
    template_affine = None
    template_header = None
    reference_path = None

    for index, config in enumerate(subject_configs):
        subject_id = config["subject_id"]
        simulation_name = config["simulation_name"]

        data, img, filepath = load_subject_nifti_ti_toolbox(
            subject_id, simulation_name, nifti_file_pattern, dtype=dtype
        )

        if reference is not None:
            _check_same_grid(subject_id, filepath, data.shape, img.affine, *reference)

        if data_4d is None:
            template_affine = img.affine.copy()
            template_header = img.header.copy()
            reference_path = filepath
            # Preallocate: np.stack + astype on a list of N volumes peaks at
            # ~3x the final array; writing each subject straight into the
            # output keeps it at 1x + one volume.
            data_4d = np.empty(data.shape + (len(subject_configs),), dtype=dtype)
        else:
            _check_same_grid(
                subject_id,
                filepath,
                data.shape,
                img.affine,
                data_4d.shape[:-1],
                template_affine,
                reference_path,
            )

        data_4d[..., index] = data
        subject_ids.append(subject_id)

        # Clear the image object to free memory
        del img, data

    # Recreate minimal template image
    template_img = nib.Nifti1Image(data_4d[..., 0], template_affine, template_header)

    gc.collect()

    return data_4d, template_img, subject_ids