From ede5147830c8d2355bc85d1fc0576d37595534a8 Mon Sep 17 00:00:00 2001 From: urme-b <148221893+urme-b@users.noreply.github.com> Date: Tue, 7 Jul 2026 17:12:00 +0600 Subject: [PATCH] fdr correction on group corr --- group/Group Correlation Matrix.ipynb | 74 +++++++++++++++++++++++++++- 1 file changed, 72 insertions(+), 2 deletions(-) diff --git a/group/Group Correlation Matrix.ipynb b/group/Group Correlation Matrix.ipynb index a4dfe51..9b00ca6 100755 --- a/group/Group Correlation Matrix.ipynb +++ b/group/Group Correlation Matrix.ipynb @@ -19,7 +19,77 @@ { "cell_type": "code", "id": "panyy1gp0ur", - "source": "import pandas as pd\nimport numpy as np\nfrom scipy import stats\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nOUTPUT = '../data/group_results'\n\n# load data\nhrv_df = pd.read_csv(f'{OUTPUT}/HRV_SDNN.csv').set_index('Participant')\npupil_df = pd.read_csv(f'{OUTPUT}/Pupil_Dilation_STD.csv').set_index('Participant')\nduration_df = pd.read_csv(f'{OUTPUT}/Psychometric_Test_Duration_STD.csv').set_index('Participant')\n\n# combine all\ncombined = pd.concat([\n hrv_df.add_suffix('_HRV_SDNN'),\n pupil_df.add_suffix('_Pupil_STD'),\n duration_df.add_suffix('_Duration_STD'),\n], axis=1)\n\ncols = combined.columns\nn = len(cols)\n\n# correlation with p-values\nr_matrix = np.zeros((n, n))\np_matrix = np.zeros((n, n))\n\nfor i in range(n):\n for j in range(n):\n r, p = stats.pearsonr(combined.iloc[:, i], combined.iloc[:, j])\n r_matrix[i, j] = r\n p_matrix[i, j] = p\n\nr_df = pd.DataFrame(r_matrix, index=cols, columns=cols)\np_df = pd.DataFrame(p_matrix, index=cols, columns=cols)\n\n# annotation with significance stars\nannot = np.empty_like(r_matrix, dtype=object)\nfor i in range(n):\n for j in range(n):\n stars = \"***\" if p_df.iloc[i, j] < 0.001 else \"**\" if p_df.iloc[i, j] < 0.01 else \"*\" if p_df.iloc[i, j] < 0.05 else \"\"\n annot[i, j] = f\"{r_df.iloc[i, j]:.2f}{stars}\"\n\n# plot heatmap\nplt.figure(figsize=(14, 11))\nsns.heatmap(r_df, annot=annot, fmt='', cmap='coolwarm', vmin=-1, vmax=1, center=0)\nplt.title('Correlation Matrix with Significance (* p<.05, ** p<.01, *** p<.001)')\nplt.tight_layout()\nplt.show()\nplt.close()", + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "from scipy import stats\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "\n", + "OUTPUT = '../data/group_results'\n", + "\n", + "# load data\n", + "hrv_df = pd.read_csv(f'{OUTPUT}/HRV_SDNN.csv').set_index('Participant')\n", + "pupil_df = pd.read_csv(f'{OUTPUT}/Pupil_Dilation_STD.csv').set_index('Participant')\n", + "duration_df = pd.read_csv(f'{OUTPUT}/Psychometric_Test_Duration_STD.csv').set_index('Participant')\n", + "\n", + "# combine all\n", + "combined = pd.concat([\n", + " hrv_df.add_suffix('_HRV_SDNN'),\n", + " pupil_df.add_suffix('_Pupil_STD'),\n", + " duration_df.add_suffix('_Duration_STD'),\n", + "], axis=1)\n", + "\n", + "cols = combined.columns\n", + "n = len(cols)\n", + "n_obs = len(combined)\n", + "\n", + "# pairwise correlations + raw p-values\n", + "r_matrix = np.zeros((n, n))\n", + "p_matrix = np.ones((n, n))\n", + "for i in range(n):\n", + " for j in range(n):\n", + " r, p = stats.pearsonr(combined.iloc[:, i], combined.iloc[:, j])\n", + " r_matrix[i, j] = r\n", + " p_matrix[i, j] = p\n", + "\n", + "r_df = pd.DataFrame(r_matrix, index=cols, columns=cols)\n", + "\n", + "# Benjamini-Hochberg FDR across the unique off-diagonal pairs only.\n", + "# With n_obs = 10 participants and dozens of pairwise tests, uncorrected\n", + "# stars overstate significance; we control the false-discovery rate instead.\n", + "iu = np.triu_indices(n, k=1)\n", + "raw_p = p_matrix[iu]\n", + "order = np.argsort(raw_p)\n", + "ranks = np.empty_like(order)\n", + "ranks[order] = np.arange(1, len(raw_p) + 1)\n", + "q = raw_p * len(raw_p) / ranks\n", + "# enforce monotonicity of BH q-values\n", + "q_sorted = np.minimum.accumulate(q[order][::-1])[::-1]\n", + "q_adj = np.empty_like(q)\n", + "q_adj[order] = np.clip(q_sorted, 0, 1)\n", + "\n", + "q_matrix = np.ones((n, n))\n", + "q_matrix[iu] = q_adj\n", + "q_matrix[(iu[1], iu[0])] = q_adj\n", + "\n", + "def star(qv):\n", + " return \"***\" if qv < 0.001 else \"**\" if qv < 0.01 else \"*\" if qv < 0.05 else \"\"\n", + "\n", + "annot = np.empty((n, n), dtype=object)\n", + "for i in range(n):\n", + " for j in range(n):\n", + " s = \"\" if i == j else star(q_matrix[i, j])\n", + " annot[i, j] = f\"{r_df.iloc[i, j]:.2f}{s}\"\n", + "\n", + "plt.figure(figsize=(14, 11))\n", + "sns.heatmap(r_df, annot=annot, fmt='', cmap='coolwarm', vmin=-1, vmax=1, center=0)\n", + "plt.title(f'Correlation Matrix (Pearson r, n={n_obs}); '\n", + " f'stars = Benjamini-Hochberg FDR q-value (* q<.05, ** q<.01, *** q<.001)')\n", + "plt.tight_layout()\n", + "plt.show()\n", + "plt.close()\n" + ], "metadata": {}, "execution_count": null, "outputs": [] @@ -46,4 +116,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +}