{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":81933,"databundleVersionId":9643020,"sourceType":"competition"}],"dockerImageVersionId":30761,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# EDA which makes sense for the *Child Mind Institute — Problematic Internet Use* competition\n\nThis notebook shows\n- a first analysis of the data\n- how to cross-validate a model\n- that regression models are better than classification models in this competition, and\n- how to tune the thresholds for rounding the regression output.\n\nThe notebook uses [polars DataFrames](https://pola.rs/). If you are more fluent with pandas than with polars, this is an opportunity to get to know polars, which is often more efficient than pandas.\n\nReference:\n- [Competition homepage](https://www.kaggle.com/competitions/child-mind-institute-problematic-internet-use)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import polars as pl\nimport polars.selectors as cs\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom matplotlib.ticker import MaxNLocator, FormatStrFormatter, PercentFormatter\nimport numpy as np\nimport seaborn as sns\nimport lightgbm\nfrom colorama import Fore, Style\nfrom scipy.optimize import minimize\n\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.metrics import cohen_kappa_score, ConfusionMatrixDisplay\n\ntarget_labels = ['None', 'Mild', 'Moderate', 'Severe']","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:30.157442Z","iopub.execute_input":"2024-09-24T19:36:30.158627Z","iopub.status.idle":"2024-09-24T19:36:34.995669Z","shell.execute_reply.started":"2024-09-24T19:36:30.158554Z","shell.execute_reply":"2024-09-24T19:36:34.994351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"season_dtype = pl.Enum(['Spring', 'Summer', 'Fall', 'Winter'])\n\ntrain = (\n    pl.read_csv('/kaggle/input/child-mind-institute-problematic-internet-use/train.csv')\n    .with_columns(pl.col('^.*Season$').cast(season_dtype))\n)\n\ntest = (\n    pl.read_csv('/kaggle/input/child-mind-institute-problematic-internet-use/test.csv')\n    .with_columns(pl.col('^.*Season$').cast(season_dtype))\n)\n\ntrain","metadata":{"execution":{"iopub.status.busy":"2024-09-24T19:36:34.998189Z","iopub.execute_input":"2024-09-24T19:36:34.999091Z","iopub.status.idle":"2024-09-24T19:36:35.237668Z","shell.execute_reply.started":"2024-09-24T19:36:34.999029Z","shell.execute_reply":"2024-09-24T19:36:35.236200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The training dataset has 3960 samples (children who participate in the study) and 80 features (not counting the `id` column and the target `sii`). According to the documentation, the full test set comprises about 3800 instances, of which 1400 are public and 2400 private.\n\n# Missing values\n\nAll columns have a substantial proportion of missing values, except `id` (not surprisingly) and the three basic demographic columns for sex, age and season of enrollment. Even the target `sii` has missing values:","metadata":{}},{"cell_type":"code","source":"missing_count = (\n    train\n    .null_count()\n    .transpose(include_header=True,\n               header_name='feature',\n               column_names=['null_count'])\n    .sort('null_count', descending=True)\n    .with_columns((pl.col('null_count') / len(train)).alias('null_ratio'))\n)\nplt.figure(figsize=(6, 15))\nplt.title('Missing values over the whole training dataset')\nplt.barh(np.arange(len(missing_count)), missing_count.get_column('null_ratio'), color='coral', label='missing')\nplt.barh(np.arange(len(missing_count)), \n         1 - missing_count.get_column('null_ratio'),\n         left=missing_count.get_column('null_ratio'),\n         color='darkseagreen', label='available')\nplt.yticks(np.arange(len(missing_count)), missing_count.get_column('feature'))\nplt.gca().xaxis.set_major_formatter(PercentFormatter(xmax=1, decimals=0))\nplt.xlim(0, 1)\nplt.legend()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:35.240548Z","iopub.execute_input":"2024-09-24T19:36:35.241169Z","iopub.status.idle":"2024-09-24T19:36:36.911516Z","shell.execute_reply.started":"2024-09-24T19:36:35.241109Z","shell.execute_reply":"2024-09-24T19:36:36.909875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If we count the missing feature values only for the usable part of the training dataset (samples are usable for supervised training if the target is known), the chart looks slightly different:","metadata":{}},{"cell_type":"code","source":"supervised_usable = (\n    train\n    .filter(pl.col('sii').is_not_null())\n)\n\nmissing_count = (\n    supervised_usable\n    .null_count()\n    .transpose(include_header=True,\n               header_name='feature',\n               column_names=['null_count'])\n    .sort('null_count', descending=True)\n    .with_columns((pl.col('null_count') / len(supervised_usable)).alias('null_ratio'))\n)\nplt.figure(figsize=(6, 15))\nplt.title(f'Missing values over the {len(supervised_usable)} samples which have a target')\nplt.barh(np.arange(len(missing_count)), missing_count.get_column('null_ratio'), color='coral', label='missing')\nplt.barh(np.arange(len(missing_count)), \n         1 - missing_count.get_column('null_ratio'),\n         left=missing_count.get_column('null_ratio'),\n         color='darkseagreen', label='available')\nplt.yticks(np.arange(len(missing_count)), missing_count.get_column('feature'))\nplt.gca().xaxis.set_major_formatter(PercentFormatter(xmax=1, decimals=0))\nplt.xlim(0, 1)\nplt.legend()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:36.914939Z","iopub.execute_input":"2024-09-24T19:36:36.915409Z","iopub.status.idle":"2024-09-24T19:36:38.541946Z","shell.execute_reply.started":"2024-09-24T19:36:36.915364Z","shell.execute_reply":"2024-09-24T19:36:38.540700Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The target\n\nThe target `sii` is available exactly for those participants for whom we have results of the Parent-Child Internet Addiction Test (PCIAT), and it is a function of the PCIAT total score:","metadata":{}},{"cell_type":"code","source":"# print(train.select(pl.col('PCIAT-PCIAT_Total').is_null() == pl.col('sii').is_null()).to_series().mean())\n\n(train\n .select(pl.col('PCIAT-PCIAT_Total'))\n .group_by(train.get_column('sii'))\n .agg(pl.col('PCIAT-PCIAT_Total').min().alias('PCIAT-PCIAT_Total min'),\n      pl.col('PCIAT-PCIAT_Total').max().alias('PCIAT-PCIAT_Total max'),\n      pl.col('PCIAT-PCIAT_Total').len().alias('count'))\n .sort('sii')\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:38.543948Z","iopub.execute_input":"2024-09-24T19:36:38.544494Z","iopub.status.idle":"2024-09-24T19:36:38.569757Z","shell.execute_reply.started":"2024-09-24T19:36:38.544439Z","shell.execute_reply":"2024-09-24T19:36:38.568490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The test dataset doesn't have any PCIAT columns (otherwise predictions would be trivial).\n\n**Insight:**\n1. We should focus on predicting the target from all other features except the PCIAT results.\n2. We know the target only for two thirds of the samples. The samples without target can perhaps be used for semi-supervised learning.\n3. We can directly predict `sii` (this is the value we have to submit), or we can predict `PCIAT-PCIAT_Total` and then transform this prediction to a `sii` prediction for submission. As `PCIAT-PCIAT_Total` is more granular and informative than `sii`, training to predict `PCIAT-PCIAT_Total` has the potential to produce a better model.","metadata":{}},{"cell_type":"code","source":"print('Columns missing in test:')\nprint([f for f in train.columns if f not in test.columns])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:38.571315Z","iopub.execute_input":"2024-09-24T19:36:38.571825Z","iopub.status.idle":"2024-09-24T19:36:38.581657Z","shell.execute_reply.started":"2024-09-24T19:36:38.571767Z","shell.execute_reply":"2024-09-24T19:36:38.579883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Demographics\n\nThe study participants are between 5 and 22 years old. There are twice as many boys as girls.","metadata":{}},{"cell_type":"code","source":"_, axs = plt.subplots(2, 1, sharex=True)\nfor sex in range(2):\n    ax = axs.ravel()[sex]\n    vc = train.filter(pl.col('Basic_Demos-Sex') == sex).get_column('Basic_Demos-Age').value_counts()\n    ax.bar(vc.get_column('Basic_Demos-Age'),\n           vc.get_column('count'),\n           color=['lightblue', 'coral'][sex],\n           label=['boys', 'girls'][sex])\n    ax.xaxis.set_major_locator(MaxNLocator(integer=True))\n    ax.set_ylabel('count')\n    ax.legend()\nplt.suptitle('Age distribution')\naxs.ravel()[1].set_xlabel('years')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:38.582875Z","iopub.execute_input":"2024-09-24T19:36:38.583748Z","iopub.status.idle":"2024-09-24T19:36:39.153064Z","shell.execute_reply.started":"2024-09-24T19:36:38.583683Z","shell.execute_reply":"2024-09-24T19:36:39.151817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The four seasons of enrollment have similar frequencies:","metadata":{}},{"cell_type":"code","source":"vc = train.get_column('Basic_Demos-Enroll_Season').value_counts()\nplt.pie(vc.get_column('count'), labels=vc.get_column('Basic_Demos-Enroll_Season'))\nplt.title('Season of enrollment')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-09-24T19:36:39.154848Z","iopub.execute_input":"2024-09-24T19:36:39.155345Z","iopub.status.idle":"2024-09-24T19:36:39.319087Z","shell.execute_reply.started":"2024-09-24T19:36:39.155290Z","shell.execute_reply":"2024-09-24T19:36:39.317416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Boys have a slightly higher risk of internet addiction than girls:","metadata":{}},{"cell_type":"code","source":"_, axs = plt.subplots(2, 1, sharex=True, sharey=True)\nfor sex in range(2):\n    ax = axs.ravel()[sex]\n    vc = train.filter(pl.col('Basic_Demos-Sex') == sex).get_column('sii').value_counts()\n    ax.bar(vc.get_column('sii'),\n           vc.get_column('count') / vc.get_column('count').sum(),\n           color=['lightblue', 'coral'][sex],\n           label=['boys', 'girls'][sex])\n    ax.set_xticks(np.arange(4), target_labels)\n    ax.yaxis.set_major_formatter(PercentFormatter(xmax=1, decimals=0))\n    ax.set_ylabel('count')\n    ax.legend()\nplt.suptitle('Target distribution')\naxs.ravel()[1].set_xlabel('Severity Impairment Index (sii)')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:39.321809Z","iopub.execute_input":"2024-09-24T19:36:39.322753Z","iopub.status.idle":"2024-09-24T19:36:39.927179Z","shell.execute_reply.started":"2024-09-24T19:36:39.322668Z","shell.execute_reply":"2024-09-24T19:36:39.925684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# A look at selected other features\n\nSome children had their body–mass index measured twice in the study. Comparing the two reported values gives an impression of data quality:\n","metadata":{}},{"cell_type":"code","source":"bmi_ratio = train.select(pl.col('Physical-BMI') / pl.col('BIA-BIA_BMI')).to_series()\ncolor = (bmi_ratio < 0.7) | (bmi_ratio > 1.3) # red if difference > 30 %\n\nplt.scatter(train.get_column('Physical-BMI'),\n            train.get_column('BIA-BIA_BMI'),\n            s=6,\n            cmap='coolwarm',\n            c=color)\nplt.gca().set_aspect('equal')\nplt.xlabel('Physical-BMI')\nplt.ylabel('BIA-BIA_BMI')\nplt.title('How much can the body–mass index change in a year?')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:39.931903Z","iopub.execute_input":"2024-09-24T19:36:39.932382Z","iopub.status.idle":"2024-09-24T19:36:40.213252Z","shell.execute_reply.started":"2024-09-24T19:36:39.932337Z","shell.execute_reply":"2024-09-24T19:36:40.211995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Systolic blood pressure should always be higher than diastolic:","metadata":{}},{"cell_type":"code","source":"color = train.get_column('Physical-Systolic_BP') <= train.get_column('Physical-Diastolic_BP')\n\nplt.scatter(train.get_column('Physical-Diastolic_BP'),\n            train.get_column('Physical-Systolic_BP'),\n            s=6,\n            cmap='coolwarm',\n            c=color)\nplt.gca().set_aspect('equal')\nplt.plot([0, 200], [0, 200], color='gray')\nplt.xlabel('Diastolic blood pressure')\nplt.ylabel('Systolic blood pressure')\nplt.title('What blood pressure is realistic?')\nplt.xticks(np.linspace(0, 200, 5))\nplt.yticks(np.linspace(0, 200, 5))\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:40.214779Z","iopub.execute_input":"2024-09-24T19:36:40.215152Z","iopub.status.idle":"2024-09-24T19:36:40.542487Z","shell.execute_reply.started":"2024-09-24T19:36:40.215112Z","shell.execute_reply":"2024-09-24T19:36:40.541280Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vc = train.get_column('Physical-HeartRate').value_counts()\ncolor = np.where(vc.get_column('Physical-HeartRate') < 50, 'r', 'b')\nplt.figure(figsize=(8, 2))\nplt.title('Histogram of Physical-HeartRate with outliers at the low end')\nplt.bar(vc.get_column('Physical-HeartRate'), vc.get_column('count'), color=color)\nplt.xlabel('Physical-HeartRate')\nplt.ylabel('count')\nplt.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:40.543840Z","iopub.execute_input":"2024-09-24T19:36:40.544198Z","iopub.status.idle":"2024-09-24T19:36:40.981142Z","shell.execute_reply.started":"2024-09-24T19:36:40.544159Z","shell.execute_reply":"2024-09-24T19:36:40.979928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The sleep disturbance scale questionnaire gives a raw score between 0 and 100. The raw score is then converted to a t score. It looks like you can drop one of the features without losing information.\n\nThe sleep disturbance scale t score is defined so that the average is 50 and the standard deviation is 10. We apparently have 29 children whose sleep is five standard deviations worse than the average of the general population. Is that plausible?","metadata":{}},{"cell_type":"code","source":"vc = train.get_column('SDS-SDS_Total_Raw').value_counts()\nplt.figure(figsize=(6, 2))\nplt.title('Sleep disturbance scale')\nplt.bar(vc.get_column('SDS-SDS_Total_Raw'), vc.get_column('count'), color='brown')\nplt.xlabel('SDS-SDS_Total_Raw')\nplt.ylabel('count')\nplt.show()\n\nplt.title('Sleep disturbance scale: conversion from raw to t score')\nplt.scatter(train.get_column('SDS-SDS_Total_Raw'),\n            train.get_column('SDS-SDS_Total_T'),\n            color='brown')\nplt.xlabel('SDS-SDS_Total_Raw')\nplt.ylabel('SDS-SDS_Total_T')\nplt.show()\n\nvc = train.get_column('SDS-SDS_Total_T').value_counts()\nplt.figure(figsize=(6, 2))\nplt.title('Sleep disturbance scale')\nplt.bar(vc.get_column('SDS-SDS_Total_T'), vc.get_column('count'), color='brown')\nplt.xlabel('SDS-SDS_Total_T')\nplt.ylabel('count')\nplt.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:40.982736Z","iopub.execute_input":"2024-09-24T19:36:40.983214Z","iopub.status.idle":"2024-09-24T19:36:42.019811Z","shell.execute_reply.started":"2024-09-24T19:36:40.983161Z","shell.execute_reply":"2024-09-24T19:36:42.018667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Insight:** There are many outliers. We have to decide whether to keep them, modify them ([winsorizing](https://en.wikipedia.org/wiki/Winsorizing)) or drop them.","metadata":{}},{"cell_type":"markdown","source":"# Actigraphy files (time series)\n\n> [Actigraphy](https://en.wikipedia.org/wiki/Actigraphy) is a non-invasive method of monitoring human rest/activity cycles. A small actigraph unit, also called an actimetry sensor, is worn for a week or more to measure gross motor activity. The unit is usually in a wristwatch-like package worn on the wrist. The movements the actigraph unit undergoes are continually recorded and some units also measure light exposure. (Wikipedia)\n\nWe have actigraphy files for a quarter of the participants (996 to be precise). The file name is always `part-0.parquet`. \n\nLooking at the file of participant `id=0417c91e`, a six-year old right-handed girl, we see that this participant started to use the accelerometer on a Tuesday (weekday=2) of the second quarter of the year at second 44100 of the day (12:15 PM), 5 days after the PCIAT test. She gave the accelerometer back on the 53rd day after the PCIAT test, a Monday of the third quarter, at 9:08 AM.\n\nThe competition data page says that `time_of_day` is in format `%H:%M:%S.%9f`. This is obviously not true. `time_of_day` is measured in nanoseconds since midnight.","metadata":{}},{"cell_type":"code","source":"actigraphy = pl.read_parquet('/kaggle/input/child-mind-institute-problematic-internet-use/series_train.parquet/id=0417c91e/part-0.parquet')\nactigraphy","metadata":{"execution":{"iopub.status.busy":"2024-09-24T19:36:42.021440Z","iopub.execute_input":"2024-09-24T19:36:42.021887Z","iopub.status.idle":"2024-09-24T19:36:42.131636Z","shell.execute_reply.started":"2024-09-24T19:36:42.021844Z","shell.execute_reply":"2024-09-24T19:36:42.130421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can plot diagrams of the time series in this file. What can we see?\n1. We clearly see a daily pattern.\n2. We see that the girl wore the accelerometer for 31 days and then took it off.\n3. The dataset has a non-wear_flag column, but that flag is always zero for this participant. \n4. The girl is in an environment where the illuminance exceeds 2500 [lux](https://en.wikipedia.org/wiki/Lux) every day (the device cannot measure more than 2500 lux). Such a high illuminance means that she is outdoors or in a room with huge windows.\n5. The girl moves a lot: she has enmo values above 2 almost every day.\n6. The time series usually contain measurements every 5 seconds, but some time steps are missing. It is not documented under what conditions time steps are skipped.\n\n**Insight:** \n- Use the ENMO and light columns and don't trust the non-wear_flag!\n- The ENMO and light columns offer themselves for analysis with a one-dimensional convolutional neural network, but if we want to start simple, we can use some basic aggregations (mean, variance, ...) of the time series as features for a gradient-boosting model.\n","metadata":{}},{"cell_type":"code","source":"def analyze_actigraphy(id, only_one_week=False, small=False):\n    actigraphy = pl.read_parquet(f'/kaggle/input/child-mind-institute-problematic-internet-use/series_train.parquet/id={id}/part-0.parquet')\n    day = actigraphy.get_column('relative_date_PCIAT') + actigraphy.get_column('time_of_day') / 86400e9\n    sample = train.filter(pl.col('id') == id)\n    age = sample.get_column('Basic_Demos-Age').item()\n    sex = ['boy', 'girl'][sample.get_column('Basic_Demos-Sex').item()]\n    actigraphy = (\n        actigraphy\n        .with_columns(\n            (day.diff() * 86400).alias('diff_seconds'),\n            (np.sqrt(np.square(pl.col('X')) + np.square(pl.col('Y')) + np.square(pl.col('Z'))).alias('norm'))\n        )\n    )\n\n    if only_one_week:\n        start = np.ceil(day.min())\n        mask = (start <= day.to_numpy()) & (day.to_numpy() <= start + 7*3)\n        mask &= ~ actigraphy.get_column('non-wear_flag').cast(bool).to_numpy()\n    else:\n        mask = np.full(len(day), True)\n        \n    if small:\n        timelines = [\n            ('enmo', 'forestgreen'),\n            ('light', 'orange'),\n        ]\n    else:\n        timelines = [\n            ('X', 'm'),\n            ('Y', 'm'),\n            ('Z', 'm'),\n#             ('norm', 'c'),\n            ('enmo', 'forestgreen'),\n            ('anglez', 'lightblue'),\n            ('light', 'orange'),\n            ('non-wear_flag', 'chocolate')\n    #         ('diff_seconds', 'k'),\n        ]\n        \n    _, axs = plt.subplots(len(timelines), 1, sharex=True, figsize=(12, len(timelines) * 1.1 + 0.5))\n    for ax, (feature, color) in zip(axs, timelines):\n        ax.set_facecolor('#eeeeee')\n        ax.scatter(day.to_numpy()[mask],\n                   actigraphy.get_column(feature).to_numpy()[mask],\n                   color=color, label=feature, s=1)\n        ax.legend(loc='upper left', facecolor='#eeeeee')\n        if feature == 'diff_seconds':\n            ax.set_ylim(-0.5, 20.5)\n    axs[-1].set_xlabel('day')\n    axs[-1].xaxis.set_major_locator(MaxNLocator(integer=True))\n    plt.tight_layout()\n    axs[0].set_title(f'id={id}, {sex}, age={age}')\n    plt.show()\n\nanalyze_actigraphy('0417c91e', only_one_week=False)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:46:05.432206Z","iopub.execute_input":"2024-09-24T19:46:05.432716Z","iopub.status.idle":"2024-09-24T19:46:08.339477Z","shell.execute_reply.started":"2024-09-24T19:46:05.432668Z","shell.execute_reply":"2024-09-24T19:46:08.338197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's look at a few more time series. The next time series covers only one day, I don't think it helps predict problematic Internet use:","metadata":{}},{"cell_type":"code","source":"analyze_actigraphy('5f9dddb4', small=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:46:11.722518Z","iopub.execute_input":"2024-09-24T19:46:11.723382Z","iopub.status.idle":"2024-09-24T19:46:12.322766Z","shell.execute_reply.started":"2024-09-24T19:46:11.723332Z","shell.execute_reply":"2024-09-24T19:46:12.321600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our next participant, a 15-year old left-handed boy, never moves much (enmo < 0.5), and he saw daylight only once in a whole month:","metadata":{}},{"cell_type":"code","source":"analyze_actigraphy('22375702', small=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:46:14.057346Z","iopub.execute_input":"2024-09-24T19:46:14.057824Z","iopub.status.idle":"2024-09-24T19:46:14.710881Z","shell.execute_reply.started":"2024-09-24T19:46:14.057777Z","shell.execute_reply":"2024-09-24T19:46:14.709577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next one: This poor nine-year old boy (id 0668373f) had a quiet life for two and a half weeks, then he had an accident with an acceleration of 12 g. After the accident, the device immediately stopped recording data; let's hope that the boy survived!","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"analyze_actigraphy('0668373f', small=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:46:16.169873Z","iopub.execute_input":"2024-09-24T19:46:16.170338Z","iopub.status.idle":"2024-09-24T19:46:16.980778Z","shell.execute_reply.started":"2024-09-24T19:46:16.170292Z","shell.execute_reply":"2024-09-24T19:46:16.979580Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our last time series shows some strange ramps in the illuminance measurement. It looks like there was missing data and somebody filled the blanks by linear interpolation. I'd prefer to get the raw data without undocumented preprocessing applied.","metadata":{}},{"cell_type":"code","source":"analyze_actigraphy('bc4eaf77', small=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:46:17.784737Z","iopub.execute_input":"2024-09-24T19:46:17.785211Z","iopub.status.idle":"2024-09-24T19:46:18.802553Z","shell.execute_reply.started":"2024-09-24T19:46:17.785165Z","shell.execute_reply":"2024-09-24T19:46:18.801265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Insight:** We need a lot of data cleaning before we can use the actigraphy data.\n\nDid you notice that for some participants I mentioned whether they were right- or left-handed? Did you wonder how I came to that conclusion? It's simple. The actigraphy device is fixed to the wrist of the non-dominant arm. Depending on the whether you choose the left or right wrist, the X coordinate of the accelerometer will return a positive or negative average.","metadata":{}},{"cell_type":"markdown","source":"# A simple classification model\n\nWe can model the competition as a multiclass classification task with the four classes 'None', 'Mild', 'Moderate', 'Severe'. \n\nIt is important that we select only the samples where the target `sii` is known and that we drop all PCIAT columns. We don't use the accelerometer data for this simple model.\n\nWe use the [scikit-learn implementation of the quadratic kappa score](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.cohen_kappa_score.html) for evaluation. ","metadata":{}},{"cell_type":"code","source":"y = supervised_usable.get_column('sii')\nX = supervised_usable.drop('id', 'sii', '^PCIAT.*$').to_pandas()\n\nkf = StratifiedKFold(shuffle=True, random_state=1)\noof = np.zeros(len(y), dtype=int)\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(X, y)):\n    X_tr = X.iloc[idx_tr]\n    X_va = X.iloc[idx_va]\n    y_tr = y[idx_tr]\n    y_va = y[idx_va]\n    \n    model = lightgbm.LGBMClassifier(verbose=-1)\n    model.fit(X_tr, y_tr)\n    y_pred = model.predict(X_va)\n    score = cohen_kappa_score(y_va, y_pred, weights='quadratic')\n    print(f\"# Fold {fold}: {score=:.3f}\")\n    oof[idx_va] = y_pred\n    \nscore = cohen_kappa_score(y, oof, weights='quadratic')\nprint(f\"{Fore.GREEN}{Style.BRIGHT}# Overall: {score=:.3f} (classification with LightGBM){Style.RESET_ALL}\")\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T19:36:48.363696Z","iopub.execute_input":"2024-09-24T19:36:48.364194Z","iopub.status.idle":"2024-09-24T19:36:59.535517Z","shell.execute_reply.started":"2024-09-24T19:36:48.364146Z","shell.execute_reply":"2024-09-24T19:36:59.534299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The [confusion matrix](https://en.wikipedia.org/wiki/Confusion_matrix) helps understand the output of the model.\n\n> Each row of the matrix represents the instances in an actual class while each column represents the instances in a predicted class \\[...\\]. The diagonal of the matrix therefore represents all instances that are correctly predicted. The name stems from the fact that it makes it easy to see whether the system is confusing two classes (i.e. commonly mislabeling one as another). (Wikipedia)\n\nThe most common mislabeling of our classification model is indicated by the highest non-diagonal entry of the confusion matrix: It classifies the mildly problematic Internet use of 454 children as non-problematic.","metadata":{}},{"cell_type":"code","source":"ConfusionMatrixDisplay.from_predictions(y, oof)\nplt.title('Confusion matrix for the simple classification model')\nplt.xticks(np.arange(4), target_labels)\nplt.yticks(np.arange(4), target_labels)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:36:59.537145Z","iopub.execute_input":"2024-09-24T19:36:59.537609Z","iopub.status.idle":"2024-09-24T19:36:59.865162Z","shell.execute_reply.started":"2024-09-24T19:36:59.537545Z","shell.execute_reply":"2024-09-24T19:36:59.863889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# A simple regression model is better\n\nWe can use a regression model as well, and it turns out that the regression model is better than the classification model (it gets the higher cv score). We interpret the target as a number from 0 to 3, and we round the predictions of the regression model to the nearest integer.","metadata":{}},{"cell_type":"code","source":"y = supervised_usable.get_column('sii')\nX = supervised_usable.drop('id', 'sii', '^PCIAT.*$').to_pandas()\n\nkf = StratifiedKFold(shuffle=True, random_state=1)\noof_raw = np.zeros(len(y), dtype=float) # oof predictions, before rounding\noof = np.zeros(len(y), dtype=int) # oof predictions, rounded\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(X, y)):\n    X_tr = X.iloc[idx_tr]\n    X_va = X.iloc[idx_va]\n    y_tr = y[idx_tr]\n    y_va = y[idx_va]\n\n    model = lightgbm.LGBMRegressor(verbose=-1)\n    model.fit(X_tr, y_tr.to_numpy())\n    y_pred = model.predict(X_va)\n    oof_raw[idx_va] = y_pred\n    y_pred = y_pred.round(0).astype(int)\n    score = cohen_kappa_score(y_va, y_pred, weights='quadratic')\n    print(f\"# Fold {fold}: {score=:.3f}\")\n    oof[idx_va] = y_pred\n\nscore = cohen_kappa_score(y, oof, weights='quadratic')\nprint(f\"{Fore.GREEN}{Style.BRIGHT}# Overall: {score=:.3f} (regression with LightGBM){Style.RESET_ALL}\")\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T19:36:59.866703Z","iopub.execute_input":"2024-09-24T19:36:59.867184Z","iopub.status.idle":"2024-09-24T19:37:02.759742Z","shell.execute_reply.started":"2024-09-24T19:36:59.867139Z","shell.execute_reply":"2024-09-24T19:37:02.758483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The confusion matrix shows that our regression model never predicts the class 'Severe' (the rightmost column of the matrix is zero):","metadata":{}},{"cell_type":"code","source":"regression_labels = [f\"{i} = {target_labels[i]}\" for i in range(4)]\n\nConfusionMatrixDisplay.from_predictions(y, oof)\nplt.title('Confusion matrix for the simple regression model')\nplt.xticks(np.arange(4), regression_labels)\nplt.yticks(np.arange(4), regression_labels)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-09-24T19:37:02.761089Z","iopub.execute_input":"2024-09-24T19:37:02.761574Z","iopub.status.idle":"2024-09-24T19:37:03.129879Z","shell.execute_reply.started":"2024-09-24T19:37:02.761515Z","shell.execute_reply":"2024-09-24T19:37:03.128633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Tuning the thresholds\n\nOur regression model predicts float values, and in the previous section of this notebook, we rounded these float values to integers because the Kaggle competition expects integer predictions.\n\nThere is no reason why rounding the float values to the integers 0, 1, 2 and 3 at the thresholds 0.5, 1.5 and 2.5 should give the best score. We can easily find thresholds which lead to a higher score, and the improvement is substantial:\n","metadata":{}},{"cell_type":"code","source":"def round_with_thresholds(raw_preds, thresholds):\n    \"\"\"Round the raw predictions using specified thresholds\n    \n    Parameters\n    ----------\n    raw_preds: raw predictions of the regressor, array of n_samples float values\n    thresholds: 3-element float array\n    \n    Returns\n    -------\n    rounded_preds: rounded predictions, array of n_samples int values in range 0..3\n    \"\"\"\n    return np.where(raw_preds < thresholds[0], 0,\n                    np.where(raw_preds < thresholds[1], 1,\n                             np.where(raw_preds < thresholds[2], 2, 3)))\n\n\ndef fun(thresholds, y_true, raw_preds):\n    \"\"\"Function to be minimized: negative quadratic kappa score\n    \n    Parameters:\n    thresholds: ndarray of shape (3, )\n    y_true: ndarray of shape (n_samples, )\n    raw_preds: ndarray of shape (n_samples, )\n    \n    Returns:\n    negative quadratic kappa score for the predictions rounded at the specified thresholds\n    \"\"\"\n    rounded_preds = round_with_thresholds(raw_preds, thresholds)\n    return - cohen_kappa_score(y_true, rounded_preds, weights='quadratic')\n\n# Determine the thresholds which give the highest quadratic kappa score\nres = minimize(fun, x0=[0.5, 1.5, 2.5], args=(y, oof_raw), method='Nelder-Mead')\nassert res.success\noof_tuned = round_with_thresholds(oof_raw, res.x)\nprint(f\"# Optimized thresholds: {res.x.round(2)}\")\nprint(f\"# Score with default rounding:     {cohen_kappa_score(y, oof, weights='quadratic'):.3f}\")\nprint(f\"# Score with optimized thresholds: {cohen_kappa_score(y, oof_tuned, weights='quadratic'):.3f}\")\n\nConfusionMatrixDisplay.from_predictions(y, oof_tuned)\nplt.title('Confusion matrix with tuned thresholds')\nplt.xticks(np.arange(4), target_labels)\nplt.yticks(np.arange(4), target_labels)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-09-24T19:37:03.131388Z","iopub.execute_input":"2024-09-24T19:37:03.131869Z","iopub.status.idle":"2024-09-24T19:37:03.685269Z","shell.execute_reply.started":"2024-09-24T19:37:03.131812Z","shell.execute_reply":"2024-09-24T19:37:03.684188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Conclusion and next steps\n\nWe have seen how to cross-validate a model for problematic Internet use and how to tune the rounding thresholds.\n\nThe logical next steps are:\n1. Connect the elements of the notebook so that it produces a submission file for the competition.\n2. Integrate the actigraphy time series into the model.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}