simnibs_analyze.meta

  1"""
  2Meta-analysis: compare intra-ROI mean e-field across spaces and ROI methods.
  3
  4Inputs : one or more  all_features_space-<X>.csv  files produced by run.py.
  5Outputs: saved in  <results_dir>/meta_analysis/
  6
  7Usage (CLI):
  8    python meta.py --results-dir /path/to/results-V3 --metric mean
  9
 10Two comparisons are implemented:
 11  1. space_comparison   — same condition, mni vs native
 12  2. roi_method_comparison — same target zone, different ROI definition
 13                             (requires explicit pairs in --roi-pairs)
 14"""
 15
 16from __future__ import annotations
 17
 18import argparse
 19from pathlib import Path
 20from typing import List, Optional
 21
 22import matplotlib.pyplot as plt
 23import pandas as pd
 24
 25from ._pipeline_io import get_features_csv_path, save_dataframe, save_figure
 26from ._logging import get_logger
 27
 28logger = get_logger(__name__)
 29
 30SPACES = ("mni", "native")
 31
 32
 33# ---------------------------------------------------------------------------
 34# Loading
 35# ---------------------------------------------------------------------------
 36
 37
 38def load_all_features(results_dir: Path, spaces: tuple = SPACES) -> pd.DataFrame:
 39    """Load and concatenate all_features CSVs found in *results_dir*/analysis/."""
 40    frames = []
 41    for space in spaces:
 42        path = get_features_csv_path(results_dir, space)
 43        if path.exists():
 44            df = pd.read_csv(path)
 45            df["space"] = space  # ensure column present even if missing
 46            frames.append(df)
 47            logger.info(f"Loaded {len(df)} rows from {path.name}")
 48        else:
 49            logger.info(f"Not found (skipped): {path.name}")
 50    if not frames:
 51        raise FileNotFoundError(
 52            f"No all_features_*.csv found in {results_dir}/analysis/"
 53        )
 54    return pd.concat(frames, ignore_index=True)
 55
 56
 57# ---------------------------------------------------------------------------
 58# Comparison 1 — space (mni vs native)
 59# ---------------------------------------------------------------------------
 60
 61
 62def compare_spaces(
 63    df: pd.DataFrame,
 64    metric: str = "mean",
 65    condition_col: str = "condition",
 66) -> pd.DataFrame:
 67    """
 68    For each condition, compare *metric* between mni and native space.
 69
 70    Returns a wide DataFrame:
 71        condition | mean_mni | std_mni | n_mni | mean_native | std_native | n_native | delta_mean
 72    """
 73    if "space" not in df.columns:
 74        raise KeyError("Column 'space' missing — load CSVs from multiple spaces.")
 75
 76    rows = []
 77    for cond, grp in df.groupby(condition_col):
 78        row: dict = {"condition": cond}
 79        for space in ("mni", "native"):
 80            sub = grp[grp["space"] == space][metric].dropna()
 81            row[f"mean_{space}"] = sub.mean()
 82            row[f"std_{space}"] = sub.std()
 83            row[f"n_{space}"] = len(sub)
 84        row["delta_mean"] = row["mean_mni"] - row["mean_native"]
 85        rows.append(row)
 86    return pd.DataFrame(rows)
 87
 88
 89def plot_space_comparison(
 90    summary: pd.DataFrame,
 91    metric: str = "mean",
 92    out_path: Optional[Path] = None,
 93) -> None:
 94    """Bar chart: mni vs native mean e-field per condition."""
 95    conditions = summary["condition"].tolist()
 96    x = range(len(conditions))
 97    width = 0.35
 98
 99    fig, ax = plt.subplots(figsize=(max(6, len(conditions) * 1.5), 5))
100    ax.bar(
101        [i - width / 2 for i in x],
102        summary["mean_mni"],
103        width,
104        yerr=summary["std_mni"],
105        label="MNI",
106        capsize=4,
107        alpha=0.8,
108    )
109    ax.bar(
110        [i + width / 2 for i in x],
111        summary["mean_native"],
112        width,
113        yerr=summary["std_native"],
114        label="Native",
115        capsize=4,
116        alpha=0.8,
117    )
118    ax.set_xticks(list(x))
119    ax.set_xticklabels(conditions, rotation=30, ha="right")
120    ax.set_ylabel(f"Intra-ROI {metric} e-field (V/m)")
121    ax.set_title("Space comparison: MNI vs Native")
122    ax.legend()
123    ax.grid(axis="y", alpha=0.3)
124    plt.tight_layout()
125    if out_path:
126        save_figure(out_path, dpi=150, bbox_inches="tight")
127    plt.close(fig)
128
129
130# ---------------------------------------------------------------------------
131# Comparison 2 — ROI method
132# ---------------------------------------------------------------------------
133
134
135def compare_roi_methods(
136    df: pd.DataFrame,
137    roi_pairs: List[tuple],
138    metric: str = "mean",
139    condition_col: str = "condition",
140) -> pd.DataFrame:
141    """
142    Compare *metric* between pairs of conditions targeting the same brain zone
143    but using different ROI definitions.
144
145    Parameters
146    ----------
147    roi_pairs : list of (cond_a, cond_b) tuples
148        e.g. [("fef_simulation", "HA-fef_simulation"),
149               ("fef_simulation", "AAL-fef_simulation")]
150
151    Returns a DataFrame with one row per pair:
152        cond_a | cond_b | mean_a | std_a | n_a | mean_b | std_b | n_b | delta_mean
153    """
154    rows = []
155    for cond_a, cond_b in roi_pairs:
156        sub_a = df[df[condition_col] == cond_a][metric].dropna()
157        sub_b = df[df[condition_col] == cond_b][metric].dropna()
158        rows.append(
159            {
160                "cond_a": cond_a,
161                "cond_b": cond_b,
162                "mean_a": sub_a.mean(),
163                "std_a": sub_a.std(),
164                "n_a": len(sub_a),
165                "mean_b": sub_b.mean(),
166                "std_b": sub_b.std(),
167                "n_b": len(sub_b),
168                "delta_mean": sub_a.mean() - sub_b.mean(),
169            }
170        )
171    return pd.DataFrame(rows)
172
173
174def plot_roi_method_comparison(
175    summary: pd.DataFrame,
176    metric: str = "mean",
177    out_path: Optional[Path] = None,
178) -> None:
179    """Bar chart: condition A vs B for each pair."""
180    n = len(summary)
181    fig, ax = plt.subplots(figsize=(max(6, n * 2), 5))
182    x = range(n)
183    width = 0.35
184
185    ax.bar(
186        [i - width / 2 for i in x],
187        summary["mean_a"],
188        width,
189        yerr=summary["std_a"],
190        label="ROI A",
191        capsize=4,
192        alpha=0.8,
193    )
194    ax.bar(
195        [i + width / 2 for i in x],
196        summary["mean_b"],
197        width,
198        yerr=summary["std_b"],
199        label="ROI B",
200        capsize=4,
201        alpha=0.8,
202    )
203    labels = [f"{r.cond_a}\nvs\n{r.cond_b}" for _, r in summary.iterrows()]
204    ax.set_xticks(list(x))
205    ax.set_xticklabels(labels, fontsize=8)
206    ax.set_ylabel(f"Intra-ROI {metric} e-field (V/m)")
207    ax.set_title("ROI method comparison")
208    ax.legend()
209    ax.grid(axis="y", alpha=0.3)
210    plt.tight_layout()
211    if out_path:
212        save_figure(out_path, dpi=150, bbox_inches="tight")
213    plt.close(fig)
214
215
216# ---------------------------------------------------------------------------
217# Entry point
218# ---------------------------------------------------------------------------
219
220
221def run(
222    results_dir: Path,
223    metric: str = "mean",
224    roi_pairs: Optional[List[tuple]] = None,
225    spaces: tuple = SPACES,
226) -> None:
227    """Run all meta-analyses and save outputs to <results_dir>/meta_analysis/."""
228    out_dir = results_dir / "meta_analysis"
229    out_dir.mkdir(parents=True, exist_ok=True)
230
231    df = load_all_features(results_dir, spaces=spaces)
232
233    # ── Space comparison ────────────────────────────────────────────────────
234    n_spaces = df["space"].nunique()
235    if n_spaces >= 2:
236        space_summary = compare_spaces(df, metric=metric)
237        save_dataframe(
238            space_summary, out_dir / f"space_comparison_{metric}.csv", index=False
239        )
240        plot_space_comparison(
241            space_summary,
242            metric=metric,
243            out_path=out_dir / f"space_comparison_{metric}.png",
244        )
245        logger.info(f"Space comparison saved → {out_dir}")
246    else:
247        logger.info(f"Space comparison skipped — only {df['space'].unique()} found.")
248
249    # ── ROI method comparison ───────────────────────────────────────────────
250    if roi_pairs:
251        roi_summary = compare_roi_methods(df, roi_pairs=roi_pairs, metric=metric)
252        save_dataframe(
253            roi_summary, out_dir / f"roi_method_comparison_{metric}.csv", index=False
254        )
255        plot_roi_method_comparison(
256            roi_summary,
257            metric=metric,
258            out_path=out_dir / f"roi_method_comparison_{metric}.png",
259        )
260        logger.info(f"ROI method comparison saved → {out_dir}")
261    else:
262        logger.info("ROI method comparison skipped — no --roi-pairs provided.")
263
264
265def _parse_args(argv=None):
266    parser = argparse.ArgumentParser(
267        description="Meta-analysis across spaces and ROI methods."
268    )
269    parser.add_argument(
270        "--results-dir",
271        type=Path,
272        required=True,
273        help="Root results directory (contains analysis/)",
274    )
275    parser.add_argument(
276        "--metric", default="mean", help="Feature column to compare (default: mean)"
277    )
278    parser.add_argument(
279        "--spaces",
280        nargs="+",
281        default=list(SPACES),
282        help="Spaces to load (default: mni native)",
283    )
284    parser.add_argument(
285        "--roi-pairs",
286        nargs="+",
287        default=[],
288        metavar="A:B",
289        help="Pairs of conditions to compare as ROI methods, e.g. "
290        "fef_simulation:HA-fef_simulation fef_simulation:AAL-fef_simulation",
291    )
292    return parser.parse_args(argv)
293
294
295if __name__ == "__main__":
296    args = _parse_args()
297    pairs = [tuple(p.split(":")) for p in args.roi_pairs]
298    run(
299        results_dir=args.results_dir,
300        metric=args.metric,
301        spaces=tuple(args.spaces),
302        roi_pairs=pairs or None,
303    )