{"cells":[{"metadata":{},"cell_type":"markdown","source":"# LGBM Quantile\n## Nested Cross-Validation for Estimating Generalization Error\nGiven that the public leaderboard forms only a small proportion of the total leaderboard, we might expect a 'shakeup' where the best performing models on the private dataset arent necessarily the best performing on the overall test set.\n\nCross-validation, where we train and validate on different subsets of the data, can either allow us to look for the best hyperparameters (i.e. those that perform the best on average over all of the subsets). However it does not allow us to estimate the generalization to the test set. This is because we select the best performing hyperparameters for the validation set so taking the cross-validation average will lead to an overly optimistic estimate.\n\nTo get around this problem, nested cross-validation has folds within folds. Each 'inner' fold estimates the best hyperparameters (or for neural networks the convergence), while each 'outer' fold takes the optimized model from the inner fold and estimates a test error on the test data for that fold.\n\nHaving set up a robust cross-validation method, we can compare different methods in a rigourous way (and a way that hopefully will translate to the private test set).","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"http://ethen8181.github.io/machine-learning/ab_tests/quantile_regression/quantile_regression.html","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n#Load the dependanciesfrom tqdm.notebook import tqdm\nimport cv2\nimport copy\nfrom pathlib import Path\nfrom sklearn.model_selection import KFold\nfrom skimage.segmentation import clear_border\nfrom skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing\nfrom skimage.measure import label, regionprops\nfrom skimage.segmentation import clear_border\nfrom skimage.filters import roberts, sobel\nfrom scipy import ndimage as ndi\nfrom skimage import measure, morphology\nfrom scipy.stats import kurtosis\nimport seaborn as sns\nimport scipy\nfrom lightgbm import LGBMRegressor\n#tuning hyperparameters\nfrom bayes_opt import BayesianOptimization\nfrom skopt  import BayesSearchCV \n#building models\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom sklearn.model_selection import train_test_split, StratifiedKFold, cross_val_score\nfrom tqdm.notebook import tqdm\nfrom sklearn.ensemble import GradientBoostingRegressor\nimport copy\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport random\nimport matplotlib.pyplot as plt\nplt.style.use('ggplot')\nimport pydicom.pixel_data_handlers.gdcm_handler as gdcm_handler \n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\nimport pydicom\nimport os\nfrom torch.utils.data import DataLoader, Dataset\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_all(seed=42):\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = False\n    \nseed_all()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef calc_metric_loss(pred_fvc,sigma,true_fvc):\n    true_fvc=np.reshape(true_fvc,pred_fvc.shape)\n    sigma[sigma<70]=70\n    delta=np.abs(pred_fvc-true_fvc)\n    delta[delta>1000]=1000\n    metric=-(np.sqrt(2)*delta/sigma)-np.log(np.sqrt(2)*sigma)\n    return -metric\n\ndef calc_fvc_loss(pred_fvc,true_fvc):\n    true_fvc=np.reshape(true_fvc,pred_fvc.shape)\n    fvc_err=np.abs(pred_fvc-true_fvc)\n    return fvc_err","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_training_loss(train, val,title='loss'):\n    plt.figure()\n    plt.plot(train, label='Train')\n    plt.plot(val, label='Val')\n    if title=='loss':\n        plt.title('Model Training Loss')\n    else:\n        plt.title('Model Metric Loss')\n    plt.xlabel('Epoch')\n    plt.ylabel('Loss')\n    plt.legend()\n    plt.savefig('training_loss')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Helper functions for preprocessing image data into 'histogram' features","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load the scans in given folder path\ndef load_scan(path):\n\n    #slices = [pydicom.read_file(path / s) for s in os.listdir(path)]\n    slices = [pydicom.read_file(path / s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n    if slice_thickness==0:\n        slice_thickness=slices[0].SliceThickness\n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices\n\ndef get_pixels_hu(slices):\n    image = np.stack([np.array(s.pixel_array,dtype=np.int16) for s in slices])\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n\n    # Set outside-of-scan pixels to 0\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    for slice_number in range(len(slices)):\n        \n        intercept = slices[slice_number].RescaleIntercept\n        slope = slices[slice_number].RescaleSlope\n        \n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slice_number] = image[slice_number].astype(np.int16)\n            \n        image[slice_number] += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)\n\ndef resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    #spacing = np.array([scan[0].SliceThickness] + scan[0].PixelSpacing, dtype=np.float32)\n    spacing = np.array([scan[0].SliceThickness] + list(scan[0].PixelSpacing), dtype=np.float32)\n    resize_factor = spacing / new_spacing\n    new_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n    return image, new_spacing\n\ndef get_segmented_lungs(im, plot=False):\n    \n    '''\n    This funtion segments the lungs from the given 2D slice.\n    '''\n    if plot == True:\n        f, plots = plt.subplots(8, 1, figsize=(5, 40))\n    '''\n    Step 1: Convert into a binary image. \n    '''\n    binary = im < -200\n    if plot == True:\n        plots[0].axis('off')\n        plots[0].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 2: Remove the blobs connected to the border of the image.\n    '''\n    cleared = clear_border(binary)\n    if plot == True:\n        plots[1].axis('off')\n        plots[1].imshow(cleared, cmap=plt.cm.bone) \n    '''\n    Step 3: Label the image.\n    '''\n    label_image = label(cleared)\n    if plot == True:\n        plots[2].axis('off')\n        plots[2].imshow(label_image, cmap=plt.cm.bone) \n    '''\n    Step 4: Keep the labels with 2 largest areas.\n    '''\n    areas = [r.area for r in regionprops(label_image)]\n    areas.sort()\n    if len(areas) > 2:\n        for region in regionprops(label_image):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       label_image[coordinates[0], coordinates[1]] = 0\n    binary = label_image > 0\n    if plot == True:\n        plots[3].axis('off')\n        plots[3].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 5: Erosion operation with a disk of radius 2. This operation is \n    seperate the lung nodules attached to the blood vessels.\n    '''\n    selem = disk(2)\n    binary = binary_erosion(binary, selem)\n    if plot == True:\n        plots[4].axis('off')\n        plots[4].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 6: Closure operation with a disk of radius 10. This operation is \n    to keep nodules attached to the lung wall.\n    '''\n    selem = disk(10)\n    binary = binary_closing(binary, selem)\n    if plot == True:\n        plots[5].axis('off')\n        plots[5].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 7: Fill in the small holes inside the binary mask of lungs.\n    '''\n    edges = roberts(binary)\n    binary = ndi.binary_fill_holes(edges)\n    if plot == True:\n        plots[6].axis('off')\n        plots[6].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 8: Superimpose the binary mask on the input image.\n    '''\n    get_high_vals = binary == 0\n    im[get_high_vals] = 0\n    if plot == True:\n        plots[7].axis('off')\n        plots[7].imshow(im, cmap=plt.cm.bone) \n        \n    return im\n\ndef get_kurtosis_stats(ids,ctscans_dir):\n    kurt=[]\n    std=[]\n    mean=[]\n    median=[]\n    for i in ids:\n        print(i)\n        #try:\n        patient_path= ctscans_dir / i\n        scan = load_scan(patient_path)\n        image=get_pixels_hu(scan)\n        image, new_spacing = resample(image, scan, new_spacing=[2,2,2])\n        image=np.asarray([get_segmented_lungs(slice) for slice in image])\n        kurt_i=kurtosis(image.ravel()[image.ravel() < -200])\n        std_i=image.ravel()[image.ravel() < -200].std()\n        mean_i=image.ravel()[image.ravel() < -200].mean()\n        median_i=np.median(image.ravel()[image.ravel() < -200])\n        print('Kurtosis: ', kurt_i)\n        print('Standard Deviation: ', std_i)\n        kurt.append(kurt_i)\n        std.append(std_i)\n        mean.append(mean_i)\n        median.append(median_i)\n        ax=sns.kdeplot(image.ravel()[(image.ravel() < 0)&(image.ravel() > -1200)], bw=0.5)\n        ax.set(xlabel='HU', ylabel='% voxels',title='Histogram of voxel characteristics')\n        plt.show()\n        plt.imshow(image[round(image.shape[0]/2),:,:])\n        plt.show()\n        #except:\n            #print('error')\n            #kurt.append(np.nan)\n            #std.append(np.nan)\n            #mean.append(np.nan)\n            #median.append(np.nan)\n    return kurt,std,mean,median","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Helper functions for preprocessing data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_and_prepare_data(add_pixel_stats=True):\n    train=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\n    test=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n    submission=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')\n    \n    #Prepare Train Data\n    train['base_Weeks']=train.groupby(['Patient'])['Weeks'].transform('min')\n    base=train[train.Weeks==train.base_Weeks]\n    base = base.rename(columns={'FVC': 'base_FVC','Percent': 'base_Percent'})\n    base.drop_duplicates(subset=['Patient', 'Weeks'], keep='first',inplace=True)\n    train=train.merge(base[['Patient','base_FVC','base_Percent']],on='Patient',how='left')\n    train['Week_passed'] = train['Weeks'] - train['base_Weeks']\n    \n    test = test.rename(columns={'Weeks': 'base_Weeks', 'FVC': 'base_FVC','Percent': 'base_Percent'})\n    # Adding Sample Submission\n    submission = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/sample_submission.csv\")\n    # In submisison file, format: ID_'week', using lambda to split the ID\n    submission['Patient'] = submission['Patient_Week'].apply(lambda x:x.split('_')[0])\n    # In submisison file, format: ID_'week', using lambda to split the Week\n    submission['Weeks'] = submission['Patient_Week'].apply(lambda x:x.split('_')[1]).astype(int)\n    test = submission.drop(columns = [\"FVC\", \"Confidence\"]).merge(test, on = 'Patient')\n    test['Week_passed'] = test['Weeks'] - test['base_Weeks']\n    test=test[train.columns.drop(['FVC','Percent'])]\n    \n    if add_pixel_stats:\n        pixel_stats=pd.read_csv('../input/osic-histogram-features/train_pixel_stats.csv')\n        train=train.merge(pixel_stats[['Patient','kurtosis','std','mean','median']],how='left',on='Patient')\n        test_ids=test.Patient.unique()\n        root_dir = Path('/kaggle/input/osic-pulmonary-fibrosis-progression')\n        ct_scans_dir=root_dir/'test'\n        pixel_stats_test=test.copy()\n        pixel_stats_test.drop_duplicates(subset=['Patient'],inplace=True)\n        k,s,m,me=get_kurtosis_stats(test_ids,ct_scans_dir)\n        pixel_stats_test['kurtosis']=np.array(k)\n        pixel_stats_test['std']=np.array(s)\n        pixel_stats_test['mean']=np.array(m)\n        pixel_stats_test['median']=np.array(me)\n        test=test.merge(pixel_stats_test[['Patient','kurtosis','std','mean','median']],how='left',on='Patient')\n    return train, test\n\ndef OH_encode(train,test):\n    #OH Encoding of categorical variables (https://www.kaggle.com/ulrich07/osic-keras-starter-with-custom-metrics)\n    COLS = ['Sex','SmokingStatus']\n    for col in COLS:\n        for mod in train[col].unique():\n            train[mod] = (train[col] == mod).astype(int)\n            test[mod] = (test[col] == mod).astype(int)\n        train.drop(col,axis=1,inplace=True)\n        test.drop(col,axis=1,inplace=True)\n    return train, test\n\ndef Scale(train):\n    from sklearn import preprocessing\n    robust_scaler = preprocessing.RobustScaler()\n    train.loc[:,train.columns.difference(['Patient','FVC','Percent','Weeks','base_Weeks'])]=robust_scaler.fit_transform(train.loc[:,train.columns.difference(['Patient','FVC','Percent','Weeks','base_Weeks'])])\n    return robust_scaler","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Pytorch Dataset class\nNeatened up since some of my other notebooks!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class OSIC(Dataset):\n    def __init__(self,patient_ids,df,scaler=None,train=True,impute_vals=None):\n        root_dir = Path('/kaggle/input/osic-pulmonary-fibrosis-progression')\n        self.df=df.copy()\n        self.df=self.df.loc[self.df.Patient.isin(patient_ids),:]\n        if not train:\n            ct_scans_dir=root_dir/'test'\n        else:\n            ctscans_dir=root_dir/'train'\n        self.df.loc[:,self.df.columns.difference(['Patient','FVC','Percent','Weeks','base_Weeks'])]=scaler.transform(self.df.loc[:,self.df.columns.difference(['Patient','FVC','Percent','Weeks','base_Weeks'])])\n        self.data=self.df.loc[:,self.df.columns.difference(['FVC','Patient','Percent'])].values\n        if train:\n            self.impute_vals=np.nanmean(self.data, axis=0)\n        else:\n            self.impute_vals=impute_vals\n        inds = np.where(np.isnan(self.data))\n        self.data[inds] = np.take(self.impute_vals, inds[1])\n        self.patients=self.df['Patient'].values\n        self.train=train\n        if self.train:\n            self.fvc=self.df['FVC'].values\n    \n    def __len__(self):\n        return len(self.df)\n\n    def __getitem__(self, idx):\n        if self.train:\n            data = {'fvc': self.fvc[idx],\n                   'data': self.data[idx]}\n        else:\n            \n            data = {'data': self.data[idx]}\n        return data","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Function to perform bayesian hyperparameter search for light gbm","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"https://www.kaggle.com/somang1418/tuning-hyperparameters-under-10-minutes-lgbm","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def bayes_parameter_opt_lgb(X, y, alpha=0.5,init_round=15, opt_round=25, n_folds=3, random_seed=6,n_estimators=10000, output_process=False):\n    # prepare data\n    train_data = lgb.Dataset(data=X, label=y, free_raw_data=False)\n    # parameters\n    def lgb_eval(learning_rate,num_leaves, feature_fraction, bagging_fraction, max_depth, max_bin, min_data_in_leaf,min_sum_hessian_in_leaf,subsample,alpha=alpha):\n        params = {'objective':'quantile','alpha':alpha,'boosting_type': 'gbdt'}\n        params['learning_rate'] = max(min(learning_rate, 1), 0)\n        params[\"num_leaves\"] = int(round(num_leaves))\n        params['feature_fraction'] = max(min(feature_fraction, 1), 0)\n        params['bagging_fraction'] = max(min(bagging_fraction, 1), 0)\n        params['max_depth'] = int(round(max_depth))\n        params['max_bin'] = int(round(max_depth))\n        params['min_data_in_leaf'] = int(round(min_data_in_leaf))\n        params['min_sum_hessian_in_leaf'] = min_sum_hessian_in_leaf\n        params['subsample'] = max(min(subsample, 1), 0)\n        \n        cv_result = lgb.cv(params, train_data, nfold=n_folds, seed=random_seed, stratified=False, verbose_eval =200,metrics='mae')\n        return -min(cv_result['l1-mean'])\n     \n    lgbBO = BayesianOptimization(lgb_eval, {'learning_rate': (0.01, 1.0),\n                                            'num_leaves': (24, 80),\n                                            'feature_fraction': (0.1, 0.9),\n                                            'bagging_fraction': (0.8, 1),\n                                            'max_depth': (5, 30),\n                                            'max_bin':(20,90),\n                                            'min_data_in_leaf': (20, 80),\n                                            'min_sum_hessian_in_leaf':(0,100),\n                                           'subsample': (0.01, 1.0)}, random_state=200,verbose=0)\n\n    \n    #n_iter: How many steps of bayesian optimization you want to perform. The more steps the more likely to find a good maximum you are.\n    #init_points: How many steps of random exploration you want to perform. Random exploration can help by diversifying the exploration space.\n    \n    lgbBO.maximize(init_points=init_round, n_iter=opt_round)\n    \n    models=[]\n    for model in range(len( lgbBO.res)):\n        models.append(lgbBO.res[model]['target'])\n    \n    \n    # return best parameters\n    return lgbBO.res[pd.Series(models).idxmax()]['target'],lgbBO.res[pd.Series(models).idxmax()]['params']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 'Inner' Fold training loop\nPerforms bayesian hyperparameter optimization on the inner train and validation set","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def train_model(ids,train,quantiles):\n    \n    np.random.shuffle(ids)\n    train_ids,val_ids=np.split(ids, [int(round(0.9 * len(ids), 0))])\n    \n    scaler=Scale(train.loc[train.Patient.isin(train_ids),:])\n    train_dataset = OSIC(train_ids,train,scaler=scaler)  \n    val_dataset = OSIC(val_ids,train,scaler=scaler,impute_vals=train_dataset.impute_vals)  \n    \n    lgb_quantile_alphas = {}\n    for quantile_alpha in quantiles:\n        # to train a quantile regression, we change the objective parameter and\n        # specify the quantile value we're interested in\n        num_round = 15000\n        trn_data = lgb.Dataset(train_dataset.data, train_dataset.fvc)\n        val_data = lgb.Dataset(val_dataset.data, val_dataset.fvc)\n        opt_params = bayes_parameter_opt_lgb(train_dataset.data, train_dataset.fvc, init_round=5, opt_round=10, n_folds=3, random_seed=6,n_estimators=10000,alpha=quantile_alpha)\n        opt_params[1][\"num_leaves\"] = int(round(opt_params[1][\"num_leaves\"]))\n        opt_params[1]['max_depth'] = int(round(opt_params[1]['max_depth']))\n        opt_params[1]['min_data_in_leaf'] = int(round(opt_params[1]['min_data_in_leaf']))\n        opt_params[1]['max_bin'] = int(round(opt_params[1]['max_bin']))\n        opt_params[1]['objective']='quantile'\n        opt_params[1]['alpha']=quantile_alpha\n        opt_params=opt_params[1]\n        clf=lgb.train(opt_params, trn_data, num_round, valid_sets = [trn_data, val_data], verbose_eval=500, early_stopping_rounds = 250)\n        lgb_quantile_alphas[quantile_alpha] = clf\n    \n    return lgb_quantile_alphas, scaler, train_dataset.impute_vals","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model validation function\nTakes the model optimised in the inner loop and applies it to 'out-of-fold' data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def val_model(lgb_quantile_alphas,val_dataset,quantiles):\n    fvc_pred = lgb_quantile_alphas[quantiles[1]].predict(val_dataset.data)\n    sigma_pred = lgb_quantile_alphas[quantiles[2]].predict(val_dataset.data)-lgb_quantile_alphas[quantiles[0]].predict(val_dataset.data)\n    fvc_loss=calc_fvc_loss(fvc_pred,val_dataset.fvc).mean()\n    metric_loss=calc_metric_loss(fvc_pred,sigma_pred,val_dataset.fvc).mean()\n    return fvc_loss, metric_loss","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Nested CV\nPerforms the inner-outer loop k times. We can then average the results of each outer fold to get a generalization estimate","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def nested_CV(train,test,quantiles,k_folds=5):\n    ids=train.Patient.unique()\n    np.random.shuffle(ids)\n    kf = KFold(n_splits=k_folds)\n    test_losses=[]\n    test_metrics=[]\n    models=[]\n    scalers=[]\n    all_impute_vals=[]\n    fold_ids=[]\n    fold=1\n    for train_index, test_index in kf.split(ids):\n        print(\"Training fold \",fold)\n        fold+=1\n        # Train\n        model, scaler, impute_vals=train_model(ids[train_index],train,quantiles)\n        # Validate\n        val_dataset = OSIC(ids[test_index],train,scaler=scaler)  \n        loss, metric = val_model(model,val_dataset,quantiles)\n        print(\"Validation fvc loss: \",loss)\n        print(\"Validation metric: \",metric)\n        test_losses.append(loss)\n        test_metrics.append(metric)\n        models.append(model)\n        scalers.append(scaler)\n        all_impute_vals.append(impute_vals)\n        fold_ids.append(test_index)\n    return models,scalers,all_impute_vals, test_losses, test_metrics, fold_ids","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Run the nested CV to obtain generalization estimates","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#add_pixel_stats=False is faster since doesn't need to preprocess the test data. Processed train data available in the attached dataset.\ntrain,test=load_and_prepare_data(add_pixel_stats=False)\ntrain,test=OH_encode(train,test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"quantiles=(0.2,0.5,0.8)\nall_models,all_scalers,all_impute_vals, test_losses, test_metrics,kf_splits=nested_CV(train,test,quantiles,k_folds=5)\nprint(\"expected generalization metric: \", np.array(test_metrics).mean(), \" std: \", np.array(test_metrics).std())\nprint(\"expected generalization loss: \", np.array(test_losses).mean(), \" std: \", np.array(test_losses).std())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"quantiles=(0.1,0.5,0.9)\nall_models,all_scalers,all_impute_vals, test_losses, test_metrics,kf_splits=nested_CV(train,test,quantiles,k_folds=5)\nprint(\"expected generalization metric: \", np.array(test_metrics).mean(), \" std: \", np.array(test_metrics).std())\nprint(\"expected generalization loss: \", np.array(test_losses).mean(), \" std: \", np.array(test_losses).std())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Test Data","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Re-run the modelling procedure for the whole train data and make predictions on the test data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Re run the training process on all of the data\nall_train_ids=train.Patient.unique()\nlgb_quantile_alphas, scaler, impute_vals=train_model(all_train_ids,train,quantiles)\n\ntest_ids=test.Patient.unique()\ntest_dataset = OSIC(test_ids,test,scaler=scaler,impute_vals=impute_vals,train=False) \n\nfvc_pred = lgb_quantile_alphas[quantiles[1]].predict(test_dataset.data)\nsigma_pred = lgb_quantile_alphas[quantiles[2]].predict(test_dataset.data)-lgb_quantile_alphas[quantiles[0]].predict(test_dataset.data)\n\ntest['FVC']=fvc_pred\ntest['Confidence']=sigma_pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test['Patient_Week']=test[\"Patient\"] + '_' + test['Weeks'].apply(str)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission=submission[['Patient_Week']].merge(test[['Patient_Week','FVC','Confidence']],on='Patient_Week')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission.csv', index=False, float_format='%.1f')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Post-Match Analysis","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Test Predictions","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(submission['FVC'],submission['Confidence'])\nplt.title('Test')\nplt.xlabel('FVC')\nplt.ylabel('Confidence')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_dataset = OSIC(all_train_ids,train,scaler=scaler,impute_vals=impute_vals,train=True) \n\n\nfvc_pred_train = lgb_quantile_alphas[quantiles[1]].predict(train_dataset.data)\nsigma_pred_train = lgb_quantile_alphas[quantiles[2]].predict(train_dataset.data)-lgb_quantile_alphas[quantiles[0]].predict(train_dataset.data)\n\nprint('train metric', calc_metric_loss(fvc_pred_train,sigma_pred_train,train_dataset.fvc).mean())\n\nplt.scatter(fvc_pred_train,sigma_pred_train)\nplt.title('Train')\nplt.xlabel('FVC')\nplt.ylabel('Confidence')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(train_dataset.fvc,fvc_pred_train)\nplt.title('Train: predicted FVC vs true FVC')\nplt.xlabel('True FVC')\nplt.ylabel('Predicted FVC')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## All","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(submission['FVC'], alpha=0.5,label='test')\nplt.hist(fvc_pred_train, alpha=0.5,label='train')\nplt.legend()\nplt.title('Histogram of FVC predictions')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(submission['Confidence'], alpha=0.5,label='test')\nplt.hist(sigma_pred_train, alpha=0.5,label='train')\nplt.legend()\nplt.title('Histogram of Confidence predictions')","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}