stats

Custom nipype interfaces for FreeSurfer stats commands.

Classes:

• AparcStats –

  Custom nipype interface for FreeSurfer aparcstats2table command.

• AparcStatsInputSpec –

  Input specification for aparcstats2table command.

• AparcStatsOutputSpec –

  Output specification for aparcstats2table command.

• AsegStats –

  Custom nipype interface for FreeSurfer asegstats2table command.

• AsegStatsInputSpec –

  Input specification for asegstats2table command.

• AsegStatsOutputSpec –

  Output specification for asegstats2table command.

Functions:

• check_metrics –

  Check metrics from stats files.

• compare_group_metrics –

  Summarize FreeSurfer metrics for named subject groups.

• gen_group_comparison_plots –

  Generate Plotly box plots comparing metrics across groups.

• gen_metric_plots –

  Generate plots from FreeSurfer stats files.

• get_stats –

  Get aseg and aparc stats from subjects.

• summarize_outlier_subjects –

  Aggregate outlier findings by subject for quick reference.

AparcStats

Bases: FSCommand

Custom nipype interface for FreeSurfer aparcstats2table command.

Methods:

• run –

  Run aparcstats2table command.

run

run(**inputs: Any) -> dict

Run aparcstats2table command.

Source code in src/pyfsviz/stats.py

def run(self, **inputs: Any) -> dict:
    """Run aparcstats2table command."""
    return super().run(**inputs)

AparcStatsInputSpec

Bases: FSTraitedSpec

Input specification for aparcstats2table command.

AparcStatsOutputSpec

Bases: TraitedSpec

Output specification for aparcstats2table command.

AsegStats

Bases: FSCommand

Custom nipype interface for FreeSurfer asegstats2table command.

Methods:

• run –

  Run asegstats2table command.

run

run(**inputs: Any) -> dict

Run asegstats2table command.

Source code in src/pyfsviz/stats.py

def run(self, **inputs: Any) -> dict:
    """Run asegstats2table command."""
    return super().run(**inputs)

AsegStatsInputSpec

Bases: FSTraitedSpec

Input specification for asegstats2table command.

AsegStatsOutputSpec

Bases: TraitedSpec

Output specification for asegstats2table command.

check_metrics

check_metrics(stats_files: list[Path], sd_threshold: float = 3.0) -> dict

Check metrics from stats files.

Parameters:

• stats_files

  (list) –

  List of paths to stats files

• sd_threshold

  (float, default: 3.0 ) –

  Standard deviation threshold, by default 3.0

Source code in src/pyfsviz/stats.py

def check_metrics(stats_files: list[Path], sd_threshold: float = 3.0) -> dict:
    """Check metrics from stats files.

    Parameters
    ----------
    stats_files : list
        List of paths to stats files
    sd_threshold : float, optional
        Standard deviation threshold, by default 3.0
    """
    metrics = _load_metrics(stats_files)
    metric_summary: dict[str, dict[str, dict[str, Any]]] = {}

    for metric, data in metrics.items():
        # Initialize metric_summary for this metric
        metric_summary[metric] = {}

        # Get column names - skip first column (subject_id) and last few columns (typically metadata)
        region_cols = [
            col
            for col in data.columns[1:]
            if col
            not in [
                "ID",
                "Measure:volume",
                "lh.aparc.a2009s_thickness",
                "rh.aparc.a2009s_thickness",
            ]
        ]
        id_col = data.columns[0]

        for region in region_cols:
            values = data[region].dropna()
            if len(values) == 0:
                metric_summary[metric][region] = {
                    "status": "no_data",
                    "message": "No data available",
                }
            else:
                values = values.astype(float)
                mean = values.mean()
                std_val = values.std()
                upper_bound = mean + sd_threshold * std_val
                lower_bound = mean - sd_threshold * std_val
                outliers = values[(values > upper_bound) | (values < lower_bound)]

                if len(outliers) > 0:
                    outlier_percentage = (len(outliers) / len(values)) * 100
                    outlier_subjects = []
                    for outlier_val in outliers:
                        outlier_rows = data[data[region] == outlier_val]
                        for _, row in outlier_rows.iterrows():
                            subject_id = row[id_col]
                            outlier_subjects.append(
                                {
                                    "subject_id": str(subject_id),
                                    "value": float(outlier_val),
                                },
                            )

                    unique_outliers = []
                    seen = set()
                    for outlier in outlier_subjects:
                        key = (outlier["subject_id"], outlier["value"])
                        if key not in seen:
                            unique_outliers.append(outlier)
                            seen.add(key)

                    metric_summary[metric][region] = {
                        "status": "outliers_detected",
                        "message": f"Found {len(outliers)} outliers ({outlier_percentage:.1f}%) beyond {sd_threshold} SD",
                        "outlier_count": len(outliers),
                        "outlier_percentage": outlier_percentage,
                        "outlier_subjects": unique_outliers,
                        "mean": mean,
                        "std": std_val,
                        "sd_threshold": sd_threshold,
                        "upper_bound": upper_bound,
                        "lower_bound": lower_bound,
                    }
                else:
                    metric_summary[metric][region] = {
                        "status": "passed",
                        "message": f"No outliers detected (mean: {mean:.2f}, ±{sd_threshold} SD: {lower_bound:.2f} to {upper_bound:.2f})",
                        "outlier_count": 0,
                        "outlier_percentage": 0.0,
                        "outlier_subjects": [],
                        "mean": mean,
                        "std": std_val,
                        "sd_threshold": sd_threshold,
                        "upper_bound": upper_bound,
                        "lower_bound": lower_bound,
                    }
    return metric_summary

compare_group_metrics

compare_group_metrics(stats_files: list[Path], groups: dict[str, list[str]]) -> dict[str, dict[str, dict[str, Any]]]

Summarize FreeSurfer metrics for named subject groups.

Parameters:

• stats_files

  (list[Path]) –

  Paths to stats CSV files from :func:get_stats.

• groups

  (dict[str, list[str]]) –

  Mapping of group name to subject IDs, e.g. {"control": [...], "patient": [...]}.

Returns:

• dict –

  Nested summaries keyed by metric file and brain region, with per-group n, mean, and std.

Source code in src/pyfsviz/stats.py

def compare_group_metrics(
    stats_files: list[Path],
    groups: dict[str, list[str]],
) -> dict[str, dict[str, dict[str, Any]]]:
    """Summarize FreeSurfer metrics for named subject groups.

    Parameters
    ----------
    stats_files
        Paths to stats CSV files from :func:`get_stats`.
    groups
        Mapping of group name to subject IDs, e.g. ``{"control": [...], "patient": [...]}``.

    Returns
    -------
    dict
        Nested summaries keyed by metric file and brain region, with per-group
        ``n``, ``mean``, and ``std``.
    """
    if len(groups) < _MIN_GROUPS:
        msg = "At least two groups are required for comparison"
        raise ValueError(msg)

    comparison: dict[str, dict[str, dict[str, Any]]] = {}
    metrics = _load_metrics(stats_files)
    group_names = list(groups)

    for metric_name, data in metrics.items():
        comparison[metric_name] = {}
        for region in _region_columns(data):  # skip id column
            grouped_values = _group_values(data, region, groups)
            comparison[metric_name][region] = {
                group_name: {
                    "n": len(grouped_values[group_name]),
                    "mean": float(grouped_values[group_name].mean()) if len(grouped_values[group_name]) else None,
                    "std": float(grouped_values[group_name].std(ddof=1))
                    if len(grouped_values[group_name]) > 1
                    else None,
                }
                for group_name in group_names
            }

    return comparison

gen_group_comparison_plots

gen_group_comparison_plots(stats_files: list[Path], groups: dict[str, list[str]]) -> list[Figure]

Generate Plotly box plots comparing metrics across groups.

Parameters:

• stats_files

  (list[Path]) –

  Paths to stats CSV files from :func:get_stats.

• groups

  (dict[str, list[str]]) –

  Mapping of group name to subject IDs.

Returns:

• list[Figure] –

  Plotly figures with one plot per metric region. Each figure stores metric and label in layout.meta for report grouping.

Source code in src/pyfsviz/stats.py

def gen_group_comparison_plots(
    stats_files: list[Path],
    groups: dict[str, list[str]],
) -> list[go.Figure]:
    """Generate Plotly box plots comparing metrics across groups.

    Parameters
    ----------
    stats_files
        Paths to stats CSV files from :func:`get_stats`.
    groups
        Mapping of group name to subject IDs.

    Returns
    -------
    list[go.Figure]
        Plotly figures with one plot per metric region. Each figure stores
        ``metric`` and ``label`` in ``layout.meta`` for report grouping.
    """
    plots: list[go.Figure] = []
    metrics = _load_metrics(stats_files)
    group_names = list(groups)

    for metric_name, data in metrics.items():
        id_col = data.columns[0]
        subject_groups = _subject_group_map(groups)
        label = _comparison_metric_label(metric_name)

        for region in _region_columns(data):
            plot_rows = []
            for _, row in data.iterrows():
                subject_id = str(row[id_col])
                group_name = subject_groups.get(subject_id)
                value = row[region]
                if group_name is None or pd.isna(value):
                    continue
                try:
                    numeric = float(value)
                except (TypeError, ValueError):
                    continue
                plot_rows.append(
                    {
                        "group": group_name,
                        "value": numeric,
                        "subject_id": subject_id,
                    },
                )

            if not plot_rows:
                continue

            plot_data = pd.DataFrame(plot_rows)
            fig = go.Figure()
            for group_name in group_names:
                group_data = plot_data[plot_data["group"] == group_name]
                if group_data.empty:
                    continue
                fig.add_trace(
                    go.Box(
                        y=_plotly_values(group_data["value"].astype(float)),
                        name=group_name,
                        text=_plotly_values(group_data["subject_id"].astype(str)),
                        boxpoints="all",
                    ),
                )

            fig.update_layout(
                autosize=True,
                height=420,
                boxmode="group",
                title={"text": region},
                yaxis={"title": {"text": region}},
                xaxis={"title": {"text": "Group"}},
                meta={"metric": metric_name, "label": label},
            )
            plots.append(fig)

    return plots

gen_metric_plots

gen_metric_plots(stats_files: list[Path]) -> list

Generate plots from FreeSurfer stats files.

Parameters:

• stats_files

  (list[Path]) –

  List of paths to stats files

Returns:

• list –

  List of plotly figure objects

Source code in src/pyfsviz/stats.py

def gen_metric_plots(stats_files: list[Path]) -> list:
    """Generate plots from FreeSurfer stats files.

    Parameters
    ----------
    stats_files: list
        List of paths to stats files

    Returns
    -------
    list
        List of plotly figure objects
    """
    plots = []
    metrics = _load_metrics(stats_files)

    for metric, data in metrics.items():
        idx_col = data.columns[0]
        if "hemi" in data.columns:
            for c in data[1:]:
                fig = go.Figure()
                fig.add_trace(
                    go.Box(
                        y=_plotly_values(data[data["hemi"] == "lh"][c]),
                        boxpoints="suspectedoutliers",
                        marker={
                            "outliercolor": "rgb(0,0,0)",
                            "line": {"outlierwidth": 1, "outliercolor": "rgb(0,0,0)"},
                        },
                        name="lh",
                        text=_plotly_values(data[data["hemi"] == "lh"][idx_col]),
                    ),
                )
                fig.add_trace(
                    go.Box(
                        y=_plotly_values(data[data["hemi"] == "rh"][c]),
                        boxpoints="suspectedoutliers",
                        marker={
                            "outliercolor": "rgb(0,0,0)",
                            "line": {"outlierwidth": 1, "outliercolor": "rgb(0,0,0)"},
                        },
                        name="rh",
                        text=_plotly_values(data[data["hemi"] == "rh"][idx_col]),
                    ),
                )
                fig.update_layout(
                    boxmode="group",
                    yaxis={"title": {"text": c}},
                    xaxis={"title": {"text": "hemisphere"}},
                    title={"text": c},
                )
                _stamp_plot_meta(fig, metric)
                plots.append(fig)
        elif any("Left-" in c for c in data.columns):
            region_groups: dict[str, dict[str, str]] = {}
            for region in data[1:]:
                # Extract base region name (remove hemisphere prefix if present)
                if region.startswith("Left-"):
                    base_region = region[5:]  # Remove 'Left-' prefix
                    hemisphere = "Left"
                elif region.startswith("Right-"):
                    base_region = region[6:]  # Remove 'Right-' prefix
                    hemisphere = "Right"
                elif region.startswith(("lh", "rh")):
                    base_region = region[2:]  # Remove 'lh' or 'rh' prefix
                    hemisphere = "Left" if region.startswith("lh") else "Right"
                else:
                    # No hemisphere prefix, treat as bilateral
                    base_region = region
                    hemisphere = "Bilateral"

                if base_region not in region_groups:
                    region_groups[base_region] = {}
                region_groups[base_region][hemisphere] = region

            for base_region, hemispheres in region_groups.items():
                fig = go.Figure()
                if len(hemispheres) > 1 and "Bilateral" not in hemispheres:
                    # Multiple hemispheres found, create combined plot
                    combined_data = []
                    for hemisphere, region_col in hemispheres.items():
                        region_data = data[[idx_col, region_col]].copy()
                        region_data = region_data.rename(columns={region_col: "value"})
                        region_data["hemisphere"] = hemisphere
                        combined_data.append(region_data)

                    if combined_data:
                        # Concatenate data from both hemispheres
                        plot_data = pd.concat(combined_data, ignore_index=True)

                        # Create box plot comparing hemispheres
                        fig.add_trace(
                            go.Box(
                                y=_plotly_values(
                                    plot_data[plot_data["hemisphere"] == "Left"]["value"],
                                ),
                                boxpoints="suspectedoutliers",
                                text=_plotly_values(
                                    plot_data[plot_data["hemisphere"] == "Left"][idx_col],
                                ),
                                name="left",
                                marker={
                                    "outliercolor": "rgb(0,0,0)",
                                    "line": {
                                        "outlierwidth": 1,
                                        "outliercolor": "rgb(0,0,0)",
                                    },
                                },
                            ),
                        )
                        fig.add_trace(
                            go.Box(
                                y=_plotly_values(
                                    plot_data[plot_data["hemisphere"] == "Right"]["value"],
                                ),
                                boxpoints="suspectedoutliers",
                                text=_plotly_values(
                                    plot_data[plot_data["hemisphere"] == "Right"][idx_col],
                                ),
                                name="right",
                                marker={
                                    "outliercolor": "rgb(0,0,0)",
                                    "line": {
                                        "outlierwidth": 1,
                                        "outliercolor": "rgb(0,0,0)",
                                    },
                                },
                            ),
                        )
                        fig.update_layout(
                            boxmode="group",
                            yaxis={"title": {"text": base_region}},
                            xaxis={"title": {"text": "hemisphere"}},
                            title={"text": base_region},
                        )
                        _stamp_plot_meta(fig, metric)
                        plots.append(fig)
                else:
                    region_col = next(iter(hemispheres.values()))
                    fig.add_trace(
                        go.Box(
                            y=_plotly_values(data[region_col]),
                            boxpoints="suspectedoutliers",
                            text=_plotly_values(data[idx_col]),
                            name=base_region,
                            marker={
                                "outliercolor": "rgb(0,0,0)",
                                "line": {
                                    "outlierwidth": 1,
                                    "outliercolor": "rgb(0,0,0)",
                                },
                            },
                        ),
                    )
                    fig.update_layout(
                        yaxis={"title": {"text": base_region}},
                        title={"text": base_region},
                    )
                    _stamp_plot_meta(fig, metric)
                    plots.append(fig)
        else:
            for region in data[1:]:
                fig = go.Figure()
                fig.add_trace(
                    go.Box(
                        y=_plotly_values(data[region]),
                        boxpoints="suspectedoutliers",
                        text=_plotly_values(data[idx_col]),
                        name=region,
                        marker={
                            "outliercolor": "rgb(0,0,0)",
                            "line": {"outlierwidth": 1, "outliercolor": "rgb(0,0,0)"},
                        },
                    ),
                )
                fig.update_layout(
                    yaxis={"title": {"text": region}},
                    title={"text": region},
                )
                _stamp_plot_meta(fig, metric)
                plots.append(fig)

    return plots

get_stats

get_stats(subjects: list[str], output_dir: str, measures: list[str] | None = None, hemis: list[str] | None = None) -> dict[str, Path | list[Path]]

Get aseg and aparc stats from subjects.

Parameters:

• subjects

  (list) –

  List of subject IDs

• output_dir

  (str) –

• measures

  (list, default: None ) –

  List of measures to get, by default None

• hemis

  (list, default: None ) –

  List of hemispheres to get, by default None

Source code in src/pyfsviz/stats.py

def get_stats(
    subjects: list[str],
    output_dir: str,
    measures: list[str] | None = None,
    hemis: list[str] | None = None,
) -> dict[str, Path | list[Path]]:
    """Get aseg and aparc stats from subjects.

    Parameters
    ----------
    subjects : list
        List of subject IDs
    output_dir : str
    measures : list, optional
        List of measures to get, by default None
    hemis : list, optional
        List of hemispheres to get, by default None
    """
    stats: dict[str, Path | list[Path]] = {}
    stats["aseg"] = _get_aseg_stats(subjects, "aseg.csv", output_dir=output_dir)
    stats["aparc"] = _get_aparc_stats(
        subjects,
        "aparc.csv",
        output_dir=output_dir,
        measures=measures,
        hemis=hemis,
    )
    return stats

summarize_outlier_subjects

summarize_outlier_subjects(quality_summary: dict[str, dict[str, dict[str, Any]]]) -> list[dict[str, Any]]

Aggregate outlier findings by subject for quick reference.

Parameters:

• quality_summary

  (dict[str, dict[str, dict[str, Any]]]) –

  Output from :func:check_metrics.

Returns:

• list[dict] –

  Sorted list of subjects with outlier counts and findings.

Source code in src/pyfsviz/stats.py

def summarize_outlier_subjects(
    quality_summary: dict[str, dict[str, dict[str, Any]]],
) -> list[dict[str, Any]]:
    """Aggregate outlier findings by subject for quick reference.

    Parameters
    ----------
    quality_summary
        Output from :func:`check_metrics`.

    Returns
    -------
    list[dict]
        Sorted list of subjects with outlier counts and findings.
    """
    subject_findings: dict[str, list[str]] = defaultdict(list)

    for metric_name, metric_data in quality_summary.items():
        for region, result in metric_data.items():
            if result.get("status") != "outliers_detected":
                continue
            for outlier in result.get("outlier_subjects", []):
                subject_id = str(outlier["subject_id"])
                finding = f"{metric_name}/{region}: {float(outlier['value']):.2f}"
                if finding not in subject_findings[subject_id]:
                    subject_findings[subject_id].append(finding)

    return sorted(
        [
            {
                "subject_id": subject_id,
                "outlier_count": len(findings),
                "findings": findings,
            }
            for subject_id, findings in subject_findings.items()
        ],
        key=lambda item: (-item["outlier_count"], item["subject_id"]),
    )

Feedback