{"cells":[{"metadata":{},"cell_type":"markdown","source":"Did you know that sklearn's BayesianRidge predict method has ```return_std=True``` where the standard deviation of the prediction is returned?\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Libraries","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import cv2\nimport pydicom\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport os\nimport math\nimport random\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\nfrom PIL import Image\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn import decomposition\nfrom typing import List, NoReturn, Union, Tuple, Optional, Text, Generic, Callable, Dict\nfrom sklearn.preprocessing import StandardScaler, MinMaxScaler, OneHotEncoder, QuantileTransformer\nfrom sklearn.model_selection import KFold, StratifiedKFold, TimeSeriesSplit\nimport itertools\nfrom sklearn.linear_model import BayesianRidge\n\n# visualize\nimport matplotlib.pyplot as plt\nimport matplotlib.style as style\nimport seaborn as sns\nfrom matplotlib_venn import venn2\nfrom matplotlib import pyplot\nfrom matplotlib.ticker import ScalarFormatter\nsns.set_context(\"talk\")\nstyle.use('fivethirtyeight')\npd.options.display.max_columns = None\n\nimport warnings\nwarnings.filterwarnings('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Config","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# CONFIG\nINPUT_DIR = \"../input/osic-pulmonary-fibrosis-progression\"\nSEED = 42\nNFOLD = 10\nSCALER = 'MinMax'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_tabular():\n    train = pd.read_csv(INPUT_DIR + '/train.csv')\n    test = pd.read_csv(INPUT_DIR + '/test.csv')\n    sub = pd.read_csv(INPUT_DIR + '/sample_submission.csv')\n    return train, test, sub\ntrain, test, sub = read_tabular()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train.shape)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(test.shape)\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nprint(sub.shape)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Feature engineering","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.merge(sub[['Patient','Weeks','Confidence','Patient_Week']], test.drop(columns=['Weeks']), on='Patient')\ntrain['where'] = 'train'\ntest['where'] = 'test'\nsub['where'] = 'sub'\ndata = pd.concat([train, test, sub], ignore_index=True)\nprint(data.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# construct train input\ndef fe(data):\n    data['min_week'] = data['Weeks']\n    data.loc[data['where'] == 'test','min_week'] = np.nan\n    data['min_week'] = data.groupby('Patient')['min_week'].transform('min')\n    \n    base = data.loc[data.Weeks == data.min_week]\n    base = base[['Patient','FVC', 'Percent']].copy()\n    base.columns = ['Patient','base_FVC', 'base_Percent']\n    base['nb'] = 1\n    base['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\n    base = base[base.nb==1]\n    base.drop('nb', axis=1, inplace=True)\n\n    data = data.merge(base, on='Patient', how='left')\n    data['base_week'] = data['Weeks'] - data['min_week']\n    del base\n    \n    train = data.loc[data['where'] == 'train', :].reset_index(drop=True)\n    test = data.loc[data['where'] == 'test', :].reset_index(drop=True)\n    sub = data.loc[data['where'] == 'sub', :].reset_index(drop=True)\n\n    return train, test, sub\ntrain, test, sub = fe(data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train.shape)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = sub\nprint(test.shape)\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# EDA","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"venn2([set(train['Patient'].values.tolist()), set(test['Patient'].values.tolist())])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot weeks vs Percent (and FVC)\ndef plot_weeks_vs(patient : str):\n    fig, ax = plt.subplots(1, 1, figsize=(7, 4))\n    colors = sns.color_palette('deep', 3)\n    weeks = train.loc[train['Patient'] == patient, 'Weeks'].values\n    \n    ax.plot(weeks, train.loc[train['Patient'] == patient, 'Percent'].values, '-o', color=colors[0], alpha=0.4)\n    ax.set_ylabel('Percent', color=colors[0])\n    ax.set_xlabel('Weeks')\n    ax.tick_params(axis='y', labelcolor=colors[0])\n    ax.set_title(patient)\n    \n    ax2 = ax.twinx()\n    ax2.plot(weeks+1, train.loc[train['Patient'] == patient, 'FVC'].values, '-^', color=colors[1], alpha=0.4)\n    ax2.plot(weeks[-3:]+1, train.loc[train['Patient'] == patient, 'FVC'].values[-3:], '-s', color=colors[2])\n    ax2.set_ylabel('FVC', color=colors[1])\n    ax2.tick_params(axis='y', labelcolor=colors[1])\n    fig.tight_layout()  # otherwise the right y-label is slightly clipped\n    plt.show()\n    \nplot_weeks_vs(train['Patient'].unique()[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_weeks_vs(train['Patient'].unique()[1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_weeks_vs(test['Patient'].unique()[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_weeks_vs(test['Patient'].unique()[2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Given there is a super-high positive correlation between FVC and Percent as a function of Weeks, I would simply use features related to weeks, FVC, and Percent.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Fitting","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Here is the groupkfold with shuffle=True, which cannot be done in the native sklearn groupkfold method.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import random\nfrom collections import Counter, defaultdict\nfrom sklearn import model_selection\n\n# ---- GroupKFold ----\nclass GroupKFold(object):\n    \"\"\"\n    GroupKFold with random shuffle with a sklearn-like structure\n    \"\"\"\n\n    def __init__(self, n_splits=4, shuffle=True, random_state=42):\n        self.n_splits = n_splits\n        self.shuffle = shuffle\n        self.random_state = random_state\n\n    def get_n_splits(self, X=None, y=None, group=None):\n        return self.n_splits\n\n    def split(self, X, y, group):\n        kf = model_selection.KFold(n_splits=self.n_splits, shuffle=self.shuffle, random_state=self.random_state)\n        unique_ids = X[group].unique()\n        for fold, (tr_group_idx, va_group_idx) in enumerate(kf.split(unique_ids)):\n            # split group\n            tr_group, va_group = unique_ids[tr_group_idx], unique_ids[va_group_idx]\n            train_idx = np.where(X[group].isin(tr_group))[0]\n            val_idx = np.where(X[group].isin(va_group))[0]\n            yield train_idx, val_idx","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# to normal\nID = 'Patient_Week'\ntarget = 'FVC'\ngroup = 'Patient'\nfeatures = ['Weeks', 'Percent', 'min_week', 'base_FVC', 'base_Percent', 'base_week']\n\n# starndarditest\nif SCALER == \"MinMax\":\n    scaler = MinMaxScaler()\nelif SCALER == \"Standard\":\n    scaler = StandardScaler()\ntrain[features] = scaler.fit_transform(train[features])\ntest[features] = scaler.transform(test[features])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(len(features))\nfeatures","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fitting","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def lb_metric(data, oof):\n    data['FVC_pred'] = oof[:, 0]\n    data['Confidence'] = oof[:, -1]\n    data['sigma_clipped'] = data['Confidence'].apply(lambda x: max(x, 70))\n    data['diff'] = abs(data['FVC'] - data['FVC_pred'])\n    data['delta'] = data['diff'].apply(lambda x: min(x, 1000))\n    data['score'] = -math.sqrt(2)*data['delta']/data['sigma_clipped'] - np.log(math.sqrt(2)*data['sigma_clipped'])\n    score = data['score'].mean()\n    return score","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nypred = np.zeros((test.shape[0], 2))\noof = np.zeros((train.shape[0], 2))\nkf = GroupKFold(n_splits=NFOLD, shuffle=True, random_state=SEED)\nkf = kf.split(train, train[target], group)\n\nfor cnt, (tr_idx, val_idx) in enumerate(kf):\n    print(f\"FOLD {cnt}\")\n    \n    # fit\n    model = BayesianRidge()\n    model.fit(train[features].values[tr_idx, :], train[target].values[tr_idx])\n    \n    # evaluate\n    yme, ystd = model.predict(train[features].values[val_idx, :], return_std=True) # return_std=True!\n    oof[val_idx, 0] = yme \n    oof[val_idx, 1] = 2 * ystd\n    yme, ystd = model.predict(test[features].values, return_std=True)\n    ypred[:, 0] += yme / NFOLD \n    ypred[:, 1] += 2 * ystd / NFOLD\n    print(r'Fold {}: score = {}'.format(cnt, lb_metric(train.iloc[val_idx], oof[val_idx, :])))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Results","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"score = lb_metric(train, oof)\nprint(f'Overall CV = {score}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(oof[:, -1])\nplt.title(\"uncertainty in prediction\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prediction","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv(INPUT_DIR + '/sample_submission.csv')\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test['FVC_pred'] = ypred[:, 0]\ntest['Confidence'] = ypred[:, 1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test[[ID, 'FVC_pred', 'Confidence']].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test[['FVC_pred', 'Confidence']].describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = submission.drop(columns=['FVC', 'Confidence']).merge(test[['Patient_Week', 'FVC_pred', 'Confidence']], \n                                                           on='Patient_Week')\nsub.columns = submission.columns\nsub.to_csv('submission.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}