MSIT_Data_Analysis.py

import os
import re
import math
import pandas as pd
import matplotlib.pyplot as plt

try:
    from scipy import stats
except ImportError:
    stats = None


OUTPUT_DIR = "msit_plots"
FILE_KEYWORD = "control_task"


def find_msit_files(base_dir="."):
    matches = []
    for root, _, files in os.walk(base_dir):
        for file_name in files:
            if file_name.lower().endswith(".csv") and FILE_KEYWORD in file_name.lower():
                matches.append(os.path.join(root, file_name))
    return sorted(matches)


def extract_session_value(df):
    for col in ("Session", "session"):
        if col in df.columns:
            values = df[col].dropna()
            if not values.empty:
                return str(values.iloc[0]).strip()
    return "unknown"


def extract_base_id(file_path):
    file_name = os.path.basename(file_path)
    return "_".join(file_name.split("_")[0:2]).lower()


def simplify_label(full_id):
    matches = re.findall(r"([a-b0-9]+_[a-b0-9]+)", str(full_id).lower())
    if matches:
        return matches[-1]
    return str(full_id).lower()


def build_run_label(file_path, df):
    base_id = extract_base_id(file_path)
    session_value = extract_session_value(df)
    full_id = f"{base_id}_session_{session_value}"
    short_label = simplify_label(full_id)
    return full_id, short_label


def load_rt_data(msit_files):
    rt_rows = []
    summary_rows = []

    for file_path in msit_files:
        df = pd.read_csv(file_path)
        full_id, run_label = build_run_label(file_path, df)

        if "first_resp.corr" not in df.columns or "first_resp.rt" not in df.columns:
            print(f"Skipping {file_path}: missing congruent columns")
            continue

        if "second_response.corr" not in df.columns or "second_response.rt" not in df.columns:
            print(f"Skipping {file_path}: missing incongruent columns")
            continue

        congruent_correct = df.loc[df["first_resp.corr"] == 1, "first_resp.rt"].dropna().astype(float)
        congruent_all_corr = df["first_resp.corr"].dropna()

        incongruent_correct = df.loc[df["second_response.corr"] == 1, "second_response.rt"].dropna().astype(float)
        incongruent_all_corr = df["second_response.corr"].dropna()

        for value in congruent_correct:
            rt_rows.append(
                {
                    "full_id": full_id,
                    "run_label": run_label,
                    "condition": "Congruent",
                    "rt": float(value),
                }
            )

        for value in incongruent_correct:
            rt_rows.append(
                {
                    "full_id": full_id,
                    "run_label": run_label,
                    "condition": "Incongruent",
                    "rt": float(value),
                }
            )

        congruent_accuracy = (congruent_all_corr == 1).sum() / len(congruent_all_corr) if len(congruent_all_corr) else float("nan")
        incongruent_accuracy = (incongruent_all_corr == 1).sum() / len(incongruent_all_corr) if len(incongruent_all_corr) else float("nan")

        summary_rows.append(
            {
                "full_id": full_id,
                "run_label": run_label,
                "congruent_accuracy": congruent_accuracy,
                "congruent_mean": congruent_correct.mean(),
                "congruent_var": congruent_correct.var(),
                "congruent_std": congruent_correct.std(),
                "congruent_n": len(congruent_correct),
                "incongruent_accuracy": incongruent_accuracy,
                "incongruent_mean": incongruent_correct.mean(),
                "incongruent_var": incongruent_correct.var(),
                "incongruent_std": incongruent_correct.std(),
                "incongruent_n": len(incongruent_correct),
                "interference_mean": incongruent_correct.mean() - congruent_correct.mean(),
            }
        )

    rt_df = pd.DataFrame(rt_rows)
    summary_df = pd.DataFrame(summary_rows).sort_values("run_label").reset_index(drop=True)
    return rt_df, summary_df


def ensure_output_dir():
    os.makedirs(OUTPUT_DIR, exist_ok=True)


def make_violin_plot(rt_df, summary_df, condition, output_name):
    condition_rt = rt_df[rt_df["condition"] == condition].copy()
    condition_summary = summary_df.copy()

    labels = condition_summary["run_label"].tolist()
    metric_prefix = condition.lower()
    data = [
        condition_rt.loc[condition_rt["run_label"] == label, "rt"].dropna().tolist()
        for label in labels
    ]

    positions = list(range(1, len(labels) + 1))

    fig, ax = plt.subplots(figsize=(max(10, len(labels) * 0.7), 6))

    violin_data = [d for d in data if len(d) > 0]
    violin_positions = [positions[i] for i, d in enumerate(data) if len(d) > 0]

    if violin_data:
        parts = ax.violinplot(
            violin_data,
            positions=violin_positions,
            widths=0.8,
            showmeans=False,
            showmedians=False,
            showextrema=False,
        )
        for body in parts["bodies"]:
            body.set_facecolor("#87aade" if condition == "Congruent" else "#d98b8b")
            body.set_edgecolor("black")
            body.set_alpha(0.5)

    means = condition_summary[f"{metric_prefix}_mean"].tolist()
    stds = condition_summary[f"{metric_prefix}_std"].fillna(0).tolist()

    ax.errorbar(
        positions,
        means,
        yerr=stds,
        fmt="o",
        color="black",
        ecolor="black",
        elinewidth=1.5,
        capsize=4,
        label="Mean ± SD",
        zorder=3,
    )

    ax.set_title(f"{condition} Reaction Time by Run")
    ax.set_ylabel("Reaction time (s)")
    ax.set_xticks(positions)
    ax.set_xticklabels(labels, rotation=45, ha="right")
    ax.grid(axis="y", alpha=0.3)
    ax.legend()
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, output_name), dpi=300, bbox_inches="tight")
    plt.close()


def make_accuracy_plot(summary_df):
    labels = summary_df["run_label"].tolist()
    x = list(range(len(labels)))

    fig, ax = plt.subplots(figsize=(max(10, len(labels) * 0.7), 5))
    ax.plot(x, summary_df["congruent_accuracy"], marker="o", label="Congruent", color="steelblue")
    ax.plot(x, summary_df["incongruent_accuracy"], marker="o", label="Incongruent", color="indianred")
    ax.set_title("Accuracy by Run")
    ax.set_ylabel("Proportion correct")
    ax.set_xticks(x)
    ax.set_xticklabels(labels, rotation=45, ha="right")
    ax.set_ylim(0.0, 1.05)
    ax.grid(axis="y", alpha=0.3)
    ax.legend()
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, "accuracy_by_run.png"), dpi=300, bbox_inches="tight")
    plt.close()


def make_interference_plot(summary_df):
    labels = summary_df["run_label"].tolist()
    x = list(range(len(labels)))

    fig, ax = plt.subplots(figsize=(max(10, len(labels) * 0.7), 5))
    ax.bar(x, summary_df["interference_mean"], color="darkslategray", alpha=0.85)
    ax.axhline(0, color="black", linewidth=1)
    ax.set_title("Interference Effect by Run")
    ax.set_ylabel("Mean incongruent RT - mean congruent RT (s)")
    ax.set_xticks(x)
    ax.set_xticklabels(labels, rotation=45, ha="right")
    ax.grid(axis="y", alpha=0.3)
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, "interference_by_run.png"), dpi=300, bbox_inches="tight")
    plt.close()


def run_stats(rt_df, summary_df):
    report_lines = []

    report_lines.append("MSIT between-run statistics")
    report_lines.append("=" * 40)
    report_lines.append("")

    if stats is None:
        report_lines.append("scipy is not installed, so inferential statistics were not run.")
        return "\n".join(report_lines)

    for condition in ["Congruent", "Incongruent"]:
        report_lines.append(condition)
        report_lines.append("-" * len(condition))

        subsets = []
        run_names = []

        for run_label in summary_df["run_label"]:
            values = rt_df.loc[
                (rt_df["run_label"] == run_label) & (rt_df["condition"] == condition),
                "rt"
            ].dropna().values

            if len(values) > 1:
                subsets.append(values)
                run_names.append(run_label)

        if len(subsets) < 2:
            report_lines.append("Not enough runs with data for comparison.")
            report_lines.append("")
            continue

        try:
            f_stat, p_anova = stats.f_oneway(*subsets)
            report_lines.append(f"One-way ANOVA across runs: F={f_stat:.4f}, p={p_anova:.6g}")
        except Exception as exc:
            report_lines.append(f"One-way ANOVA failed: {exc}")

        try:
            h_stat, p_kw = stats.kruskal(*subsets)
            report_lines.append(f"Kruskal-Wallis across runs: H={h_stat:.4f}, p={p_kw:.6g}")
        except Exception as exc:
            report_lines.append(f"Kruskal-Wallis failed: {exc}")

        try:
            lev_stat, p_lev = stats.levene(*subsets)
            report_lines.append(f"Levene variance test across runs: W={lev_stat:.4f}, p={p_lev:.6g}")
        except Exception as exc:
            report_lines.append(f"Levene test failed: {exc}")

        report_lines.append("")

    return "\n".join(report_lines)


def main():
    ensure_output_dir()

    msit_files = find_msit_files(".")
    if not msit_files:
        raise FileNotFoundError("No MSIT CSV files found.")

    rt_df, summary_df = load_rt_data(msit_files)

    if rt_df.empty or summary_df.empty:
        raise ValueError("No usable MSIT data found in the discovered CSV files.")

    summary_df.to_csv(os.path.join(OUTPUT_DIR, "msit_run_summary.csv"), index=False)

    make_violin_plot(rt_df, summary_df, "Congruent", "congruent_rt_violin.png")
    make_violin_plot(rt_df, summary_df, "Incongruent", "incongruent_rt_violin.png")
    make_accuracy_plot(summary_df)
    make_interference_plot(summary_df)

    stats_report = run_stats(rt_df, summary_df)
    with open(os.path.join(OUTPUT_DIR, "msit_stats_report.txt"), "w", encoding="utf-8") as f:
        f.write(stats_report)

    print(f"Saved plots and report in ./{OUTPUT_DIR}")
    print(stats_report)


if __name__ == "__main__":
    main()