From ffa14693b351f3275bdcd73ba48d8aff50db648e Mon Sep 17 00:00:00 2001 From: Tyler Date: Wed, 15 Jul 2026 16:24:13 -0700 Subject: [PATCH] start to figure out stats --- flares.py | 841 +++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 617 insertions(+), 224 deletions(-) diff --git a/flares.py b/flares.py index 5f708ec..520d81a 100644 --- a/flares.py +++ b/flares.py @@ -24,6 +24,7 @@ from concurrent.futures import ProcessPoolExecutor, as_completed from queue import Empty import time import multiprocessing as mp +import itertools # External library imports import matplotlib.pyplot as plt @@ -1488,38 +1489,62 @@ def make_design_matrix(raw_haemo, short_chans): except: pass - # 2) Create design matrix + design_matrix = make_first_level_design_matrix( + raw=raw_haemo, + stim_dur=STIM_DUR, + hrf_model=HRF_MODEL, + drift_model=DRIFT_MODEL, + high_pass=HIGH_PASS, + drift_order=DRIFT_ORDER, + fir_delays=FIR_DELAYS, + min_onset=MIN_ONSET, + oversampling=OVERSAMPLING + ) + + # 3) Average and Append Short Channels if SHORT_CHANNEL_REGRESSION and not FOLDING_BYP: - design_matrix = make_first_level_design_matrix( - raw=raw_haemo, - stim_dur=STIM_DUR, - hrf_model=HRF_MODEL, - drift_model=DRIFT_MODEL, - high_pass=HIGH_PASS, - drift_order=DRIFT_ORDER, - fir_delays=FIR_DELAYS, - add_regs=short_chans.get_data().T, - add_reg_names=short_chans.ch_names, - min_onset=MIN_ONSET, - oversampling=OVERSAMPLING - ) - else: - design_matrix = make_first_level_design_matrix( - raw=raw_haemo, - stim_dur=STIM_DUR, - hrf_model=HRF_MODEL, - drift_model=DRIFT_MODEL, - high_pass=HIGH_PASS, - drift_order=DRIFT_ORDER, - fir_delays=FIR_DELAYS, - min_onset=MIN_ONSET, - oversampling=OVERSAMPLING - ) + if short_chans is not None and len(short_chans.ch_names) > 0: + ch_types = short_chans.get_channel_types() + + # Scenario A: Short channels are already converted to Hemoglobin (hbo/hbr) + if "hbo" in ch_types or "hbr" in ch_types: + hbo_chs = [ch for ch, t in zip(short_chans.ch_names, ch_types) if t == "hbo"] + hbr_chs = [ch for ch, t in zip(short_chans.ch_names, ch_types) if t == "hbr"] + + if hbo_chs: + hbo_data = short_chans.copy().pick(hbo_chs).get_data() + design_matrix["ShortHbO"] = np.mean(hbo_data, axis=0) + if hbr_chs: + hbr_data = short_chans.copy().pick(hbr_chs).get_data() + design_matrix["ShortHbR"] = np.mean(hbr_data, axis=0) + print(f"Successfully added averaged ShortHbO ({len(hbo_chs)} chs) and ShortHbR ({len(hbr_chs)} chs) to the matrix.") + + # Scenario B: Short channels are raw wavelengths (760nm, 850nm, etc.) + else: + wavelength_groups = {} + for ch_name in short_chans.ch_names: + # Look for the wavelength number (digits) at the end of the channel name + match = re.search(r'(\d+)$', ch_name) + if match: + wl = match.group(1) + wavelength_groups.setdefault(wl, []).append(ch_name) + + if wavelength_groups: + for wl, chs in wavelength_groups.items(): + wl_data = short_chans.copy().pick(chs).get_data() + col_name = f"Short_{wl}" + design_matrix[col_name] = np.mean(wl_data, axis=0) + print(f"Successfully added averaged short channels by wavelength: {list(wavelength_groups.keys())}") + else: + # Emergency fallback: if names have no digits, average all of them together + design_matrix["Short_Avg"] = np.mean(short_chans.get_data(), axis=0) + print("Could not detect wavelengths. Averaged all short channels into 'Short_Avg'.") + else: + print("Warning: SHORT_CHANNEL_REGRESSION is True, but no short channels were found.") print(design_matrix.head()) print(design_matrix.columns) - fig, ax1 = plt.subplots(figsize=(10, 6), constrained_layout=True) _ = plot_design_matrix(design_matrix, axes=ax1) @@ -1933,6 +1958,8 @@ def fold_channels(raw: BaseRaw, p_name: str, progress_queue=None) -> dict[str, l return channel_results + + def individual_significance(raw_haemo, glm_est): fig_individual_significances = [] # List to store figures @@ -2090,170 +2117,170 @@ def individual_significance(raw_haemo, glm_est): return fig_individual_significances # TODO: Hardcoded -def group_significance( - raw_haemo, - all_cha: pd.DataFrame, - condition: str, - correction: str = "fdr_bh" -) -> plt.Figure: - """ - Compute group-level significance using weighted Stouffer's method and plot results. +# def group_significance( +# raw_haemo, +# all_cha: pd.DataFrame, +# condition: str, +# correction: str = "fdr_bh" +# ) -> plt.Figure: +# """ +# Compute group-level significance using weighted Stouffer's method and plot results. - Args: - raw_haemo: Raw haemoglobin MNE object (used for optode positions) - all_cha: DataFrame with columns including 'ID', 'Condition', 'p_value', 'theta', 'df', 'ch_name', 'Chroma' - condition: condition prefix, e.g., 'Activity' - correction: p-value correction method ('fdr_bh' or 'bonferroni') +# Args: +# raw_haemo: Raw haemoglobin MNE object (used for optode positions) +# all_cha: DataFrame with columns including 'ID', 'Condition', 'p_value', 'theta', 'df', 'ch_name', 'Chroma' +# condition: condition prefix, e.g., 'Activity' +# correction: p-value correction method ('fdr_bh' or 'bonferroni') - Returns: - Matplotlib Figure with group-level theta values and significance. - """ +# Returns: +# Matplotlib Figure with group-level theta values and significance. +# """ - assert "ID" in all_cha.columns, "'ID' column missing in input data" - assert len(raw_haemo) >= 1, "At least one raw haemoglobin object is required" +# assert "ID" in all_cha.columns, "'ID' column missing in input data" +# assert len(raw_haemo) >= 1, "At least one raw haemoglobin object is required" - condition_prefix = f"{condition}_delay" +# condition_prefix = f"{condition}_delay" - # Filter relevant data - ch_summary = all_cha.query( - "Condition.str.startswith(@condition_prefix) and Chroma == 'hbo'", - engine='python' - ).copy() +# # Filter relevant data +# ch_summary = all_cha.query( +# "Condition.str.startswith(@condition_prefix) and Chroma == 'hbo'", +# engine='python' +# ).copy() - logger.info("=== ch_summary head ===") - logger.info(ch_summary.head()) +# logger.info("=== ch_summary head ===") +# logger.info(ch_summary.head()) - logger.info("\nSummary stats:") - logger.info(f"Total rows: {len(ch_summary)}") - logger.info(f"Unique subjects: {ch_summary['ID'].nunique() if 'ID' in ch_summary.columns else 'ID column missing'}") - logger.info(f"Unique conditions: {ch_summary['Condition'].unique()}") - logger.info(f"Unique channels (Source-Detector pairs): {ch_summary.groupby(['Source', 'Detector']).ngroups}") +# logger.info("\nSummary stats:") +# logger.info(f"Total rows: {len(ch_summary)}") +# logger.info(f"Unique subjects: {ch_summary['ID'].nunique() if 'ID' in ch_summary.columns else 'ID column missing'}") +# logger.info(f"Unique conditions: {ch_summary['Condition'].unique()}") +# logger.info(f"Unique channels (Source-Detector pairs): {ch_summary.groupby(['Source', 'Detector']).ngroups}") - logger.info("\nSample p_values:") - logger.info(ch_summary['p_value'].describe()) +# logger.info("\nSample p_values:") +# logger.info(ch_summary['p_value'].describe()) - if ch_summary.empty: - raise ValueError(f"No data found for condition prefix: {condition_prefix}") +# if ch_summary.empty: +# raise ValueError(f"No data found for condition prefix: {condition_prefix}") - # --- For debugging - logger.info(f"Total rows after filtering for condition '{condition_prefix}': {len(ch_summary)}") - logger.info(f"Unique channels: {ch_summary['ch_name'].nunique()}") - logger.info(f"Participants: {ch_summary['ID'].nunique()}") +# # --- For debugging +# logger.info(f"Total rows after filtering for condition '{condition_prefix}': {len(ch_summary)}") +# logger.info(f"Unique channels: {ch_summary['ch_name'].nunique()}") +# logger.info(f"Participants: {ch_summary['ID'].nunique()}") - # Step 1: Select the peak regressor (~6s after stimulus onset) - peak_regressor = f"{condition}_delay_6" - peak_data = ch_summary[ch_summary["Condition"] == peak_regressor].copy() +# # Step 1: Select the peak regressor (~6s after stimulus onset) +# peak_regressor = f"{condition}_delay_6" +# peak_data = ch_summary[ch_summary["Condition"] == peak_regressor].copy() - logger.info(f"\n=== Logging all values for {peak_regressor} ===") - for row in peak_data.itertuples(index=False): - logger.info( - f"Subject: {row.ID}, " - f"Channel: {row.ch_name}, " - f"Source: {row.Source}, Detector: {row.Detector}, " - f"theta: {row.theta:.4f}, " - f"p_value: {row.p_value:.6f}, " - f"df: {row.df}" - ) +# logger.info(f"\n=== Logging all values for {peak_regressor} ===") +# for row in peak_data.itertuples(index=False): +# logger.info( +# f"Subject: {row.ID}, " +# f"Channel: {row.ch_name}, " +# f"Source: {row.Source}, Detector: {row.Detector}, " +# f"theta: {row.theta:.4f}, " +# f"p_value: {row.p_value:.6f}, " +# f"df: {row.df}" +# ) - if peak_data.empty: - raise ValueError(f"No data found for peak regressor: {peak_regressor}") +# if peak_data.empty: +# raise ValueError(f"No data found for peak regressor: {peak_regressor}") - # Step 2: Combine per-channel stats across subjects - group_results = [] +# # Step 2: Combine per-channel stats across subjects +# group_results = [] - for (src, det), group in peak_data.groupby(["Source", "Detector"]): - pvals = group["p_value"].values - thetas = group["theta"].values - dfs = group["df"].values +# for (src, det), group in peak_data.groupby(["Source", "Detector"]): +# pvals = group["p_value"].values +# thetas = group["theta"].values +# dfs = group["df"].values - # Weighted Stouffer's method - weights = np.sqrt(dfs) - z_scores = norm.isf(pvals) - combined_z = np.sum(weights * z_scores) / np.sqrt(np.sum(weights**2)) - combined_p = norm.sf(combined_z) +# # Weighted Stouffer's method +# weights = np.sqrt(dfs) +# z_scores = norm.isf(pvals) +# combined_z = np.sum(weights * z_scores) / np.sqrt(np.sum(weights**2)) +# combined_p = norm.sf(combined_z) - theta_avg = np.average(thetas, weights=weights) +# theta_avg = np.average(thetas, weights=weights) - group_results.append({ - "Source": src, - "Detector": det, - "theta_avg": theta_avg, - "combined_p": combined_p - }) +# group_results.append({ +# "Source": src, +# "Detector": det, +# "theta_avg": theta_avg, +# "combined_p": combined_p +# }) - # Step 3: Create combined_df - combined_df = pd.DataFrame(group_results) +# # Step 3: Create combined_df +# combined_df = pd.DataFrame(group_results) - # Step 4: Multiple comparisons correction - _, pvals_corr, _, significant = multipletests( - combined_df["combined_p"], alpha=0.05, method=correction - ) +# # Step 4: Multiple comparisons correction +# _, pvals_corr, _, significant = multipletests( +# combined_df["combined_p"], alpha=0.05, method=correction +# ) - combined_df["pval_corr"] = pvals_corr - combined_df["significant"] = significant +# combined_df["pval_corr"] = pvals_corr +# combined_df["significant"] = significant - logger.info(f"Used peak regressor: {peak_regressor}") - logger.info(f"Channels tested: {len(combined_df)}") - logger.info(f"Significant channels after correction: {combined_df['significant'].sum()}") - # Get optode positions from the first raw file - raw = raw_haemo - src_pos, det_pos = {}, {} - for ch in raw.info["chs"]: - ch_name = ch["ch_name"] - if not ch_name or not ch["loc"].any(): - continue - parts = ch_name.split()[0] - src_str, det_str = parts.split("_") - src_num = int(src_str[1:]) - det_num = int(det_str[1:]) - src_pos[src_num] = ch["loc"][3:5] - det_pos[det_num] = ch["loc"][6:8] +# logger.info(f"Used peak regressor: {peak_regressor}") +# logger.info(f"Channels tested: {len(combined_df)}") +# logger.info(f"Significant channels after correction: {combined_df['significant'].sum()}") +# # Get optode positions from the first raw file +# raw = raw_haemo +# src_pos, det_pos = {}, {} +# for ch in raw.info["chs"]: +# ch_name = ch["ch_name"] +# if not ch_name or not ch["loc"].any(): +# continue +# parts = ch_name.split()[0] +# src_str, det_str = parts.split("_") +# src_num = int(src_str[1:]) +# det_num = int(det_str[1:]) +# src_pos[src_num] = ch["loc"][3:5] +# det_pos[det_num] = ch["loc"][6:8] - # Plotting parameters - ABS_SIGNIFICANCE_THETA_VALUE = 1 - P_THRESHOLD = 0.05 - cmap = plt.get_cmap("seismic") - norm = mcolors.Normalize(vmin=-ABS_SIGNIFICANCE_THETA_VALUE, vmax=ABS_SIGNIFICANCE_THETA_VALUE) +# # Plotting parameters +# ABS_SIGNIFICANCE_THETA_VALUE = 1 +# P_THRESHOLD = 0.05 +# cmap = plt.get_cmap("seismic") +# norm = mcolors.Normalize(vmin=-ABS_SIGNIFICANCE_THETA_VALUE, vmax=ABS_SIGNIFICANCE_THETA_VALUE) - fig, ax = plt.subplots(figsize=(8, 6)) +# fig, ax = plt.subplots(figsize=(8, 6)) - # Plot optodes - for pos in src_pos.values(): - ax.scatter(*pos, s=120, c="k", marker="o", edgecolors="white", linewidths=1, zorder=3) - for pos in det_pos.values(): - ax.scatter(*pos, s=120, c="k", marker="s", edgecolors="white", linewidths=1, zorder=3) +# # Plot optodes +# for pos in src_pos.values(): +# ax.scatter(*pos, s=120, c="k", marker="o", edgecolors="white", linewidths=1, zorder=3) +# for pos in det_pos.values(): +# ax.scatter(*pos, s=120, c="k", marker="s", edgecolors="white", linewidths=1, zorder=3) - # Plot connections colored by average theta, solid if significant - for row in combined_df.itertuples(): - src, det = int(row.Source), int(row.Detector) - tval, pval = row.theta_avg, row.pval_corr - if src in src_pos and det in det_pos: - x = [src_pos[src][0], det_pos[det][0]] - y = [src_pos[src][1], det_pos[det][1]] - linestyle = "-" if pval <= P_THRESHOLD else "--" - ax.plot(x, y, linestyle=linestyle, color=cmap(norm(tval)), linewidth=4, alpha=0.9, zorder=2) +# # Plot connections colored by average theta, solid if significant +# for row in combined_df.itertuples(): +# src, det = int(row.Source), int(row.Detector) +# tval, pval = row.theta_avg, row.pval_corr +# if src in src_pos and det in det_pos: +# x = [src_pos[src][0], det_pos[det][0]] +# y = [src_pos[src][1], det_pos[det][1]] +# linestyle = "-" if pval <= P_THRESHOLD else "--" +# ax.plot(x, y, linestyle=linestyle, color=cmap(norm(tval)), linewidth=4, alpha=0.9, zorder=2) - # Colorbar - sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) - sm.set_array([]) - cbar = plt.colorbar(sm, ax=ax, shrink=0.85) - cbar.set_label(f"Average {condition_prefix.rstrip('_')} θ-value (HbO)", fontsize=11) +# # Colorbar +# sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) +# sm.set_array([]) +# cbar = plt.colorbar(sm, ax=ax, shrink=0.85) +# cbar.set_label(f"Average {condition_prefix.rstrip('_')} θ-value (HbO)", fontsize=11) - # Format axes - ax.set_aspect("equal") - ax.set_title(f"Group-level θ-values for {condition_prefix.rstrip('_')} (HbO)", fontsize=14) - ax.set_xlabel("X position (m)", fontsize=11) - ax.set_ylabel("Y position (m)", fontsize=11) - ax.grid(True, alpha=0.3) +# # Format axes +# ax.set_aspect("equal") +# ax.set_title(f"Group-level θ-values for {condition_prefix.rstrip('_')} (HbO)", fontsize=14) +# ax.set_xlabel("X position (m)", fontsize=11) +# ax.set_ylabel("Y position (m)", fontsize=11) +# ax.grid(True, alpha=0.3) - all_x = [p[0] for p in src_pos.values()] + [p[0] for p in det_pos.values()] - all_y = [p[1] for p in src_pos.values()] + [p[1] for p in det_pos.values()] - ax.set_xlim(min(all_x) - 0.01, max(all_x) + 0.01) - ax.set_ylim(min(all_y) - 0.01, max(all_y) + 0.01) +# all_x = [p[0] for p in src_pos.values()] + [p[0] for p in det_pos.values()] +# all_y = [p[1] for p in src_pos.values()] + [p[1] for p in det_pos.values()] +# ax.set_xlim(min(all_x) - 0.01, max(all_x) + 0.01) +# ax.set_ylim(min(all_y) - 0.01, max(all_y) + 0.01) - fig.tight_layout() - fig.show() +# fig.tight_layout() +# fig.show() @@ -3046,42 +3073,92 @@ def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]: return raw -def run_second_level_analysis(df_contrasts, raw, p, bounds): + + + + + + + +# Data science requires good statistics. This section is requiring work. Am working on figuring this out. Will likely add this to a new viewer. + +from scipy import stats +import matplotlib.cm as cm + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors +from scipy import stats +from statsmodels.stats.multitest import multipletests +import logging + +logger = logging.getLogger(__name__) + + +def run_second_level_analysis(df_contrasts, raw, p, bounds, min_subjects=5, + correction_method='fdr_bh', color_by='t_val'): """ - Perform second-level analysis using contrast data from multiple participants. + Perform one-sample second-level analysis using contrast data from + multiple participants (tests each channel's group mean against zero). Parameters ---------- df_contrasts : pd.DataFrame Combined contrast results from multiple participants. Must include: ['ch_name', 'effect', 'ID'] - - Returns - ------- - pd.DataFrame - Group-level t-values, p-values, and mean effect per channel. + raw : mne.io.Raw + Raw object containing sensor geometry info. + p : float + P-value significance threshold (e.g., 0.05) for solid vs. dashed lines. + Applied to the *corrected* p-value. + bounds : float + Symmetric colormap limits (e.g., t-value limits like 4.0, or beta + value limits like 1e-6). Must match `color_by`. + min_subjects : int, default 5 + Minimum number of subjects required per channel to run a t-test. + n=2 gives 1 degree of freedom and is not a reliable estimate. + correction_method : str or None, default 'fdr_bh' + Multiple comparisons correction passed to + statsmodels.stats.multitest.multipletests. None skips correction. + color_by : {'t_val', 'mean_beta'}, default 't_val' + Which statistic drives the colormap. Match `bounds` to whichever is + chosen (t-values ~3-5; beta/concentration values are much smaller, + e.g. ~1e-6). """ if not all(col in df_contrasts.columns for col in ['ch_name', 'effect', 'ID']): raise ValueError("Input DataFrame must include 'ch_name', 'effect', and 'ID' columns.") + if color_by not in ('t_val', 'mean_beta'): + raise ValueError("color_by must be 't_val' or 'mean_beta'.") + + n_before = len(df_contrasts) + df_contrasts = df_contrasts.dropna(subset=['effect']) + n_dropped = n_before - len(df_contrasts) + if n_dropped: + logger.warning(f"Dropped {n_dropped} rows with NaN 'effect' values.") + channels = df_contrasts['ch_name'].unique() group_results = [] for ch in channels: ch_data = df_contrasts[df_contrasts['ch_name'] == ch] - if ch_data['ID'].nunique() < 2: - logger.warning(f"Skipping channel {ch} — not enough subjects.") + # Collapse to one value per subject first (averages multiple + # runs/sessions per subject) to preserve the t-test's independence + # assumption. + ch_data = ch_data.groupby('ID', as_index=False)['effect'].mean() + + if ch_data['ID'].nunique() < min_subjects: + logger.warning( + f"Skipping channel {ch} — only {ch_data['ID'].nunique()} subject(s), " + f"need at least {min_subjects}." + ) continue Y = ch_data['effect'].values - design_matrix = np.ones((len(Y), 1)) # intercept-only - model = OLSModel(design_matrix) - result = model.fit(Y) - - t_val = result.t(0).item() - p_val = 2 * t.sf(np.abs(t_val), df=result.df_model) + t_val, p_val = stats.ttest_1samp(Y, 0) mean_beta = np.mean(Y) group_results.append({ @@ -3091,118 +3168,358 @@ def run_second_level_analysis(df_contrasts, raw, p, bounds): 'mean_beta': mean_beta, 'n_subjects': len(Y) }) + if not group_results: - # Create a "Warning" figure instead of a map fig, ax = plt.subplots(figsize=(8, 4)) - ax.text(0.5, 0.5, + ax.text(0.5, 0.5, f"Second-Level Analysis Aborted\n\n" f"Reason: All {len(channels)} channels skipped.\n" - f"Requirement: At least 2 subjects (IDs) per channel.\n" - f"Current Subject Count: {df_contrasts['ID'].nunique()}", + f"Requirement: At least {min_subjects} subjects (IDs) per channel.\n" + f"Current Subject Count: {df_contrasts['ID'].nunique()}", ha='center', va='center', fontsize=12, color='darkred', bbox=dict(facecolor='white', alpha=0.5, edgecolor='red')) ax.set_axis_off() plt.show() return pd.DataFrame() - + df_group = pd.DataFrame(group_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: + logger.warning("No multiple comparisons correction applied — " + "raw p-values will be used for thresholding.") + df_group['p_corrected'] = df_group['p_val'] + df_group['significant'] = df_group['p_val'] <= p + logger.info("Second-level results:\n%s", df_group) - - # Extract the cource and detector positions from raw + # Extract source/detector positions (ignoring auxiliary/stim channels) src_pos: dict[int, tuple[float, float]] = {} det_pos: dict[int, tuple[float, float]] = {} for ch in getattr(raw, "info")["chs"]: ch_name = ch['ch_name'] - if not ch_name or not ch['loc'].any(): + if not ch_name or not ch['loc'].any() or '_' not in ch_name: continue - parts = ch_name.split()[0] - src_str, det_str = parts.split('_') - src_num = int(src_str[1:]) - det_num = int(det_str[1:]) - src_pos[src_num] = ch['loc'][3:5] - det_pos[det_num] = ch['loc'][6:8] + try: + parts = ch_name.split()[0] + src_str, det_str = parts.split('_') + src_num = int(src_str[1:]) + det_num = int(det_str[1:]) + src_pos[src_num] = ch['loc'][3:5] + det_pos[det_num] = ch['loc'][6:8] + except Exception as e: + logger.debug(f"Skipping non-fNIRS channel geometry for {ch_name}: {e}") - # Set up the plot - fig, ax = plt.subplots(figsize=(8, 6)) # type: ignore + fig, ax = plt.subplots(figsize=(8, 6)) - # Plot the sources for pos in src_pos.values(): - ax.scatter(pos[0], pos[1], s=120, c='k', marker='o', edgecolors='white', linewidths=1, zorder=3) # type: ignore - - # Plot the detectors + ax.scatter(pos[0], pos[1], s=120, c='k', marker='o', edgecolors='white', linewidths=1, zorder=3) for pos in det_pos.values(): - ax.scatter(pos[0], pos[1], s=120, c='k', marker='s', edgecolors='white', linewidths=1, zorder=3) # type: ignore + ax.scatter(pos[0], pos[1], s=120, c='k', marker='s', edgecolors='white', linewidths=1, zorder=3) - # Ensure that the colors stay within the boundaries even if they are over or under the max/min values norm = mcolors.Normalize(vmin=-bounds, vmax=bounds) + cmap = plt.get_cmap('seismic') - cmap: mcolors.Colormap = plt.get_cmap('seismic') - - # Plot connections with avg t-values for _, row in df_group.iterrows(): ch = row['ch_name'] - pval = row['p_val'] tval = row['t_val'] + mean_val = row['mean_beta'] + is_sig = row['significant'] if '_' not in ch: - logger.info(f"Skipping channel with unexpected format (no underscore): {ch}") continue src_str, det_str = ch.split('_') det_parts = det_str.split() - detector_id = det_parts[0] # e.g. "D1" + detector_id = det_parts[0] hemo_type = det_parts[1].lower() if len(det_parts) > 1 else '' - logger.info(f"Parsing channel: {ch} -> src_str: {src_str}, det_str: {detector_id}, hemo_type: {hemo_type}") - if hemo_type != 'hbo': - logger.info(f"Skipping channel {ch} because hemo_type is not HbO: {hemo_type}") continue try: src = int(src_str[1:]) det = int(detector_id[1:]) - logger.info(f"Parsed src: {src}, det: {det}") - - except Exception as e: - logger.info(f"Error parsing source/detector from channel '{ch}': {e}") + except Exception: continue if src in src_pos and det in det_pos: x = [src_pos[src][0], det_pos[det][0]] y = [src_pos[src][1], det_pos[det][1]] - style = '-' if pval <= p else '--' - color = cmap(norm(tval)) - logger.info(f"Plotting {ch}: t={tval:.2f}, p={pval:.3f}, color={color}, style={style}") + style = '-' if is_sig else '--' + color_val = tval if color_by == 't_val' else mean_val + color = cmap(norm(color_val)) ax.plot(x, y, linestyle=style, color=color, linewidth=4, alpha=0.9, zorder=2) - - # Format the Colorbar sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) sm.set_array([]) - cbar = plt.colorbar(sm, ax=ax, shrink=0.85) # type: ignore - cbar.set_label(f'Average value (hbo)', fontsize=11) # type: ignore + cbar = plt.colorbar(sm, ax=ax, shrink=0.85) + cbar.set_label( + 'Group t-value (HbO)' if color_by == 't_val' else 'Group mean beta (HbO)', + fontsize=11 + ) - # Formatting the subplots ax.set_aspect('equal') - ax.set_title(f"Average values (HbO)", fontsize=14) # type: ignore - ax.set_xlabel('X position (m)', fontsize=11) # type: ignore - ax.set_ylabel('Y position (m)', fontsize=11) # type: ignore - ax.grid(True, alpha=0.3) # type: ignore + correction_label = correction_method if correction_method else 'uncorrected' + ax.set_title(f"Group-Level Activation Map (HbO), p<{p} ({correction_label})", fontsize=14) + ax.set_xlabel('X position (m)', fontsize=11) + ax.set_ylabel('Y position (m)', fontsize=11) + ax.grid(True, alpha=0.3) - # Set axis limits to be 1cm more than the optode positions all_x = [pos[0] for pos in src_pos.values()] + [pos[0] for pos in det_pos.values()] all_y = [pos[1] for pos in src_pos.values()] + [pos[1] for pos in det_pos.values()] ax.set_xlim(min(all_x)-0.01, max(all_x)+0.01) ax.set_ylim(min(all_y)-0.01, max(all_y)+0.01) fig.tight_layout() - plt.show() # type: ignore + plt.show() return df_group + + +def run_second_level_analysis_two_sample(df_contrasts, raw, p, bounds, + group_a_ids, group_b_ids, + min_subjects=5, + correction_method='fdr_bh', + color_by='t_val', + group_a_label='Group A', + group_b_label='Group B'): + """ + Compare the same contrast between two independent groups of participants + (e.g. two age ranges, patients vs controls) using Welch's t-test per + channel. This is NOT the same statistical test as run_second_level_analysis + — that function asks "is this contrast non-zero within one group?"; this + one asks "does this contrast differ between two separate groups of + people?" Group sizes and variances are not assumed equal (Welch's + correction), since two independently recruited groups will rarely match + on either. + + Parameters + ---------- + df_contrasts : pd.DataFrame + Combined contrast results from multiple participants across BOTH + groups. Must include: ['ch_name', 'effect', 'ID'] + raw : mne.io.Raw + Raw object containing sensor geometry info. + p : float + P-value significance threshold (e.g., 0.05) for solid vs. dashed lines. + Applied to the *corrected* p-value. + bounds : float + Symmetric colormap limits. Must match `color_by`: + - 't_val': typically ~3-5 + - 'mean_diff': match your effect units (e.g. ~1e-6) + group_a_ids, group_b_ids : list + The 'ID' values (e.g. file paths) belonging to each group. IDs not + present in either list are excluded from the comparison — this lets + you pass a df_contrasts containing more groups than the two you're + currently comparing. + min_subjects : int, default 5 + Minimum number of subjects required in EACH group for a channel to + be tested. A channel is skipped if either group falls short, not + just the combined total — an imbalanced 8-vs-2 split is not the same + as a balanced 5-vs-5 split even though both total 10. + correction_method : str or None, default 'fdr_bh' + Multiple comparisons correction passed to + statsmodels.stats.multitest.multipletests. None skips correction. + color_by : {'t_val', 'mean_diff'}, default 't_val' + Which statistic drives the colormap. + group_a_label, group_b_label : str + Display labels for the two groups (used in the title only). + """ + + if not all(col in df_contrasts.columns for col in ['ch_name', 'effect', 'ID']): + raise ValueError("Input DataFrame must include 'ch_name', 'effect', and 'ID' columns.") + + if color_by not in ('t_val', 'mean_diff'): + raise ValueError("color_by must be 't_val' or 'mean_diff'.") + + overlap = set(group_a_ids) & set(group_b_ids) + if overlap: + raise ValueError( + f"{len(overlap)} ID(s) appear in both groups — a participant can't be " + f"in both the A and B samples of an independent two-sample test: {overlap}" + ) + + n_before = len(df_contrasts) + df_contrasts = df_contrasts.dropna(subset=['effect']) + n_dropped = n_before - len(df_contrasts) + if n_dropped: + logger.warning(f"Dropped {n_dropped} rows with NaN 'effect' values.") + + # Keep only rows belonging to one of the two groups being compared. + keep_ids = set(group_a_ids) | set(group_b_ids) + df_contrasts = df_contrasts[df_contrasts['ID'].isin(keep_ids)] + + channels = df_contrasts['ch_name'].unique() + group_results = [] + + for ch in channels: + ch_data = df_contrasts[df_contrasts['ch_name'] == ch] + + # Collapse to one value per subject first (averages multiple + # runs/sessions per subject). + ch_data = ch_data.groupby('ID', as_index=False)['effect'].mean() + + Y_a = ch_data[ch_data['ID'].isin(group_a_ids)]['effect'].values + Y_b = ch_data[ch_data['ID'].isin(group_b_ids)]['effect'].values + + if len(Y_a) < min_subjects or len(Y_b) < min_subjects: + logger.warning( + f"Skipping channel {ch} — group sizes are {len(Y_a)} (A) and " + f"{len(Y_b)} (B), need at least {min_subjects} in EACH group." + ) + continue + + # Welch's t-test: does not assume equal variance or equal N between + # the two independent groups. + t_val, p_val = stats.ttest_ind(Y_a, Y_b, equal_var=False) + mean_a = np.mean(Y_a) + mean_b = np.mean(Y_b) + mean_diff = mean_a - mean_b + + group_results.append({ + 'ch_name': ch, + 't_val': t_val, + 'p_val': p_val, + 'mean_a': mean_a, + 'mean_b': mean_b, + 'mean_diff': mean_diff, + 'n_a': len(Y_a), + 'n_b': len(Y_b), + }) + + if not group_results: + fig, ax = plt.subplots(figsize=(8, 4)) + ax.text(0.5, 0.5, + f"Two-Sample Second-Level Analysis Aborted\n\n" + f"Reason: All {len(channels)} channels skipped.\n" + f"Requirement: At least {min_subjects} subjects per group, per channel.\n" + f"Group A total IDs: {len(set(group_a_ids))} | " + f"Group B total IDs: {len(set(group_b_ids))}", + ha='center', va='center', fontsize=12, color='darkred', + bbox=dict(facecolor='white', alpha=0.5, edgecolor='red')) + ax.set_axis_off() + plt.show() + return pd.DataFrame() + + df_group = pd.DataFrame(group_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: + logger.warning("No multiple comparisons correction applied — " + "raw p-values will be used for thresholding.") + df_group['p_corrected'] = df_group['p_val'] + df_group['significant'] = df_group['p_val'] <= p + + logger.info("Two-sample second-level results:\n%s", df_group) + + # Extract source/detector positions + src_pos: dict[int, tuple[float, float]] = {} + det_pos: dict[int, tuple[float, float]] = {} + for ch in getattr(raw, "info")["chs"]: + ch_name = ch['ch_name'] + if not ch_name or not ch['loc'].any() or '_' not in ch_name: + continue + try: + parts = ch_name.split()[0] + src_str, det_str = parts.split('_') + src_num = int(src_str[1:]) + det_num = int(det_str[1:]) + src_pos[src_num] = ch['loc'][3:5] + det_pos[det_num] = ch['loc'][6:8] + except Exception as e: + logger.debug(f"Skipping non-fNIRS channel geometry for {ch_name}: {e}") + + fig, ax = plt.subplots(figsize=(8, 6)) + + for pos in src_pos.values(): + ax.scatter(pos[0], pos[1], s=120, c='k', marker='o', edgecolors='white', linewidths=1, zorder=3) + for pos in det_pos.values(): + ax.scatter(pos[0], pos[1], s=120, c='k', marker='s', edgecolors='white', linewidths=1, zorder=3) + + norm = mcolors.Normalize(vmin=-bounds, vmax=bounds) + cmap = plt.get_cmap('seismic') + + for _, row in df_group.iterrows(): + ch = row['ch_name'] + tval = row['t_val'] + mean_diff = row['mean_diff'] + is_sig = row['significant'] + + if '_' not in ch: + continue + + src_str, det_str = ch.split('_') + det_parts = det_str.split() + detector_id = det_parts[0] + hemo_type = det_parts[1].lower() if len(det_parts) > 1 else '' + + if hemo_type != 'hbo': + continue + + try: + src = int(src_str[1:]) + det = int(detector_id[1:]) + except Exception: + continue + + if src in src_pos and det in det_pos: + x = [src_pos[src][0], det_pos[det][0]] + y = [src_pos[src][1], det_pos[det][1]] + style = '-' if is_sig else '--' + color_val = tval if color_by == 't_val' else mean_diff + color = cmap(norm(color_val)) + ax.plot(x, y, linestyle=style, color=color, linewidth=4, alpha=0.9, zorder=2) + + sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + cbar = plt.colorbar(sm, ax=ax, shrink=0.85) + cbar.set_label( + f'{group_a_label} vs {group_b_label}: t-value (HbO)' if color_by == 't_val' + else f'{group_a_label} vs {group_b_label}: mean difference (HbO)', + fontsize=11 + ) + + ax.set_aspect('equal') + correction_label = correction_method if correction_method else 'uncorrected' + ax.set_title( + f"Group Difference Map (HbO): {group_a_label} vs {group_b_label}, " + f"p<{p} ({correction_label})", + fontsize=14 + ) + ax.set_xlabel('X position (m)', fontsize=11) + ax.set_ylabel('Y position (m)', fontsize=11) + ax.grid(True, alpha=0.3) + + all_x = [pos[0] for pos in src_pos.values()] + [pos[0] for pos in det_pos.values()] + all_y = [pos[1] for pos in src_pos.values()] + [pos[1] for pos in det_pos.values()] + ax.set_xlim(min(all_x)-0.01, max(all_x)+0.01) + ax.set_ylim(min(all_y)-0.01, max(all_y)+0.01) + + fig.tight_layout() + plt.show() + + return df_group + + + + + + + + def calculate_dpf(file_path): # order is hbo / hbr with h5py.File(file_path, 'r') as f: @@ -4163,6 +4480,22 @@ def process_participant(file_path, progress_callback=None): # Step 13: Haemoglobin Concentration raw_haemo = beer_lambert_law(raw_od, ppf=calculate_dpf(file_path)) + + # Temporary test + if SHORT_CHANNEL and not FOLDING_BYP and short_chans is not None: + try: + logger.info("Converting raw short channels to hemoglobin concentration...") + # 1. Convert raw short channel intensity to optical density + short_od = optical_density(short_chans) + + # 2. Convert short channel optical density to hemoglobin concentration + # (using the same DPF calculation function as your main data) + short_chans = beer_lambert_law(short_od, ppf=calculate_dpf(file_path)) + + logger.info("Successfully converted short channels to HbO/HbR concentration.") + except Exception as e: + logger.error(f"Failed to convert short channels to hemoglobin: {e}") + fig_raw_haemo_bll = raw_haemo.plot(duration=raw_haemo.times[-1], n_channels=raw_haemo.info['nchan'], title="HbO and HbR Signals", show=False) fig_individual["BLL"] = fig_raw_haemo_bll if progress_callback: progress_callback(17) @@ -4299,6 +4632,13 @@ def process_participant(file_path, progress_callback=None): # Step 23: Compute Contrast Results contrast_results = {} + # 0 is NOT a baseline. + # AI Explaination: + # When you regress a single condition (e.g., "Tapping_Right") against zero, the GLM asks: + # "Is the signal during Tapping_Right significantly higher than the average signal across the entire run?" + # Because of systemic physiology (the global blood pressure rise that happens during almost any active task), + # the answer is almost always "Yes, the whole head is higher than the average." + if HRF_MODEL == "fir": for cond, contrast_vector in contrast_dict.items(): @@ -4307,6 +4647,59 @@ def process_participant(file_path, progress_callback=None): df["ID"] = file_path contrast_results[cond] = df + else: + exclude_list = ["drift", "constant", "short"] + task_cols = [c for c in design_matrix.columns if not any(ex in c.lower() for ex in exclude_list)] + + # Dictionary to hold all our contrast results + contrast_results = {} + + # 2. Loop through every task to create "Simple" contrasts (Task vs 0) + for cond in task_cols: + # Create a vector of zeros + vec = np.zeros(len(design_matrix.columns)) + # Set the index of our condition to 1 + vec[list(design_matrix.columns).index(cond)] = 1 + + contrast = glm_est.compute_contrast(vec) + df = contrast.to_dataframe() + df["ID"] = file_path + contrast_results[f"{cond}_vs_Zero"] = df + + # 3. Loop through all combinations to create "Differential" contrasts (Task A vs Task B) + # This uses itertools.combinations to get every unique pair + for cond_a, cond_b in itertools.combinations(task_cols, 2): + vec = np.zeros(len(design_matrix.columns)) + vec[list(design_matrix.columns).index(cond_a)] = 1 + vec[list(design_matrix.columns).index(cond_b)] = -1 + + contrast = glm_est.compute_contrast(vec) + df = contrast.to_dataframe() + df["ID"] = file_path + + # Name the key nicely, e.g., "Tapping_Right_vs_Tapping_Left" + contrast_results[f"{cond_a}_vs_{cond_b}"] = df + + #NOTE: Temporary + export_dir = os.path.join(os.path.dirname(file_path), "exported_contrasts") + os.makedirs(export_dir, exist_ok=True) + base = os.path.splitext(os.path.basename(file_path))[0] + + combined_contrasts = [] + for contrast_name, df in contrast_results.items(): + df = df.copy() + df["contrast_name"] = contrast_name + combined_contrasts.append(df) + + if combined_contrasts: + pd.concat(combined_contrasts, ignore_index=True).to_csv( + os.path.join(export_dir, f"{base}_contrasts.csv"), index=False + ) + + df_ind.to_csv(os.path.join(export_dir, f"{base}_df_ind.csv"), index=False) + design_matrix.to_csv(os.path.join(export_dir, f"{base}_design_matrix.csv")) + + logger.info(f"Exported contrast/ROI/design-matrix CSVs to {export_dir}") cha["ID"] = file_path if progress_callback: progress_callback(27)