{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Lifestyle Archetypes and Problematic Internet Use\n\n**Question:** Given patterns of sleep, physical activity and body composition, do distinct lifestyle profiles emerge among children and adolescents, and how do these profiles relate to problematic internet use (PIU) and mental health screening scores?\n\n**Important analytical boundary (Responsible Analysis Requirement, Section 8):**\n- This is an **unsupervised, exploratory segmentation** task, not a predictive or diagnostic modelling task.\n- PIU and mental health scores here are **screening instrument results, not clinical diagnoses**.\n- Clusters describe **population-level lifestyle patterns**. They must never be used to label, diagnose, or make judgements about any individual child.\n- Any outcome association reported below (Section 10) is **preliminary and subset-derived** (train-set labels only) and is reported as an association, not a validated clinical finding.\n\n## Notebook structure\n1. Data dictionary review\n2. Missingness documentation\n3. Actigraphy (accelerometer) feature extraction\n4. Feature engineering & merge\n5. Exploratory data analysis\n6. Dimensionality reduction (UMAP)\n7. Clustering comparison (K-means, Agglomerative, DBSCAN)\n8. Cluster validation (silhouette, Davies-Bouldin, Calinski-Harabasz)\n9. Stability analysis (subsampling + noise sensitivity)\n10. Archetype characterisation & relation to PIU / mental health\n11. Visual summary (UMAP plot, radar charts, heatmap)\n12. Optional stretch: robustness of cluster-PIU association\n13. Written summary & limitations\n","metadata":{}},{"cell_type":"code","source":"import os\n\nos.listdir(\"/kaggle/input\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:15:17.329895Z","iopub.execute_input":"2026-07-15T12:15:17.330213Z","iopub.status.idle":"2026-07-15T12:15:17.341159Z","shell.execute_reply.started":"2026-07-15T12:15:17.330183Z","shell.execute_reply":"2026-07-15T12:15:17.340132Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\n\nos.listdir(\"/kaggle/input/competitions\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:17:02.926947Z","iopub.execute_input":"2026-07-15T12:17:02.927284Z","iopub.status.idle":"2026-07-15T12:17:02.933791Z","shell.execute_reply.started":"2026-07-15T12:17:02.927254Z","shell.execute_reply":"2026-07-15T12:17:02.932932Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pathlib import Path\n\nDATA_DIR = Path(\"/kaggle/input/competitions/child-mind-institute-problematic-internet-use\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:18:52.082073Z","iopub.execute_input":"2026-07-15T12:18:52.082812Z","iopub.status.idle":"2026-07-15T12:18:52.086938Z","shell.execute_reply.started":"2026-07-15T12:18:52.082779Z","shell.execute_reply":"2026-07-15T12:18:52.086054Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\n\nos.listdir(DATA_DIR)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:19:31.031798Z","iopub.execute_input":"2026-07-15T12:19:31.032133Z","iopub.status.idle":"2026-07-15T12:19:31.039228Z","shell.execute_reply.started":"2026-07-15T12:19:31.032102Z","shell.execute_reply":"2026-07-15T12:19:31.038434Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"TRAIN_CSV = DATA_DIR / \"train.csv\"\nTEST_CSV = DATA_DIR / \"test.csv\"\nSAMPLE_SUBMISSION = DATA_DIR / \"sample_submission.csv\"\nDATA_DICTIONARY = DATA_DIR / \"data_dictionary.csv\"\n\nSERIES_TRAIN_DIR = DATA_DIR / \"series_train.parquet\"\nSERIES_TEST_DIR = DATA_DIR / \"series_test.parquet\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:22:17.403169Z","iopub.execute_input":"2026-07-15T12:22:17.403607Z","iopub.status.idle":"2026-07-15T12:22:17.408918Z","shell.execute_reply.started":"2026-07-15T12:22:17.403576Z","shell.execute_reply":"2026-07-15T12:22:17.407999Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip -q install umap-learn pyarrow scikit-learn-extra missingno","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:26:19.885868Z","iopub.execute_input":"2026-07-15T12:26:19.886703Z","iopub.status.idle":"2026-07-15T12:27:21.134969Z","shell.execute_reply.started":"2026-07-15T12:26:19.886666Z","shell.execute_reply":"2026-07-15T12:27:21.133898Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport glob\nimport warnings\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom scipy import stats\nfrom scipy.stats import f_oneway, chi2_contingency, kruskal\n\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.preprocessing import StandardScaler\n\nfrom sklearn.decomposition import PCA\n\nfrom sklearn.cluster import (\n    KMeans,\n    AgglomerativeClustering,\n    DBSCAN\n)\n\nfrom sklearn.metrics import (\n    silhouette_score,\n    davies_bouldin_score,\n    calinski_harabasz_score\n)\n\nimport umap.umap_ as umap\n\nimport pyarrow.parquet as pq","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:31:06.752447Z","iopub.execute_input":"2026-07-15T12:31:06.752771Z","iopub.status.idle":"2026-07-15T12:31:36.398737Z","shell.execute_reply.started":"2026-07-15T12:31:06.752743Z","shell.execute_reply":"2026-07-15T12:31:36.397715Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"DATA_DIR = Path(\n    \"/kaggle/input/competitions/child-mind-institute-problematic-internet-use\"\n)\n\nTRAIN_CSV = DATA_DIR / \"train.csv\"\nTEST_CSV = DATA_DIR / \"test.csv\"\nDATA_DICTIONARY = DATA_DIR / \"data_dictionary.csv\"\nSAMPLE_SUBMISSION = DATA_DIR / \"sample_submission.csv\"\n\nSERIES_DIR = DATA_DIR / \"series_train.parquet\"\nSERIES_TEST_DIR = DATA_DIR / \"series_test.parquet\"\n\nCACHE_PATH = Path(\"/kaggle/working/actigraphy_summary.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:32:04.642533Z","iopub.execute_input":"2026-07-15T12:32:04.643453Z","iopub.status.idle":"2026-07-15T12:32:04.648309Z","shell.execute_reply.started":"2026-07-15T12:32:04.643410Z","shell.execute_reply":"2026-07-15T12:32:04.647415Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train = pd.read_csv(TRAIN_CSV)\ntest = pd.read_csv(TEST_CSV)\ndictionary = pd.read_csv(DATA_DICTIONARY)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:32:08.577339Z","iopub.execute_input":"2026-07-15T12:32:08.577703Z","iopub.status.idle":"2026-07-15T12:32:08.652616Z","shell.execute_reply.started":"2026-07-15T12:32:08.577674Z","shell.execute_reply":"2026-07-15T12:32:08.651895Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 1. Understanding the Dataset (Data Dictionary Review)\n\nBefore creating any new features or running any analysis, it is important to understand what each group of variables represents. The **Child Mind Institute (CMI) Healthy Brain Network** dataset organises variables into different sections based on the questionnaire, assessment, or measurement from which they were collected.\n\nThe main groups of variables are:\n\n| Variable Group | Full Name | What It Contains |\n|----------------|-----------|------------------|\n| **Basic_Demos** | Basic Demographics | Participant information such as age, sex, and basic background characteristics. |\n| **CGAS** | Children's Global Assessment Scale | A clinician-rated measure of a child's overall psychological, social, and school functioning. Higher scores generally indicate better functioning. |\n| **Physical** | Physical Measurements | Height, weight, Body Mass Index (BMI), waist circumference, blood pressure, heart rate, and other physical health measurements. |\n| **FGC** | FitnessGram Components | Measures of physical fitness such as grip strength, flexibility, muscular endurance, and cardiovascular fitness. |\n| **BIA** | Bioelectrical Impedance Analysis | Body composition measurements including body fat percentage, lean muscle mass, total body water, and related metrics. |\n| **PAQ_A** | Physical Activity Questionnaire for Adolescents | A self-reported questionnaire measuring the physical activity levels of adolescents (typically ages 14–20) over the previous seven days. |\n| **PAQ_C** | Physical Activity Questionnaire for Children | A child-friendly version of the physical activity questionnaire designed for younger participants (typically ages 8–14). |\n| **PCIAT** | Parent–Child Internet Addiction Test | A questionnaire assessing problematic internet use (PIU). It measures behaviours such as excessive internet use, difficulty controlling internet use, and its impact on daily life. This is the primary outcome used later to compare lifestyle archetypes—it is **not** used to create the clusters themselves. |\n| **SDS** | Sleep Disturbance Scale for Children | A questionnaire measuring different aspects of children's sleep, including sleep quality, sleep disorders, and disturbances. |\n| **PreInt_EduHx** | Prenatal, Intervention and Educational History | Information about prenatal history, developmental milestones, educational support, previous interventions, and learning history. |\n\n### Accelerometer (Wearable) Data\n\nIn addition to the questionnaire and clinical measurements, the dataset includes wearable accelerometer data stored separately as:\n\n- `series_train.parquet`\n- `series_test.parquet`\n\nThese files contain continuous movement data collected from wearable devices (actigraphy). Unlike the questionnaire data, they record participants' activity throughout the day and night.\n\nFrom these time-series data, we can engineer lifestyle features such as:\n\n- Average nightly sleep duration\n- Night-to-night sleep variability\n- Physical activity intensity\n- Daily activity consistency\n- Circadian (day–night) activity patterns\n\nThese engineered features provide objective measurements of participants' lifestyles and form the basis of the clustering analysis used to identify different lifestyle archetypes.","metadata":{}},{"cell_type":"code","source":"print(\"train shape:\", train.shape)\nprint(\"test shape :\", test.shape)\nprint(\"dictionary shape:\", dictionary.shape)\n\ndisplay(dictionary.head(20))\n\n# Column families present in train, grouped by instrument prefix\nfamilies = sorted(set(c.split(\"-\")[0] for c in train.columns if \"-\" in c))\nprint(\"\\nInstrument families in train.csv:\")\nfor fam in families:\n    cols = [c for c in train.columns if c.startswith(fam + \"-\")]\n    print(f\"  {fam:15s} -> {len(cols)} columns\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:41:27.607595Z","iopub.execute_input":"2026-07-15T12:41:27.608292Z","iopub.status.idle":"2026-07-15T12:41:27.651296Z","shell.execute_reply.started":"2026-07-15T12:41:27.608259Z","shell.execute_reply":"2026-07-15T12:41:27.650314Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 2. Missingness Documentation\n\nSeveral instruments (FitnessGram, BIA, PAQ) were only administered to subsets of participants (e.g. `PAQ_A` for adolescents, `PAQ_C` for children), so missingness is expected to be structured rather than random. We document it explicitly before deciding how to handle it, per the responsible-analysis requirement to make every data decision transparent.","metadata":{}},{"cell_type":"code","source":"missing = (\n    train.isna().mean().sort_values(ascending=False).mul(100).round(1)\n    .rename(\"pct_missing\")\n    .to_frame()\n)\nmissing[\"n_missing\"] = train.isna().sum().reindex(missing.index)\ndisplay(missing.head(30))\n\nplt.figure(figsize=(8, 10))\ntop_missing = missing.head(30)\nsns.barplot(x=top_missing[\"pct_missing\"], y=top_missing.index, color=\"#4C72B0\")\nplt.xlabel(\"% missing\")\nplt.title(\"Top 30 columns by missingness (train.csv)\")\nplt.tight_layout()\nplt.show()\n\ntry:\n    import missingno as msno\n    msno.matrix(train.sample(min(300, len(train)), random_state=42))\n    plt.title(\"Missingness pattern (random sample of 300 rows)\")\n    plt.show()\nexcept ImportError:\n    print(\"missingno not available - skipping matrix plot\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:42:02.739123Z","iopub.execute_input":"2026-07-15T12:42:02.739450Z","iopub.status.idle":"2026-07-15T12:42:03.667894Z","shell.execute_reply.started":"2026-07-15T12:42:02.739424Z","shell.execute_reply":"2026-07-15T12:42:03.666955Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Observation to carry forward:** questionnaire-instrument missingness (PAQ, FGC, BIA) is largely structural (age-gated instruments, or a participant simply not completing a station), while the accelerometer series is present for only a subset of participants (`series_train.parquet` has one folder per child who wore the device; not every `id` in `train.csv` has a matching recording). We handle this explicitly in Section 3-4: children without actigraphy retain missing wearable features and are imputed at the model-input stage, with imputation applied only to clustering *inputs*, never to the outcome variables used later in Section 10.","metadata":{}},{"cell_type":"markdown","source":"## 3. Actigraphy (Accelerometer) Feature Extraction\n\n`series_train.parquet` stores one sub-folder per child (`id=<participant_id>/part-0.parquet`) with 5-second epoch-level accelerometer data: `X`, `Y`, `Z`, `enmo` (Euclidean Norm Minus One, a standard activity-intensity metric), `anglez`, `light`, `battery_voltage`, `non-wear_flag`, `time_of_day` (nanoseconds since midnight), `weekday` (1=Monday...7=Sunday) and `relative_date_PCIAT` (day index relative to the PCIAT assessment).\n\nWe reduce each child's raw time series to a small number of **lifestyle-relevant summary features**: how much they sleep, how variable that sleep is, how physically active they are, how variable that activity is, and how different their weekday vs weekend activity looks. This mirrors the \"reduce dimensionality\" step suggested for the wearable data before it ever reaches UMAP.\n\n**Method transparency (documented simplifications):**\n- Epochs are 5 seconds long (`epoch_seconds = 5`).\n- Non-wear epochs (`non-wear_flag == 1`) are excluded from all calculations.\n- A day is only used if it has at least 16 hours of valid wear time (`min_valid_hours_per_day = 16`), to avoid days dominated by device removal biasing sleep/activity estimates.\n- \"Night-time\" is defined as clock hours 22:00-08:00, derived from `time_of_day`. This is a fixed-window heuristic, **not** a validated polysomnography-grade sleep-staging algorithm.\n- An epoch is flagged \"low activity\" if `enmo < 0.02` (g), a threshold consistent with common wrist-worn accelerometer sleep-detection conventions in the actigraphy literature. Sleep duration for a night is the count of low-activity epochs within the night window, converted to hours.\n- \"Weekday\" = days 1-5, \"weekend\" = days 6-7 (`weekday` column).\n\nThis is intentionally a pragmatic, well-documented simplification appropriate for population-level segmentation - not a clinical sleep score.","metadata":{}},{"cell_type":"code","source":"EPOCH_SECONDS = 5\nLOW_ACTIVITY_ENMO_THRESHOLD = 0.02\nMIN_VALID_HOURS_PER_DAY = 16\nNIGHT_START_HOUR = 22\nNIGHT_END_HOUR = 8\n\n\ndef _hour_of_day(time_of_day_ns):\n    \"\"\"Convert the 'time_of_day' nanoseconds-since-midnight column to a\n    fractional clock hour in [0, 24).\"\"\"\n    return (time_of_day_ns / 1e9) / 3600.0\n\n\ndef extract_actigraphy_features(participant_dir):\n    \"\"\"Summarise one child's accelerometer recording into lifestyle features.\n\n    Returns a dict of features, or None if the recording has no usable\n    (sufficiently-worn) days.\n    \"\"\"\n    parquet_files = glob.glob(os.path.join(str(participant_dir), \"*.parquet\"))\n    if not parquet_files:\n        return None\n\n    df = pq.read_table(parquet_files[0]).to_pandas()\n    if \"non-wear_flag\" in df.columns:\n        df = df[df[\"non-wear_flag\"] == 0]\n    if df.empty:\n        return None\n\n    df[\"hour\"] = _hour_of_day(df[\"time_of_day\"].to_numpy())\n    df[\"is_night\"] = (df[\"hour\"] >= NIGHT_START_HOUR) | (df[\"hour\"] < NIGHT_END_HOUR)\n    df[\"is_weekend\"] = df[\"weekday\"].isin([6, 7])\n\n    # keep only days with sufficient valid wear time\n    epochs_per_day = df.groupby(\"relative_date_PCIAT\").size()\n    valid_hours_per_day = epochs_per_day * EPOCH_SECONDS / 3600.0\n    valid_days = valid_hours_per_day[valid_hours_per_day >= MIN_VALID_HOURS_PER_DAY].index\n    df = df[df[\"relative_date_PCIAT\"].isin(valid_days)]\n    if df.empty:\n        return None\n\n    # --- sleep duration per night ---\n    night_df = df[df[\"is_night\"]]\n    sleep_epochs = night_df[night_df[\"enmo\"] < LOW_ACTIVITY_ENMO_THRESHOLD]\n    sleep_hours_per_day = (\n        sleep_epochs.groupby(\"relative_date_PCIAT\").size() * EPOCH_SECONDS / 3600.0\n    )\n    sleep_duration_mean = sleep_hours_per_day.mean() if len(sleep_hours_per_day) else np.nan\n    sleep_duration_std = sleep_hours_per_day.std() if len(sleep_hours_per_day) > 1 else np.nan\n\n    # --- daytime activity intensity ---\n    day_df = df[~df[\"is_night\"]]\n    daily_mean_enmo = day_df.groupby(\"relative_date_PCIAT\")[\"enmo\"].mean()\n    activity_intensity_mean = daily_mean_enmo.mean() if len(daily_mean_enmo) else np.nan\n    activity_intensity_std = daily_mean_enmo.std() if len(daily_mean_enmo) > 1 else np.nan\n\n    # --- weekday vs weekend consistency ---\n    weekday_mean = day_df.loc[~day_df[\"is_weekend\"], \"enmo\"].mean()\n    weekend_mean = day_df.loc[day_df[\"is_weekend\"], \"enmo\"].mean()\n    if pd.notna(weekday_mean) and weekday_mean > 0 and pd.notna(weekend_mean):\n        weekday_weekend_ratio = weekend_mean / weekday_mean\n    else:\n        weekday_weekend_ratio = np.nan\n\n    n_valid_days = len(valid_days)\n\n    return {\n        \"sleep_duration_mean\": sleep_duration_mean,\n        \"sleep_duration_std\": sleep_duration_std,\n        \"activity_intensity_mean\": activity_intensity_mean,\n        \"activity_intensity_std\": activity_intensity_std,\n        \"weekday_weekend_activity_ratio\": weekday_weekend_ratio,\n        \"n_valid_actigraphy_days\": n_valid_days,\n    }\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:42:45.733697Z","iopub.execute_input":"2026-07-15T12:42:45.734076Z","iopub.status.idle":"2026-07-15T12:42:45.747012Z","shell.execute_reply.started":"2026-07-15T12:42:45.734047Z","shell.execute_reply":"2026-07-15T12:42:45.746042Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 4. Feature Engineering & Merge\n\nWe now merge the actigraphy summary onto the tabular data and assemble the final **lifestyle feature set** used for clustering. Following the brief, demographics (age, sex) and outcome measures (PIU / mental health scores) are deliberately **excluded from the clustering inputs** - they are held back for Section 10 (characterisation) so that the segmentation reflects lifestyle only, not the outcomes we later relate it to.\n\n**Clustering input features:**\n\n| Feature | Source | Represents |\n|---|---|---|\n| `sleep_duration_mean`, `sleep_duration_std` | Actigraphy | Sleep amount & regularity |\n| `activity_intensity_mean`, `activity_intensity_std` | Actigraphy | Physical activity level & consistency |\n| `weekday_weekend_activity_ratio` | Actigraphy | Routine consistency across the week |\n| `BMI` | Physical exam | Body composition |\n| `BIA-BIA_Fat`, `BIA-BIA_FFMI`, `BIA-BIA_SMM` | Bioelectrical impedance analysis | Body composition (fat mass, fat-free mass index, skeletal muscle mass) |\n| `PAQ_Total` | Physical Activity Questionnaire (self/parent-report) | Self-reported activity level, complementing the objective wearable measure |\n\n**Held back for characterisation only (Section 10):** `Basic_Demos-Age`, `Basic_Demos-Sex`, `PCIAT-PCIAT_Total`, `sii`, `CGAS-CGAS_Score`, `SDS-SDS_Total_T`, `PreInt_EduHx-computerinternet_hoursday`.","metadata":{}},{"cell_type":"code","source":"if CACHE_PATH.exists():\n    actigraphy_summary = pd.read_csv(CACHE_PATH)\n    print(f\"Loaded cached actigraphy summary: {actigraphy_summary.shape}\")\nelse:\n    participant_dirs = sorted(glob.glob(str(SERIES_DIR / \"id=*\")))\n    print(f\"Found {len(participant_dirs)} participants with actigraphy recordings\")\n\n    try:\n        from tqdm.notebook import tqdm\n    except ImportError:\n        tqdm = lambda x, **kw: x\n\n    records = []\n    for pdir in tqdm(participant_dirs, desc=\"Extracting actigraphy features\"):\n        participant_id = os.path.basename(pdir).replace(\"id=\", \"\")\n        with warnings.catch_warnings():\n            warnings.simplefilter(\"ignore\")\n            feats = extract_actigraphy_features(pdir)\n        if feats is not None:\n            feats[\"id\"] = participant_id\n            records.append(feats)\n\n    actigraphy_summary = pd.DataFrame.from_records(records)\n    CACHE_PATH.parent.mkdir(parents=True, exist_ok=True)\n    actigraphy_summary.to_csv(CACHE_PATH, index=False)\n    print(f\"Extracted and cached actigraphy summary: {actigraphy_summary.shape}\")\n\ndisplay(actigraphy_summary.head())\nprint(\"\\nCoverage: {}/{} train participants have usable actigraphy features\".format(\n    train[\"id\"].isin(actigraphy_summary[\"id\"]).sum(), len(train)\n))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:43:19.961302Z","iopub.execute_input":"2026-07-15T12:43:19.962152Z","iopub.status.idle":"2026-07-15T12:45:36.845829Z","shell.execute_reply.started":"2026-07-15T12:43:19.962120Z","shell.execute_reply":"2026-07-15T12:45:36.845058Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df = train.merge(actigraphy_summary, on=\"id\", how=\"left\")\n\n# harmonise PAQ_A / PAQ_C (age-gated instruments) into a single lifestyle feature\ndf[\"PAQ_Total\"] = df[\"PAQ_A-PAQ_A_Total\"].combine_first(df[\"PAQ_C-PAQ_C_Total\"])\n\n# prefer the direct physical-exam BMI, falling back to BIA-derived BMI\ndf[\"BMI\"] = df[\"Physical-BMI\"].combine_first(df[\"BIA-BIA_BMI\"])\n\nCLUSTER_FEATURES = [\n    \"sleep_duration_mean\",\n    \"sleep_duration_std\",\n    \"activity_intensity_mean\",\n    \"activity_intensity_std\",\n    \"weekday_weekend_activity_ratio\",\n    \"BMI\",\n    \"BIA-BIA_Fat\",\n    \"BIA-BIA_FFMI\",\n    \"BIA-BIA_SMM\",\n    \"PAQ_Total\",\n]\n\nCHARACTERISATION_VARS = [\n    \"Basic_Demos-Age\",\n    \"Basic_Demos-Sex\",\n    \"PCIAT-PCIAT_Total\",\n    \"sii\",\n    \"CGAS-CGAS_Score\",\n    \"SDS-SDS_Total_T\",\n    \"PreInt_EduHx-computerinternet_hoursday\",\n]\n\nprint(\"Missingness of clustering input features:\")\ndisplay(df[CLUSTER_FEATURES].isna().mean().mul(100).round(1).rename(\"pct_missing\").to_frame())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:46:26.660849Z","iopub.execute_input":"2026-07-15T12:46:26.661298Z","iopub.status.idle":"2026-07-15T12:46:26.689854Z","shell.execute_reply.started":"2026-07-15T12:46:26.661245Z","shell.execute_reply":"2026-07-15T12:46:26.688835Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We require at least the actigraphy-derived features to be present (a child must have worn the device to contribute a lifestyle profile), then median-impute any remaining sporadic gaps in the body-composition / questionnaire columns. Median imputation is chosen over mean because several body-composition variables are right-skewed; it is applied only to the small residual missingness left after dropping actigraphy-less rows, and is documented here as a modelling choice rather than silently applied.","metadata":{}},{"cell_type":"code","source":"has_actigraphy = df[\"sleep_duration_mean\"].notna() & df[\"activity_intensity_mean\"].notna()\nprint(f\"Rows with usable actigraphy: {has_actigraphy.sum()} / {len(df)}\")\n\nmodel_df = df.loc[has_actigraphy, [\"id\"] + CLUSTER_FEATURES + CHARACTERISATION_VARS].reset_index(drop=True)\n\nimputer = SimpleImputer(strategy=\"median\")\nX_imputed = imputer.fit_transform(model_df[CLUSTER_FEATURES])\nX_imputed = pd.DataFrame(X_imputed, columns=CLUSTER_FEATURES, index=model_df.index)\n\nprint(f\"\\nFinal modelling sample: {model_df.shape[0]} children, {len(CLUSTER_FEATURES)} lifestyle features\")\ndisplay(X_imputed.describe().T)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:46:52.222063Z","iopub.execute_input":"2026-07-15T12:46:52.222508Z","iopub.status.idle":"2026-07-15T12:46:52.273693Z","shell.execute_reply.started":"2026-07-15T12:46:52.222463Z","shell.execute_reply":"2026-07-15T12:46:52.272478Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 5. Exploratory Data Analysis\n\nDistributions and pairwise correlations of the engineered lifestyle features, before standardisation.","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(2, 5, figsize=(20, 7))\nfor ax, col in zip(axes.flat, CLUSTER_FEATURES):\n    sns.histplot(X_imputed[col], kde=True, ax=ax, color=\"#4C72B0\")\n    ax.set_title(col, fontsize=10)\nplt.tight_layout()\nplt.show()\n\nplt.figure(figsize=(9, 7))\ncorr = X_imputed.corr()\nsns.heatmap(corr, annot=True, fmt=\".2f\", cmap=\"vlag\", center=0, square=True)\nplt.title(\"Pairwise correlation of lifestyle features\")\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:47:17.294269Z","iopub.execute_input":"2026-07-15T12:47:17.295264Z","iopub.status.idle":"2026-07-15T12:47:21.737355Z","shell.execute_reply.started":"2026-07-15T12:47:17.295217Z","shell.execute_reply":"2026-07-15T12:47:21.736532Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6. Dimensionality Reduction (UMAP)\n\nFeatures are standardised first (z-scored), since they are on very different scales (hours, g, kg/m^2, questionnaire totals). We then explore how UMAP's `n_neighbors` (local vs global structure) and `min_dist` (how tightly points are allowed to pack) affect the shape of the embedding before committing to final parameters, as required by the brief.","metadata":{}},{"cell_type":"code","source":"scaler = StandardScaler()\nX_scaled = scaler.fit_transform(X_imputed)\nX_scaled = pd.DataFrame(X_scaled, columns=CLUSTER_FEATURES, index=X_imputed.index)\n\n# quick PCA sanity check - how much variance do the raw features already carry\npca = PCA(random_state=42).fit(X_scaled)\nplt.figure(figsize=(6, 4))\nplt.plot(np.cumsum(pca.explained_variance_ratio_), marker=\"o\")\nplt.xlabel(\"Number of components\")\nplt.ylabel(\"Cumulative explained variance\")\nplt.title(\"PCA on standardised lifestyle features\")\nplt.grid(alpha=0.3)\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:48:38.665105Z","iopub.execute_input":"2026-07-15T12:48:38.666200Z","iopub.status.idle":"2026-07-15T12:48:38.846858Z","shell.execute_reply.started":"2026-07-15T12:48:38.666159Z","shell.execute_reply":"2026-07-15T12:48:38.846024Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"neighbor_grid = [5, 15, 50]\nmin_dist_grid = [0.0, 0.5]\n\nfig, axes = plt.subplots(len(neighbor_grid), len(min_dist_grid), figsize=(10, 14))\nfor i, n_neighbors in enumerate(neighbor_grid):\n    for j, min_dist in enumerate(min_dist_grid):\n        reducer = umap.UMAP(\n            n_neighbors=n_neighbors,\n            min_dist=min_dist,\n            n_components=2,\n            random_state=42,\n        )\n        emb = reducer.fit_transform(X_scaled)\n        ax = axes[i, j]\n        ax.scatter(emb[:, 0], emb[:, 1], s=6, alpha=0.5, color=\"#4C72B0\")\n        ax.set_title(f\"n_neighbors={n_neighbors}, min_dist={min_dist}\", fontsize=9)\n        ax.set_xticks([]); ax.set_yticks([])\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:49:06.539840Z","iopub.execute_input":"2026-07-15T12:49:06.540174Z","iopub.status.idle":"2026-07-15T12:49:29.268818Z","shell.execute_reply.started":"2026-07-15T12:49:06.540147Z","shell.execute_reply":"2026-07-15T12:49:29.267676Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Parameter choice rationale:** a moderate `n_neighbors` (15) balances local detail against the global structure needed for population-level archetypes (very small `n_neighbors` fragments the embedding into noisy local clusters; very large `n_neighbors` over-smooths genuine subgroup structure). A small but non-zero `min_dist` (0.1) keeps groups visually separable without artificially collapsing genuinely continuous variation into artefactual tight blobs. We fit both a 2D embedding (for visualisation) and a higher-dimensional embedding (for clustering itself, which benefits from retaining more structure than 2D allows).","metadata":{}},{"cell_type":"code","source":"UMAP_N_NEIGHBORS = 15\nUMAP_MIN_DIST = 0.1\n\nreducer_2d = umap.UMAP(\n    n_neighbors=UMAP_N_NEIGHBORS, min_dist=UMAP_MIN_DIST, n_components=2, random_state=42\n)\nembedding_2d = reducer_2d.fit_transform(X_scaled)\n\nreducer_5d = umap.UMAP(\n    n_neighbors=UMAP_N_NEIGHBORS, min_dist=UMAP_MIN_DIST, n_components=5, random_state=42\n)\nembedding_5d = reducer_5d.fit_transform(X_scaled)\n\nprint(\"2D embedding shape:\", embedding_2d.shape)\nprint(\"5D embedding shape (used for clustering):\", embedding_5d.shape)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:50:04.564498Z","iopub.execute_input":"2026-07-15T12:50:04.565009Z","iopub.status.idle":"2026-07-15T12:50:08.630446Z","shell.execute_reply.started":"2026-07-15T12:50:04.564947Z","shell.execute_reply":"2026-07-15T12:50:08.629443Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 7. Clustering Comparison\n\nWe compare three families of clustering algorithm on the 5D UMAP embedding, each with a different bias: K-means (spherical, evenly-sized clusters), Agglomerative/hierarchical (nested, variable-shape clusters), and DBSCAN (density-based, can flag noise/outliers rather than forcing every child into a cluster). Testing more than one method guards against reporting an artefact of a single algorithm's assumptions.","metadata":{}},{"cell_type":"code","source":"# --- K-means: elbow + silhouette across k ---\nk_range = range(2, 9)\ninertias, sil_scores_km = [], []\nfor k in k_range:\n    km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(embedding_5d)\n    inertias.append(km.inertia_)\n    sil_scores_km.append(silhouette_score(embedding_5d, km.labels_))\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 4))\naxes[0].plot(list(k_range), inertias, marker=\"o\")\naxes[0].set_xlabel(\"k\"); axes[0].set_ylabel(\"Inertia\"); axes[0].set_title(\"K-means elbow\")\naxes[1].plot(list(k_range), sil_scores_km, marker=\"o\", color=\"darkorange\")\naxes[1].set_xlabel(\"k\"); axes[1].set_ylabel(\"Silhouette score\"); axes[1].set_title(\"K-means silhouette vs k\")\nplt.tight_layout()\nplt.show()\n\nbest_k = list(k_range)[int(np.argmax(sil_scores_km))]\nprint(f\"Silhouette-selected k for K-means: {best_k}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:50:27.224516Z","iopub.execute_input":"2026-07-15T12:50:27.224876Z","iopub.status.idle":"2026-07-15T12:50:28.239000Z","shell.execute_reply.started":"2026-07-15T12:50:27.224847Z","shell.execute_reply":"2026-07-15T12:50:28.238091Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"kmeans_final = KMeans(n_clusters=best_k, n_init=10, random_state=42).fit(embedding_5d)\nlabels_kmeans = kmeans_final.labels_","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:51:18.904042Z","iopub.execute_input":"2026-07-15T12:51:18.904403Z","iopub.status.idle":"2026-07-15T12:51:18.922372Z","shell.execute_reply.started":"2026-07-15T12:51:18.904351Z","shell.execute_reply":"2026-07-15T12:51:18.921313Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Agglomerative (Ward linkage) ---\nfrom scipy.cluster.hierarchy import dendrogram, linkage\n\nZ = linkage(embedding_5d, method=\"ward\")\nplt.figure(figsize=(12, 4))\ndendrogram(Z, truncate_mode=\"lastp\", p=30, show_leaf_counts=True)\nplt.title(\"Agglomerative clustering (Ward linkage) - truncated dendrogram\")\nplt.xlabel(\"Cluster size\"); plt.ylabel(\"Distance\")\nplt.show()\n\nsil_scores_agg = []\nfor k in k_range:\n    agg = AgglomerativeClustering(n_clusters=k, linkage=\"ward\").fit(embedding_5d)\n    sil_scores_agg.append(silhouette_score(embedding_5d, agg.labels_))\n\nplt.figure(figsize=(6, 4))\nplt.plot(list(k_range), sil_scores_agg, marker=\"o\", color=\"seagreen\")\nplt.xlabel(\"k\"); plt.ylabel(\"Silhouette score\"); plt.title(\"Agglomerative silhouette vs k\")\nplt.show()\n\nbest_k_agg = list(k_range)[int(np.argmax(sil_scores_agg))]\nprint(f\"Silhouette-selected k for Agglomerative: {best_k_agg}\")\nagg_final = AgglomerativeClustering(n_clusters=best_k_agg, linkage=\"ward\").fit(embedding_5d)\nlabels_agg = agg_final.labels_","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:51:38.578863Z","iopub.execute_input":"2026-07-15T12:51:38.579682Z","iopub.status.idle":"2026-07-15T12:51:39.318428Z","shell.execute_reply.started":"2026-07-15T12:51:38.579636Z","shell.execute_reply":"2026-07-15T12:51:39.317199Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- DBSCAN: eps tuning via k-distance graph ---\nfrom sklearn.neighbors import NearestNeighbors\n\nMIN_SAMPLES = 10\nneighbors = NearestNeighbors(n_neighbors=MIN_SAMPLES).fit(embedding_5d)\ndistances, _ = neighbors.kneighbors(embedding_5d)\nk_distances = np.sort(distances[:, -1])\n\nplt.figure(figsize=(6, 4))\nplt.plot(k_distances)\nplt.xlabel(\"Points sorted by distance\"); plt.ylabel(f\"{MIN_SAMPLES}-NN distance\")\nplt.title(\"DBSCAN eps selection (k-distance graph - look for the 'elbow')\")\nplt.grid(alpha=0.3)\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:51:57.662635Z","iopub.execute_input":"2026-07-15T12:51:57.662965Z","iopub.status.idle":"2026-07-15T12:51:57.831260Z","shell.execute_reply.started":"2026-07-15T12:51:57.662938Z","shell.execute_reply":"2026-07-15T12:51:57.830277Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Elbow read from the plot above - adjust EPS if the elbow sits elsewhere for your run\nEPS = float(np.percentile(k_distances, 90))\nprint(f\"Selected eps (90th percentile of k-distance, as an elbow proxy): {EPS:.3f}\")\n\ndbscan_final = DBSCAN(eps=EPS, min_samples=MIN_SAMPLES).fit(embedding_5d)\nlabels_dbscan = dbscan_final.labels_\nn_clusters_dbscan = len(set(labels_dbscan)) - (1 if -1 in labels_dbscan else 0)\nn_noise = int((labels_dbscan == -1).sum())\nprint(f\"DBSCAN found {n_clusters_dbscan} clusters and flagged {n_noise} points as noise ({n_noise/len(labels_dbscan):.1%})\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:52:17.141331Z","iopub.execute_input":"2026-07-15T12:52:17.142295Z","iopub.status.idle":"2026-07-15T12:52:17.158348Z","shell.execute_reply.started":"2026-07-15T12:52:17.142260Z","shell.execute_reply":"2026-07-15T12:52:17.157396Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 8. Cluster Validation\n\nThree complementary internal-validation metrics, each measuring something different:\n- **Silhouette score** (-1 to 1, higher better): how well-separated and internally cohesive clusters are, on average, per point.\n- **Davies-Bouldin index** (>=0, lower better): average similarity between each cluster and its most-similar other cluster - penalises clusters that are close together relative to their own spread.\n- **Calinski-Harabasz index** (higher better): ratio of between-cluster to within-cluster dispersion - sensitive to the number of clusters, so it's read alongside the other two rather than alone.\n\nWe compare all three algorithms' chosen solutions on the same footing, and exclude DBSCAN's noise points (-1) from its own metric computation since those points are explicitly \"no cluster\", not a nominal cluster.","metadata":{}},{"cell_type":"code","source":"def validation_row(name, X, labels):\n    mask = labels != -1  # exclude DBSCAN noise points from validation of the clustered points\n    n_clusters = len(set(labels[mask]))\n    if n_clusters < 2:\n        return {\"method\": name, \"n_clusters\": n_clusters, \"silhouette\": np.nan,\n                \"davies_bouldin\": np.nan, \"calinski_harabasz\": np.nan}\n    return {\n        \"method\": name,\n        \"n_clusters\": n_clusters,\n        \"silhouette\": silhouette_score(X[mask], labels[mask]),\n        \"davies_bouldin\": davies_bouldin_score(X[mask], labels[mask]),\n        \"calinski_harabasz\": calinski_harabasz_score(X[mask], labels[mask]),\n    }\n\nvalidation_table = pd.DataFrame([\n    validation_row(\"K-means\", embedding_5d, labels_kmeans),\n    validation_row(\"Agglomerative\", embedding_5d, labels_agg),\n    validation_row(\"DBSCAN\", embedding_5d, labels_dbscan),\n]).round(3)\n\ndisplay(validation_table)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:52:37.311004Z","iopub.execute_input":"2026-07-15T12:52:37.311354Z","iopub.status.idle":"2026-07-15T12:52:37.393828Z","shell.execute_reply.started":"2026-07-15T12:52:37.311326Z","shell.execute_reply":"2026-07-15T12:52:37.393126Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 9. Stability Analysis\n\nA clustering solution that looks clean on one run but collapses under mild perturbation is not a reliable finding. We test the K-means solution (our leading candidate, pending Section 8's results) two ways:\n\n1. **Subsampling stability:** repeatedly re-cluster random 80% subsamples and measure agreement (Adjusted Rand Index) with cluster labels on the shared, overlapping points.\n2. **Noise sensitivity:** add small Gaussian noise to the standardised input features, re-embed and re-cluster, and measure agreement with the original labelling.\n\nARI ranges from ~0 (agreement no better than chance) to 1 (perfect agreement).","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import adjusted_rand_score\n\nN_TRIALS = 10\nSUBSAMPLE_FRAC = 0.8\n\nsubsample_aris = []\nrng = np.random.RandomState(42)\nfor trial in range(N_TRIALS):\n    idx = rng.choice(len(embedding_5d), size=int(len(embedding_5d) * SUBSAMPLE_FRAC), replace=False)\n    idx_sorted = np.sort(idx)\n    km_sub = KMeans(n_clusters=best_k, n_init=10, random_state=trial).fit(embedding_5d[idx_sorted])\n    ari = adjusted_rand_score(labels_kmeans[idx_sorted], km_sub.labels_)\n    subsample_aris.append(ari)\n\nprint(f\"Subsampling stability (k={best_k}, {N_TRIALS} trials, {SUBSAMPLE_FRAC:.0%} subsamples):\")\nprint(f\"  mean ARI = {np.mean(subsample_aris):.3f}  (std = {np.std(subsample_aris):.3f})\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:53:24.447220Z","iopub.execute_input":"2026-07-15T12:53:24.447565Z","iopub.status.idle":"2026-07-15T12:53:24.576729Z","shell.execute_reply.started":"2026-07-15T12:53:24.447538Z","shell.execute_reply":"2026-07-15T12:53:24.575896Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"NOISE_SIGMA = 0.1  # in standardised-feature units\n\nnoise_aris = []\nfor trial in range(N_TRIALS):\n    rs = np.random.RandomState(100 + trial)\n    X_noisy = X_scaled.to_numpy() + rs.normal(0, NOISE_SIGMA, X_scaled.shape)\n    reducer_noisy = umap.UMAP(\n        n_neighbors=UMAP_N_NEIGHBORS, min_dist=UMAP_MIN_DIST, n_components=5, random_state=trial\n    )\n    emb_noisy = reducer_noisy.fit_transform(X_noisy)\n    km_noisy = KMeans(n_clusters=best_k, n_init=10, random_state=trial).fit(emb_noisy)\n    ari = adjusted_rand_score(labels_kmeans, km_noisy.labels_)\n    noise_aris.append(ari)\n\nprint(f\"Noise sensitivity (sigma={NOISE_SIGMA}, {N_TRIALS} trials):\")\nprint(f\"  mean ARI = {np.mean(noise_aris):.3f}  (std = {np.std(noise_aris):.3f})\")\n\nprint(\n    \"\\nInterpretation guide: ARI > 0.6 suggests a reasonably stable solution; \"\n    \"0.3-0.6 suggests moderate stability worth flagging as a limitation; \"\n    \"< 0.3 suggests the solution should be treated as provisional/exploratory only.\"\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:53:45.478613Z","iopub.execute_input":"2026-07-15T12:53:45.478972Z","iopub.status.idle":"2026-07-15T12:54:06.783306Z","shell.execute_reply.started":"2026-07-15T12:53:45.478944Z","shell.execute_reply":"2026-07-15T12:54:06.782372Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 10. Archetype Characterisation\n\nWe adopt the K-means solution (pending the validation/stability results above being the strongest of the three) as the final archetype assignment, attach it back to the modelling dataframe, and describe each archetype in terms of its lifestyle features, demographics, and - as a **preliminary, subset-derived association only** - its relationship to PIU and mental health screening scores.","metadata":{}},{"cell_type":"code","source":"FINAL_LABELS = labels_kmeans  # swap to labels_agg / labels_dbscan here if validation favours another method\n\nmodel_df = model_df.copy()\nmodel_df[\"archetype\"] = FINAL_LABELS\nmodel_df[\"archetype\"] = model_df[\"archetype\"].astype(\"category\")\n\nprint(\"Archetype sizes:\")\ndisplay(model_df[\"archetype\"].value_counts().sort_index())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:54:44.125822Z","iopub.execute_input":"2026-07-15T12:54:44.126373Z","iopub.status.idle":"2026-07-15T12:54:44.138080Z","shell.execute_reply.started":"2026-07-15T12:54:44.126346Z","shell.execute_reply":"2026-07-15T12:54:44.137309Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- lifestyle-feature profile per archetype ---\nprofile = model_df.groupby(\"archetype\")[CLUSTER_FEATURES].agg([\"mean\", \"std\"])\ndisplay(profile.round(2))\n\n# --- ANOVA / Kruskal-Wallis: do lifestyle features differ significantly by archetype? ---\nprint(\"\\nSignificance of lifestyle-feature differences across archetypes:\")\nfor feat in CLUSTER_FEATURES:\n    groups = [g[feat].dropna().values for _, g in model_df.groupby(\"archetype\")]\n    f_stat, p_anova = f_oneway(*groups)\n    h_stat, p_kw = kruskal(*groups)\n    print(f\"  {feat:32s}  ANOVA p={p_anova:.4f}   Kruskal-Wallis p={p_kw:.4f}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:55:09.167660Z","iopub.execute_input":"2026-07-15T12:55:09.167994Z","iopub.status.idle":"2026-07-15T12:55:09.244400Z","shell.execute_reply.started":"2026-07-15T12:55:09.167967Z","shell.execute_reply":"2026-07-15T12:55:09.243603Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- sex distribution across archetypes (chi-square) ---\nsex_table = pd.crosstab(model_df[\"archetype\"], model_df[\"Basic_Demos-Sex\"])\nchi2, p_chi2, dof, _ = chi2_contingency(sex_table)\nprint(\"Sex distribution by archetype:\")\ndisplay(sex_table)\nprint(f\"Chi-square test: chi2={chi2:.2f}, dof={dof}, p={p_chi2:.4f}\")\n\nprint(\"\\nAge distribution by archetype:\")\ndisplay(model_df.groupby(\"archetype\")[\"Basic_Demos-Age\"].describe().round(2))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:55:30.054315Z","iopub.execute_input":"2026-07-15T12:55:30.054671Z","iopub.status.idle":"2026-07-15T12:55:30.107094Z","shell.execute_reply.started":"2026-07-15T12:55:30.054642Z","shell.execute_reply":"2026-07-15T12:55:30.106356Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Relating archetypes to PIU and mental health screening scores\n\n**Caution:** `PCIAT-PCIAT_Total`, `sii`, `CGAS-CGAS_Score` and `SDS-SDS_Total_T` were **not** used to form the archetypes above; they are examined here only to describe how the lifestyle-based segmentation relates to these screening outcomes. This is an **association from train-set labels on a subset of participants with usable actigraphy**, and should be read as hypothesis-generating, not as a validated clinical finding.","metadata":{}},{"cell_type":"code","source":"outcome_vars = [\"PCIAT-PCIAT_Total\", \"sii\", \"CGAS-CGAS_Score\", \"SDS-SDS_Total_T\",\n                \"PreInt_EduHx-computerinternet_hoursday\"]\n\nfig, axes = plt.subplots(1, len(outcome_vars), figsize=(22, 4))\nfor ax, col in zip(axes, outcome_vars):\n    sns.boxplot(data=model_df, x=\"archetype\", y=col, ax=ax, palette=\"Set2\")\n    ax.set_title(col, fontsize=9)\nplt.tight_layout()\nplt.show()\n\nprint(\"Significance of outcome differences across archetypes (Kruskal-Wallis, robust to non-normal screening scores):\")\nfor col in outcome_vars:\n    groups = [g[col].dropna().values for _, g in model_df.groupby(\"archetype\")]\n    groups = [g for g in groups if len(g) > 0]\n    if len(groups) < 2:\n        continue\n    h_stat, p_val = kruskal(*groups)\n    print(f\"  {col:42s}  Kruskal-Wallis p={p_val:.4f}\")\n\nprint(\"\\nsii (Severity Impairment Index category) distribution by archetype:\")\ndisplay(pd.crosstab(model_df[\"archetype\"], model_df[\"sii\"], normalize=\"index\").round(2))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:55:50.472677Z","iopub.execute_input":"2026-07-15T12:55:50.473645Z","iopub.status.idle":"2026-07-15T12:55:51.247056Z","shell.execute_reply.started":"2026-07-15T12:55:50.473609Z","shell.execute_reply":"2026-07-15T12:55:51.246234Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 11. Visual Summary\n\nUMAP scatter coloured by final archetype, plus one radar chart per archetype and a heatmap of standardised feature means - the requested \"one clear summary visualisation per archetype\".","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(7, 6))\nscatter = plt.scatter(\n    embedding_2d[:, 0], embedding_2d[:, 1],\n    c=FINAL_LABELS, cmap=\"tab10\", s=12, alpha=0.7\n)\nplt.legend(*scatter.legend_elements(), title=\"Archetype\", loc=\"best\")\nplt.title(\"UMAP projection coloured by lifestyle archetype\")\nplt.xlabel(\"UMAP-1\"); plt.ylabel(\"UMAP-2\")\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:56:16.454337Z","iopub.execute_input":"2026-07-15T12:56:16.454709Z","iopub.status.idle":"2026-07-15T12:56:16.759704Z","shell.execute_reply.started":"2026-07-15T12:56:16.454680Z","shell.execute_reply":"2026-07-15T12:56:16.758932Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# standardised feature means per archetype, for the radar charts and heatmap\narchetype_means_scaled = (\n    pd.DataFrame(X_scaled.values, columns=CLUSTER_FEATURES, index=model_df.index)\n    .groupby(model_df[\"archetype\"])\n    .mean()\n)\n\ndef radar_chart(ax, values, labels, title):\n    n = len(labels)\n    angles = np.linspace(0, 2 * np.pi, n, endpoint=False).tolist()\n    values = list(values) + [values[0]]\n    angles = angles + [angles[0]]\n    ax.plot(angles, values, color=\"#4C72B0\", linewidth=2)\n    ax.fill(angles, values, color=\"#4C72B0\", alpha=0.25)\n    ax.set_xticks(angles[:-1])\n    ax.set_xticklabels(labels, fontsize=7)\n    ax.set_title(title, fontsize=11, pad=15)\n    ax.set_yticklabels([])\n\nn_arch = archetype_means_scaled.shape[0]\nfig, axes = plt.subplots(1, n_arch, subplot_kw=dict(polar=True), figsize=(5 * n_arch, 5))\nif n_arch == 1:\n    axes = [axes]\nfor ax, (arch_id, row) in zip(axes, archetype_means_scaled.iterrows()):\n    n_members = (model_df[\"archetype\"] == arch_id).sum()\n    radar_chart(ax, row.values, CLUSTER_FEATURES, f\"Archetype {arch_id}  (n={n_members})\")\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:56:47.761174Z","iopub.execute_input":"2026-07-15T12:56:47.761587Z","iopub.status.idle":"2026-07-15T12:56:48.178629Z","shell.execute_reply.started":"2026-07-15T12:56:47.761555Z","shell.execute_reply":"2026-07-15T12:56:48.177824Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(9, max(3, 0.6 * n_arch + 2)))\nsns.heatmap(\n    archetype_means_scaled, annot=True, fmt=\".2f\", cmap=\"vlag\", center=0,\n    cbar_kws={\"label\": \"Standardised mean (z-score)\"}\n)\nplt.title(\"Archetype profile heatmap (standardised lifestyle feature means)\")\nplt.ylabel(\"Archetype\")\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:57:05.757524Z","iopub.execute_input":"2026-07-15T12:57:05.757870Z","iopub.status.idle":"2026-07-15T12:57:06.021207Z","shell.execute_reply.started":"2026-07-15T12:57:05.757842Z","shell.execute_reply":"2026-07-15T12:57:06.020226Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 12. Optional Stretch: Robustness of the Archetype-PIU Association\n\nAs a robustness check, we re-cluster with a different `k` and a different UMAP `n_neighbors`, and see whether the direction of the archetype-PIU association (which archetypes have higher/lower `PCIAT-PCIAT_Total`) holds up, rather than being an artefact of one specific parameter choice.","metadata":{}},{"cell_type":"code","source":"alt_k = max(2, best_k - 1) if best_k > 2 else best_k + 1\nalt_reducer = umap.UMAP(n_neighbors=30, min_dist=UMAP_MIN_DIST, n_components=5, random_state=7)\nalt_embedding = alt_reducer.fit_transform(X_scaled)\nalt_labels = KMeans(n_clusters=alt_k, n_init=10, random_state=7).fit_predict(alt_embedding)\n\nalt_df = model_df.copy()\nalt_df[\"alt_archetype\"] = alt_labels\n\nprint(f\"Alternative run: k={alt_k}, n_neighbors=30 (vs. original k={best_k}, n_neighbors={UMAP_N_NEIGHBORS})\")\nprint(\"\\nPCIAT-PCIAT_Total mean by alternative archetype:\")\ndisplay(alt_df.groupby(\"alt_archetype\")[\"PCIAT-PCIAT_Total\"].mean().round(1))\n\ngroups = [g[\"PCIAT-PCIAT_Total\"].dropna().values for _, g in alt_df.groupby(\"alt_archetype\")]\nh_stat, p_val = kruskal(*[g for g in groups if len(g) > 0])\nprint(f\"\\nKruskal-Wallis p={p_val:.4f} (alternative parameterisation)\")\nprint(\n    \"Robustness read: if archetypes with the lowest sleep / highest activity-variability \"\n    \"profile also show the highest PCIAT scores under this alternative parameterisation, \"\n    \"that supports the original finding being a genuine pattern rather than a clustering artefact.\"\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-15T12:57:36.110295Z","iopub.execute_input":"2026-07-15T12:57:36.110691Z","iopub.status.idle":"2026-07-15T12:57:38.716408Z","shell.execute_reply.started":"2026-07-15T12:57:36.110662Z","shell.execute_reply":"2026-07-15T12:57:38.715688Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}