{ "cells": [ { "cell_type": "markdown", "id": "5ede6850", "metadata": {}, "source": [ "Effect on Time Binding" ] }, { "cell_type": "code", "execution_count": 4, "id": "4034081f", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "from scipy.stats import mannwhitneyu\n", "\n", "# -----------------------------------------------------------\n", "# 1. Load Data\n", "# -----------------------------------------------------------\n", "df_scr = pd.read_csv('TTB_SCR.csv') # Self-Control Reflection Group\n", "df_slr = pd.read_csv('TTB_SLR.csv') # Selfless Reflection Group\n", "df_ctr = pd.read_csv('TTB_CTR.csv') # Control Group\n", "\n", "# -----------------------------------------------------------\n", "# 2. Helper Function to Print Mean ± SD (×1000)\n", "# -----------------------------------------------------------\n", "def summarize_variable(df_list, labels, column):\n", " for df, label in zip(df_list, labels):\n", " mean_val = df[column].mean() * 1000\n", " sd_val = df[column].std() * 1000\n", "\n", "# -----------------------------------------------------------\n", "# 3. Descriptive Statistics (Scaled to Match Table)\n", "# -----------------------------------------------------------\n", "groups = [df_scr, df_slr, df_ctr]\n", "labels = [\"Group SCR\", \"Group SLR\", \"Group CTR\"]\n", "\n", "summarize_variable(groups, labels, 'BindingScore_action')\n", "summarize_variable(groups, labels, 'BindingScore_tone')\n", "summarize_variable(groups, labels, 'TTB')\n", "\n", "# -----------------------------------------------------------\n", "# 4. Inferential Statistics (Mann–Whitney U Tests)\n", "# -----------------------------------------------------------\n", "ttb_scr = df_scr['TTB']\n", "ttb_slr = df_slr['TTB']\n", "ttb_ctr = df_ctr['TTB']\n", "\n", "mw_scr_slr = mannwhitneyu(ttb_scr, ttb_slr, alternative='greater')\n", "mw_slr_ctr = mannwhitneyu(ttb_slr, ttb_ctr, alternative='less')\n", "mw_scr_ctr = mannwhitneyu(ttb_scr, ttb_ctr, alternative='greater')\n" ] }, { "cell_type": "markdown", "id": "4db2ffc1", "metadata": {}, "source": [ "Kruskwal wallis for TTB" ] }, { "cell_type": "code", "execution_count": 5, "id": "745cb376", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "from scipy.stats import kruskal\n", "import numpy as np\n", "\n", "# Read CSVs\n", "df_slr = pd.read_csv('df_corr_SLR_AF.csv')\n", "df_ctr = pd.read_csv('df_corr_CTR_AF.csv')\n", "df_scr = pd.read_csv('df_corr_SCR_AF.csv')\n", "\n", "# Get TTB columns\n", "ttb_slr = df_slr['TTB']\n", "ttb_ctr = df_ctr['TTB']\n", "ttb_scr = df_scr['TTB']\n", "\n", "# Kruskal-Wallis test\n", "kw_stat, kw_p = kruskal(ttb_slr, ttb_ctr, ttb_scr)\n", "\n", "# Effect size (eta squared)\n", "groups = [ttb_slr, ttb_ctr, ttb_scr]\n", "all_data = np.concatenate(groups)\n", "group_labels = np.concatenate([[i]*len(g) for i, g in enumerate(groups)])\n", "ss_between = sum([len(g) * (np.mean(g) - np.mean(all_data))**2 for g in groups])\n", "ss_total = np.sum((all_data - np.mean(all_data))**2)\n", "eta_sq = ss_between / ss_total\n" ] }, { "cell_type": "markdown", "id": "d609b840", "metadata": {}, "source": [ "EEG Features and stastical analysis\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "928c884c", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import scipy.stats as stats\n", "from sklearn.decomposition import KernelPCA, PCA\n", "from gtda.time_series import TakensEmbedding, PermutationEntropy\n", "from gtda.homology import VietorisRipsPersistence\n", "from gtda.diagrams import BettiCurve, Amplitude\n", "from sklearn.preprocessing import StandardScaler\n", "from scipy.signal import resample\n", "\n", "\n", "# --- 2. Load Topological Data Analysis (TDA) Features ---\n", "# Data files contain dictionaries of TDA features ('amp', 'pe', 'dg') \n", "# SCR = Self-Centered Reflection, SLR = Selfless Reflection, CTR = Control\n", "# 'af' and 'tp' are two distinct EEG channels.\n", "\n", "# Load 'af' channel data\n", "scr_tda_af = np.load('SCR_tda_af.npy', allow_pickle=True).item()\n", "slr_tda_af = np.load('SLR_tda_af.npy', allow_pickle=True).item()\n", "ctr_tda_af = np.load('ctr_tda_af.npy', allow_pickle=True).item()\n", "\n", "# Load 'tp' channel data\n", "scr_tda_tp = np.load('SCR_tda_tp.npy', allow_pickle=True).item()\n", "slr_tda_tp = np.load('SLR_tda_tp.npy', allow_pickle=True).item()\n", "ctr_tda_tp = np.load('ctr_tda_tp.npy', allow_pickle=True).item()\n", "\n", "\n", "# --- 3. Betti Curve Area Feature Calculation Function ---\n", "def get_betti_areas(diagrams):\n", " \"\"\"Computes Betti Curves and calculates the area under each curve.\"\"\"\n", " BC = BettiCurve(n_bins=100)\n", " bcs = BC.fit_transform(diagrams)\n", " X = BC.samplings_\n", " ba_sw = []\n", " \n", " for bc in bcs:\n", " temp_areas = []\n", " for bci, bc_hg in enumerate(bc):\n", " # Calculate area using trapezoidal rule\n", " temp_areas.append(np.trapz(bc_hg, X[bci]))\n", " ba_sw.append(temp_areas)\n", " \n", " return ba_sw\n", "\n", "# Define common column names (Feature_Dimension_Channel)\n", "columns=['A0_tp','A1_tp','A2_tp','A0_af','A1_af','A2_af',\n", " 'PE0_tp','PE1_tp','PE2_tp','PE0_af','PE1_af','PE2_af']\n", "\n", "# --- 4. Feature Extraction and DataFrame Creation (SLR Group) ---\n", "\n", "# SLR 'af' channel features: Amplitude and Permutation Entropy for dimensions 0, 1, 2\n", "A0_af, A1_af, A2_af = slr_tda_af['amp'][:,0], slr_tda_af['amp'][:,1], slr_tda_af['amp'][:,2]\n", "PE0_af, PE1_af, PE2_af = slr_tda_af['pe'][:,0], slr_tda_af['pe'][:,1], slr_tda_af['pe'][:,2]\n", "# Calculate Betti Area (BA) - Note: BA is computed but not included in final 'arr_slr'\n", "BA_af = np.asarray(get_betti_areas(slr_tda_af['dg'])) \n", "BA0_af, BA1_af, BA2_af = BA_af[:,0], BA_af[:,1], BA_af[:,2]\n", "\n", "# SLR 'tp' channel features\n", "A0_tp, A1_tp, A2_tp = slr_tda_tp['amp'][:,0], slr_tda_tp['amp'][:,1], slr_tda_tp['amp'][:,2]\n", "PE0_tp, PE1_tp, PE2_tp = slr_tda_tp['pe'][:,0], slr_tda_tp['pe'][:,1], slr_tda_tp['pe'][:,2]\n", "BA_tp = np.asarray(get_betti_areas(slr_tda_tp['dg']))\n", "BA0_tp, BA1_tp, BA2_tp = BA_tp[:,0], BA_tp[:,1], BA_tp[:,2]\n", "\n", "arr_slr=[A0_tp,A1_tp,A2_tp,A0_af,A1_af,A2_af,\n", " PE0_tp,PE1_tp,PE2_tp,PE0_af,PE1_af,PE2_af]\n", "df_slr = pd.DataFrame(np.asarray(arr_slr).T, columns=columns)\n", "\n", "\n", "# --- 5. Feature Extraction and DataFrame Creation (CTR Group) ---\n", "\n", "# CTR 'af' features\n", "A0_af, A1_af, A2_af = ctr_tda_af['amp'][:,0], ctr_tda_af['amp'][:,1], ctr_tda_af['amp'][:,2]\n", "PE0_af, PE1_af, PE2_af = ctr_tda_af['pe'][:,0], ctr_tda_af['pe'][:,1], ctr_tda_af['pe'][:,2]\n", "# CTR 'tp' features\n", "A0_tp, A1_tp, A2_tp = ctr_tda_tp['amp'][:,0], ctr_tda_tp['amp'][:,1], ctr_tda_tp['amp'][:,2]\n", "PE0_tp, PE1_tp, PE2_tp = ctr_tda_tp['pe'][:,0], ctr_tda_tp['pe'][:,1], ctr_tda_tp['pe'][:,2]\n", "\n", "# Combine A and PE features for CTR\n", "arr_ctr=[A0_tp,A1_tp,A2_tp,A0_af,A1_af,A2_af,\n", " PE0_tp,PE1_tp,PE2_tp,PE0_af,PE1_af,PE2_af]\n", "df_ctr = pd.DataFrame(np.asarray(arr_ctr).T, columns=columns)\n", "\n", "\n", "# --- 6. Feature Extraction and DataFrame Creation (SCR Group) ---\n", "\n", "# SCR 'af' features\n", "A0_af, A1_af, A2_af = scr_tda_af['amp'][:,0], scr_tda_af['amp'][:,1], scr_tda_af['amp'][:,2]\n", "PE0_af, PE1_af, PE2_af = scr_tda_af['pe'][:,0], scr_tda_af['pe'][:,1], scr_tda_af['pe'][:,2]\n", "# SCR 'tp' features\n", "A0_tp, A1_tp, A2_tp = scr_tda_tp['amp'][:,0], scr_tda_tp['amp'][:,1], scr_tda_tp['amp'][:,2]\n", "PE0_tp, PE1_tp, PE2_tp = scr_tda_tp['pe'][:,0], scr_tda_tp['pe'][:,1], scr_tda_tp['pe'][:,2]\n", "\n", "# Combine A and PE features for SCR\n", "arr_scr=[A0_tp,A1_tp,A2_tp,A0_af,A1_af,A2_af,\n", " PE0_tp,PE1_tp,PE2_tp,PE0_af,PE1_af,PE2_af]\n", "df_scr = pd.DataFrame(np.asarray(arr_scr).T, columns=columns)\n", "\n", "\n", "# --- 7. Betti Curve Computation for Visualization ---\n", "BC = BettiCurve()\n", "\n", "# Compute Betti curves and filtration values for visualization\n", "BC_ctr_af = BC.fit_transform(ctr_tda_af['dg'])\n", "x_ctr_af = BC._samplings\n", "BC_ctr_tp = BC.fit_transform(ctr_tda_tp['dg'])\n", "x_ctr_tp = BC._samplings\n", "BC_scr_af = BC.fit_transform(scr_tda_af['dg'])\n", "x_scr_af = BC._samplings\n", "BC_scr_tp = BC.fit_transform(scr_tda_tp['dg'])\n", "x_scr_tp = BC._samplings\n", "BC_slr_af = BC.fit_transform(slr_tda_af['dg'])\n", "x_slr_af = BC._samplings\n", "BC_slr_tp = BC.fit_transform(slr_tda_tp['dg'])\n", "x_slr_tp = BC._samplings\n", "\n", "\n", "columns=['A0_tp','A1_tp','A2_tp','A0_af','A1_af','A2_af',\n", " 'PE0_tp','PE1_tp','PE2_tp','PE0_af','PE1_af','PE2_af']\n", "\n", "\n", "\n", "mean_SCResults = {}\n", "for col in columns:\n", " mean_SCResults[col] = [\n", " np.round(df_scr[col].dropna().mean(),2), \n", " np.round(df_slr[col].dropna().mean(),2),\n", " np.round(df_ctr[col].dropna().mean(),2)]\n", " \n", "kruskal_SCResults = {}\n", "for col in columns:\n", " stat, p_value = stats.kruskal(df_slr[col].dropna(), df_scr[col].dropna(), df_ctr[col].dropna())\n", " kruskal_SCResults[col] = {'statistic': stat, 'p_value': p_value} # Keep p_value internally for a full record\n", " \n", "\n", "# Create a DataFrame from the mean_SCResults dictionary\n", "df_means = pd.DataFrame(mean_SCResults, index=['SCR (Self-Centered)', 'SLR (Selfless)', 'CTR (Control)']).T\n", "\n", "\n", "# ----------------------------------------------------------------------\n", "\n", "## Kruskal-Wallis H-Test Statistic (SLR vs SCR vs CTR) 📊\n", "\n", "# Create a DataFrame containing only the Statistic\n", "kruskal_statistic_data = {\n", " 'Statistic': [v['statistic'] for v in kruskal_SCResults.values()],\n", "}\n", "df_kruskal_stat = pd.DataFrame(kruskal_statistic_data, index=columns)\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "ac4c3a34", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "from scipy.stats import kruskal as scipy_kruskal\n", "\n", "# List of (measure name, H statistic) tuples from a previous analysis\n", "kruskal_SCResults = [\n", " (\"A0_tp\", 36.22630681692385),\n", " (\"A1_tp\", 46.08434264322767),\n", " (\"A2_tp\", 45.20543796724746),\n", " (\"A0_af\", 32.0750728776909),\n", " (\"A1_af\", 29.350138467612737),\n", " (\"A2_af\", 22.72526992420717),\n", " (\"PE0_tp\", 36.02875117157811),\n", " (\"PE1_tp\", 45.20366676260687),\n", " (\"PE2_tp\", 5.859694854650513),\n", " (\"PE0_af\", 32.0750728776909),\n", " (\"PE1_af\", 28.670451140581974),\n", " (\"PE2_af\", 21.612283673183185),\n", "]\n", "\n", "# Constants for effect size calculation\n", "n = 31 + 31 + 30 # Total sample size (n=92)\n", "k = 3 # Number of groups \n", "\n", "\n", "# Epsilon-squared ($\\epsilon^2$) formula for Kruskal-Wallis\n", "# $\\epsilon^2 = (H - k + 1) / (n - k)$\n", "def epsilon_squared(H, n, k):\n", " \"\"\"Calculates epsilon-squared effect size for Kruskal-Wallis test.\"\"\"\n", " return (H - k + 1) / (n - k)\n", "\n", "\n", "for name, H in kruskal_SCResults:\n", " eps2 = epsilon_squared(H, n, k)\n", " # Output formatted to three decimal places for clarity\n", "\n", "\n", "\n", "# ==============================================================================\n", "# SECTION 2: Kruskal-Wallis Tests for AF Channel Data\n", "# This section loads data for the 'AF' channel and runs Kruskal-Wallis tests\n", "# for 'vn0', 'vn1', and 'vn2' measures across the three groups.\n", "# ==============================================================================\n", "\n", "\n", "df_SLR_af = pd.read_csv('df_corr_SLR_AF.csv')\n", "df_ctr_af = pd.read_csv('df_corr_CTR_AF.csv')\n", "df_SCR_af = pd.read_csv('df_corr_SCR_AF.csv')\n", "\n", "# --- Helper Function for Kruskal-Wallis ---\n", "# Renamed from the original to clearly indicate it uses the imported scipy function\n", "def run_kruskal_with_eps2(group1, group2, group3):\n", " \"\"\"\n", " Runs Kruskal-Wallis test and calculates H, p-value, and epsilon-squared.\n", " \"\"\"\n", " # Perform the Kruskal-Wallis H-test\n", " H, p = scipy_kruskal(group1, group2, group3)\n", "\n", " # Re-calculate constants: n and k are determined by the input data\n", " n = len(group1) + len(group2) + len(group3)\n", " k = 3 # Number of groups\n", "\n", " # Calculate Epsilon-squared ($\\epsilon^2$) effect size\n", " eps2 = (H - k + 1) / (n - k)\n", " return H, p, eps2\n", "# ------------------------------------------\n", "\n", "# Extracting 'vn0', 'vn1', 'vn2' columns for the three groups\n", "# 'vn' stands for von Neumann or Hodge spectral entropy measure.\n", "\n", "# vn0 data extraction\n", "vn0_SLR = df_SLR_af['vn0']\n", "vn0_ctr = df_ctr_af['vn0']\n", "vn0_SCR = df_SCR_af['vn0']\n", "\n", "# vn1 data extraction\n", "vn1_SLR = df_SLR_af['vn1']\n", "vn1_ctr = df_ctr_af['vn1']\n", "vn1_SCR = df_SCR_af['vn1']\n", "\n", "# vn2 data extraction\n", "vn2_SLR = df_SLR_af['vn2']\n", "vn2_ctr = df_ctr_af['vn2']\n", "vn2_SCR = df_SCR_af['vn2']\n", "\n", "# Run Kruskal-Wallis tests and print results\n", "H_vn0, p_vn0, epsilon_squared_vn0 = run_kruskal_with_eps2(vn0_SLR, vn0_ctr, vn0_SCR)\n", "\n", "H_vn1, p_vn1, epsilon_squared_vn1 = run_kruskal_with_eps2(vn1_SLR, vn1_ctr, vn1_SCR)\n", "\n", "H_vn2, p_vn2, epsilon_squared_vn2 = run_kruskal_with_eps2(vn2_SLR, vn2_ctr, vn2_SCR)\n", "\n", "\n", "# ==============================================================================\n", "# SECTION 3: Kruskal-Wallis Tests for TP Channel Data\n", "# This section repeats the analysis for the 'TP' channel data.\n", "# ==============================================================================\n", "\n", "\n", "# Read CSVs for the TP channel data\n", "df_SLR_tp = pd.read_csv('df_corr_SLR_TP.csv')\n", "df_ctr_tp = pd.read_csv('df_corr_CTR_TP.csv')\n", "df_SCR_tp = pd.read_csv('df_corr_SCR_TP.csv')\n", "\n", "# Extracting 'vn0', 'vn1', 'vn2' columns for the three groups\n", "\n", "# vn0 data extraction\n", "vn0_SLR = df_SLR_tp['vn0']\n", "vn0_ctr = df_ctr_tp['vn0']\n", "vn0_SCR = df_SCR_tp['vn0']\n", "\n", "# vn1 data extraction\n", "vn1_SLR = df_SLR_tp['vn1']\n", "vn1_ctr = df_ctr_tp['vn1']\n", "vn1_SCR = df_SCR_tp['vn1']\n", "\n", "# vn2 data extraction\n", "vn2_SLR = df_SLR_tp['vn2']\n", "vn2_ctr = df_ctr_tp['vn2']\n", "vn2_SCR = df_SCR_tp['vn2']\n", "\n", "H_vn0, p_vn0, epsilon_squared_vn0 = run_kruskal_with_eps2(vn0_SLR, vn0_ctr, vn0_SCR)\n", "\n", "H_vn1, p_vn1, epsilon_squared_vn1 = run_kruskal_with_eps2(vn1_SLR, vn1_ctr, vn1_SCR)\n", "\n", "H_vn2, p_vn2, epsilon_squared_vn2 = run_kruskal_with_eps2(vn2_SLR, vn2_ctr, vn2_SCR)\n" ] }, { "cell_type": "markdown", "id": "dcb2700c", "metadata": {}, "source": [ "Correlation of TTB with Hodge spectral entropy" ] }, { "cell_type": "code", "execution_count": 8, "id": "c4c38957", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "from scipy.stats import pearsonr\n", "\n", "# Load all CSV files\n", "df_corr_SLR_AF = pd.read_csv('df_corr_SLR_AF.csv')\n", "df_corr_CTR_AF = pd.read_csv('df_corr_CTR_AF.csv')\n", "df_corr_SCR_AF = pd.read_csv('df_corr_SCR_AF.csv')\n", "\n", "df_corr_SLR_TP = pd.read_csv('df_corr_SLR_TP.csv')\n", "df_corr_CTR_TP = pd.read_csv('df_corr_CTR_TP.csv')\n", "df_corr_SCR_TP = pd.read_csv('df_corr_SCR_TP.csv')\n", "\n", "# List of dataframes and their labels\n", "dfs = [\n", " df_corr_SCR_TP,\n", " df_corr_SCR_AF,\n", " df_corr_SLR_TP,\n", " df_corr_SLR_AF,\n", " df_corr_CTR_TP,\n", " df_corr_CTR_AF\n", "]\n", "\n", "dfs_names = ['SCR_TP', 'SCR_AF', 'SLR_TP', 'SLR_AF', 'CTR_TP', 'CTR_AF']\n", "\n", "# Loop through each dataframe and compute correlation + p-value for TTB and vn2\n", "for name, df in zip(dfs_names, dfs):\n", " # Remove rows with NaN in the two columns\n", " subset = df[['TTB', 'vn2']].dropna()\n", "\n", " # Calculate Pearson correlation and p-value\n", " corr, pval = pearsonr(subset['TTB'], subset['vn2'])\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.16" } }, "nbformat": 4, "nbformat_minor": 5 }