diff --git a/data_processing.ipynb b/data_processing.ipynb
index 2ffc250..69f8608 100644
--- a/data_processing.ipynb
+++ b/data_processing.ipynb
@@ -89,7 +89,9 @@
     "import nibabel as nib\n",
     "import pandas as pd\n",
     "from scipy.interpolate import interp1d\n",
-    "from scipy.ndimage import uniform_filter1d"
+    "from scipy.ndimage import uniform_filter1d\n",
+    "from scipy.stats import f_oneway\n",
+    "from statsmodels.stats.multicomp import pairwise_tukeyhsd"
    ]
   },
   {
@@ -325,7 +327,8 @@
    "source": [
     "# Make figure of B1+ values along the spinal cord across shim methods\n",
     "\n",
-    "plt.rcParams['font.family'] = 'Arial'\n",
+    "# Go back to root data folder\n",
+    "os.chdir(os.path.join(path_data))\n",
     "\n",
     "def smooth_data(data, window_size=20):\n",
     "    \"\"\" Apply a simple moving average to smooth the data. \"\"\"\n",
@@ -406,6 +409,53 @@
     "plt.tight_layout()\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a66be786",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create tables and perform statistics\n",
+    "\n",
+    "# Go back to root data folder\n",
+    "os.chdir(os.path.join(path_data))\n",
+    "\n",
+    "# Data storage\n",
+    "data_summary = []\n",
+    "\n",
+    "# Compute mean and SD of B1+ for each subject and each shim mode\n",
+    "for subject in subjects:\n",
+    "    for shim_mode in shim_modes:\n",
+    "        all_data = []\n",
+    "        for file_path in glob.glob(os.path.join(path_data, subject, \"fmap\", f\"TB1map_*{shim_mode}*.csv\")):\n",
+    "            df = pd.read_csv(file_path)\n",
+    "            wa_data = df['WA()']\n",
+    "            all_data.extend(wa_data)\n",
+    "\n",
+    "        if all_data:\n",
+    "            A = np.mean(all_data)\n",
+    "            SD = np.std(all_data)\n",
+    "            data_summary.append([subject, shim_mode, A, SD])\n",
+    "\n",
+    "# Convert to DataFrame and save to CSV\n",
+    "df_summary = pd.DataFrame(data_summary, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation'])\n",
+    "df_summary.to_csv('subject_shim_mode_summary.csv', index=False)\n",
+    "\n",
+    "# Step 3: Compute statistics across subjects\n",
+    "df_grouped = df_summary.groupby('Shim_Mode').agg({'Average': 'mean', 'Standard_Deviation': 'mean'}).reset_index()\n",
+    "df_grouped.to_csv('average_across_subjects.csv', index=False)\n",
+    "\n",
+    "# Step 4: ANOVA and Posthoc Tests\n",
+    "anova_result = f_oneway(*[group[\"Average\"].values for name, group in df_summary.groupby(\"Shim_Mode\")])\n",
+    "print(\"ANOVA Result:\", anova_result)\n",
+    "\n",
+    "# Perform posthoc test if ANOVA is significant\n",
+    "if anova_result.pvalue < 0.05:\n",
+    "    posthoc_result = pairwise_tukeyhsd(df_summary['Average'], df_summary['Shim_Mode'])\n",
+    "    print(\"Posthoc Tukey HSD Result:\\n\", posthoc_result)"
+   ]
   }
  ],
  "metadata": {