more stats
This commit is contained in:
@@ -3324,6 +3324,345 @@ def run_cross_group_second_level_analysis(df_roi_all, file_paths_a, file_paths_b
|
||||
|
||||
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import scipy.stats as stats
|
||||
import matplotlib.pyplot as plt
|
||||
import seaborn as sns
|
||||
from statsmodels.stats.multitest import multipletests
|
||||
import logging
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
def run_roi_paired_contrast_analysis(df_roi_all, roi_pairs, condition,
|
||||
target_chroma='hbo', min_subjects=5,
|
||||
p_threshold=0.05, correction_method=None,
|
||||
roi_a_label=None, roi_b_label=None):
|
||||
"""
|
||||
Paired within-subject ROI contrast (e.g. contralateral minus ipsilateral
|
||||
motor ROI), as a companion to run_roi_second_level_analysis rather than a
|
||||
replacement for it. Where run_roi_second_level_analysis tests each ROI's
|
||||
theta against zero independently (still contaminated by systemic/global
|
||||
physiology shared across the whole head), this function computes, per
|
||||
subject, (ROI_A theta - ROI_B theta) for a single condition and tests
|
||||
THAT difference against zero. Any systemic component that's roughly equal
|
||||
in both ROIs cancels out in the subtraction itself, rather than being
|
||||
inferred afterwards by comparing two separate p-values.
|
||||
|
||||
This is the more powerful, more directly interpretable test whenever you
|
||||
already have a specific hypothesis about which two ROIs should differ
|
||||
(e.g. laterality) — use run_roi_second_level_analysis for open-ended
|
||||
per-ROI screening, and this function for a pre-specified paired
|
||||
comparison you want to report as a single confirmatory statistic.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
df_roi_all : pd.DataFrame
|
||||
Combined individual-level ROI results across subjects.
|
||||
Must include: ['ROI', 'Condition', 'Chroma', 'theta', 'ID']
|
||||
roi_pairs : tuple(str, str) or list of tuple(str, str)
|
||||
One (roi_a, roi_b) pair, or several. Each pair is tested
|
||||
independently as (roi_a - roi_b). Passing several pairs lets you
|
||||
e.g. test left-hand-tap laterality and right-hand-tap laterality
|
||||
(different `condition` values) in one call/figure.
|
||||
condition : str or list of str
|
||||
The 'Condition' value to filter to for the paired test. If
|
||||
`roi_pairs` has multiple pairs and you want a different condition
|
||||
per pair, pass a list of the same length as `roi_pairs`; otherwise
|
||||
a single value is used for every pair.
|
||||
target_chroma : str, default 'hbo'
|
||||
Chromophore to test. HbO and HbR should never be tested together.
|
||||
min_subjects : int, default 5
|
||||
Minimum number of subjects with BOTH ROI_A and ROI_B present (after
|
||||
dropping NaNs) required to run the test. Below this, the pair is
|
||||
skipped with a warning rather than silently reported.
|
||||
p_threshold : float, default 0.05
|
||||
Significance threshold applied to the (optionally corrected) p-value.
|
||||
correction_method : str or None, default None
|
||||
Multiple comparisons correction across the pairs tested in this call
|
||||
(statsmodels.stats.multitest.multipletests method name, e.g.
|
||||
'fdr_bh'). Left off by default since a single pre-specified paired
|
||||
contrast typically doesn't need correction — turn it on if you're
|
||||
testing several pairs in the same call and want to control for that.
|
||||
roi_a_label, roi_b_label : str or list of str, optional
|
||||
Display labels for each pair's ROI_A/ROI_B (defaults to the raw ROI
|
||||
names). If testing multiple pairs, pass lists matching `roi_pairs`.
|
||||
|
||||
Returns
|
||||
-------
|
||||
pd.DataFrame with one row per tested pair:
|
||||
['roi_a', 'roi_b', 'condition', 't_val', 'p_val', 'p_corrected',
|
||||
'significant', 'mean_diff', 'n_subjects']
|
||||
"""
|
||||
|
||||
required_cols = ['ROI', 'Condition', 'Chroma', 'theta', 'ID']
|
||||
if not all(col in df_roi_all.columns for col in required_cols):
|
||||
raise ValueError(f"Input ROI DataFrame must include: {required_cols}")
|
||||
|
||||
# Normalize inputs to lists so single-pair and multi-pair calls share code.
|
||||
if isinstance(roi_pairs, tuple):
|
||||
roi_pairs = [roi_pairs]
|
||||
n_pairs = len(roi_pairs)
|
||||
|
||||
if isinstance(condition, str):
|
||||
conditions = [condition] * n_pairs
|
||||
else:
|
||||
if len(condition) != n_pairs:
|
||||
raise ValueError("If passing a list of conditions, it must match len(roi_pairs).")
|
||||
conditions = list(condition)
|
||||
|
||||
def _expand_labels(labels, default_from):
|
||||
if labels is None:
|
||||
return [None] * n_pairs
|
||||
if isinstance(labels, str):
|
||||
return [labels] * n_pairs
|
||||
if len(labels) != n_pairs:
|
||||
raise ValueError("Label list length must match len(roi_pairs).")
|
||||
return list(labels)
|
||||
|
||||
roi_a_labels = _expand_labels(roi_a_label, roi_pairs)
|
||||
roi_b_labels = _expand_labels(roi_b_label, roi_pairs)
|
||||
|
||||
df_chroma = df_roi_all[df_roi_all['Chroma'] == target_chroma].copy()
|
||||
df_chroma = df_chroma.dropna(subset=['theta'])
|
||||
|
||||
results = []
|
||||
diff_data_for_plot = [] # keep per-subject diffs around for plotting
|
||||
|
||||
for (roi_a, roi_b), cond, lbl_a, lbl_b in zip(roi_pairs, conditions, roi_a_labels, roi_b_labels):
|
||||
df_cond = df_chroma[df_chroma['Condition'] == cond]
|
||||
|
||||
a_vals = df_cond[df_cond['ROI'] == roi_a].groupby('ID', as_index=False)['theta'].mean()
|
||||
b_vals = df_cond[df_cond['ROI'] == roi_b].groupby('ID', as_index=False)['theta'].mean()
|
||||
|
||||
# Inner join on ID: only subjects with BOTH ROIs present for this
|
||||
# condition contribute to the paired test.
|
||||
merged = a_vals.merge(b_vals, on='ID', suffixes=('_a', '_b'))
|
||||
merged['diff'] = merged['theta_a'] - merged['theta_b']
|
||||
|
||||
n_subs = merged['ID'].nunique()
|
||||
if n_subs < min_subjects:
|
||||
logger.warning(
|
||||
f"Skipping pair ({roi_a} - {roi_b}) for condition '{cond}' — "
|
||||
f"only {n_subs} subject(s) have both ROIs, need at least {min_subjects}."
|
||||
)
|
||||
continue
|
||||
|
||||
Y = merged['diff'].values
|
||||
t_val, p_val = stats.ttest_1samp(Y, 0)
|
||||
mean_diff = np.mean(Y)
|
||||
|
||||
results.append({
|
||||
'roi_a': roi_a,
|
||||
'roi_b': roi_b,
|
||||
'label_a': lbl_a or roi_a,
|
||||
'label_b': lbl_b or roi_b,
|
||||
'condition': cond,
|
||||
't_val': t_val,
|
||||
'p_val': p_val,
|
||||
'mean_diff': mean_diff,
|
||||
'n_subjects': n_subs,
|
||||
})
|
||||
diff_data_for_plot.append(merged.assign(pair=f"{lbl_a or roi_a} - {lbl_b or roi_b}\n({cond})"))
|
||||
|
||||
if not results:
|
||||
print("\n[ERROR] No ROI pairs met the minimum subject threshold.\n")
|
||||
return pd.DataFrame()
|
||||
|
||||
df_group = pd.DataFrame(results)
|
||||
|
||||
if correction_method is not None:
|
||||
reject, p_corrected, _, _ = multipletests(df_group['p_val'].values, method=correction_method)
|
||||
df_group['p_corrected'] = p_corrected
|
||||
df_group['significant'] = reject
|
||||
else:
|
||||
df_group['p_corrected'] = df_group['p_val']
|
||||
df_group['significant'] = df_group['p_val'] <= p_threshold
|
||||
|
||||
# --- Print report ---
|
||||
print("\n" + "=" * 70)
|
||||
print(f" PAIRED ROI CONTRAST RESULTS ({target_chroma.upper()})")
|
||||
print("=" * 70)
|
||||
df_print = df_group.copy()
|
||||
df_print['mean_diff'] = df_print['mean_diff'].apply(lambda x: f"{x:.4f}")
|
||||
df_print['t_val'] = df_print['t_val'].apply(lambda x: f"{x:.3f}")
|
||||
df_print['p_val'] = df_print['p_val'].apply(lambda x: f"{x:.4f}")
|
||||
df_print['p_corrected'] = df_print['p_corrected'].apply(lambda x: f"{x:.4f}")
|
||||
print(df_print[['label_a', 'label_b', 'condition', 'mean_diff', 't_val',
|
||||
'p_val', 'p_corrected', 'significant', 'n_subjects']].to_string(index=False))
|
||||
print("=" * 70 + "\n")
|
||||
|
||||
# --- Plot: one bar per pair, individual subject differences overlaid ---
|
||||
sns.set_theme(style="whitegrid")
|
||||
fig, ax = plt.subplots(figsize=(max(6, 2.2 * len(results)), 6))
|
||||
|
||||
plot_df = pd.concat(diff_data_for_plot, ignore_index=True)
|
||||
|
||||
sns.barplot(
|
||||
data=plot_df, x='pair', y='diff', ax=ax,
|
||||
errorbar=('ci', 95), capsize=0.1,
|
||||
color='lightgray', edgecolor='black', linewidth=1.5, zorder=1
|
||||
)
|
||||
sns.swarmplot(
|
||||
data=plot_df, x='pair', y='diff', ax=ax,
|
||||
color='darkblue', size=8, alpha=0.7, zorder=2
|
||||
)
|
||||
|
||||
ax.axhline(0, color='black', linewidth=1, linestyle='--')
|
||||
|
||||
global_max = plot_df['diff'].max()
|
||||
for i, row in df_group.iterrows():
|
||||
pair_label = f"{row['label_a']} - {row['label_b']}\n({row['condition']})"
|
||||
pair_points = plot_df[plot_df['pair'] == pair_label]['diff']
|
||||
max_y = pair_points.max() if len(pair_points) > 0 else 0
|
||||
text_y = max_y + (abs(global_max) * 0.08 if global_max else 0.1)
|
||||
|
||||
p_val_corr = row['p_corrected']
|
||||
if p_val_corr < 0.001:
|
||||
sig_symbol = "***"
|
||||
elif p_val_corr < 0.01:
|
||||
sig_symbol = "**"
|
||||
elif p_val_corr < p_threshold:
|
||||
sig_symbol = "*"
|
||||
else:
|
||||
sig_symbol = "n.s."
|
||||
|
||||
ax.text(
|
||||
i, text_y, f"{sig_symbol}\np = {p_val_corr:.3f}",
|
||||
ha='center', va='bottom', fontsize=11, fontweight='bold',
|
||||
color='red' if p_val_corr < p_threshold else 'gray'
|
||||
)
|
||||
|
||||
ax.set_ylabel(r'Paired ROI Difference ($\Delta$ HbO, A $-$ B)', fontsize=12)
|
||||
ax.set_xlabel('')
|
||||
correction_lbl = f"({correction_method} corrected)" if correction_method else "(uncorrected — single pre-specified contrast)"
|
||||
ax.set_title(
|
||||
f"Paired ROI Contrast ({target_chroma.upper()})\n"
|
||||
f"Significance threshold: p < {p_threshold} {correction_lbl}",
|
||||
fontsize=13, fontweight='bold', pad=15
|
||||
)
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
return df_group
|
||||
|
||||
|
||||
|
||||
def aggregate_channel_contrasts_to_roi(df_contrasts, roi_json_path, weighted=True):
|
||||
"""
|
||||
Combine already-computed per-channel CONTRAST results (e.g. your
|
||||
'2.0_vs_3.0' rows from contrasts.csv / contrast_results) into
|
||||
per-subject, per-ROI values — so a joint-fit task contrast can be tested
|
||||
at the ROI level using the same one-sample machinery as
|
||||
run_roi_second_level_analysis / run_roi_paired_contrast_analysis.
|
||||
|
||||
This exists because mne_nirs.statistics.RegressionResults has a built-in
|
||||
.to_dataframe_region_of_interest() that does inverse-variance-weighted
|
||||
channel combination, but the ContrastResults object returned by
|
||||
glm_est.compute_contrast() does NOT have that method. This function
|
||||
replicates the same weighting logic (weight each channel by the inverse
|
||||
of its GLM fit's variance) manually, on the already-exported contrast
|
||||
dataframe, rather than requiring you to go back and recompute anything
|
||||
from raw GLM objects.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
df_contrasts : pd.DataFrame
|
||||
Combined per-channel contrast results across subjects/contrasts
|
||||
(i.e. your contrasts.csv format). Must include:
|
||||
['ch_name', 'effect', 'stat', 'Chroma', 'contrast_name', 'ID']
|
||||
`stat` must be the t-statistic (ContrastType == 't'), since standard
|
||||
error is recovered as effect / stat.
|
||||
roi_json_path : str
|
||||
Path to the same regions.json used elsewhere in the pipeline, with
|
||||
the structure: {"regions_of_interest": [{"name": ..., "channels": [...]}]}
|
||||
`channels` entries should be bare source-detector names (e.g. "S1_D1"),
|
||||
matching the convention already used for the GLM-level ROI loading.
|
||||
weighted : bool, default True
|
||||
If True, combine channels within an ROI using inverse-variance
|
||||
weighting (weight = 1 / se^2), matching MNE-NIRS's own default
|
||||
behavior for to_dataframe_region_of_interest. If False, channels are
|
||||
weighted equally (a plain mean).
|
||||
|
||||
Returns
|
||||
-------
|
||||
pd.DataFrame with columns ['ROI', 'Condition', 'Chroma', 'theta', 'ID'],
|
||||
directly usable as `df_roi_all` in run_roi_second_level_analysis or
|
||||
run_roi_paired_contrast_analysis. 'Condition' holds the contrast name
|
||||
(e.g. '2.0_vs_3.0'), and 'theta' holds the ROI-combined contrast effect.
|
||||
"""
|
||||
|
||||
required_cols = ['ch_name', 'effect', 'stat', 'Chroma', 'contrast_name', 'ID']
|
||||
print(df_contrasts.columns)
|
||||
if not all(col in df_contrasts.columns for col in required_cols):
|
||||
raise ValueError(f"Input contrast DataFrame must include: {required_cols}")
|
||||
|
||||
# --- Load ROI definitions and build a channel-base -> ROI lookup ---
|
||||
# Channel base names (e.g. "S1_D1") map to both hbo/hbr rows via the
|
||||
# ch_name column ("S1_D1 hbo" / "S1_D1 hbr"), so we key on the base name.
|
||||
with open(roi_json_path, 'r') as f:
|
||||
roi_data = json.load(f)
|
||||
|
||||
ch_base_to_roi = {}
|
||||
for region in roi_data.get("regions_of_interest", []):
|
||||
roi_name = region["name"]
|
||||
for ch_base in region["channels"]:
|
||||
if ch_base in ch_base_to_roi:
|
||||
logger.warning(
|
||||
f"Channel '{ch_base}' assigned to multiple ROIs "
|
||||
f"('{ch_base_to_roi[ch_base]}' and '{roi_name}') — "
|
||||
f"using '{roi_name}' (last one wins)."
|
||||
)
|
||||
ch_base_to_roi[ch_base] = roi_name
|
||||
|
||||
df = df_contrasts.copy()
|
||||
df['ch_base'] = df['ch_name'].str.split().str[0] # "S1_D1 hbo" -> "S1_D1"
|
||||
df['ROI'] = df['ch_base'].map(ch_base_to_roi)
|
||||
|
||||
n_unassigned = df['ROI'].isna().sum()
|
||||
if n_unassigned:
|
||||
logger.warning(
|
||||
f"{n_unassigned} channel-rows did not match any ROI in "
|
||||
f"'{roi_json_path}' and will be excluded."
|
||||
)
|
||||
df = df.dropna(subset=['ROI'])
|
||||
|
||||
# Recover standard error from the t-statistic: t = effect / se -> se = effect / t
|
||||
with np.errstate(divide='ignore', invalid='ignore'):
|
||||
df['se'] = df['effect'] / df['stat']
|
||||
# A zero or near-zero t-stat gives an undefined/huge se; drop those rows
|
||||
# from the weighting rather than let them explode the ROI average.
|
||||
bad_se = ~np.isfinite(df['se']) | (df['se'] == 0)
|
||||
if bad_se.any():
|
||||
logger.warning(f"Dropping {bad_se.sum()} channel-rows with non-finite "
|
||||
f"standard error (t-stat ~ 0) from ROI aggregation.")
|
||||
df = df[~bad_se]
|
||||
|
||||
if weighted:
|
||||
df['weight'] = 1.0 / (df['se'] ** 2)
|
||||
else:
|
||||
df['weight'] = 1.0
|
||||
|
||||
group_cols = ['ROI', 'contrast_name', 'Chroma', 'ID']
|
||||
|
||||
def _weighted_mean(g):
|
||||
return np.average(g['effect'], weights=g['weight'])
|
||||
|
||||
roi_theta = (
|
||||
df.groupby(group_cols, group_keys=False)
|
||||
.apply(lambda g: pd.Series({'theta': _weighted_mean(g)}))
|
||||
.reset_index()
|
||||
)
|
||||
|
||||
roi_theta = roi_theta.rename(columns={'contrast_name': 'Condition'})
|
||||
return roi_theta[['ROI', 'Condition', 'Chroma', 'theta', 'ID']]
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user