{ "cells": [ { "cell_type": "code", "execution_count": 2, "id": "8180413b-33aa-4136-9aea-6a7c7251a020", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The autoreload extension is already loaded. To reload it, use:\n", " %reload_ext autoreload\n" ] } ], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", "%matplotlib inline\n", "\n", "import sys\n", "sys.path.append('../waves')\n", "from IPython.display import display\n", "from dbfread import DBF, FieldParser\n", "import matplotlib.pyplot as plt\n", "from utils_functions import *\n", "from build_df_functions import *\n", "from dbfread import DBF\n", "from dbfread.field_parser import FieldParser\n", "import pandas as pd\n", "import seaborn as sns\n", "import scipy.stats as stats\n", "from scipy.stats import wilcoxon\n", "from scipy.stats import t\n", "from scipy.stats import shapiro\n", "from collections import defaultdict\n", "import re\n", "from sklearn.cluster import KMeans\n", "from sklearn.preprocessing import StandardScaler\n", "import statsmodels.formula.api as smf\n", "import statsmodels.api as sm\n", "from sklearn.ensemble import RandomForestClassifier\n", "from sklearn.datasets import make_classification\n", "from sklearn.model_selection import train_test_split\n", "from sklearn.metrics import classification_report\n", "from sklearn.linear_model import Ridge\n", "from pandas.api.types import (is_numeric_dtype, is_categorical_dtype, is_object_dtype, is_integer_dtype, is_float_dtype)\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib.ticker import MaxNLocator\n", "from pandas.api.types import is_numeric_dtype, is_bool_dtype\n", "import scipy.stats as sts\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore') \n", "from statsmodels.stats.multitest import multipletests\n", "\n", "pd.options.mode.chained_assignment = None # default='warn'\n", "\n", "\n", "def fmt_small_conditional(x):\n", " if isinstance(x, (float, int)) and x != 0 and abs(x) < 1e-3:\n", " return f\"{x:.2e}\"\n", " elif isinstance(x, (float, int)):\n", " return f\"{x:.3f}\"\n", " else:\n", " return x\n", " \n", "fmt_small = '{:.2f}'\n", "fmt_integer = '{:.0f}'\n", "fmt_boolean = lambda x: 'True' if x else 'False'" ] }, { "cell_type": "code", "execution_count": 4, "id": "2371c501-a335-4baa-bef2-2d3e4026d421", "metadata": {}, "outputs": [], "source": [ "feature_df = pd.read_csv('./data/FERTILI_features.csv', sep='\\t')\n", "#feature_df = feature_df.rename(columns={'Avg_Peak-to-Peak_Distance': 'avg_peak_peak_distance', 'Intra-monthly_Variability': 'monthly_Variability', 'Sawtooth Up_r2': 'Sawtooth_up_r2', 'Sawtooth Up_period': 'Sawtooth_up_period', 'Sawtooth Up_phase': 'Sawtooth_up_phase'})\n", "feature_df.columns = feature_df.columns.str.replace('-', '_').str.replace(' ', '_')" ] }, { "cell_type": "code", "execution_count": 5, "id": "1d9edaa8-f33e-49ee-b3a0-eb2ef92fb164", "metadata": {}, "outputs": [], "source": [ "# Add age group\n", "feature_df = add_age_group(feature_df, age_bins=[18, 35, 42])\n", "category_mapping = {'18-35': 1, '35-42': 2}\n", "feature_df['Age_Group'] = feature_df['Age_Group'].replace(category_mapping)" ] }, { "cell_type": "markdown", "id": "a4d78692-0971-4e85-b5ca-391240d8cb05", "metadata": { "tags": [] }, "source": [ "#### Some stats" ] }, { "cell_type": "code", "execution_count": 6, "id": "ef768f01-e0d1-4be1-a1bd-5cdadbd09c79", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "num of cycles 6175\n", "num of participants 857\n" ] } ], "source": [ "print('num of cycles', len(feature_df))\n", "print('num of participants', len(feature_df.CODICE.unique()))" ] }, { "cell_type": "code", "execution_count": 9, "id": "2fd611d5-ad52-44c5-9e1f-a3617806603e", "metadata": {}, "outputs": [], "source": [ "features = feature_df.columns.drop(['CODICE', 'QUALIFI'])" ] }, { "cell_type": "code", "execution_count": 10, "id": "b32f902a-8043-4909-975c-8e55059de35f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "7.20536756126021" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "count_cycles = feature_df.groupby(\"CODICE\")[\"Age_at_Menses\"].count()\n", "count_cycles.mean()" ] }, { "cell_type": "code", "execution_count": 11, "id": "f346a0e1-c08c-41c4-afa8-a0fce983f79e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Median = 4 cycles (IQR 2–9)\n", "Mean = 7.2 ± 7.9 cycles\n", "Range = 1–63 cycles\n" ] } ], "source": [ "\n", "med = count_cycles.median()\n", "q25 = count_cycles.quantile(0.25)\n", "q75 = count_cycles.quantile(0.75)\n", "mean = count_cycles.mean()\n", "sd = count_cycles.std()\n", "rng = (count_cycles.min(), count_cycles.max())\n", "mad = sts.median_abs_deviation(count_cycles)\n", "\n", "print(f\"Median = {med:.0f} cycles (IQR {q25:.0f}–{q75:.0f})\")\n", "print(f\"Mean = {mean:.1f} ± {sd:.1f} cycles\")\n", "print(f\"Range = {rng[0]}–{rng[1]} cycles\")" ] }, { "cell_type": "code", "execution_count": 41, "id": "f888cec3-37dc-45da-8fbb-fd7f1be89398", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Saved: violin_metrics_raw_points.png\n", "Saved: violin_metrics_standardized_points.png\n" ] } ], "source": [ "\n", "\n", "\n", "plt.rcParams.update({\n", " 'font.size': 15, # default font size\n", " 'font.weight': 'bold', # make all text bold\n", " 'axes.titleweight': 'bold',\n", " 'axes.labelweight': 'bold'})\n", "\n", "# helpers \n", "def zscore(series):\n", " s = pd.to_numeric(series, errors='coerce').astype(float)\n", " mu = np.nanmean(s)\n", " sd = np.nanstd(s, ddof=0)\n", " return (s - mu) / sd if sd and not np.isnan(sd) else s * np.nan\n", "\n", "def violin_with_points(ax, values, title=None, ylab=None, max_points=5000):\n", " # Force to float to avoid boolean dtype issues inside numpy.percentile\n", " v = pd.to_numeric(pd.Series(values), errors='coerce').astype(float).dropna().values\n", " if v.size == 0:\n", " ax.set_title(title or \"\", fontsize=9)\n", " ax.set_xticks([])\n", " ax.set_ylabel(ylab or \"\")\n", " ax.grid(True, axis='y', alpha=0.2)\n", " return\n", "\n", " # If all values are identical, violinplot can still draw; it’ll degenerate to a line.\n", " ax.violinplot([v], showmeans=False, showmedians=True, showextrema=False)\n", "\n", " # Jittered raw points (cap to keep files small)\n", " n_plot = min(len(v), max_points)\n", " vv = v[np.random.choice(len(v), size=n_plot, replace=False)] if len(v) > n_plot else v\n", " x = np.random.normal(loc=1.0, scale=0.03, size=len(vv))\n", " ax.plot(x, vv, 'o', markersize=2, alpha=0.5)\n", "\n", " ax.set_xticks([])\n", " ax.set_ylabel(ylab or \"\")\n", " ax.set_title(title or \"\", fontsize=9)\n", " ax.yaxis.set_major_locator(MaxNLocator(nbins=6))\n", " ax.grid(True, axis='y', alpha=0.2)\n", "\n", "def make_violin_grid(df, features, zscore_flag=False, filename=\"violin_metrics.png\", n_cols=5):\n", " # Keep numeric columns; drop booleans explicitly (they’re “numeric” but cause issues)\n", " numeric_feats = []\n", " for f in features:\n", " if f in df.columns and is_numeric_dtype(df[f]) and not is_bool_dtype(df[f]):\n", " numeric_feats.append(f)\n", " if not numeric_feats:\n", " raise ValueError(\"No numeric (non-boolean) metrics available to plot.\")\n", "\n", " # Prepare data (z-score if requested)\n", " data = df.copy()\n", " if zscore_flag:\n", " for f in numeric_feats:\n", " data[f] = zscore(data[f])\n", "\n", " n = len(numeric_feats)\n", " n_rows = int(np.ceil(n / n_cols))\n", " fig, axes = plt.subplots(n_rows, n_cols, figsize=(3.2*n_cols, 2.8*n_rows), squeeze=False)\n", " axes = axes.flatten()\n", "\n", " last_i = -1\n", " for i, feat in enumerate(numeric_feats):\n", " last_i = i\n", " ax = axes[i]\n", " vals = pd.to_numeric(data[feat], errors='coerce').astype(float) # <- cast to float here too\n", " violin_with_points(ax, vals, title=feat + (\" (z)\" if zscore_flag else \"\"), \n", " ylab=\"Z-score\" if zscore_flag else \"Raw\")\n", "\n", " # Remove any extra axes\n", " for j in range(last_i + 1, len(axes)):\n", " fig.delaxes(axes[j])\n", "\n", " fig.tight_layout()\n", " fig.savefig(filename, dpi=300, bbox_inches=\"tight\")\n", " plt.close(fig)\n", " print(f\"Saved: {filename}\")\n", "\n", "# ---------- run (two figures as requested) ----------\n", "\n", "make_violin_grid(feature_df, features_to_test_reduced, zscore_flag=False, filename=\"violin_metrics_raw_points.png\")\n", "make_violin_grid(feature_df, features_to_test_reduced, zscore_flag=True, filename=\"violin_metrics_standardized_points.png\")" ] }, { "cell_type": "markdown", "id": "8389bbec-5268-4d4a-9354-136b46de31c1", "metadata": { "tags": [] }, "source": [ "#### Removal of pregnancies cycles" ] }, { "cell_type": "code", "execution_count": 58, "id": "91143cbc-d85c-49f6-9f7c-4c1f858506f4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "num of cycles 5674\n", "num of participants 753\n" ] }, { "data": { "text/plain": [ "Age_Group\n", "1 4654\n", "2 1020\n", "Name: count, dtype: int64" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_non_conceptive_subset = feature_df.loc[~feature_df['conceptive']]\n", "\n", "print('num of cycles', len(df_non_conceptive_subset))\n", "print('num of participants', len(df_non_conceptive_subset.CODICE.unique()))\n", "\n", "df_non_conceptive_subset.Age_Group.value_counts()" ] }, { "cell_type": "markdown", "id": "da1b03f2-baf3-423e-9dc4-b1e8149a8f21", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "#### Age group comparison" ] }, { "cell_type": "code", "execution_count": 34, "id": "8dce8461-9032-478b-a4e9-fa140bfedf5d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Index(['Mean_Temperature', 'Standard_Deviation', 'Skewness', 'Kurtosis',\n", " 'Average_Slope', 'Period', 'Amplitude', 'Avg_Peak_to_Peak_Distance',\n", " 'M5', 'L5', 'Relative_Amplitude', 'Intra_monthly_Variability',\n", " 'Cosine_r2', 'Square_r2', 'Zigzag_r2', 'Sawtooth_Up_r2', 'Half_Sine_r2',\n", " 'Gaussian_r2', 'dur_fol', 'num_T_vals_fol', 'mean_fol', 'std_fol',\n", " 'amp_raw_fol', 'amp_smooth_fol', 'smoothness_fol', 'dur_lut',\n", " 'num_T_vals_lut', 'mean_lut', 'std_lut', 'amp_raw_lut',\n", " 'amp_smooth_lut', 'smoothness_lut', 'acrophase_day', 'nadir_day',\n", " 'acrophase_angle_lut', 'nadir_angle_fol', 'slope_m_to_nadir',\n", " 'slope_nadir_to_acrophase', 'slope_acrophase_to_m'],\n", " dtype='object')" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "features_to_test = features[~features.str.contains('_phase|_period')]\n", "features_to_test = features_to_test.drop(['Age_at_Menses', 'Age_Group', 'conceptive', 'fertile_intercourse'])\n", "features_to_test" ] }, { "cell_type": "code", "execution_count": 35, "id": "5bac8c6e-5821-460a-b018-f911eec8c215", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAAsTAAALEwEAmpwYAAARCklEQVR4nO3db4xc113G8e+DW9oobhuHtCsTBxyBQSQ1DeoqVCov1k0hFolIQAS5SpEjBZkXqdSKIOr0DX8kCwvRCqRQhGmqGgKsrJQqVqOCgttVqURI41JwnTSK1ZjISWSrbZJ2qyrI4ceLvRFTZ9Y79s6f3TPfj2TNnTNn5p795e6zJ2fu3ElVIUlqyw9NegCSpOEz3CWpQYa7JDXIcJekBhnuktSg1016AABXXHFFbd26dWz7+973vsell146tv2tJ9amP+uyPGvT3zjqcvTo0W9W1Vv7PbYmwn3r1q089thjY9vfwsICc3NzY9vfemJt+rMuy7M2/Y2jLkn+e7nHXJaRpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGrYlPqEqTsHXvQwP1O7n/phGPRBo+Z+6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQQOGe5GSSY0m+muSxru3yJA8neaq73dTT/54kJ5I8meTGUQ1ektTfhczcd1TVdVU1293fCxypqm3Ake4+Sa4BdgHXAjuBjyfZMMQxS5JWsJplmVuAg932QeDWnvb5qnq5qp4GTgDXr2I/kqQLlKpauVPyNPACUMBfVdWBJC9W1WU9fV6oqk1J7gUeqar7u/b7gM9V1QPnvOYeYA/AzMzMO+fn54f1M61ocXGRjRs3jm1/68k01ebYsy8N1G/7lW+ZqrpcKGvT3zjqsmPHjqM9qyk/4HUDvsa7q+q5JG8DHk7y9fP0TZ+21/wFqaoDwAGA2dnZmpubG3Aoq7ewsMA497eeTFNt7tj70ED9Tt4+N1V1uVDWpr9J12WgZZmqeq67PQN8hqVlltNJNgN0t2e67qeAq3qevgV4blgDliStbMVwT3Jpkje9ug38EvA14DCwu+u2G3iw2z4M7EryhiRXA9uAR4c9cEnS8gZZlpkBPpPk1f5/X1X/lOTLwKEkdwLPALcBVNXxJIeAx4GzwF1V9cpIRi9J6mvFcK+qbwDv6NP+LeCGZZ6zD9i36tFJki6Kn1CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGvm/QApGHauvehSQ9BWhOcuUtSgwx3SWqQ4S5JDXLNXeuCa+nShRl45p5kQ5L/SPLZ7v7lSR5O8lR3u6mn7z1JTiR5MsmNoxi4JGl5F7Is80HgiZ77e4EjVbUNONLdJ8k1wC7gWmAn8PEkG4YzXEnSIAYK9yRbgJuAT/Q03wIc7LYPArf2tM9X1ctV9TRwArh+KKOVJA0kVbVyp+QB4I+BNwG/W1U3J3mxqi7r6fNCVW1Kci/wSFXd37XfB3yuqh445zX3AHsAZmZm3jk/Pz+sn2lFi4uLbNy4cWz7W0/Wam2OPfvSxPa9/cq3rNm6rAXWpr9x1GXHjh1Hq2q232MrvqGa5GbgTFUdTTI3wP7Sp+01f0Gq6gBwAGB2drbm5gZ56eFYWFhgnPtbT9Zqbe6Y4BuqJ2+fW7N1WQusTX+TrssgZ8u8G/iVJL8MvBF4c5L7gdNJNlfV80k2A2e6/qeAq3qevwV4bpiDliSd34pr7lV1T1VtqaqtLL1R+vmqej9wGNjdddsNPNhtHwZ2JXlDkquBbcCjQx+5JGlZqznPfT9wKMmdwDPAbQBVdTzJIeBx4CxwV1W9suqRSpIGdkHhXlULwEK3/S3ghmX67QP2rXJskqSL5OUHJKlBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkF/WIa1g696HuHv72RWvb3Ny/01jGpG0MsNdGpJBvy3KPwIaB5dlJKlBhrskNchwl6QGGe6S1CDDXZIaZLhLUoM8FVITNejpg5IujDN3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQSuGe5I3Jnk0yX8mOZ7kD7v2y5M8nOSp7nZTz3PuSXIiyZNJbhzlDyBJeq1BZu4vA++pqncA1wE7k7wL2AscqaptwJHuPkmuAXYB1wI7gY8n2TCCsUuSlrFiuNeSxe7u67t/BdwCHOzaDwK3dtu3APNV9XJVPQ2cAK4f5qAlSeeXqlq509LM+yjwk8BfVNWHk7xYVZf19HmhqjYluRd4pKru79rvAz5XVQ+c85p7gD0AMzMz75yfnx/Wz7SixcVFNm7cOLb9rSfjrs2xZ18a275WY+YSOP394bzW9ivfMpwXWiP8fepvHHXZsWPH0aqa7ffYQNdzr6pXgOuSXAZ8Jsnbz9M9/V6iz2seAA4AzM7O1tzc3CBDGYqFhQXGub/1ZNy1uWOdXM/97u1n+eix4Xz9wcnb54byOmuFv0/9TbouF3S2TFW9CCywtJZ+OslmgO72TNftFHBVz9O2AM+tdqCSpMENcrbMW7sZO0kuAd4LfB04DOzuuu0GHuy2DwO7krwhydXANuDRIY9bknQeg/x/5mbgYLfu/kPAoar6bJJ/Aw4luRN4BrgNoKqOJzkEPA6cBe7qlnUkMfhXC57cf9OIR6KWrRjuVfVfwM/1af8WcMMyz9kH7Fv16CRJF8VPqEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBg3ne8MkDZ3XfddqGO4aiUGDSdJouCwjSQ0y3CWpQYa7JDXIcJekBhnuktQgz5aRdNG27n2Iu7ef5Y4Vzo7ydM3xc+YuSQ1y5q4L4vnr0vrgzF2SGmS4S1KDDHdJapDhLkkNWjHck1yV5AtJnkhyPMkHu/bLkzyc5KnudlPPc+5JciLJk0luHOUPIEl6rUFm7meBu6vqZ4B3AXcluQbYCxypqm3Ake4+3WO7gGuBncDHk2wYxeAlSf2tGO5V9XxVfaXb/i7wBHAlcAtwsOt2ELi1274FmK+ql6vqaeAEcP2Qxy1JOo9U1eCdk63AF4G3A89U1WU9j71QVZuS3As8UlX3d+33AZ+rqgfOea09wB6AmZmZd87Pz6/yRxnc4uIiGzduHNv+1pOVanPs2ZfGOJq1Y+YSOP39SY+iv+1XvmVi+z727EsD1WaSY5yUceTMjh07jlbVbL/HBv4QU5KNwKeBD1XVd5Is27VP22v+glTVAeAAwOzsbM3NzQ06lFVbWFhgnPtbT1aqzUofM2/V3dvP8tFja/Mzfydvn5vYvu/oLj+wUm0mOcZJmXTODHS0Jnk9S8H+d1X1j13z6SSbq+r5JJuBM137KeCqnqdvAZ4b1oAljZ6fRF7/Vgz3LE3R7wOeqKqP9Tx0GNgN7O9uH+xp//skHwN+FNgGPDrMQUtaX/w+2PEbZOb+buA3gWNJvtq1fYSlUD+U5E7gGeA2gKo6nuQQ8DhLZ9rcVVWvDHvgkqTlrRjuVfUl+q+jA9ywzHP2AftWMS5JU8gZ/vCszXeIJA3MQFQ/Xn5AkhrkzF2aEp4BM12cuUtSg5y5N+5CZ2uDfB+mpLXPmbskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDvCqkpHXHb59amTN3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoNWDPckn0xyJsnXetouT/Jwkqe62009j92T5ESSJ5PcOKqBS5KWN8gnVD8F3Av8TU/bXuBIVe1Psre7/+Ek1wC7gGuBHwX+JclPVdUrwx32dBv003mSpteKM/eq+iLw7XOabwEOdtsHgVt72uer6uWqeho4AVw/nKFKkgZ1sdeWmamq5wGq6vkkb+varwQe6el3qmt7jSR7gD0AMzMzLCwsXORQLtzi4uJY9zdsd28/O7LXnrlktK+/XlmX5a3l2kzy93zSOTPsC4elT1v161hVB4ADALOzszU3NzfkoSxvYWGBce5v2O4Y4bLM3dvP8tFjXk/uXNZleWu5Nidvnxuo34UsdQ56MbJJ58zF/hc5nWRzN2vfDJzp2k8BV/X02wI8t5oBStJaMugfgk/tvHTEIzm/iz0V8jCwu9veDTzY074ryRuSXA1sAx5d3RAlSRdqxZl7kn8A5oArkpwCfh/YDxxKcifwDHAbQFUdT3IIeBw4C9zlmTKSNH4rhntVvW+Zh25Ypv8+YN9qBiVJWp21+S6IJA3BNH8mxMsPSFKDDHdJapDhLkkNMtwlqUG+obqGTPObP5KGy5m7JDXIcJekBhnuktQgw12SGmS4S1KDPFtmDDwLRtK4OXOXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDfJUyFXwFEdJa5Uzd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatDIPqGaZCfw58AG4BNVtX9U+7rQT4revf0sd/jpUkkNG0m4J9kA/AXwi8Ap4MtJDlfV46PYnyStNceefWmgSeTJ/TeNZP+jWpa5HjhRVd+oqv8B5oFbRrQvSdI5UlXDf9Hk14GdVfVb3f3fBH6+qj7Q02cPsKe7+9PAk0MfyPKuAL45xv2tJ9amP+uyPGvT3zjq8uNV9dZ+D4xqzT192n7gr0hVHQAOjGj/55XksaqancS+1zpr0591WZ616W/SdRnVsswp4Kqe+1uA50a0L0nSOUYV7l8GtiW5OskPA7uAwyPalyTpHCNZlqmqs0k+APwzS6dCfrKqjo9iXxdpIstB64S16c+6LM/a9DfRuozkDVVJ0mT5CVVJapDhLkkNaj7ck1yV5AtJnkhyPMkHu/bLkzyc5KnudtOkxzpO56nLHyR5NslXu3+/POmxjlOSNyZ5NMl/dnX5w659qo8XOG9tpvqYeVWSDUn+I8lnu/sTPWaaX3NPshnYXFVfSfIm4ChwK3AH8O2q2p9kL7Cpqj48uZGO13nq8hvAYlX96STHNylJAlxaVYtJXg98Cfgg8GtM8fEC563NTqb4mHlVkt8BZoE3V9XNSf6ECR4zzc/cq+r5qvpKt/1d4AngSpYuh3Cw63aQpWCbGuepy1SrJYvd3dd3/4opP17gvLWZekm2ADcBn+hpnugx03y490qyFfg54N+Bmap6HpaCDnjbBIc2UefUBeADSf4rySendPlhQ5KvAmeAh6vK46WzTG1gyo8Z4M+A3wP+t6dtosfM1IR7ko3Ap4EPVdV3Jj2etaJPXf4S+AngOuB54KOTG91kVNUrVXUdS5+svj7J2yc8pDVjmdpM9TGT5GbgTFUdnfRYek1FuHfrg58G/q6q/rFrPt2tO7+6/nxmUuOblH51qarT3S/w/wJ/zdIVPqdSVb0ILLC0pjz1x0uv3tp4zPBu4FeSnGTpCrjvSXI/Ez5mmg/37k2g+4AnqupjPQ8dBnZ327uBB8c9tklari6vHoydXwW+Nu6xTVKStya5rNu+BHgv8HWm/HiB5Wsz7cdMVd1TVVuqaitLl1r5fFW9nwkfM9NwtswvAP8KHOP/18M+wtL68iHgx4BngNuq6tsTGeQEnKcu72Ppf68LOAn89qvrhtMgyc+y9ObXBpYmP4eq6o+S/AhTfLzAeWvzt0zxMdMryRzwu93ZMhM9ZpoPd0maRs0vy0jSNDLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoP+D3DsagxmAdH9AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "df_non_conceptive_subset['Age_at_Menses'].hist(bins=30)" ] }, { "cell_type": "code", "execution_count": 36, "id": "5c1d775c-284d-4b8b-ac84-b1c5262466d3", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Mean_TemperatureStandard_DeviationSkewnessKurtosisAverage_SlopePeriodAmplitudeAvg_Peak_to_Peak_DistanceM5L5Relative_AmplitudeIntra_monthly_VariabilityCosine_r2Cosine_periodCosine_phaseSquare_r2Square_periodSquare_phaseZigzag_r2Zigzag_periodZigzag_phaseSawtooth_Up_r2Sawtooth_Up_periodSawtooth_Up_phaseHalf_Sine_r2Half_Sine_periodHalf_Sine_phaseGaussian_r2Gaussian_periodGaussian_phasedur_folnum_T_vals_folmean_folstd_folamp_raw_folamp_smooth_folsmoothness_foldur_lutnum_T_vals_lutmean_lutstd_lutamp_raw_lutamp_smooth_lutsmoothness_lutacrophase_daynadir_dayacrophase_angle_lutnadir_angle_folslope_m_to_nadirslope_nadir_to_acrophaseslope_acrophase_to_mCODICEAge_at_MensesQUALIFI
meanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstdmeanstd
Age_Group
136.6930.1730.1970.0480.190.498-0.8680.6150.0740.02325.9893.4460.5230.1474.3301.45536.9470.18236.4580.1850.0070.0020.6660.4190.6700.19528.9345.23812.8929.2140.6790.16725.7224.23210.3218.6450.6510.19028.3184.50314.878.3700.6440.14931.2445.21617.7678.5680.7140.17529.6025.02913.7468.7290.6370.1632.7724.26916.68211.65016.2053.28211.8844.01036.5590.1770.1240.0560.3980.2030.2490.1460.1140.0549.7663.5109.3813.41636.8620.1740.1320.0570.3930.1950.3460.1480.1060.05323.4743.60817.6793.671138.73739.167196.99619.5600.0020.0170.060.039-0.0290.0303.173348e+072.330931e+0729.0803.0181.00.0
236.7480.1660.1910.0460.140.513-0.9100.5780.0700.02224.6973.7160.5040.1404.2591.44336.9860.17436.5210.1720.0060.0020.6090.3650.7090.16627.9995.15112.0438.7650.7170.14924.9053.6149.2087.5540.6850.16327.6714.33014.188.1420.6810.13531.4795.22618.3038.6390.7470.15029.2385.10213.3199.2300.6670.1532.7834.28515.34011.39115.1712.97011.2383.61636.6140.1720.1170.0550.3740.1950.2330.1460.1070.0539.5523.6829.3053.61936.9110.1620.1220.0570.3560.1910.3210.1450.1000.05221.9893.47216.6863.441134.16040.444198.75123.9380.0040.0190.060.040-0.0260.0273.639942e+072.220592e+0736.9821.4381.00.0
\n", "
" ], "text/plain": [ " Mean_Temperature Standard_Deviation Skewness Kurtosis Average_Slope Period Amplitude Avg_Peak_to_Peak_Distance M5 L5 Relative_Amplitude Intra_monthly_Variability Cosine_r2 Cosine_period Cosine_phase Square_r2 Square_period Square_phase Zigzag_r2 Zigzag_period Zigzag_phase Sawtooth_Up_r2 Sawtooth_Up_period Sawtooth_Up_phase Half_Sine_r2 Half_Sine_period Half_Sine_phase Gaussian_r2 Gaussian_period Gaussian_phase dur_fol num_T_vals_fol mean_fol std_fol amp_raw_fol amp_smooth_fol smoothness_fol dur_lut num_T_vals_lut mean_lut std_lut amp_raw_lut amp_smooth_lut smoothness_lut acrophase_day nadir_day acrophase_angle_lut nadir_angle_fol slope_m_to_nadir slope_nadir_to_acrophase slope_acrophase_to_m CODICE Age_at_Menses QUALIFI \n", " mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std mean std\n", "Age_Group \n", "1 36.693 0.173 0.197 0.048 0.19 0.498 -0.868 0.615 0.074 0.023 25.989 3.446 0.523 0.147 4.330 1.455 36.947 0.182 36.458 0.185 0.007 0.002 0.666 0.419 0.670 0.195 28.934 5.238 12.892 9.214 0.679 0.167 25.722 4.232 10.321 8.645 0.651 0.190 28.318 4.503 14.87 8.370 0.644 0.149 31.244 5.216 17.767 8.568 0.714 0.175 29.602 5.029 13.746 8.729 0.637 0.16 32.772 4.269 16.682 11.650 16.205 3.282 11.884 4.010 36.559 0.177 0.124 0.056 0.398 0.203 0.249 0.146 0.114 0.054 9.766 3.510 9.381 3.416 36.862 0.174 0.132 0.057 0.393 0.195 0.346 0.148 0.106 0.053 23.474 3.608 17.679 3.671 138.737 39.167 196.996 19.560 0.002 0.017 0.06 0.039 -0.029 0.030 3.173348e+07 2.330931e+07 29.080 3.018 1.0 0.0\n", "2 36.748 0.166 0.191 0.046 0.14 0.513 -0.910 0.578 0.070 0.022 24.697 3.716 0.504 0.140 4.259 1.443 36.986 0.174 36.521 0.172 0.006 0.002 0.609 0.365 0.709 0.166 27.999 5.151 12.043 8.765 0.717 0.149 24.905 3.614 9.208 7.554 0.685 0.163 27.671 4.330 14.18 8.142 0.681 0.135 31.479 5.226 18.303 8.639 0.747 0.150 29.238 5.102 13.319 9.230 0.667 0.15 32.783 4.285 15.340 11.391 15.171 2.970 11.238 3.616 36.614 0.172 0.117 0.055 0.374 0.195 0.233 0.146 0.107 0.053 9.552 3.682 9.305 3.619 36.911 0.162 0.122 0.057 0.356 0.191 0.321 0.145 0.100 0.052 21.989 3.472 16.686 3.441 134.160 40.444 198.751 23.938 0.004 0.019 0.06 0.040 -0.026 0.027 3.639942e+07 2.220592e+07 36.982 1.438 1.0 0.0" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd.set_option('display.max_columns', None) # Show all columns\n", "pd.set_option('display.expand_frame_repr', False) # Prevent column wrapping\n", "\n", "summary_stats = df_non_conceptive_subset.groupby(\"Age_Group\").describe().round(3)\n", "summary_stats = summary_stats.loc[:, (slice(None), ['mean', 'std'])]\n", "summary_stats" ] }, { "cell_type": "code", "execution_count": 37, "id": "4342aa39-6107-4ba8-ab98-f346fb1a37d1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "bonferroni threshold value 0.0013\n", "number of comparisons: 39\n" ] } ], "source": [ "from scipy.stats import chi2_contingency\n", "p_values = {}\n", "# Loop through each feature\n", "for feature in features:\n", " if is_categorical_dtype(df_non_conceptive_subset[feature]):\n", " df_non_conceptive_subset[feature] = df_non_conceptive_subset[feature].astype(\"category\")\n", " ct = pd.crosstab(feature_df[\"Age_Group\"], feature_df[feature])\n", "\n", " # 3) Run Chi‐square test of association\n", " chi2, p_value, dof, expected = chi2_contingency(ct)\n", " test_used = \"Chi2\"\n", "\n", " else:\n", " \n", " groups = [group.dropna().values for _, group in df_non_conceptive_subset.groupby(\"Age_Group\")[feature]]\n", " # Test for normality (Shapiro-Wilk test on the first group as an approximation)\n", " stat, p_normality = stats.shapiro(groups[0]) if len(groups[0]) >= 3 else (None, 1.0)\n", " if p_normality > 0.05: \n", " stat, p_value = stats.f_oneway(*groups)\n", " test_used = \"ANOVA\"\n", " else:\n", " stat, p_value = stats.kruskal(*groups)\n", " test_used = \"Kruskal-Wallis\"\n", " \n", " p_values[feature] = {\"Test\": test_used, \"p-value\": p_value}\n", "\n", "# Convert to DataFrame\n", "p_values_df = pd.DataFrame(p_values).T\n", "\n", "#Positive result ?\n", "alpha=0.05\n", "#p_values_df['positive_test'] = p_values_df['p-value']\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
 p-valuepositive_test_bonferronimean_AgeGroup_1mean_AgeGroup_2std_AgeGroup_1std_AgeGroup_2
Feature      
Period7.17e-31True25.9924.703.453.72
acrophase_day3.73e-29True23.4721.993.613.47
dur_fol3.32e-20True16.2015.173.282.97
Mean_Temperature6.94e-20True36.6936.750.170.17
L51.90e-19True36.4636.520.180.17
mean_lut3.07e-16True36.8636.910.170.16
nadir_day3.44e-16True17.6816.693.673.44
mean_fol1.49e-15True36.5636.610.180.17
Sawtooth_Up_r26.87e-14True0.640.680.150.14
Square_r21.46e-12True0.680.720.170.15
M58.45e-11True36.9536.990.180.17
Relative_Amplitude7.39e-09True0.010.010.000.00
nadir_angle_fol1.98e-08True197.00198.7519.5623.94
amp_raw_lut2.21e-08True0.390.360.200.19
Average_Slope2.93e-08True0.070.070.020.02
std_lut4.44e-08True0.130.120.060.06
Half_Sine_r26.91e-08True0.710.750.170.15
Cosine_r21.07e-07True0.670.710.200.17
Gaussian_r21.74e-07True0.640.670.160.15
amp_smooth_lut5.37e-07True0.350.320.150.14
Zigzag_r23.09e-06True0.650.690.190.16
Amplitude8.30e-06True0.520.500.150.14
smoothness_fol2.98e-05True0.110.110.050.05
std_fol1.33e-04True0.120.120.060.06
amp_smooth_fol1.57e-04True0.250.230.150.15
num_T_vals_fol1.77e-04True11.8811.244.013.62
amp_raw_fol2.84e-04True0.400.370.200.20
Intra_monthly_Variability4.20e-04True0.670.610.420.36
Skewness5.00e-04True0.190.140.500.51
smoothness_lut9.44e-04True0.110.100.050.05
Standard_Deviation0.001True0.200.190.050.05
acrophase_angle_lut0.002False138.74134.1639.1740.44
slope_acrophase_to_m0.018False-0.03-0.030.030.03
Kurtosis0.049False-0.87-0.910.610.58
slope_m_to_nadir0.107False0.000.000.020.02
dur_lut0.149False9.779.553.513.68
Avg_Peak_to_Peak_Distance0.188False4.334.261.461.44
num_T_vals_lut0.722False9.389.303.423.62
slope_nadir_to_acrophase0.929False0.060.060.040.04
\n" ], "text/plain": [ "" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "# Apply formatting styles\n", "styled = final_summary.sort_values(by='p-value', ascending=True).style.format({\n", " 'Test': lambda x: x,\n", " 'p-value': fmt_small_conditional,\n", " 'positive_test_bonferroni': fmt_boolean,\n", " 'count_AgeGroup_1': fmt_integer,\n", " 'mean_AgeGroup_1': fmt_small,\n", " 'std_AgeGroup_1': fmt_small,\n", " 'count_AgeGroup_2': fmt_integer,\n", " 'mean_AgeGroup_2': fmt_small,\n", " 'std_AgeGroup_2': fmt_small,})\n", "\n", "# 3) display it\n", "styled\n" ] }, { "cell_type": "markdown", "id": "b88931ac-0fd3-4bc5-a69b-f47bfeff7344", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "#### LMM for age at menses for every feature" ] }, { "cell_type": "code", "execution_count": 52, "id": "969bf619-af65-47ba-9183-a17ffbad7bf3", "metadata": {}, "outputs": [], "source": [ "features_to_test_reduced = ['L5', 'M5', 'mean_fol', 'mean_lut', 'Mean_Temperature',\\\n", "'acrophase_angle_lut', 'acrophase_day', 'dur_fol', 'dur_lut', 'nadir_angle_fol', 'nadir_day', 'Period', \\\n", "'amp_smooth_fol', 'amp_smooth_lut', 'Amplitude', 'Average_Slope', 'Kurtosis', 'Skewness', 'slope_acrophase_to_m', \\\n", "'slope_m_to_nadir', 'slope_nadir_to_acrophase', 'smoothness_fol', 'smoothness_lut', 'Standard_Deviation', 'std_fol', 'std_lut', \\\n", "'Cosine_r2', 'Gaussian_r2', 'Half_Sine_r2', 'Sawtooth_Up_r2', 'Square_r2', 'Zigzag_r2']" ] }, { "cell_type": "code", "execution_count": 53, "id": "796a580d-5a7a-4fe4-8f09-f1095613651e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "753" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(df_non_conceptive_subset.groupby(['CODICE'])['Age_at_Menses'].count())" ] }, { "cell_type": "code", "execution_count": 60, "id": "1aae2970-e98e-4960-ad70-508f4e1ee90d", "metadata": {}, "outputs": [], "source": [ "all_features = list(set(df_non_conceptive_subset.columns.to_list()) - set(['Age_at_Menses', 'Age_Group', 'Window_Start', 'Window_End', 'DONNA', 'DATA', 'QUALIFI', 'conceptive_group']))\n" ] }, { "cell_type": "code", "execution_count": 61, "id": "b66ebb5e-0df7-4042-b91d-2e9ede6ea7a1", "metadata": {}, "outputs": [], "source": [ "group_keys = ['CODICE', 'Age_at_Menses']\n", "feature_cols = [c for c in all_features if c not in group_keys]\n", "\n", "feature_df_grouped = (df_non_conceptive_subset.groupby(group_keys)[feature_cols].mean(numeric_only=True).reset_index())" ] }, { "cell_type": "code", "execution_count": 62, "id": "e1f874b7-8adf-4c3a-ba5a-a50fdbd3b591", "metadata": {}, "outputs": [], "source": [ "three_way_table = feature_df_grouped.pivot(index='CODICE', columns='Age_at_Menses', values=all_features)\n", "\n", "# Convert to a NumPy tensor (Patients x Age_Years x Features)\n", "three_way_tensor = np.array([three_way_table[feature].values for feature in all_features])\n", "three_way_tensor = np.moveaxis(three_way_tensor, 0, -1) " ] }, { "cell_type": "markdown", "id": "9d5a158b-3af6-4597-9514-e901ecb45ed9", "metadata": {}, "source": [ "##### Recentering age to mean" ] }, { "cell_type": "code", "execution_count": 63, "id": "cd343c15-18be-49b5-9c13-8836bb8a64b7", "metadata": {}, "outputs": [], "source": [ "# ------------------------------------------------------------------\n", "# 0 · Prep containers\n", "# ------------------------------------------------------------------\n", "coeffs, pvals, ci_low, ci_high, icc = {}, {}, {}, {}, {}\n", "\n", "# ------------------------------------------------------------------\n", "# 1 · Loop over every feature\n", "# ------------------------------------------------------------------\n", "for feature in features_to_test_reduced:\n", " feature_data = three_way_table[feature]\n", "\n", " df_long = (feature_data.reset_index()\n", " .melt(id_vars=\"CODICE\",\n", " var_name=\"Age_at_Menses\",\n", " value_name=\"Value\")\n", " .merge(site_lookup, on=\"CODICE\", how=\"left\") # <- adds 'site'\n", " .dropna(subset=[\"Value\"]) )\n", "\n", " # make sure it's treated as categorical (creates dummies automatically later)\n", " df_long[\"site\"] = df_long[\"site\"].astype(\"category\")\n", "\n", "\n", " # ----- treat age as numeric & centre at its mean ------------------\n", " df_long['Age_at_Menses'] = (pd.to_numeric(df_long['Age_at_Menses'], errors='coerce') )\n", " mu_age = df_long['Age_at_Menses'].mean()\n", " df_long['Age_c'] = df_long['Age_at_Menses'] - mu_age # centred\n", "\n", " # ----- fit random-intercept model; skip feature if it fails -------\n", " try:\n", " #Choose model to run\n", "\n", " ### Model using participants as random effects:\n", " model = smf.mixedlm(\"Value ~ Age_c\", df_long, groups=df_long[\"CODICE\"])\n", " res = model.fit(reml=False, warn_convergence=False)\n", "\n", " coeffs[feature] = res.params['Age_c']\n", " pvals[feature] = res.pvalues['Age_c']\n", " ci_low[feature], ci_high[feature] = res.conf_int().loc['Age_c']\n", "\n", " # 2) extract variance components\n", " var_between = res.cov_re.iloc[0, 0] # random‐intercept variance\n", " var_within = res.scale # residual variance\n", "\n", " # 3) compute and store ICC\n", " icc[feature] = var_between / (var_between + var_within)\n", "\n", " except Exception as err:\n", " warnings.warn(f\"{feature}: model failed - {err}\")\n", " coeffs[feature] = pvals[feature] = ci_low[feature] = ci_high[feature] = np.nan\n" ] }, { "cell_type": "code", "execution_count": 64, "id": "b2a3d10b-9ce5-498f-b7d9-3f7940c73878", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
 CoefficientCI_low_95%CI_high_95%ICCP_valueBonferroniTest
Period-0.194-0.239-0.1490.6132.11e-17True
acrophase_day-0.192-0.238-0.1460.4412.52e-16True
Gaussian_r20.0050.0030.0070.3102.50e-07True
dur_fol-0.118-0.163-0.0720.5593.84e-07True
Sawtooth_Up_r20.0040.0020.0060.3252.47e-06True
nadir_day-0.109-0.158-0.0610.4608.66e-06True
L50.0040.0020.0060.8094.59e-04True
Square_r20.0030.0010.0050.3155.87e-04True
Zigzag_r20.0040.0010.0060.2730.001True
amp_smooth_lut-0.003-0.004-0.0010.2530.001True
Half_Sine_r20.0030.0010.0050.3100.001True
Cosine_r20.0040.0010.0060.2920.001True
slope_m_to_nadir2.56e-048.47e-054.28e-040.1530.003False
dur_lut-0.067-0.115-0.0200.5490.005False
std_lut-8.52e-04-0.001-2.50e-040.1840.006False
Average_Slope-3.68e-04-6.48e-04-8.86e-050.3500.010False
smoothness_fol-8.07e-04-0.001-1.78e-040.2990.012False
Amplitude-0.002-0.004-4.78e-040.2970.012False
amp_smooth_fol-0.002-0.004-2.38e-040.2590.025False
acrophase_angle_lut-0.471-0.895-0.0470.2240.030False
smoothness_lut-6.00e-04-0.001-2.66e-050.2150.040False
nadir_angle_fol0.215-0.0020.4320.1980.052False
Standard_Deviation-4.98e-04-0.0015.72e-050.2920.079False
std_fol-5.27e-04-0.0011.08e-040.2500.104False
slope_acrophase_to_m2.32e-04-7.71e-055.42e-040.1060.141False
Kurtosis-0.004-0.0100.0020.1370.219False
M5-0.001-0.0040.0010.7220.301False
Skewness-0.002-0.0080.0030.1970.352False
mean_fol9.28e-04-0.0010.0030.8050.432False
Mean_Temperature8.56e-04-0.0010.0030.8270.439False
slope_nadir_to_acrophase1.19e-04-2.79e-045.17e-040.1590.558False
mean_lut-3.57e-05-0.0020.0020.7760.976False
\n" ], "text/plain": [ "" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# ------------------------------------------------------------------\n", "# 2 · Put everything in a tidy frame\n", "# ------------------------------------------------------------------\n", "mixed_model_summary_table6 = (\n", " pd.DataFrame({\n", " 'Coefficient' : coeffs,\n", " 'CI_low_95%' : ci_low,\n", " 'CI_high_95%' : ci_high,\n", " 'ICC' : icc,\n", " 'P_value' : pvals,})\n", " .dropna(subset=['P_value']) # drop failed models\n", " .sort_values('P_value'))\n", "\n", "# ------------------------------------------------------------------\n", "# 3 · Multiple-testing correction\n", "# ------------------------------------------------------------------\n", "alpha = 0.05\n", "# --- Bonferroni\n", "mixed_model_summary_table6['BonferroniTest'] = (mixed_model_summary_table6['P_value'] < alpha / len(mixed_model_summary_table6))\n", "# --- Benjamini–Hochberg FDR\n", "#_, qvals, _, _ = multipletests(mixed_model_summary['P_value'], method='fdr_bh', alpha=alpha)\n", "#mixed_model_summary['Q_value'] = qvals\n", "#mixed_model_summary['sig_FDR_0.05'] = qvals < alpha\n", "\n", "# 2) apply different formats per column\n", "styled_table_section_6 = mixed_model_summary_table6.style.format({\n", " 'Coefficient': fmt_small_conditional, \n", " 'CI_low_95%': fmt_small_conditional, \n", " 'CI_high_95%': fmt_small_conditional, \n", " 'ICC': fmt_small_conditional,\n", " 'P_value': fmt_small_conditional, \n", " 'BonferroniTest': fmt_boolean,\n", " #'Q_value': fmt_small\n", "})\n", "\n", "# 3) display it\n", "styled_table_section_6" ] }, { "cell_type": "markdown", "id": "12350ad6-9f39-45ec-a112-8c4dc4ad562c", "metadata": { "tags": [] }, "source": [ "#### intra subject variability" ] }, { "cell_type": "code", "execution_count": 65, "id": "e382592b-5470-4ef3-a5f7-831c3cd3996b", "metadata": {}, "outputs": [], "source": [ "WINDOW_CYCLES = 6 \n", "\n", "# what statistic do we want?\n", "USE_CV = True # True ⇒ coefficient of variation; False ⇒ plain std\n", "MIN_VALID = 4 # require full window to call the std “valid”\n" ] }, { "cell_type": "code", "execution_count": 66, "id": "43ecf291-cb02-46d0-b362-abf11e77ad5d", "metadata": {}, "outputs": [], "source": [ "import statsmodels.formula.api as smf\n", "### Table of STDs\n", "\n", "variability_tables = {} # feature ➜ wide table (CODICE × Age) of within–window SD/CV\n", "\n", "for feature in features_to_test_reduced:\n", " # A. long table (same as before)\n", " df_long = (three_way_table[feature].reset_index() \n", " .melt(id_vars='CODICE', var_name='Age_at_Menses', value_name='Value')\n", " .dropna(subset=['Value'])\n", " .assign(Age_at_Menses=lambda d: pd.to_numeric(d['Age_at_Menses'], errors='coerce')).sort_values(['CODICE', 'Age_at_Menses']) )\n", "\n", " # B. rolling variability within each subject\n", " def _rolling_stat(group):\n", " vals = group['Value']\n", " roll = vals.rolling(window=WINDOW_CYCLES, min_periods=MIN_VALID)\n", " if USE_CV:\n", " return roll.std() / roll.mean()\n", " else:\n", " return roll.std()\n", "\n", " df_long['RollVar'] = (df_long.groupby('CODICE', group_keys=False).apply(_rolling_stat))\n", "\n", " # keep only rows where rolling SD/CV is defined\n", " df_var = (df_long.dropna(subset=['RollVar']).pivot(index='CODICE', columns='Age_at_Menses',values='RollVar'))\n", " variability_tables[feature] = df_var\n" ] }, { "cell_type": "code", "execution_count": 67, "id": "1bbe773d-1dad-4d74-a16e-06e95105c90a", "metadata": {}, "outputs": [], "source": [ "coef_var, p_var, icc_var = {}, {}, {}\n", "\n", "for feature, table_var in variability_tables.items():\n", "\n", " df_long = (table_var\n", " .reset_index()\n", " .melt(id_vars='CODICE',var_name='Age_at_Menses',value_name='VarValue')\n", " .dropna(subset=['VarValue'])\n", " .assign(Age_at_Menses=lambda d: pd.to_numeric(d['Age_at_Menses'], errors='coerce')))\n", "\n", " # centre age\n", " mu_age = df_long['Age_at_Menses'].mean()\n", " df_long['Age_c'] = df_long['Age_at_Menses'] - mu_age\n", "\n", " try:\n", " mod = smf.mixedlm(\"VarValue ~ Age_c\", df_long, groups=df_long[\"CODICE\"])\n", " res = mod.fit(reml=False, warn_convergence=False)\n", "\n", " # 1) store your slope and p-value as before\n", " coef_var[feature] = res.params['Age_c']\n", " p_var[feature] = res.pvalues['Age_c']\n", "\n", " # 2) extract variance components\n", " var_between = res.cov_re.iloc[0, 0] # random‐intercept variance\n", " var_within = res.scale # residual variance\n", "\n", " # 3) compute and store ICC\n", " icc_var[feature] = var_between / (var_between + var_within)\n", "\n", " except Exception as e:\n", " warnings.warn(f\"{feature}: model failed ({e})\")\n", " coef_var[feature] = p_var[feature] = icc_var[feature] = np.nan\n" ] }, { "cell_type": "code", "execution_count": 70, "id": "c25afa39-c5b7-441e-aeaf-6a6160ed2ad0", "metadata": {}, "outputs": [], "source": [ "alpha = 0.05\n", "\n", "summary_var_table = (pd.DataFrame({'Coefficient_SDvsAge' : coef_var, 'P_value': p_var, 'ICC': icc_var }).dropna(subset=['P_value']).sort_values('P_value'))\n", "\n", "# Bonferroni\n", "summary_var_table['BonferroniTest'] = (summary_var_table['P_value'] < alpha / len(summary_var_table))" ] }, { "cell_type": "code", "execution_count": 71, "id": "3f3d4a21-5473-4178-8582-f49fa21e48a7", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
 Coefficient_SDvsAgeICCP_valueBonferroniTest
Period0.0040.76432.65e-18True
std_lut0.0110.52491.52e-10True
amp_smooth_lut0.0080.57841.71e-07True
smoothness_lut0.0090.51901.95e-07True
dur_lut0.0070.68114.97e-07True
Square_r2-0.0050.66381.10e-05True
Average_Slope0.0040.57165.04e-05True
Zigzag_r2-0.0060.65615.08e-05True
Cosine_r2-0.0060.66896.12e-05True
smoothness_fol0.0060.59821.15e-04True
std_fol0.0060.56671.29e-04True
Half_Sine_r2-0.0040.68792.72e-04True
amp_smooth_fol0.0060.56624.13e-04True
Sawtooth_Up_r2-0.0030.66890.007False
Gaussian_r2-0.0030.63780.008False
slope_m_to_nadir0.5940.00480.034False
Mean_Temperature-1.12e-050.53620.141False
acrophase_day5.41e-040.60060.181False
mean_fol-1.03e-050.59690.252False
slope_acrophase_to_m-0.0050.66960.262False
dur_fol-5.87e-040.54330.290False
nadir_angle_fol5.59e-040.52690.352False
Standard_Deviation7.10e-040.57340.390False
Amplitude6.42e-040.62770.515False
Skewness-0.1030.00000.586False
Kurtosis-0.0290.10890.641False
mean_lut-3.59e-060.51240.688False
M5-4.30e-060.58080.723False
nadir_day-2.28e-040.57580.726False
acrophase_angle_lut-3.71e-040.61510.745False
slope_nadir_to_acrophase-0.0090.01790.786False
L5-2.09e-060.53460.802False
\n" ], "text/plain": [ "" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "summary_var_table = summary_var_table[['Coefficient_SDvsAge', 'ICC', 'P_value', 'BonferroniTest']]\n", "\n", "# 2) apply different formats per column\n", "styled_table = summary_var_table_7.style.format({\n", " 'Coefficient_SDvsAge': fmt_small_conditional,\n", " 'P_value': fmt_small_conditional,\n", " 'ICC': '{:.4f}',\n", " 'Q_value': fmt_small_conditional})\n", "\n", "# 3) display it\n", "styled_table" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.11" } }, "nbformat": 4, "nbformat_minor": 5 }