{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Basic Data Augmentation & Feature Reduction\n---\nA common problem with small datasets is that models will identify patterns from random variance in the data rather than statistically reliable trends. This problem is exacerbated when we have a large number of features which can interact with eachother. Without oversampling, there are approximately 4194 distinct rows that can be generated with the LANL Earthquake Prediction dataset. How can we increase this number while avoiding oversampling and the consequent risk of leakage?\n\nThis kernel presents a rudimentary approach to data augmentation, so that new data can be generated which will diminish the effect of spurious feature interactions. We can then use this to identify the features in our dataset that are most vulnerable to this problem, along with the features that are more statistically robust.\n\nThe augmented data will use a proportion of features from the original data and retain their values. The remaining features will be values that have been sampled at random from each feature's overall distribution. "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport numpy.random\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_absolute_error\nfrom tqdm import tqdm_notebook as tqdm\nimport matplotlib.pyplot as plt\nimport seaborn as sns\ntrain = pd.read_csv('../input/eq-100-features/train_100.csv')\ntest = pd.read_csv('../input/eq-100-features/test_100.csv')\ny = pd.read_csv('../input/eq-100-features/y.csv').values.flatten()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"First let's quickly establish a baseline CV score using the function from [one of my other kernels](https://www.kaggle.com/bigironsphere/parameter-tuning-in-one-function-with-hyperopt), `quick_hyperopt()`."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#import required packages\nimport lightgbm as lgb\nimport xgboost as xgb\nimport catboost as cb\nimport gc\nfrom hyperopt import hp, tpe, Trials, STATUS_OK\nfrom hyperopt.fmin import fmin\nfrom hyperopt.pyll.stochastic import sample\n#optional but advised\nimport warnings\nwarnings.filterwarnings('ignore')\n\n#GLOBAL HYPEROPT PARAMETERS\nNUM_EVALS = 1000 #number of hyperopt evaluation rounds\nN_FOLDS = 5 #number of cross-validation folds on data in each evaluation round\n\n#LIGHTGBM PARAMETERS\nLGBM_MAX_LEAVES = 2**10 #maximum number of leaves per tree for LightGBM\nLGBM_MAX_DEPTH = 25 #maximum tree depth for LightGBM\nEVAL_METRIC_LGBM_REG = 'mae' #LightGBM regression metric. Note that 'rmse' is more commonly used \nEVAL_METRIC_LGBM_CLASS = 'auc'#LightGBM classification metric\n\n#XGBOOST PARAMETERS\nXGB_MAX_LEAVES = 2**12 #maximum number of leaves when using histogram splitting\nXGB_MAX_DEPTH = 25 #maximum tree depth for XGBoost\nEVAL_METRIC_XGB_REG = 'mae' #XGBoost regression metric\nEVAL_METRIC_XGB_CLASS = 'auc' #XGBoost classification metric\n\n#CATBOOST PARAMETERS\nCB_MAX_DEPTH = 8 #maximum tree depth in CatBoost\nOBJECTIVE_CB_REG = 'MAE' #CatBoost regression metric\nOBJECTIVE_CB_CLASS = 'Logloss' #CatBoost classification metric\n\n#OPTIONAL OUTPUT\nBEST_SCORE = 0\n\ndef quick_hyperopt(data, labels, package='lgbm', num_evals=NUM_EVALS, diagnostic=False, Class=False):\n    \n    #==========\n    #LightGBM\n    #==========\n    \n    if package=='lgbm':\n        \n        print('Running {} rounds of LightGBM parameter optimisation:'.format(num_evals))\n        #clear space\n        gc.collect()\n        \n        integer_params = ['max_depth',\n                         'num_leaves',\n                          'max_bin',\n                         'min_data_in_leaf',\n                         'min_data_in_bin']\n        \n        def objective(space_params):\n            \n            #cast integer params from float to int\n            for param in integer_params:\n                space_params[param] = int(space_params[param])\n            \n            #extract nested conditional parameters\n            if space_params['boosting']['boosting'] == 'goss':\n                top_rate = space_params['boosting'].get('top_rate')\n                other_rate = space_params['boosting'].get('other_rate')\n                #0 <= top_rate + other_rate <= 1\n                top_rate = max(top_rate, 0)\n                top_rate = min(top_rate, 0.5)\n                other_rate = max(other_rate, 0)\n                other_rate = min(other_rate, 0.5)\n                space_params['top_rate'] = top_rate\n                space_params['other_rate'] = other_rate\n            \n            subsample = space_params['boosting'].get('subsample', 1.0)\n            space_params['boosting'] = space_params['boosting']['boosting']\n            space_params['subsample'] = subsample\n            \n            if Class:\n                cv_results = lgb.cv(space_params, train, nfold = N_FOLDS, stratified=True,\n                                    early_stopping_rounds=100, metrics=EVAL_METRIC_LGBM_CLASS, seed=42)\n                best_loss = 1 - cv_results['auc-mean'][-1]\n                \n            else:\n                cv_results = lgb.cv(space_params, train, nfold = N_FOLDS, stratified=False,\n                                    early_stopping_rounds=100, metrics=EVAL_METRIC_LGBM_REG, seed=42)\n                best_loss = cv_results['l1-mean'][-1] #'l2-mean' for rmse\n            \n            return{'loss':best_loss, 'status': STATUS_OK }\n        \n        train = lgb.Dataset(data, labels)\n                \n        #integer and string parameters, used with hp.choice()\n        boosting_list = [{'boosting': 'gbdt',\n                          'subsample': hp.uniform('subsample', 0.5, 1)},\n                         {'boosting': 'goss',\n                          'subsample': 1.0,\n                         'top_rate': hp.uniform('top_rate', 0, 0.5),\n                         'other_rate': hp.uniform('other_rate', 0, 0.5)}] #if including 'dart', make sure to set 'n_estimators'\n        \n        if Class:\n            metric_list = ['auc'] #modify as required for other classification metrics\n            objective_list = ['binary', 'cross_entropy']\n        \n        else:\n            metric_list = ['MAE', 'RMSE'] \n            objective_list = ['huber', 'gamma', 'fair', 'tweedie']\n        \n        \n        space ={'boosting' : hp.choice('boosting', boosting_list),\n                'num_leaves' : hp.quniform('num_leaves', 2, LGBM_MAX_LEAVES, 1),\n                'max_depth': hp.quniform('max_depth', 2, LGBM_MAX_DEPTH, 1),\n                'max_bin': hp.quniform('max_bin', 32, 255, 1),\n                'min_data_in_leaf': hp.quniform('min_data_in_leaf', 1, 256, 1),\n                'min_data_in_bin': hp.quniform('min_data_in_bin', 1, 256, 1),\n                'min_gain_to_split' : hp.quniform('min_gain_to_split', 0.1, 5, 0.01),\n                'lambda_l1' : hp.uniform('lambda_l1', 0, 5),\n                'lambda_l2' : hp.uniform('lambda_l2', 0, 5),\n                'learning_rate' : hp.loguniform('learning_rate', np.log(0.005), np.log(0.2)),\n                'metric' : hp.choice('metric', metric_list),\n                'objective' : hp.choice('objective', objective_list),\n                'feature_fraction' : hp.quniform('feature_fraction', 0.5, 1, 0.01),\n                'bagging_fraction' : hp.quniform('bagging_fraction', 0.5, 1, 0.01)\n            }\n        \n        #optional: activate GPU for LightGBM\n        #follow compilation steps here:\n        #https://www.kaggle.com/vinhnguyen/gpu-acceleration-for-lightgbm/\n        #then uncomment lines below:\n        #space['device'] = 'gpu'\n        #space['gpu_platform_id'] = 0,\n        #space['gpu_device_id'] =  0\n\n        trials = Trials()\n        best = fmin(fn=objective,\n                    space=space,\n                    algo=tpe.suggest,\n                    max_evals=num_evals, \n                    trials=trials)\n                \n        #fmin() will return the index of values chosen from the lists/arrays in 'space'\n        #to obtain actual values, index values are used to subset the original lists/arrays\n        best['boosting'] = boosting_list[best['boosting']]['boosting']#nested dict, index twice\n        best['metric'] = metric_list[best['metric']]\n        best['objective'] = objective_list[best['objective']]\n                \n        #cast floats of integer params to int\n        for param in integer_params:\n            best[param] = int(best[param])\n        \n        print('{' + '\\n'.join('{}: {}'.format(k, v) for k, v in best.items()) + '}')\n        if diagnostic:\n            return(best, trials)\n        else:\n            return(best)\n    \n    #==========\n    #XGBoost\n    #==========\n    \n    if package=='xgb':\n        \n        print('Running {} rounds of XGBoost parameter optimisation:'.format(num_evals))\n        #clear space\n        gc.collect()\n        \n        integer_params = ['max_depth']\n        \n        def objective(space_params):\n            \n            for param in integer_params:\n                space_params[param] = int(space_params[param])\n                \n            #extract multiple nested tree_method conditional parameters\n            #libera te tutemet ex inferis\n            if space_params['tree_method']['tree_method'] == 'hist':\n                max_bin = space_params['tree_method'].get('max_bin')\n                space_params['max_bin'] = int(max_bin)\n                if space_params['tree_method']['grow_policy']['grow_policy']['grow_policy'] == 'depthwise':\n                    grow_policy = space_params['tree_method'].get('grow_policy').get('grow_policy').get('grow_policy')\n                    space_params['grow_policy'] = grow_policy\n                    space_params['tree_method'] = 'hist'\n                else:\n                    max_leaves = space_params['tree_method']['grow_policy']['grow_policy'].get('max_leaves')\n                    space_params['grow_policy'] = 'lossguide'\n                    space_params['max_leaves'] = int(max_leaves)\n                    space_params['tree_method'] = 'hist'\n            else:\n                space_params['tree_method'] = space_params['tree_method'].get('tree_method')\n                \n            #for classification replace EVAL_METRIC_XGB_REG with EVAL_METRIC_XGB_CLASS\n            cv_results = xgb.cv(space_params, train, nfold=N_FOLDS, metrics=[EVAL_METRIC_XGB_REG],\n                             early_stopping_rounds=100, stratified=False, seed=42)\n            \n            best_loss = cv_results['test-mae-mean'].iloc[-1] #or 'test-rmse-mean' if using RMSE\n            #for classification, comment out the line above and uncomment the line below:\n            #best_loss = 1 - cv_results['test-auc-mean'].iloc[-1]\n            #if necessary, replace 'test-auc-mean' with 'test-[your-preferred-metric]-mean'\n            return{'loss':best_loss, 'status': STATUS_OK }\n        \n        train = xgb.DMatrix(data, labels)\n        \n        #integer and string parameters, used with hp.choice()\n        boosting_list = ['gbtree', 'gblinear'] #if including 'dart', make sure to set 'n_estimators'\n        metric_list = ['MAE', 'RMSE'] \n        #for classification comment out the line above and uncomment the line below\n        #metric_list = ['auc']\n        #modify as required for other classification metrics classification\n        \n        tree_method = [{'tree_method' : 'exact'},\n               {'tree_method' : 'approx'},\n               {'tree_method' : 'hist',\n                'max_bin': hp.quniform('max_bin', 2**3, 2**7, 1),\n                'grow_policy' : {'grow_policy': {'grow_policy':'depthwise'},\n                                'grow_policy' : {'grow_policy':'lossguide',\n                                                  'max_leaves': hp.quniform('max_leaves', 32, XGB_MAX_LEAVES, 1)}}}]\n        \n        #if using GPU, replace 'exact' with 'gpu_exact' and 'hist' with\n        #'gpu_hist' in the nested dictionary above\n        \n        objective_list_reg = ['reg:linear', 'reg:gamma', 'reg:tweedie']\n        objective_list_class = ['reg:logistic', 'binary:logistic']\n        #for classification change line below to 'objective_list = objective_list_class'\n        objective_list = objective_list_reg\n        \n        space ={'boosting' : hp.choice('boosting', boosting_list),\n                'tree_method' : hp.choice('tree_method', tree_method),\n                'max_depth': hp.quniform('max_depth', 2, XGB_MAX_DEPTH, 1),\n                'reg_alpha' : hp.uniform('reg_alpha', 0, 5),\n                'reg_lambda' : hp.uniform('reg_lambda', 0, 5),\n                'min_child_weight' : hp.uniform('min_child_weight', 0, 5),\n                'gamma' : hp.uniform('gamma', 0, 5),\n                'learning_rate' : hp.loguniform('learning_rate', np.log(0.005), np.log(0.2)),\n                'eval_metric' : hp.choice('eval_metric', metric_list),\n                'objective' : hp.choice('objective', objective_list),\n                'colsample_bytree' : hp.quniform('colsample_bytree', 0.1, 1, 0.01),\n                'colsample_bynode' : hp.quniform('colsample_bynode', 0.1, 1, 0.01),\n                'colsample_bylevel' : hp.quniform('colsample_bylevel', 0.1, 1, 0.01),\n                'subsample' : hp.quniform('subsample', 0.5, 1, 0.05),\n                'nthread' : -1\n            }\n        \n        trials = Trials()\n        best = fmin(fn=objective,\n                    space=space,\n                    algo=tpe.suggest,\n                    max_evals=num_evals, \n                    trials=trials)\n        \n        best['tree_method'] = tree_method[best['tree_method']]['tree_method']\n        best['boosting'] = boosting_list[best['boosting']]\n        best['eval_metric'] = metric_list[best['eval_metric']]\n        best['objective'] = objective_list[best['objective']]\n        \n        #cast floats of integer params to int\n        for param in integer_params:\n            best[param] = int(best[param])\n        if 'max_leaves' in best:\n            best['max_leaves'] = int(best['max_leaves'])\n        if 'max_bin' in best:\n            best['max_bin'] = int(best['max_bin'])\n        \n        print('{' + '\\n'.join('{}: {}'.format(k, v) for k, v in best.items()) + '}')\n        \n        if diagnostic:\n            return(best, trials)\n        else:\n            return(best)\n    \n    #==========\n    #CatBoost\n    #==========\n    \n    if package=='cb':\n        \n        print('Running {} rounds of CatBoost parameter optimisation:'.format(num_evals))\n        \n        #clear memory \n        gc.collect()\n            \n        integer_params = ['depth',\n                          #'one_hot_max_size', #for categorical data\n                          'min_data_in_leaf',\n                          'max_bin']\n        \n        def objective(space_params):\n                        \n            #cast integer params from float to int\n            for param in integer_params:\n                space_params[param] = int(space_params[param])\n                \n            #extract nested conditional parameters\n            if space_params['bootstrap_type']['bootstrap_type'] == 'Bayesian':\n                bagging_temp = space_params['bootstrap_type'].get('bagging_temperature')\n                space_params['bagging_temperature'] = bagging_temp\n                \n            if space_params['grow_policy']['grow_policy'] == 'LossGuide':\n                max_leaves = space_params['grow_policy'].get('max_leaves')\n                space_params['max_leaves'] = int(max_leaves)\n                \n            space_params['bootstrap_type'] = space_params['bootstrap_type']['bootstrap_type']\n            space_params['grow_policy'] = space_params['grow_policy']['grow_policy']\n                           \n            #random_strength cannot be < 0\n            space_params['random_strength'] = max(space_params['random_strength'], 0)\n            #fold_len_multiplier cannot be < 1\n            space_params['fold_len_multiplier'] = max(space_params['fold_len_multiplier'], 1)\n                       \n            #for classification set stratified=True\n            cv_results = cb.cv(train, space_params, fold_count=N_FOLDS, \n                             early_stopping_rounds=25, stratified=False, partition_random_seed=42)\n           \n            best_loss = cv_results['test-MAE-mean'].iloc[-1] #'test-RMSE-mean' for RMSE\n            #for classification, comment out the line above and uncomment the line below:\n            #best_loss = cv_results['test-Logloss-mean'].iloc[-1]\n            #if necessary, replace 'test-Logloss-mean' with 'test-[your-preferred-metric]-mean'\n            \n            return{'loss':best_loss, 'status': STATUS_OK}\n        \n        train = cb.Pool(data, labels.astype('float32'))\n        \n        #integer and string parameters, used with hp.choice()\n        bootstrap_type = [{'bootstrap_type':'Poisson'}, \n                           {'bootstrap_type':'Bayesian',\n                            'bagging_temperature' : hp.loguniform('bagging_temperature', np.log(1), np.log(50))},\n                          {'bootstrap_type':'Bernoulli'}] \n        LEB = ['No', 'AnyImprovement', 'Armijo'] #remove 'Armijo' if not using GPU\n        #score_function = ['Correlation', 'L2', 'NewtonCorrelation', 'NewtonL2']\n        grow_policy = [{'grow_policy':'SymmetricTree'},\n                       {'grow_policy':'Depthwise'},\n                       {'grow_policy':'Lossguide',\n                        'max_leaves': hp.quniform('max_leaves', 2, 32, 1)}]\n        eval_metric_list_reg = ['MAE', 'RMSE', 'Poisson']\n        eval_metric_list_class = ['Logloss', 'AUC', 'F1']\n        #for classification change line below to 'eval_metric_list = eval_metric_list_class'\n        eval_metric_list = eval_metric_list_reg\n                \n        space ={'depth': hp.quniform('depth', 2, CB_MAX_DEPTH, 1),\n                'max_bin' : hp.quniform('max_bin', 1, 32, 1), #if using CPU just set this to 254\n                'l2_leaf_reg' : hp.uniform('l2_leaf_reg', 0, 5),\n                'min_data_in_leaf' : hp.quniform('min_data_in_leaf', 1, 50, 1),\n                'random_strength' : hp.loguniform('random_strength', np.log(0.005), np.log(5)),\n                #'one_hot_max_size' : hp.quniform('one_hot_max_size', 2, 16, 1), #uncomment if using categorical features\n                'bootstrap_type' : hp.choice('bootstrap_type', bootstrap_type),\n                'learning_rate' : hp.uniform('learning_rate', 0.05, 0.25),\n                'eval_metric' : hp.choice('eval_metric', eval_metric_list),\n                'objective' : OBJECTIVE_CB_REG,\n                #'score_function' : hp.choice('score_function', score_function), #crashes kernel - reason unknown\n                'leaf_estimation_backtracking' : hp.choice('leaf_estimation_backtracking', LEB),\n                'grow_policy': hp.choice('grow_policy', grow_policy),\n                #'colsample_bylevel' : hp.quniform('colsample_bylevel', 0.1, 1, 0.01),# CPU only\n                'fold_len_multiplier' : hp.loguniform('fold_len_multiplier', np.log(1.01), np.log(2.5)),\n                'od_type' : 'Iter',\n                'od_wait' : 25,\n                'task_type' : 'GPU',\n                'verbose' : 0\n            }\n        \n        #optional: run CatBoost without GPU\n        #uncomment line below\n        #space['task_type'] = 'CPU'\n            \n        trials = Trials()\n        best = fmin(fn=objective,\n                    space=space,\n                    algo=tpe.suggest,\n                    max_evals=num_evals, \n                    trials=trials)\n        \n        #unpack nested dicts first\n        best['bootstrap_type'] = bootstrap_type[best['bootstrap_type']]['bootstrap_type']\n        best['grow_policy'] = grow_policy[best['grow_policy']]['grow_policy']\n        best['eval_metric'] = eval_metric_list[best['eval_metric']]\n        \n        #best['score_function'] = score_function[best['score_function']] \n        #best['leaf_estimation_method'] = LEM[best['leaf_estimation_method']] #CPU only\n        best['leaf_estimation_backtracking'] = LEB[best['leaf_estimation_backtracking']]        \n        \n        #cast floats of integer params to int\n        for param in integer_params:\n            best[param] = int(best[param])\n        if 'max_leaves' in best:\n            best['max_leaves'] = int(best['max_leaves'])\n        \n        print('{' + '\\n'.join('{}: {}'.format(k, v) for k, v in best.items()) + '}')\n        \n        if diagnostic:\n            return(best, trials)\n        else:\n            return(best)\n    \n    else:\n        print('Package not recognised. Please use \"lgbm\" for LightGBM, \"xgb\" for XGBoost or \"cb\" for CatBoost.')              ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lgbm_params = quick_hyperopt(train, y, 'lgbm', 150)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With some parameter tuning these features can obtain a CV score of just under 2. Now we can use another quick function to isolate the feature split and gain importance scores from a shuffled KFold. "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def quick_kfold_imp(X, y, test=None, params=None, n_fold=5, random_state=42):\n    \n    MAE = 0\n    folds = KFold(n_splits=n_fold, shuffle=True, random_state=random_state)\n    test_preds = np.zeros(len(test))\n    #obtain both split and gain importance\n    imp_s = np.zeros(X.shape[1])\n    imp_g = np.zeros(X.shape[1])\n\n    if type(y) is not np.ndarray:\n        y = y.values.flatten()\n        \n    for train_idx, valid_idx in folds.split(y):\n         \n        X_train, X_valid = X.iloc[train_idx, :], X.iloc[valid_idx, :]\n        y_train, y_valid = y[train_idx], y[valid_idx]\n        \n        model = lgb.LGBMRegressor(**params, n_estimators = 50000, n_jobs = -1, eval_metric='mae', importance_type='split')\n        model.fit(X_train, y_train,eval_set=[(X_train, y_train), (X_valid, y_valid)],\n                  verbose=0, early_stopping_rounds=200)\n        val_preds = model.predict(X_valid)\n        imp_s += model.feature_importances_/n_fold\n        \n        model = lgb.LGBMRegressor(**params, n_estimators = 50000, n_jobs = -1, eval_metric='mae', importance_type='gain')\n        model.fit(X_train, y_train,eval_set=[(X_train, y_train), (X_valid, y_valid)],\n                  verbose=0, early_stopping_rounds=200)\n        val_preds = model.predict(X_valid)\n        MAE += mean_absolute_error(y_valid, val_preds)/n_fold\n        imp_g += model.feature_importances_/n_fold\n        test_preds += model.predict(test)/n_fold\n                \n    print('OOF MAE: {}'.format(MAE))\n    \n    return(imp_s, imp_g, test_preds)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We'll run it in a loop with different seeds to get a more accurate picture of the feature importance."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_features = train.columns\nimp_split = np.zeros(train.shape[1])\nimp_gain = np.zeros(train.shape[1])\ntest_preds = np.zeros(len(test))\n\nN=25\nfor i in range(0, N):\n    imp_s, imp_g, preds = quick_kfold_imp(train, y, test, lgbm_params, random_state=i)\n    imp_split += imp_s/N\n    imp_gain += imp_g/N\n    test_preds += preds/N\n\ninitial_imp = pd.DataFrame({'feature':train_features,\n                          'importance_split':imp_split,\n                          'importance_gain':imp_gain,\n                          'importance_score':np.log((imp_split*imp_gain))})\n\ninitial_imp.sort_values('importance_score', ascending=False, inplace=True)\ninitial_imp.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"initial_imp.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The importance score was calculated as the natural log of the gain score multiplied by the split score."},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_height = int(np.floor(len(initial_imp)/5))\nplt.figure(figsize=(12, plot_height));\nsns.barplot(x='importance_score', y='feature', data=initial_imp);\nplt.title('Original Feature Scores');\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we have our baseline feature scores and predictions for the original data."},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_1 = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv')\nsub_1['time_to_failure'] = test_preds\nsub_1.to_csv('sub_orginal_data.csv', index=False)\nsub_1['time_to_failure'].hist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Generate Augmented Data\n---\nNow that we have established a baseline, we can produce our augmented dataset. In this version I will be substituting 50% of the features in each row with randomly sampled values from that feature's actual distribution."},{"metadata":{"trusted":true},"cell_type":"code","source":"a = np.arange(0, train.shape[1])\n#initialise aug dataframe - remember to set dtype!\ntrain_aug = pd.DataFrame(index=train.index, columns=train.columns, dtype='float64')\n\nfor i in tqdm(range(0, len(train))):\n    #ratio of features to be randomly sampled\n    AUG_FEATURE_RATIO = 0.5\n    #to integer count\n    AUG_FEATURE_COUNT = np.floor(train.shape[1]*AUG_FEATURE_RATIO).astype('int16')\n    \n    #randomly sample half of columns that will contain random values\n    aug_feature_index = np.random.choice(train.shape[1], AUG_FEATURE_COUNT, replace=False)\n    aug_feature_index.sort()\n    \n    #obtain indices for features not in aug_feature_index\n    feature_index = np.where(np.logical_not(np.in1d(a, aug_feature_index)))[0]\n        \n    #first insert real values for features in feature_index\n    train_aug.iloc[i, feature_index] = train.iloc[i, feature_index]\n              \n    #random row index to randomly sampled values for each features\n    rand_row_index = np.random.choice(len(train), len(aug_feature_index), replace=True)\n        \n    #for each feature being randomly sampled, extract value from random row in train\n    for n, j in enumerate(aug_feature_index):\n        train_aug.iloc[i, j] = train.iloc[rand_row_index[n], j]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Please note that `pandas` will set the datatype of its columns as `'object'` unless you specify otherwise.  I mention this because the above code, which takes less than 1 minute to process 4194 rows of 100 features, will take around an hour if `dtype` isn't set to `'float64'`!\n\nComparing the first few rows of the regular data and the augmented data:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_aug.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head(3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see, for a basic sanity check, that for each row in `train_aug` that half the values are the same as in `train`, and the remaining half are seemingly random. The random values for each feature were sampled from that feature's original distribution, and the overall distributions for each variable in `train` shouldn't be very different as a result. We can do another quick sanity check for this:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train['var_0'].hist(bins=25)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_aug['var_0'].hist(bins=25)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The distributions look almost identical and now we can examine the MAE for the augmented data, along with the feature importances.\n\n# Augmented Data Evaluation\n---\n\nThe corresponding y-values for the augmented rows will be the same for the original data. We can run `quick_hyperopt()` again to get an idea of the optimal CV score."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_all = pd.concat([train, train_aug])\ny_all = np.append(y, y)\n\nprint('Original train data shape: {}'.format(train.shape))\nprint('Augmented train data shape: {}'.format(train_all.shape))\n\nparams_all = quick_hyperopt(train_all, y_all, 'lgbm', 150)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The combined data has a higher CV, indicating that LightGBM has been less eager to identify predictive feature interactions. Now we can examine the feature importances, and how they have changed relative to the original data. We'll run the augmentation process within a loop to ensure an even distribution of random features."},{"metadata":{"trusted":true},"cell_type":"code","source":"imp_split_all = np.zeros(train_aug.shape[1])\nimp_gain_all = np.zeros(train_aug.shape[1])\ntest_preds_aug = np.zeros(len(test)) \n\nN = 25\nfor i in tqdm(range(0, N)):\n    \n    a = np.arange(0, train.shape[1])\n    np.random.seed(i)\n    train_aug = pd.DataFrame(index=train.index, columns=train.columns, dtype='float64')\n\n    for i in range(0, len(train)):\n        #ratio of features to be randomly sampled\n        AUG_FEATURE_RATIO = 0.5\n        #to integer count\n        AUG_FEATURE_COUNT = np.floor(train.shape[1]*AUG_FEATURE_RATIO).astype('int16')\n    \n        #randomly sample half of columns that will contain random values\n        aug_feature_index = np.random.choice(train.shape[1], AUG_FEATURE_COUNT, replace=False)\n        aug_feature_index.sort()\n    \n        #obtain indices for features not in aug_feature_index\n        feature_index = np.where(np.logical_not(np.in1d(a, aug_feature_index)))[0]\n        \n        #first insert real values for features in feature_index\n        train_aug.iloc[i, feature_index] = train.iloc[i, feature_index]\n              \n        #random row index to randomly sampled values for each features\n        rand_row_index = np.random.choice(len(train), len(aug_feature_index), replace=True)\n        \n        #for each feature being randomly sampled, extract value from random row in train\n        for n, j in enumerate(aug_feature_index):\n            train_aug.iloc[i, j] = train.iloc[rand_row_index[n], j]\n    \n    \n    train_all = pd.concat([train, train_aug])\n        \n    imp_s, imp_g, preds = quick_kfold_imp(train_all, y_all, test, params=params_all, random_state=i)\n    imp_split_all += imp_s/N\n    imp_gain_all += imp_g/N\n    test_preds_aug += preds/N\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aug_imp = pd.DataFrame({'feature':train_features,\n                          'importance_split':imp_split_all,\n                          'importance_gain':imp_gain_all,\n                          'importance_score':np.log((imp_split_all*imp_gain_all))})\n\naug_imp.sort_values('importance_score', ascending=False, inplace=True)\naug_imp.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From inspection, the most important features appear to have changed. "},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_height = int(np.floor(len(aug_imp)/5))\nplt.figure(figsize=(12, plot_height));\nsns.barplot(x='importance_score', y='feature', data=aug_imp);\nplt.title('Augmented Data Feature Scores');\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can then prepare a submission made with the augmented data and its different feature importances for comparison."},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_aug = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv')\nsub_aug['time_to_failure'] = test_preds_aug\nsub_aug.to_csv('sub_aug_data.csv', index=False)\nsub_aug['time_to_failure'].hist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Feature Importance Change\n---\n\nSo what is the relative change in feature importance for each feature? LightGBM can be inconsistent with feature importance scores so they'll have to be scaled."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import MinMaxScaler\nscaler = MinMaxScaler()\ninitial_imp.iloc[:, 1:] = scaler.fit_transform(initial_imp.iloc[:, 1:])\n\nscaler = MinMaxScaler()\naug_imp.iloc[:, 1:] = scaler.fit_transform(aug_imp.iloc[:, 1:])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aug_cols = ['feature'] + [x + '_all' for x in aug_imp.columns if x != 'feature']\naug_imp.columns = aug_cols\nfeature_df = initial_imp.merge(aug_imp, on='feature', how='inner')\nfeature_df['score_change'] = feature_df['importance_score_all'] - feature_df['importance_score']\nfeature_df.sort_values('score_change', ascending=False, inplace=True)\nfeature_df.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_height = int(np.floor(len(feature_df)/5))\nplt.figure(figsize=(12, plot_height));\nsns.barplot(x='score_change', y='feature', data=feature_df);\nplt.title('Difference in Feature Scores');\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This has clearly identified some features that became *less* valuable when more randomness is introduced into the dataset.  This indicates that these particular features may have been attributed weight by the model due to random variance instead of a genuine relationship with the target. As a last experiment, we can isolate the features whose change in feature importance was in the bottom two quintiles and remove them. Their predictions can then be evaluated on the test set."},{"metadata":{"trusted":true},"cell_type":"code","source":"MIN_SCORE = np.percentile(feature_df.score_change.values.flatten(), 20)\nfeatures = feature_df.loc[feature_df.score_change >= MIN_SCORE, :].feature.values.flatten()\ntrain = train[features]\ntest = test[features]\n\nparams_final = quick_hyperopt(train, y, 'lgbm', 150)\n\ntest_preds_final = np.zeros(len(test)) \n\nN = 30\nfor i in tqdm(range(0, N)):\n    \n    a = np.arange(0, train.shape[1])\n    np.random.seed(i)\n    train_aug = pd.DataFrame(index=train.index, columns=train.columns, dtype='float64')\n\n    for i in range(0, len(train)):\n        #ratio of features to be randomly sampled\n        AUG_FEATURE_RATIO = 0.5\n        #to integer count\n        AUG_FEATURE_COUNT = np.floor(train.shape[1]*AUG_FEATURE_RATIO).astype('int16')\n    \n        #randomly sample half of columns that will contain random values\n        aug_feature_index = np.random.choice(train.shape[1], AUG_FEATURE_COUNT, replace=False)\n        aug_feature_index.sort()\n    \n        #obtain indices for features not in aug_feature_index\n        feature_index = np.where(np.logical_not(np.in1d(a, aug_feature_index)))[0]\n        \n        #first insert real values for features in feature_index\n        train_aug.iloc[i, feature_index] = train.iloc[i, feature_index]\n              \n        #random row index to randomly sampled values for each features\n        rand_row_index = np.random.choice(len(train), len(aug_feature_index), replace=True)\n        \n        #for each feature being randomly sampled, extract value from random row in train\n        for n, j in enumerate(aug_feature_index):\n            train_aug.iloc[i, j] = train.iloc[rand_row_index[n], j]\n    \n    \n    train_all = pd.concat([train, train_aug])\n        \n    _,  _, preds = quick_kfold_imp(train_all, y_all, test, params=params_final, random_state=i)\n    test_preds_final += preds/N\n\nsub_final = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv')\nsub_final['time_to_failure'] = test_preds_final\nsub_final.to_csv('sub_quintile_features_removed.csv', index=False)\nsub_final['time_to_failure'].hist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data Score Evaluation\n---\n\nThe original features were a random selectin of features I had generated. The augmented data used all these features, but with the augmentation process detailed above to double the number of rows. The 'refined' data was the original data, minus bottom quartile of the features isolated in the Feature Importance Change section above. The three submissions scored:\n\n* **original features** Hyperopt CV: 1.987 LB: 1.454\n* **augmented data** Hyperopt CV: 2.012 LB 1.445\n* **refined data** Hyperopt CV: 2.005 **LB 1.435**\n\n\nData augmentation, and feature selection via data augmentation can clearly yield positive results. This is a basic introduction to data augmentation and doubtless you can improve on the methods in this kernel. Hopefully you will find a way to make your models more accurate on the test set."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}