{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom sklearn import linear_model\nfrom sklearn.model_selection import KFold\n\ntrain1 = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/train.csv\")\ntest = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/test.csv\")\n\n################### Team CV scheme ########################\n\ndef mask_first(x):\n    result = np.ones_like(x)\n    result[0] = 0\n    return result\n\nclass RobustCV(KFold):\n    \n    def __init__(self, df, *args, **kwargs):\n        \n        self.df = df.sort_values(['Patient', 'Weeks'])\n        \n        self.base_features = self.df.groupby('Patient').head(1)\n        self.base_features.columns = [c if c == 'Patient' else f'base_{c}' for c in self.base_features.columns]\n\n        self.df = self.df[['Patient', 'Weeks', 'FVC', ]].set_index(['Patient', 'Weeks'])\n        \n        super().__init__(*args, **kwargs)\n    \n    def split(self, patients, val_last_three_only=False, train_drop_first=False, train_last_three_only=False):\n        for train_idx, val_idx in super().split(patients):\n            train_patient = patients[train_idx]\n            val_patient = patients[val_idx]\n            assert set(train_patient).intersection(set(val_patient)) == set()\n            \n            train = self.df.loc[train_patient].reset_index()\n            \n            # drop first row (baseline week)\n            if train_drop_first:\n                mask = train.groupby('Patient')['Patient'].transform(mask_first).astype(bool)\n                train = train[mask]\n            \n            if train_last_three_only:\n                train = train.groupby('Patient').tail(3)\n\n            \n            train = pd.merge(train, self.base_features, left_on='Patient', right_on='Patient')\n            \n            \n            \n            val = self.df.loc[val_patient].reset_index()\n            \n            # drop first row (baseline week)\n            mask = val.groupby('Patient')['Patient'].transform(mask_first).astype(bool)\n            val = val[mask]\n            \n            # merge with base features\n            val = pd.merge(val, self.base_features, left_on='Patient', right_on='Patient')\n            \n            # if last three take last 3 week only per patient\n            if val_last_three_only:\n                val = val.groupby('Patient').tail(3)\n            \n            yield train, val\nROOT = '../input/osic-pulmonary-fibrosis-progression'\ntrain_df = pd.read_csv(f'{ROOT}/train.csv')\nN_FOLDS = 5\nkf = RobustCV(train_df, N_FOLDS, shuffle=True, random_state=2020)\nfor n_fold, (train, val) in enumerate(kf.split(train_df['Patient'].unique(), val_last_three_only=False, train_drop_first=False, train_last_three_only=False)):\n    if n_fold==0:\n        val['fold'] = n_fold\n        folddf = val.copy()        \n    else:\n        val['fold'] = n_fold\n        folddf = folddf.append(val.copy()).reset_index(drop=True)\nfolddf = folddf[['Patient', 'fold']].drop_duplicates(ignore_index=True)\n\ntrain1 = pd.merge(train1, folddf, on='Patient', how='inner')\n\n########################### Flag last three records per patient for CV results ##################################\ndef mlthree(s):\n    s2 = pd.Series(0, index=s.index)\n    s2.iloc[-1] = 1\n    if len(s2)>1:\n        s2.iloc[-2] = 1\n    if len(s2)>2:\n        s2.iloc[-3] = 1\n    return s2\n\ntrain1[\"forscore\"] = train1.groupby(['Patient'])['Weeks'].apply(mlthree)\ntest['forscore'] = 0\n######################################################################\n\n\n##################### Natural \"agumentation\" use every week (other than last 3) as baseline ###########################\n\ntrain1 = pd.merge(train1,\n                 ((train1.rename(columns={\"Weeks\": \"realbasewk\"}))[['Patient', 'realbasewk']]).groupby('Patient').min(),\n                 how='left', on='Patient')\n\ntrain2 = pd.merge(train1[train1['realbasewk']!=train1['Weeks']],\n         (train1[train1['forscore']==0].rename(columns={\"Weeks\": \"basewk\", \n                                            \"FVC\": \"base\",\n                                           \"Percent\": \"basepercent\"}))[['Patient', 'basewk', 'base', 'basepercent']],\n         how='outer', on='Patient').reset_index()\n\ntrain2['original'] = (train2['realbasewk']==train2['basewk']) # Flag for original records, only use those for CV scoring\n\n######################################################################\n\n\n\n################### Old CV scheme ########################\n#np.random.seed(42)\n#\n#folds = pd.DataFrame({'Patient': np.unique(train2['Patient'])}) #train2.loc[train2['train']==1, 'patno']\n#folds['fold'] = np.int32(np.random.choice(a=np.repeat([0,1,2,3,4], \n#                                            repeats=np.ceil(len(folds)/5)), \n#                                size=len(folds), \n#                                replace=False))\n##train2 = pd.merge(train2, folds, on='Patient', how='left') #.reset_index()\n#train2.loc[np.isnan(train2.fold.values), 'fold'] = -1\n\n\n################### Add some more features ########################\ntrain2 = train2.assign(logchg = lambda x: np.log(x.FVC) - np.log(x.base),\n                       logbase = lambda x: np.log(x.base),\n                       logbasepercent = lambda x: np.log(x.basepercent),\n                       logpercent = lambda x: np.log(x.Percent),\n                       logpercentchg = lambda x: np.log(x.Percent) - np.log(x.basepercent),\n                       chg = lambda x: x.FVC - x.base,\n                       relweek = lambda x: x.Weeks - x.basewk,\n                       logFVC = lambda x: np.log(x.FVC),\n                       sexm = lambda x: np.select([x.Sex.eq('Male'), x.Sex.eq('Female')], [1, 0]),\n                       smoker = lambda x: np.select( [x.SmokingStatus.eq(\"Never smoked\"),\n                                        x.SmokingStatus.eq(\"Ex-smoker\"),\n                                        x.SmokingStatus.eq(\"Currently smokes\")],\n                                      [0, 1, 2]))\n\ntrain2['rwk_bwk'] = train2['relweek']*train2['basewk']\ntrain2['rwk_smk1'] = train2['relweek']*(train2['smoker']>0)\ntrain2['rwk_smk2'] = train2['relweek']*(train2['smoker']>1)\ntrain2['rwk_lbp'] = train2['relweek']*train2['logbasepercent']\n\ntrain2['rwk_base'] = train2['relweek']*train2['base']\ntrain2['rwk_logbase'] = train2['relweek']*train2['logbase']\ntrain2['rwk_bp'] = train2['relweek']*train2['basepercent']\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# A bunch of random number seeds to ensure repeated model training produces varying, but reproducicble results\nmyseeds = [538, 42, 23, 666, 2020, 4629, 25479,1337, 13, 649913,\n           73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173,\n    179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,\n    283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,\n    419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541,\n    547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659,\n    661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809,\n    811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,\n    947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,\n    1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,\n    1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,\n    1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,\n    1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,\n    1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,\n    1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987,\n    1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,\n    2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287,\n    2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,\n    2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617,\n    2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,\n    2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903,\n    2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,\n    3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257,\n    3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,\n    3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571,\n    3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,\n    3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907,\n    3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,\n    4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231,\n    4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,\n    4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583,\n    4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get CV splits and data ready\nfold_splits = [ [ train2[ (train2['fold']!=fff )  ].index.tolist(), \n                 train2[ (train2['fold']==fff) & (train2['original']==True)].index.tolist()] for fff in range(5) ]\n# Note: we will model logFVC, because that might be numerically more stable (at least performed a little better) \ny = np.array( train2['logFVC'] ).flatten()\nX = np.array( train2[ ['relweek', 'logbase', 'rwk_smk2', 'rwk_logbase']] )\n\n\n# Number of RANSAC models to fit and take the median prediction over\nnum_ransac = 250\n\n# Dataframe that will contain out-of-fold predictions\noofs0 = train2.copy()\n\n# Loop over folds\nfor foldid in range(5):\n    trainidx, predidx = fold_splits[foldid][0], fold_splits[foldid][1]\n    \n    mepreds = np.empty([num_ransac, len(predidx)]) # For writing predictions into\n    \n    # Fit RANSAC models for the fold \n    for ransac_model in range(num_ransac):        \n        ransac = linear_model.RANSACRegressor(random_state=myseeds[ransac_model])\n        ransac.fit(X[trainidx], y[trainidx])\n        mepreds[ransac_model, :] = ransac.predict(X[predidx]) # write predictions from individual models\n        \n    oofs0.loc[predidx, 'pred'] = np.exp( np.median(mepreds, axis=0) ) # Point estimate is median prediction of the 250 models\n    oofs0.loc[predidx, 'predsd'] = np.std( np.exp( mepreds), axis=0)  # Point estimate is median prediction of the 250 models\n\n#print(train2[['Percent', 'ransacloo1']])\n#print(train2[['FVC', 'ransacloo2']])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Looking at how the MAE behaves gives some indication on whether a constant 'Confidence'-value would be okay\ninvestigate = oofs0.loc[oofs0['pred'].notnull(), ['FVC', 'pred', 'predsd', 'forscore', 'relweek']]\ninvestigate['logMAE'] = np.log(np.absolute(investigate['FVC']-investigate['pred']))\ninvestigate['logsd'] = np.log(investigate['predsd']) \n\nimport seaborn as sns\n\n# Firstly, no hint that the variation in the model predictions is a good indicator \n#sns.regplot(data=investigate[investigate['forscore']==1], x=\"logsd\", y=\"logMAE\",  lowess=True)\n\n# It seem like error is higher the further away we are from the baseline week (makes sense)\nsns.regplot(data=investigate[investigate['forscore']==1], x=\"relweek\", y=\"logMAE\",  lowess=True)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print( np.mean(np.absolute(oofs0['FVC']-oofs0['pred'])) ) #201.48929130152206\n#print( np.std(np.absolute(oofs0['FVC']-oofs0['pred'])) ) #201.48929130152206\n\nprint( -np.mean(np.absolute(oofs0['FVC']-oofs0['pred']))*np.sqrt(2)/290 - np.log(np.sqrt(2)*290)) \n\nconfs = 235 + oofs0['relweek']/63 * 80 \n#-6.987823939981736, 235,80: -6.985560881063153, 235,85: -6.9856983175228615, 235,75:-6.985543406781326  240: -6.985843635925146\nprint( np.mean( -np.absolute(oofs0['FVC']-oofs0['pred'])*np.sqrt(2) / confs  - np.log(np.sqrt(2)* confs )) )\n\n\nprint( np.mean(np.absolute(oofs0[oofs0['forscore']==1]['FVC']-oofs0[oofs0['forscore']==1]['pred'])) ) #201.48929130152206\n#print( np.std(np.absolute(oofs0['FVC']-oofs0['pred'])) ) #201.48929130152206\n\nprint( -np.mean(np.absolute(oofs0[oofs0['forscore']==1]['FVC']-oofs0[oofs0['forscore']==1]['pred']))*np.sqrt(2)/290 - np.log(np.sqrt(2)*290)) \n\nconfs = 235 + oofs0[oofs0['forscore']==1]['relweek']/63 * 80 \n#-6.987823939981736, 235,80: -6.985560881063153, 235,85: -6.9856983175228615, 235,75:-6.985543406781326  240: -6.985843635925146\nprint( np.mean( -np.absolute(oofs0[oofs0['forscore']==1]['FVC']-oofs0[oofs0['forscore']==1]['pred'])*np.sqrt(2) / confs  - np.log(np.sqrt(2)* confs )) )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"oofs1 = oofs0.loc[oofs0['pred'].notnull(), ['Patient', 'Weeks', 'relweek', 'pred', 'fold', 'forscore', 'FVC']].rename(columns={'FVC':'trueFVC', \"pred\": \"FVC\", 'forscore': 'last3'}).reset_index()\noofs1['Confidence'] = 235 + oofs1['relweek']/63 * 80 \n\noofs1['Patient_Week'] = oofs1['Patient'].astype(str) + '_' + oofs1['Weeks'].astype(str)\n\n#oofs1[ ['Patient', 'Weeks', 'FVC', 'Confidence'] ]\noofs1 = oofs1[ ['Patient_Week', 'FVC', 'Confidence', 'fold', 'last3', 'trueFVC'] ]\n\noofs1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"oofs1.to_csv('ransacaug_oof.csv')","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}