engine

tit.stats.engine

Core statistical computation engine for cluster-based permutation testing.

Provides vectorised implementations of voxelwise t-tests, correlations, cluster-based permutation correction, connected-component cluster analysis, and MNE-style p-value computation.

All orchestration functions accept an optional logger. If None the module-level logger (tit.stats.engine) is used.

See Also

tit.stats.permutation : High-level orchestration that drives this engine. tit.stats.config : Configuration dataclasses consumed by the 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, 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,
    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._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)

    test_data = np.zeros((n_test, n_total), dtype=np.float32)
    for idx, (i, j, k) in enumerate(test_coords):
        test_data[idx, :] = all_data[i, j, k, :].astype(np.float32)
    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.randint(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,
        )
        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_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.randint(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])
        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])
        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),
    )

    # 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,
    )

pval_from_histogram

pval_from_histogram(observed_stats, null_distribution, tail=0, sampled=True)

Compute per-cluster p-values from a permutation null distribution.

Implements the MNE-Python approach based on Maris & Oostenveld (2007).

Parameters

observed_stats : array-like Observed cluster-level statistics. null_distribution : array-like Max-cluster statistics from each permutation. tail : {0, 1, -1}, optional Tail of the test: 0 for two-sided, 1 for greater, -1 for less. sampled : bool, optional True (default) when null_distribution is a random sample of m permutations rather than the exhaustive enumeration of the permutation group. The p-value is then the unbiased, valid estimator (b + 1) / (m + 1) of Phipson & Smyth (2010), where b counts the permutations at least as extreme as the observation. Adding the observed statistic to its own null is what makes the test exact: the naive b / m can return p = 0, which is not a valid p-value and makes the test anti-conservative at the very tail where cluster inference lives. Pass sampled=False only when null_distribution really is the complete enumeration of every possible relabelling, in which case b / m is exact.

Returns

numpy.ndarray P-values, one per observed cluster. With sampled=True the smallest attainable p-value is 1 / (m + 1), never zero.

References

Maris, E. & Oostenveld, R. (2007). J. Neurosci. Methods 164(1), 177-190. Phipson, B. & Smyth, G. K. (2010). Stat. Appl. Genet. Mol. Biol. 9(1), 39.

Source code in tit/stats/engine.py

def pval_from_histogram(observed_stats, null_distribution, tail=0, sampled=True):
    """Compute per-cluster p-values from a permutation null distribution.

    Implements the MNE-Python approach based on Maris & Oostenveld (2007).

    Parameters
    ----------
    observed_stats : array-like
        Observed cluster-level statistics.
    null_distribution : array-like
        Max-cluster statistics from each permutation.
    tail : {0, 1, -1}, optional
        Tail of the test: 0 for two-sided, 1 for greater, -1 for less.
    sampled : bool, optional
        ``True`` (default) when *null_distribution* is a **random sample** of
        ``m`` permutations rather than the exhaustive enumeration of the
        permutation group.  The p-value is then the unbiased, valid estimator
        ``(b + 1) / (m + 1)`` of Phipson & Smyth (2010), where ``b`` counts the
        permutations at least as extreme as the observation.  Adding the
        observed statistic to its own null is what makes the test exact: the
        naive ``b / m`` can return ``p = 0``, which is not a valid p-value and
        makes the test anti-conservative at the very tail where cluster
        inference lives.  Pass ``sampled=False`` only when *null_distribution*
        really is the complete enumeration of every possible relabelling, in
        which case ``b / m`` is exact.

    Returns
    -------
    numpy.ndarray
        P-values, one per observed cluster.  With ``sampled=True`` the smallest
        attainable p-value is ``1 / (m + 1)``, never zero.

    References
    ----------
    Maris, E. & Oostenveld, R. (2007). J. Neurosci. Methods 164(1), 177-190.
    Phipson, B. & Smyth, G. K. (2010). Stat. Appl. Genet. Mol. Biol. 9(1), 39.
    """
    observed_stats = np.atleast_1d(observed_stats)
    null_distribution = np.asarray(null_distribution)
    m = null_distribution.size

    if tail == -1:
        counts = np.array(
            [np.sum(null_distribution <= obs) for obs in observed_stats], dtype=float
        )
    elif tail == 1:
        counts = np.array(
            [np.sum(null_distribution >= obs) for obs in observed_stats], dtype=float
        )
    else:
        counts = np.array(
            [
                np.sum(np.abs(null_distribution) >= np.abs(obs))
                for obs in observed_stats
            ],
            dtype=float,
        )

    if sampled:
        return (counts + 1.0) / (m + 1.0)
    if m == 0:
        return np.ones_like(counts)
    return counts / m

tail_from_alternative

tail_from_alternative(alternative)

MNE-style tail code for an alternative string.

Source code in tit/stats/engine.py

def tail_from_alternative(alternative):
    """MNE-style tail code for an ``alternative`` string."""
    return {"greater": 1, "less": -1}.get(alternative, 0)

label_signed

label_signed(mask, t_stats, alternative)

Connected components that never merge voxels of opposite t sign.

scipy.ndimage.label on a bare significance mask happily fuses a positive and a negative supra-threshold blob that happen to touch; the signed cluster mass of the fused object is then the difference of two real effects and can be ~0. Cluster inference requires clusters to be sign-homogeneous, so positive and negative voxels are labelled separately and the negative labels offset past the positive ones.

One-sided alternatives keep only the relevant sign, matching the correlation path and MNE's tail=+/-1 behaviour.

Returns (labelled_array, n_clusters).

Source code in tit/stats/engine.py

def label_signed(mask, t_stats, alternative):
    """Connected components that never merge voxels of opposite ``t`` sign.

    ``scipy.ndimage.label`` on a bare significance mask happily fuses a
    positive and a negative supra-threshold blob that happen to touch; the
    signed cluster mass of the fused object is then the *difference* of two
    real effects and can be ~0.  Cluster inference requires clusters to be
    sign-homogeneous, so positive and negative voxels are labelled separately
    and the negative labels offset past the positive ones.

    One-sided alternatives keep only the relevant sign, matching the
    correlation path and MNE's ``tail=+/-1`` behaviour.

    Returns ``(labelled_array, n_clusters)``.
    """
    if alternative == "greater":
        return label(mask & (t_stats > 0))
    if alternative == "less":
        return label(mask & (t_stats < 0))
    pos, n_pos = label(mask & (t_stats > 0))
    neg, n_neg = label(mask & (t_stats < 0))
    if n_neg:
        neg_sel = neg > 0
        pos[neg_sel] = neg[neg_sel] + n_pos
    return pos, n_pos + n_neg

tail_statistic

tail_statistic(sizes, masses, cluster_stat, tail)

Orient the cluster statistic so that larger is more extreme.

The null distribution is a distribution of maxima, so observed and permuted statistics must be measured on the same, monotone "extremeness" scale:

size -- always non-negative; clusters are already sign-restricted, so the raw size is the extremeness. mass -- signed. tail=+1 uses the mass, tail=-1 its negation, tail=0 its absolute value.

Taking a plain max() of signed masses under a left or two-sided tail (the v2.x behaviour) selects the mass closest to zero among negative clusters, i.e. the least extreme one, and so produces a null that is far too small.

Source code in tit/stats/engine.py

def tail_statistic(sizes, masses, cluster_stat, tail):
    """Orient the cluster statistic so that *larger is more extreme*.

    The null distribution is a distribution of maxima, so observed and permuted
    statistics must be measured on the same, monotone "extremeness" scale:

    ``size``   -- always non-negative; clusters are already sign-restricted, so
                  the raw size is the extremeness.
    ``mass``   -- signed.  ``tail=+1`` uses the mass, ``tail=-1`` its negation,
                  ``tail=0`` its absolute value.

    Taking a plain ``max()`` of signed masses under a left or two-sided tail
    (the v2.x behaviour) selects the mass *closest to zero* among negative
    clusters, i.e. the least extreme one, and so produces a null that is far
    too small.
    """
    if cluster_stat == "size":
        return np.abs(np.asarray(sizes, dtype=float))
    masses = np.asarray(masses, dtype=float)
    if tail == 1:
        return masses
    if tail == -1:
        return -masses
    return np.abs(masses)

correlation

correlation(voxel_data, effect_sizes, correlation_type='pearson', weights=None, voxel_data_preranked=False)

Vectorised correlation for all voxels at once.

Returns (r_values, t_values, p_values) – each shape (n_voxels,).

Source code in tit/stats/engine.py

def correlation(
    voxel_data,
    effect_sizes,
    correlation_type="pearson",
    weights=None,
    voxel_data_preranked=False,
):
    """Vectorised correlation for all voxels at once.

    Returns (r_values, t_values, p_values) – each shape ``(n_voxels,)``.
    """
    n_voxels, n_subjects = voxel_data.shape

    if correlation_type == "spearman":
        if not voxel_data_preranked:
            voxel_data = np.apply_along_axis(rankdata, 1, voxel_data)
        effect_sizes = rankdata(effect_sizes)

    if weights is None:
        x_mean = np.mean(voxel_data, axis=1, keepdims=True)
        y_mean = np.mean(effect_sizes)
        x_centered = voxel_data - x_mean
        y_centered = effect_sizes - y_mean
        cov_xy = np.sum(x_centered * y_centered, axis=1) / (n_subjects - 1)
        std_x = np.std(voxel_data, axis=1, ddof=1)
        std_y = np.std(effect_sizes, ddof=1)
    else:
        weights = np.asarray(weights, dtype=np.float64)
        weight_sum = np.sum(weights)
        x_mean = np.sum(weights * voxel_data, axis=1, keepdims=True) / weight_sum
        y_mean = np.sum(weights * effect_sizes) / weight_sum
        x_centered = voxel_data - x_mean
        y_centered = effect_sizes - y_mean
        cov_xy = np.sum(weights * x_centered * y_centered, axis=1) / weight_sum
        var_x = np.sum(weights * (x_centered**2), axis=1) / weight_sum
        var_y = np.sum(weights * (y_centered**2)) / weight_sum
        std_x = np.sqrt(var_x)
        std_y = np.sqrt(var_y)

    denom = std_x * std_y
    valid_mask = denom > 1e-10

    r_values = np.zeros(n_voxels)
    r_values[valid_mask] = cov_xy[valid_mask] / denom[valid_mask]
    r_values = np.clip(r_values, -1.0, 1.0)

    df = n_subjects - 2
    r_squared = r_values**2
    denom_t = 1 - r_squared
    denom_t[denom_t < 1e-10] = 1e-10
    t_values = r_values * np.sqrt(df / denom_t)
    p_values = 2 * sp_stats.t.sf(np.abs(t_values), df)

    return r_values, t_values, p_values

correlation_voxelwise

correlation_voxelwise(subject_data, effect_sizes, weights=None, correlation_type='pearson', log=None)

Voxelwise correlation between E-field and continuous outcome (ACES-style).

Returns (r_values, t_statistics, p_values, valid_mask) – each (x, y, z).

Source code in tit/stats/engine.py

def correlation_voxelwise(
    subject_data,
    effect_sizes,
    weights=None,
    correlation_type="pearson",
    log=None,
):
    """Voxelwise correlation between E-field and continuous outcome (ACES-style).

    Returns (r_values, t_statistics, p_values, valid_mask) – each ``(x, y, z)``.
    """
    _log = log or logger
    n_subjects = subject_data.shape[-1]
    if len(effect_sizes) != n_subjects:
        raise ValueError(
            f"effect_sizes length ({len(effect_sizes)}) != n_subjects ({n_subjects})"
        )
    if weights is not None and len(weights) != n_subjects:
        raise ValueError(
            f"weights length ({len(weights)}) != n_subjects ({n_subjects})"
        )
    if n_subjects < 3:
        raise ValueError(f"Need >= 3 subjects for correlation, got {n_subjects}")

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

    weight_str = " (weighted)" if weights is not None else ""
    _log.info(
        "Voxelwise %s correlation%s – %d subjects",
        correlation_type.capitalize(),
        weight_str,
        n_subjects,
    )

    shape = subject_data.shape[:3]
    r_values = np.zeros(shape)
    t_statistics = np.zeros(shape)
    p_values = np.ones(shape)

    valid_mask = np.any(subject_data > 0, axis=-1)
    valid_coords = np.argwhere(valid_mask)
    n_valid = len(valid_coords)
    _log.info("Valid voxels: %d", n_valid)

    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, :]

    r_1d, t_1d, p_1d = correlation(
        voxel_data, effect_sizes, correlation_type=correlation_type, weights=weights
    )

    for idx, (i, j, k) in enumerate(valid_coords):
        r_values[i, j, k] = r_1d[idx]
        t_statistics[i, j, k] = t_1d[idx]
        p_values[i, j, k] = p_1d[idx]

    _log.info(
        "r range [%.4f, %.4f], mean |r| %.4f",
        np.min(r_values[valid_mask]),
        np.max(r_values[valid_mask]),
        np.mean(np.abs(r_values[valid_mask])),
    )
    _log.info(
        "Uncorrected p<0.05: %d,  p<0.01: %d",
        np.sum((p_values < 0.05) & valid_mask),
        np.sum((p_values < 0.01) & valid_mask),
    )

    return r_values, t_statistics, p_values, valid_mask

ttest_ind

ttest_ind(test_data, n_resp, n_non_resp, alternative='two-sided')

Vectorised independent-samples t-test. Returns (t_stats, p_values).

Source code in tit/stats/engine.py

def ttest_ind(test_data, n_resp, n_non_resp, alternative="two-sided"):
    """Vectorised independent-samples t-test. Returns (t_stats, p_values)."""
    resp_data = test_data[:, :n_resp]
    non_resp_data = test_data[:, n_resp:]

    resp_means = np.mean(resp_data, axis=1)
    non_resp_means = np.mean(non_resp_data, axis=1)

    # Pooled variance from each group's SUM OF SQUARED DEVIATIONS, not from
    # ``(n - 1) * np.var(..., ddof=1)``.
    #
    # The two are algebraically the same for n >= 2, but a group of ONE subject makes
    # ``np.var(x, ddof=1)`` a 0/0 ``nan``, and ``(n - 1) * nan`` is ``0 * nan == nan``, not the
    # 0 the pooled estimator calls for. So every voxel of a 2-vs-1 design came out ``nan``,
    # was counted as degenerate by ``ttest_voxelwise`` ("... 0 perfectly separated, N constant"),
    # emptied ``valid_mask``, and the run died with "No voxel could be tested" having written
    # nothing but its log — a design that is under-powered but perfectly well defined
    # (df = n1 + n2 - 2 = 1) reported as data with no variance in it. Summing the deviations
    # gives the singleton group its true contribution, which is exactly zero.
    resp_ss = np.sum((resp_data - resp_means[:, None]) ** 2, axis=1)
    non_resp_ss = np.sum((non_resp_data - non_resp_means[:, None]) ** 2, axis=1)

    denominator = n_resp + n_non_resp - 2
    pooled_vars = (resp_ss + non_resp_ss) / denominator

    se_diff = np.sqrt(pooled_vars * (1 / n_resp + 1 / n_non_resp))
    t_stats = _safe_t(resp_means - non_resp_means, se_diff)

    df = n_resp + n_non_resp - 2
    match alternative:
        case "two-sided":
            p_values = 2 * sp_stats.t.sf(np.abs(t_stats), df)
        case "greater":
            p_values = sp_stats.t.sf(t_stats, df)
        case "less":
            p_values = sp_stats.t.sf(-t_stats, df)
        case _:
            raise ValueError("alternative must be 'two-sided', 'greater', or 'less'")

    return t_stats, p_values

ttest_rel

ttest_rel(test_data, n_resp, alternative='two-sided')

Vectorised paired-samples t-test. Returns (t_stats, p_values).

Source code in tit/stats/engine.py

def ttest_rel(test_data, n_resp, alternative="two-sided"):
    """Vectorised paired-samples t-test. Returns (t_stats, p_values)."""
    resp_data = test_data[:, :n_resp]
    non_resp_data = test_data[:, n_resp:]

    diff = resp_data - non_resp_data
    diff_means = np.mean(diff, axis=1)
    diff_stds = np.std(diff, axis=1, ddof=1)

    se = diff_stds / np.sqrt(n_resp)
    t_stats = _safe_t(diff_means, se)

    df = n_resp - 1
    match alternative:
        case "two-sided":
            p_values = 2 * sp_stats.t.sf(np.abs(t_stats), df)
        case "greater":
            p_values = sp_stats.t.sf(t_stats, df)
        case "less":
            p_values = sp_stats.t.sf(-t_stats, df)
        case _:
            raise ValueError("alternative must be 'two-sided', 'greater', or 'less'")

    return t_stats, p_values

ttest_voxelwise

ttest_voxelwise(responders, non_responders, test_type='unpaired', alternative='two-sided', log=None)

Voxelwise t-test. Returns (p_values, t_statistics, valid_mask).

Source code in tit/stats/engine.py

def ttest_voxelwise(
    responders,
    non_responders,
    test_type="unpaired",
    alternative="two-sided",
    log=None,
):
    """Voxelwise t-test. Returns (p_values, t_statistics, valid_mask)."""
    _log = log or logger

    if test_type not in ("paired", "unpaired"):
        raise ValueError("test_type must be 'paired' or 'unpaired'")

    if test_type == "paired" and responders.shape[-1] != non_responders.shape[-1]:
        raise ValueError(
            f"Paired test requires equal sample sizes: "
            f"{responders.shape[-1]} vs {non_responders.shape[-1]}"
        )

    test_name = "Paired" if test_type == "paired" else "Unpaired"
    _log.info("Voxelwise %s t-test (alternative=%s)", test_name, alternative)

    shape = responders.shape[:3]
    p_values = np.ones(shape)
    t_statistics = np.zeros(shape)

    valid_mask = np.any(responders > 0, axis=-1) | np.any(non_responders > 0, axis=-1)
    valid_coords = np.argwhere(valid_mask)
    _log.info("Valid voxels: %d", len(valid_coords))

    n_resp = responders.shape[-1]
    n_non_resp = non_responders.shape[-1]

    idx_i, idx_j, idx_k = valid_coords[:, 0], valid_coords[:, 1], valid_coords[:, 2]
    resp_extracted = responders[idx_i, idx_j, idx_k, :].astype(np.float32)
    non_resp_extracted = non_responders[idx_i, idx_j, idx_k, :].astype(np.float32)
    voxel_data = np.concatenate([resp_extracted, non_resp_extracted], axis=1)

    if test_type == "paired":
        t_1d, p_1d = ttest_rel(voxel_data, n_resp, alternative=alternative)
    else:
        t_1d, p_1d = ttest_ind(voxel_data, n_resp, n_non_resp, alternative=alternative)

    # Voxels whose t is not finite are degenerate (zero within-group variance:
    # ``nan`` for identical constant groups, ``+/-inf`` for perfect separation).
    # They carry no usable cluster mass, so they are dropped from valid_mask
    # rather than being allowed to poison a cluster with inf/nan.
    degenerate = ~np.isfinite(t_1d)
    if np.any(degenerate):
        _log.warning(
            "Excluding %d voxel(s) with zero within-group variance "
            "(degenerate t: %d perfectly separated, %d constant)",
            int(degenerate.sum()),
            int(np.isinf(t_1d).sum()),
            int(np.isnan(t_1d).sum()),
        )
        valid_mask = valid_mask.copy()
        valid_mask[
            idx_i[degenerate], idx_j[degenerate], idx_k[degenerate]
        ] = False
        t_1d = np.where(degenerate, 0.0, t_1d)
        p_1d = np.where(degenerate, 1.0, p_1d)

    t_statistics[idx_i, idx_j, idx_k] = t_1d
    p_values[idx_i, idx_j, idx_k] = p_1d

    # An empty mask is a real, reportable state of the data, not something a caller can use.
    # Every reduction over `values[valid_mask]` downstream — the caller's own `np.min(p_values
    # [valid_mask])` first of all — raises "zero-size array to reduction operation minimum",
    # a numpy traceback that says nothing about why. Two things produce it: no voxel is
    # non-zero in either group, and (since the degenerate exclusion above) every testable voxel
    # has zero within-group variance. Say which of the two, here, where both are known.
    if not valid_mask.any():
        raise ValueError(
            "No voxel could be tested. "
            + (
                "Every voxel with data has zero within-group variance — a constant value "
                "across every subject in a group — and was excluded as degenerate."
                if len(valid_coords) > 0
                else "No voxel is non-zero in either group; the input images or the mask are empty."
            )
            + " Check that the input images are the fields you meant and that the groups vary."
        )

    return p_values, t_statistics, valid_mask

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