{"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"},{"sourceId":208161292,"sourceType":"kernelVersion"}],"dockerImageVersionId":30775,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# OverView\n- Target variable: PCIAT-PCIAT_Total\n- Hierarchical Bayesian Model with 5 features\n- NaN features are imutated by median of train features.\n- 7-fold Stratified KFold\n- CV: 0.441, Public: 0.447, Private: 0.473","metadata":{}},{"cell_type":"markdown","source":"# Libraries & Data","metadata":{}},{"cell_type":"markdown","source":"PyMC can be installed using a notebook created from dataset \"[PyMC Built Distribution](https://www.kaggle.com/datasets/carlosgaravatti/pymc-built-distribution)\".","metadata":{}},{"cell_type":"code","source":"!pip install --no-index -U --find-links=/kaggle/input/pymc-dependencies pymc==5.10.3","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:14:48.594284Z","iopub.execute_input":"2024-12-22T13:14:48.595274Z","iopub.status.idle":"2024-12-22T13:15:10.844969Z","shell.execute_reply.started":"2024-12-22T13:14:48.595209Z","shell.execute_reply":"2024-12-22T13:15:10.843518Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')\nimport pandas as pd\nimport numpy as np\nfrom pathlib import Path\nfrom copy import deepcopy\n\nimport plotly.graph_objects as go\nimport plotly.express as px\nfrom plotly.subplots import make_subplots\nimport plotly.io as pio\nfrom plotly.offline import plot\npio.renderers.default = \"kaggle\"\n\nfrom sklearn.metrics import cohen_kappa_score, confusion_matrix\nfrom sklearn.model_selection import StratifiedGroupKFold\nfrom sklearn.impute import SimpleImputer\n\nimport pymc as pm","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:15:11.008251Z","iopub.execute_input":"2024-12-22T13:15:11.008764Z","iopub.status.idle":"2024-12-22T13:15:18.729659Z","shell.execute_reply.started":"2024-12-22T13:15:11.008706Z","shell.execute_reply":"2024-12-22T13:15:18.728393Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"ROOT = Path('/kaggle/input/child-mind-institute-problematic-internet-use/')\nTRAIN_ACTIGRAPHY = ROOT / 'series_train.parquet'\nTEST_ACTIGRAPHY = ROOT / 'series_test.parquet'","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:08.436177Z","iopub.execute_input":"2024-12-22T13:25:08.436976Z","iopub.status.idle":"2024-12-22T13:25:08.442488Z","shell.execute_reply.started":"2024-12-22T13:25:08.436937Z","shell.execute_reply":"2024-12-22T13:25:08.441437Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train = pd.read_csv(ROOT / 'train.csv')\ntrain = train.dropna(subset=['sii'])\ntrain","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:08.815780Z","iopub.execute_input":"2024-12-22T13:25:08.816180Z","iopub.status.idle":"2024-12-22T13:25:08.950177Z","shell.execute_reply.started":"2024-12-22T13:25:08.816144Z","shell.execute_reply":"2024-12-22T13:25:08.949094Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test = pd.read_csv(ROOT / 'test.csv')\ntest","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:09.216805Z","iopub.execute_input":"2024-12-22T13:25:09.217190Z","iopub.status.idle":"2024-12-22T13:25:09.256333Z","shell.execute_reply.started":"2024-12-22T13:25:09.217158Z","shell.execute_reply":"2024-12-22T13:25:09.255240Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# EDA","metadata":{}},{"cell_type":"markdown","source":"QWK is calculated by \"sii\", but \"sii\" is derived from \"PCIAT-PCIAT_Total\".\nTherefore, I focused on characteristics of \"PCIAT-PCAIT_Total\", which are a distribution or a correlation with othre variables.","metadata":{}},{"cell_type":"markdown","source":"## Distribution of PCIAT-PCIAT_Total","metadata":{}},{"cell_type":"markdown","source":"It seems to be a censored normal distribution or zero-inflated poisson distribution.","metadata":{}},{"cell_type":"code","source":"# 全体の分布\nfig = px.histogram(\n    train, x='PCIAT-PCIAT_Total'\n)\nfig.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:15.270823Z","iopub.execute_input":"2024-12-22T13:25:15.271199Z","iopub.status.idle":"2024-12-22T13:25:17.071636Z","shell.execute_reply.started":"2024-12-22T13:25:15.271168Z","shell.execute_reply":"2024-12-22T13:25:17.070603Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Demographic effect on PCIAT-PCIAT_Total","metadata":{}},{"cell_type":"markdown","source":"As age increases up to around 13 years old, the median PCIAT score also rises. This is likely due to easier access to the internet and longer usage history, which may contribute to higher PCIAT scores.<br>\nFor individuals over 18 years old, the sample size is too small to capture clear distribution trends<br>\nGender differences are generally observed across age groups, with males tending to have higher scores than females.<br>\nBased on these observations, categorizing age and considering gender differences are deemed necessary for modeling.","metadata":{}},{"cell_type":"code","source":"fig = px.box(\n    train.assign(Sex=lambda df: df['Basic_Demos-Sex'].map({0: 'Male', 1: 'Female'})),\n    x='Basic_Demos-Age', y='PCIAT-PCIAT_Total', points='all', color='Sex'\n)\nfig.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:17.073522Z","iopub.execute_input":"2024-12-22T13:25:17.073976Z","iopub.status.idle":"2024-12-22T13:25:17.184443Z","shell.execute_reply.started":"2024-12-22T13:25:17.073930Z","shell.execute_reply":"2024-12-22T13:25:17.183363Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Distribution of PCIAT-PCIAT_Total by age category and gender","metadata":{}},{"cell_type":"markdown","source":"Up to the age of 13, age is categorized in 3-year intervals, while ages 14 and above are grouped into a single category.<br>\nThe younger age categories tend to exhibit distributions that appear truncated at zero.","metadata":{}},{"cell_type":"code","source":"def _categorize_age(x):\n    if x < 14:\n        return (x - 5) // 3\n    else:\n        return 3\n\ntmp = train.copy()\n# Age Category\ntmp['Age_Category'] = tmp['Basic_Demos-Age'].apply(_categorize_age)\nnum_age_categories = int(tmp['Age_Category'].nunique())\n\nfig = make_subplots(\n    rows=num_age_categories, cols=1,\n    shared_xaxes=True, vertical_spacing=0.05,\n    subplot_titles=[f'Category {age_cat}' for age_cat in range(num_age_categories)]\n)\nfor age_cat in range(num_age_categories):\n    tmp_age = tmp.query(f'Age_Category=={age_cat}')\n    fig.add_trace(\n        go.Histogram(\n            x=tmp_age[tmp_age['Basic_Demos-Sex']==0]['PCIAT-PCIAT_Total'],\n            marker_color='#636EFA',\n            opacity=0.6,\n            showlegend=True if age_cat==0 else False,\n            legendgroup=\"Male\",\n            legendgrouptitle_text=\"Male\",\n            name=\"\",\n        ),\n        row=age_cat+1, col=1\n    )\n    fig.add_trace(\n        go.Histogram(\n            x=tmp_age[tmp_age['Basic_Demos-Sex']==1]['PCIAT-PCIAT_Total'],\n            marker_color='#EF553B',\n            showlegend=True if age_cat==0 else False,\n            legendgroup=\"Female\",\n            legendgrouptitle_text=\"Female\",\n            name=\"\",\n        ),\n        row=age_cat+1, col=1\n    )\nfig.update_layout(\n    width=900, height=600,\n    barmode='overlay'\n)\nfig.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:21.484970Z","iopub.execute_input":"2024-12-22T13:25:21.485337Z","iopub.status.idle":"2024-12-22T13:25:21.658089Z","shell.execute_reply.started":"2024-12-22T13:25:21.485307Z","shell.execute_reply":"2024-12-22T13:25:21.656890Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Correlation analysis","metadata":{}},{"cell_type":"markdown","source":"The actigraph data has a high number of missing values, so it will not be used in this model.","metadata":{}},{"cell_type":"code","source":"# Variables that will be used for calculating correlation coefficient\nnum_features = list(train.select_dtypes(include=['number']).columns)\nnum_features = [feature for feature in num_features if 'PCIAT' not in feature and feature!='sii']\nnum_features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:26.437640Z","iopub.execute_input":"2024-12-22T13:25:26.438609Z","iopub.status.idle":"2024-12-22T13:25:26.448477Z","shell.execute_reply.started":"2024-12-22T13:25:26.438567Z","shell.execute_reply":"2024-12-22T13:25:26.447057Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Physical-Height showed the highest correlation coefficient, but this is likely a pseudo-correlation due to age.","metadata":{}},{"cell_type":"code","source":"corr_df = (\n    train[['PCIAT-PCIAT_Total']+num_features]\n    .corr()\n    .loc[num_features, ['PCIAT-PCIAT_Total']]\n    .reset_index()\n    .rename(columns={'index': 'feature', 'PCIAT-PCIAT_Total': 'corr'})  \n)\ncorr_df = corr_df.sort_values(['corr'])\n\nfig = px.bar(corr_df, x='feature', y='corr')\nfig.update_layout(\n    height=500,\n    yaxis_range=(-1, 1)\n)\nfig.show()\n\ncorr_df.round(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:29.600623Z","iopub.execute_input":"2024-12-22T13:25:29.601850Z","iopub.status.idle":"2024-12-22T13:25:29.708976Z","shell.execute_reply.started":"2024-12-22T13:25:29.601808Z","shell.execute_reply":"2024-12-22T13:25:29.707783Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"When calculating the partial correlation coefficient with age as a covariate, a correlation over 0.2 was observed with only PreInt_EduHx-computerinternet_hoursday and SDS-SDS_Total_T/Raw.","metadata":{}},{"cell_type":"code","source":"def partial_corr(cov_matrix):\n    cov_XY = cov_matrix[0, 1]\n    cov_XZ = cov_matrix[0, 2]\n    cov_YZ = cov_matrix[1, 2]\n    var_Z = cov_matrix[2, 2]\n    \n    numerator = cov_XY - (cov_XZ * cov_YZ / var_Z)\n    denominator = np.sqrt(\n        (cov_matrix[0, 0] - cov_XZ**2 / var_Z) *\n        (cov_matrix[1, 1] - cov_YZ**2 / var_Z)\n    )\n    return numerator / denominator\n    \npcorr_df = pd.DataFrame()\nfor feature in num_features:\n    if feature!='Basic_Demos-Age':\n        cov_matrix = train[['PCIAT-PCIAT_Total', feature, 'Basic_Demos-Age']].cov(ddof=0).values\n        \n        result = partial_corr(cov_matrix)\n        pcorr_df = pd.concat(\n            [pcorr_df, pd.DataFrame([[feature, result]], columns=['feature', 'pcorr'])], \n            ignore_index=True)\n\npcorr_df = pcorr_df.sort_values(['pcorr'])\n\nfig = px.bar(pcorr_df, x='feature', y='pcorr')\nfig.update_layout(\n    height=500,\n    yaxis_range=(-1, 1)\n)\nfig.show()\npcorr_df.round(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:35.370762Z","iopub.execute_input":"2024-12-22T13:25:35.371918Z","iopub.status.idle":"2024-12-22T13:25:35.498047Z","shell.execute_reply.started":"2024-12-22T13:25:35.371872Z","shell.execute_reply":"2024-12-22T13:25:35.496809Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"SDS variabels have strong correlation. Therefore, SDS-SDS_Total_T is only selected.","metadata":{}},{"cell_type":"code","source":"train[['SDS-SDS_Total_T', 'SDS-SDS_Total_Raw']].corr().round(4)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:49.587931Z","iopub.execute_input":"2024-12-22T13:25:49.588626Z","iopub.status.idle":"2024-12-22T13:25:49.600911Z","shell.execute_reply.started":"2024-12-22T13:25:49.588586Z","shell.execute_reply":"2024-12-22T13:25:49.599754Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig = px.scatter(\n    tmp.sort_values(['Age_Category']).astype({'Age_Category': 'str'}),\n    x='SDS-SDS_Total_T', y='PCIAT-PCIAT_Total',\n    color='Age_Category', trendline='ols'\n)\nfig.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:52.230039Z","iopub.execute_input":"2024-12-22T13:25:52.230474Z","iopub.status.idle":"2024-12-22T13:25:52.331888Z","shell.execute_reply.started":"2024-12-22T13:25:52.230430Z","shell.execute_reply":"2024-12-22T13:25:52.330736Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig = px.box(\n    tmp.sort_values(['Age_Category']),\n    x='PreInt_EduHx-computerinternet_hoursday', y='PCIAT-PCIAT_Total',\n    color='Age_Category'\n)\nfig.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:53.354759Z","iopub.execute_input":"2024-12-22T13:25:53.355213Z","iopub.status.idle":"2024-12-22T13:25:53.619900Z","shell.execute_reply.started":"2024-12-22T13:25:53.355173Z","shell.execute_reply":"2024-12-22T13:25:53.618775Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Although no correlation was observed with the CGAS score, a relationship with the target variable may still exist, so the scatter plot will be examined.<br>\nThe relationship with the target variable is weak up to a CGAS score of around 80, but a negative correlation is observed above 80. Therefore, a binary feature based on a threshold of 80 will be created for the CGAS score.","metadata":{}},{"cell_type":"code","source":"fig = px.scatter(\n    tmp.sort_values(['Age_Category']).astype({'Age_Category': 'str'}),\n    x='CGAS-CGAS_Score', y='PCIAT-PCIAT_Total',\n    color='Age_Category', trendline='ols'\n)\nfig.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:25:57.710100Z","iopub.execute_input":"2024-12-22T13:25:57.710550Z","iopub.status.idle":"2024-12-22T13:25:57.810601Z","shell.execute_reply.started":"2024-12-22T13:25:57.710510Z","shell.execute_reply":"2024-12-22T13:25:57.809481Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Cross Validation","metadata":{}},{"cell_type":"markdown","source":"**Target Variables**<br>\nPCIAT-PCIAT_Total\n\n**Model**<br>\nHierarchical Bayesian Model\n\n**Features**<br>\n- Age Category (new)\n- Basic_Demos-Sex\n- SDS-SDS_Total_T\n- PreInt_EduHx-computerinternet_hoursday\n- Binary CGAS Score (new)\n\n**Data Split**<br>\n- Stratified(Group)KFold<br>\nIn this dataset, there were no individuals with multiple data points, but in medical datasets, such cases are common. Therefore, as a reference, StratifiedGroupKFold is used.\n- N Folds=7\n\n**Thresholds for Deciding sii**<br>\nRandom Grid Search<br>\nRandomly generate three thresholds, repeatedly calculate the Quadratic Weighted Kappa (QWK) on the training data, and select the threshold that maximizes the QWK for the test data.","metadata":{}},{"cell_type":"code","source":"N_FOLDS = 7\nSEED = 42","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:26:27.120511Z","iopub.execute_input":"2024-12-22T13:26:27.121364Z","iopub.status.idle":"2024-12-22T13:26:27.125964Z","shell.execute_reply.started":"2024-12-22T13:26:27.121316Z","shell.execute_reply":"2024-12-22T13:26:27.124776Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def quadratic_weighted_kappa(y_true, y_pred):\n    return cohen_kappa_score(y_true, y_pred, weights='quadratic')\n\ndef threshold_rounder(oof_non_rounded, thresholds):\n    return np.where(oof_non_rounded < thresholds[0], 0,\n                    np.where(oof_non_rounded < thresholds[1], 1,\n                             np.where(oof_non_rounded < thresholds[2], 2, 3)))\n\ndef optimize_thresholds(train_predictions, true_classes, n_iter=1000, lower=20, upper=70):\n    np.random.seed(SEED)\n    \n    best_kappa = -np.inf\n    best_thresholds = None\n\n    for _ in range(n_iter):\n        thresholds = np.sort(np.random.uniform(lower, upper, size=3))\n        predicted_classes = threshold_rounder(train_predictions, thresholds)\n\n        kappa = quadratic_weighted_kappa(true_classes, predicted_classes)\n        if kappa > best_kappa:\n            best_kappa = kappa\n            best_thresholds = thresholds\n\n    return best_thresholds, best_kappa\n\ndef feature_engineering(df):\n    df['Age_Category'] = df['Basic_Demos-Age'].apply(_categorize_age)\n    df['Binary_CGAS_Score'] = df['CGAS-CGAS_Score'].apply(lambda x: 1 if x>=80 else 0) # NaN is changed into 0.\n\n    return df\n\ndef model_pipeline(\n    train_df,\n    test_df,\n    mode='cv',\n    fillna_features=['SDS-SDS_Total_T', 'PreInt_EduHx-computerinternet_hoursday']\n):\n    mutable = True\n    if mode=='train':\n        mutable = False\n    # Feature Engineering\n    train_df = feature_engineering(train_df)\n    test_df = feature_engineering(test_df)\n\n    # Filling NaN by median by age category\n    imputers_by_age_cat = []\n    age_categories = train_df['Age_Category'].sort_values().unique()\n    for age_cat in age_categories:\n        imputer = SimpleImputer(strategy='median')\n        train_idx_by_age_cat = train_df.query(f'Age_Category=={age_cat}').index\n        test_idx_by_age_cat = test_df.query(f'Age_Category=={age_cat}').index\n        train_df.loc[train_idx_by_age_cat, fillna_features] = imputer.fit_transform(train_df.loc[train_idx_by_age_cat, fillna_features])\n        test_df.loc[test_idx_by_age_cat, fillna_features] = imputer.transform(test_df.loc[test_idx_by_age_cat, fillna_features])\n        imputers_by_age_cat.append(imputer)\n\n    # Modeling\n    n_groups = len(age_categories)\n    with pm.Model() as model:\n        model.add_coord('data_ind', values=range(train_df.shape[0]), mutable=mutable)\n        age_category = pm.Data('age_category', train_df['Age_Category'].values, dims='data', mutable=mutable)\n        pre_int = pm.Data('pre_int', train_df['PreInt_EduHx-computerinternet_hoursday'].values, dims='data', mutable=mutable)\n        sds_total = pm.Data('sds_total', train_df['SDS-SDS_Total_T'].values, dims='data', mutable=mutable)\n        sex = pm.Data('sex', train_df['Basic_Demos-Sex'].values, dims='data', mutable=mutable)\n        cgas = pm.Data('cgas', train_df['Binary_CGAS_Score'].values, dims='data', mutable=mutable)\n        obs = pm.Data('obs', train_df['PCIAT-PCIAT_Total'].values, dims='data', mutable=mutable)\n        \n        mu_pop = pm.Normal(\"mu_pop\", mu=10, sigma=10)\n        sigma_pop = pm.HalfNormal(\"sigma_pop\", sigma=10)\n\n        mu_age = pm.Normal(\"mu_age\", mu=mu_pop, sigma=sigma_pop, shape=n_groups)\n        sigma_age = pm.HalfNormal(\"sigma_age\", sigma=5, shape=n_groups)\n\n        beta_preint_eduhx_age = pm.Normal(\"beta_preint_eduhx_age\", mu=0, sigma=10, shape=n_groups)\n        beta_sds_age = pm.Normal(\"beta_sds_age\", mu=0, sigma=10, shape=n_groups)\n        beta_cgas_age = pm.Normal(\"beta_cgas_age\", mu=0, sigma=10, shape=n_groups)\n        beta_gender = pm.Normal(\"beta_gender\", mu=0, sigma=10)\n        \n        mu = (\n            mu_age[age_category] +\n            beta_preint_eduhx_age[age_category] * pre_int +\n            beta_sds_age[age_category] * sds_total +\n            beta_cgas_age[age_category] * cgas +\n            beta_gender * sex \n        )\n\n        normal_dist = pm.Normal.dist(mu=mu, sigma=sigma_age[age_category])\n        pciat_score_obs = pm.Censored(\n            \"pciat_score_obs\",\n            normal_dist, lower=0, upper=100,\n            observed=obs\n        )\n\n        trace = pm.sample(\n            draws=1500,       \n            tune=500,\n            chains=2,\n            cores=4,\n            return_inferencedata=True,\n            progressbar=True, \n            random_seed=SEED\n        )\n\n    # Sampling posterior predictive of Train/Test Data\n    pre_trace = deepcopy(trace)\n    save_trace = deepcopy(trace)\n    if mode=='cv':\n        with model:\n            # Train\n            pm.set_data(\n                {\n                    'obs': np.zeros(train_df.shape[0])\n                }\n            )\n            pre_trace.extend(pm.sample_posterior_predictive(pre_trace, random_seed=SEED))\n            \n            # Test\n            pm.set_data(\n                {\n                    'age_category': test_df['Age_Category'].values,\n                    'pre_int': test_df['PreInt_EduHx-computerinternet_hoursday'].values,\n                    'sds_total': test_df['SDS-SDS_Total_T'].values,\n                    'sex': test_df['Basic_Demos-Sex'].values,\n                    'cgas': test_df['Binary_CGAS_Score'].values,\n                    'obs': np.zeros(test_df.shape[0])\n                }\n                , \n                coords={\n                    \"data_ind\": (np.arange(test_df.shape[0]))+train_df.shape[0]\n                }\n            )\n            trace.extend(pm.sample_posterior_predictive(trace, random_seed=SEED))\n    \n        # Deciding prediction score by aggrigating sample values\n        train_df['pred_score'] = pre_trace.posterior_predictive[\"pciat_score_obs\"].mean(dim=[\"chain\", \"draw\"]).values\n        test_df['pred_score'] = trace.posterior_predictive[\"pciat_score_obs\"].mean(dim=[\"chain\", \"draw\"]).values\n    \n        # Thresholds which change PCIAT-PCAIT_Total into sii\n        thresholds, _ = optimize_thresholds(train_df['pred_score'], train_df['sii']) \n        train_df['pred'] = threshold_rounder(train_df['pred_score'].values, thresholds)\n        test_df['pred'] = threshold_rounder(test_df['pred_score'].values, thresholds)\n\n        # CV Results\n        print('Tresholds', thresholds)\n        print('Train score:', round(quadratic_weighted_kappa(train_df['sii'], train_df['pred']), 4))\n        print('Test score:', round(quadratic_weighted_kappa(test_df['sii'], test_df['pred']), 4))\n        print(confusion_matrix(test_df['sii'], test_df['pred']))\n\n        return model, save_trace, thresholds, imputers_by_age_cat, train_df, test_df\n        \n    elif mode=='train':\n        return model, save_trace","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T14:32:36.282622Z","iopub.execute_input":"2024-12-22T14:32:36.283145Z","iopub.status.idle":"2024-12-22T14:32:36.314743Z","shell.execute_reply.started":"2024-12-22T14:32:36.283104Z","shell.execute_reply":"2024-12-22T14:32:36.313512Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Data Split\ncv = train.copy().reset_index(drop=True)\ncv['sii'] = cv['sii'].astype('int')\ncv['fold'] = -1\n\nsgkf = StratifiedGroupKFold(n_splits=N_FOLDS)\nfor i, (train_index, test_index) in enumerate(sgkf.split(cv, cv['sii'], cv['id']), start=1):\n    cv.loc[test_index, 'fold'] = i\n    \ndisplay(pd.crosstab(cv['fold'], cv['sii'], margins=True))\ndisplay(pd.crosstab(cv['fold'], cv['Basic_Demos-Age'], margins=True))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:29:52.269309Z","iopub.execute_input":"2024-12-22T13:29:52.269733Z","iopub.status.idle":"2024-12-22T13:29:53.963105Z","shell.execute_reply.started":"2024-12-22T13:29:52.269684Z","shell.execute_reply":"2024-12-22T13:29:53.962029Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Cross Validation\ncv_res = pd.DataFrame()\ntrained_models = []\nfor test_fold in range(1, N_FOLDS + 1):\n    print('Fold', test_fold)\n    cv_train = cv.query(f'fold!={test_fold}')\n    cv_test = cv.query(f'fold=={test_fold}')\n\n    model, save_trace, thresholds, imputers_by_age_cat, cv_train, cv_test = model_pipeline(cv_train, cv_test)\n\n    trained_models.append([model, save_trace, thresholds, imputers_by_age_cat])\n\n    cv_res = pd.concat(\n        [\n            cv_res,\n            cv_test\n        ],\n        ignore_index=True\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:30:00.696585Z","iopub.execute_input":"2024-12-22T13:30:00.697030Z","iopub.status.idle":"2024-12-22T13:43:43.361697Z","shell.execute_reply.started":"2024-12-22T13:30:00.696993Z","shell.execute_reply":"2024-12-22T13:43:43.360215Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# QWK\ncv_metrics_source_df = cv_res.groupby('fold').apply(\n    lambda df: pd.Series(\n        {\n            'quadratic_kappa': quadratic_weighted_kappa(df['sii'], df['pred'])\n        }\n    )\n).reset_index()\ndisplay(cv_metrics_source_df.round(4))\nprint('CV Score:', round(cv_metrics_source_df['quadratic_kappa'].mean(), 4))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:43:43.364158Z","iopub.execute_input":"2024-12-22T13:43:43.364539Z","iopub.status.idle":"2024-12-22T13:43:43.399945Z","shell.execute_reply.started":"2024-12-22T13:43:43.364502Z","shell.execute_reply":"2024-12-22T13:43:43.398744Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Prediction of test data","metadata":{}},{"cell_type":"markdown","source":"In the seven models built using cross-validation, predictions are made on the test data for each model. The final prediction is obtained by averaging and rounding the results of these predictions.<br>\nIn the submitted notebook, the seed value was not fixed when sampling the posterior distribution for the test data, resulting in slight variations in the Public and Private scores.<br>\nBelow, I provide the score obtained when the seed value was changed in three different ways and the notebook was resubmitted.<br>\n||Public|Private|\n|---|---:|---:|\n|Submission|0.447|0.473|\n|Pattern 1(this notebook)|0.443|0.476|\n|Pattern 2|0.439|0.476|\n|Pattern 3|0.442|0.473|","metadata":{}},{"cell_type":"code","source":"def prediction(\n    train_df,\n    test_df,\n    model, \n    trace, \n    imputers_by_age_cat,\n    seed,\n    fillna_features=['SDS-SDS_Total_T', 'PreInt_EduHx-computerinternet_hoursday']\n):\n    # Feature Engineering\n    test_df = feature_engineering(test_df)\n\n    # Filling NaN by median by age category\n    for age_cat, imputer in enumerate(imputers_by_age_cat):\n        test_idx_by_age_cat = test_df.query(f'Age_Category=={age_cat}').index\n        if len(test_idx_by_age_cat)>0:\n            test_df.loc[test_idx_by_age_cat, fillna_features] = imputer.transform(test_df.loc[test_idx_by_age_cat, fillna_features])\n\n    with model:\n        pm.set_data(\n            {\n                'age_category': test_df['Age_Category'].values,\n                'pre_int': test_df['PreInt_EduHx-computerinternet_hoursday'].values,\n                'sds_total': test_df['SDS-SDS_Total_T'].values,\n                'sex': test_df['Basic_Demos-Sex'].values,\n                'cgas': test_df['Binary_CGAS_Score'].values,\n                'obs': np.zeros(test_df.shape[0])\n            }\n            , \n            coords={\n                \"data_ind\": (np.arange(test_df.shape[0]))+train_df.shape[0]\n            }\n        )\n        trace.extend(pm.sample_posterior_predictive(trace, random_seed=seed))\n\n    test_df['pred_score'] = trace.posterior_predictive[\"pciat_score_obs\"].mean(dim=[\"chain\", \"draw\"]).values\n    \n    return threshold_rounder(test_df['pred_score'].values, thresholds)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:43:43.401508Z","iopub.execute_input":"2024-12-22T13:43:43.401995Z","iopub.status.idle":"2024-12-22T13:43:43.412971Z","shell.execute_reply.started":"2024-12-22T13:43:43.401935Z","shell.execute_reply":"2024-12-22T13:43:43.411758Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"np.random.seed(SEED)\ntest_pred_sii = np.zeros([test.shape[0], N_FOLDS])\nfor fold, (model, trace, thresholds, imputers_by_age_cat) in enumerate(trained_models, start=1):\n    tmp_test = test.copy()\n    cv_train = cv.query(f'fold!={fold}')\n    seed = np.random.randint(0, 1001)\n    test_pred_sii[:, fold - 1] = prediction(cv_train, tmp_test, model, trace, imputers_by_age_cat, seed)\n\nsubmission = test[['id']]\nsubmission['sii'] = np.int32(np.rint(test_pred_sii.mean(axis=1)))\nsubmission.to_csv('submission.csv', index=False)\nsubmission","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:43:43.415574Z","iopub.execute_input":"2024-12-22T13:43:43.416098Z","iopub.status.idle":"2024-12-22T13:43:47.377179Z","shell.execute_reply.started":"2024-12-22T13:43:43.416046Z","shell.execute_reply":"2024-12-22T13:43:47.376122Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Model Structure & Parameters","metadata":{}},{"cell_type":"code","source":"# Training Model by all train dataset\ntmp_test = test.copy()\nmodel, save_trace = model_pipeline(train, tmp_test, mode='train')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T14:32:54.018055Z","iopub.execute_input":"2024-12-22T14:32:54.018526Z","iopub.status.idle":"2024-12-22T14:34:26.368546Z","shell.execute_reply.started":"2024-12-22T14:32:54.018489Z","shell.execute_reply":"2024-12-22T14:34:26.367258Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pm.model_to_graphviz(\n    model,\n    \n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T14:40:03.340350Z","iopub.execute_input":"2024-12-22T14:40:03.340755Z","iopub.status.idle":"2024-12-22T14:40:03.547686Z","shell.execute_reply.started":"2024-12-22T14:40:03.340719Z","shell.execute_reply":"2024-12-22T14:40:03.546178Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import arviz as az\nimport matplotlib.pyplot as plt\naz.plot_posterior(save_trace)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:55:14.622207Z","iopub.execute_input":"2024-12-22T13:55:14.622665Z","iopub.status.idle":"2024-12-22T13:55:17.826404Z","shell.execute_reply.started":"2024-12-22T13:55:14.622632Z","shell.execute_reply":"2024-12-22T13:55:17.825341Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"az.plot_trace(save_trace, legend=True)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-22T13:59:50.064306Z","iopub.execute_input":"2024-12-22T13:59:50.064845Z","iopub.status.idle":"2024-12-22T13:59:53.103833Z","shell.execute_reply.started":"2024-12-22T13:59:50.064804Z","shell.execute_reply":"2024-12-22T13:59:53.102651Z"}},"outputs":[],"execution_count":null}]}