{"cells":[{"metadata":{"_uuid":"6493ba2b058a2768f4c3e64bceea10d5b98bdcf5"},"cell_type":"markdown","source":"

Feature Selection Methods

\n\nIn this kernel I will be trying a few different methods for feature selection in a regression problem. The dataset used has only 4194 samples, so it's possible to test a diverse set of parameters and number of features.\n\n

About the dataset

\n\nThe data for this notebook is from the [LANL Earthquake Prediction](https://www.kaggle.com/c/LANL-Earthquake-Prediction/data). The goal of this competition is to use seismic signals to predict the timing of laboratory earthquakes. We have only one original feature (acoustic data) from a single experiment and must make predictions for chunks of 150,000 data points. The most common aproach has been grouping the training data in chunks with the same size as the test set and extracting features. Details about the experiments can be found in [this discussion](https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/77526)."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import time\nimport numpy as np\nimport pandas as pd\nfrom boruta import BorutaPy\n# Plotly for Visualizations\nimport plotly.graph_objs as go\nimport plotly.figure_factory as ff\nfrom plotly.offline import init_notebook_mode, iplot\n# Sklearn\nfrom sklearn.feature_selection import SelectKBest, f_regression, mutual_info_regression, SelectFromModel\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold, train_test_split\nfrom sklearn.linear_model import LinearRegression, LassoCV\nfrom sklearn.ensemble import RandomForestRegressor\n# scipy (feature engineering)\nfrom scipy.signal import hilbert\nfrom scipy.signal import hann\nfrom scipy.signal import convolve\nfrom scipy import stats\nimport lightgbm as lgb\nimport warnings\n# Configurations\nwarnings.simplefilter(action='ignore', category=FutureWarning)\nwarnings.simplefilter(action='ignore', category=UserWarning)\nwarnings.simplefilter(action='ignore', category=RuntimeWarning)\ninit_notebook_mode(connected=True)\nRANDOM_SEED = 19\nnp.random.seed(RANDOM_SEED)\n\ndef lineplot(x_axis, y_axis, title=''):\n trace = go.Scatter(x=x_axis, y=y_axis, mode='lines+markers')\n\n layout = go.Layout(\n title=title, \n showlegend=False,\n xaxis=dict(\n title='Number of features',\n titlefont=dict(size=14, color='#7f7f7f')\n ),\n yaxis=dict(\n title='Validation score',\n titlefont=dict(size=14, color='#7f7f7f')\n )\n )\n fig = go.Figure(data=[trace], layout=layout)\n iplot(fig)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b21c4eca023b43c761d3ab6073cb5723dc836b0e"},"cell_type":"markdown","source":"Columns:\n* acoustic_data - the seismic signal\n* time_to_failure - the time (in seconds) until the next laboratory earthquake"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"data_type = {'acoustic_data': np.int16, 'time_to_failure': np.float64}\ntrain = pd.read_csv('../input/train.csv', dtype=data_type)\ntrain.head(3)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d0d60070d80cb62c950a522ae390db79672cc108"},"cell_type":"markdown","source":"

Feature Engineering

\n\nAs described in the introduction, the original dataset has approximately 600 million points, which will be grouped in 4194 chunks of 150,000 observations. Statistical features are extracted from each chunk according to [Lukyanenko's kernel](https://www.kaggle.com/artgor/earthquakes-fe-more-features-and-samples)."},{"metadata":{"trusted":true,"_uuid":"14b2f6c42ed1e259e863012a8878961c8dfcb341","_kg_hide-input":true},"cell_type":"code","source":"rows = 150_000\nsegments = int(np.floor(train.shape[0] / rows))\n\ndef add_trend_feature(arr, abs_values=False):\n idx = np.array(range(len(arr)))\n if abs_values:\n arr = np.abs(arr)\n lr = LinearRegression()\n lr.fit(idx.reshape(-1, 1), arr)\n return lr.coef_[0]\n\ndef classic_sta_lta(x, length_sta, length_lta):\n \n sta = np.cumsum(x ** 2)\n\n # Convert to float\n sta = np.require(sta, dtype=np.float)\n\n # Copy for LTA\n lta = sta.copy()\n\n # Compute the STA and the LTA\n sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n sta /= length_sta\n lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n lta /= length_lta\n\n # Pad zeros\n sta[:length_lta - 1] = 0\n\n # Avoid division by zero by setting zero values to tiny float\n dtiny = np.finfo(0.0).tiny\n idx = lta < dtiny\n lta[idx] = dtiny\n\n return sta / lta\n\nX_tr = pd.DataFrame(index=range(segments), dtype=np.float64)\n\ny_tr = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])\n\ntotal_mean = train['acoustic_data'].mean()\ntotal_std = train['acoustic_data'].std()\ntotal_max = train['acoustic_data'].max()\ntotal_min = train['acoustic_data'].min()\ntotal_sum = train['acoustic_data'].sum()\ntotal_abs_sum = np.abs(train['acoustic_data']).sum()\n\ndef calc_change_rate(x):\n change = (np.diff(x) / x[:-1]).values\n change = change[np.nonzero(change)[0]]\n change = change[~np.isnan(change)]\n change = change[change != -np.inf]\n change = change[change != np.inf]\n return np.mean(change)\n\nfor segment in range(segments):\n seg = train.iloc[segment*rows:segment*rows+rows]\n x = pd.Series(seg['acoustic_data'].values)\n y = seg['time_to_failure'].values[-1]\n \n y_tr.loc[segment, 'time_to_failure'] = y\n X_tr.loc[segment, 'mean'] = x.mean()\n X_tr.loc[segment, 'std'] = x.std()\n X_tr.loc[segment, 'max'] = x.max()\n X_tr.loc[segment, 'min'] = x.min()\n \n X_tr.loc[segment, 'mean_change_abs'] = np.mean(np.diff(x))\n X_tr.loc[segment, 'mean_change_rate'] = calc_change_rate(x)\n X_tr.loc[segment, 'abs_max'] = np.abs(x).max() \n X_tr.loc[segment, 'std_first_50000'] = x[:50000].std()\n X_tr.loc[segment, 'std_last_50000'] = x[-50000:].std()\n X_tr.loc[segment, 'std_first_10000'] = x[:10000].std()\n X_tr.loc[segment, 'std_last_10000'] = x[-10000:].std()\n \n X_tr.loc[segment, 'avg_first_50000'] = x[:50000].mean()\n X_tr.loc[segment, 'avg_last_50000'] = x[-50000:].mean()\n X_tr.loc[segment, 'avg_first_10000'] = x[:10000].mean()\n X_tr.loc[segment, 'avg_last_10000'] = x[-10000:].mean()\n \n X_tr.loc[segment, 'min_first_50000'] = x[:50000].min()\n X_tr.loc[segment, 'min_last_50000'] = x[-50000:].min()\n X_tr.loc[segment, 'min_first_10000'] = x[:10000].min()\n X_tr.loc[segment, 'min_last_10000'] = x[-10000:].min()\n \n X_tr.loc[segment, 'max_first_50000'] = x[:50000].max()\n X_tr.loc[segment, 'max_last_50000'] = x[-50000:].max()\n X_tr.loc[segment, 'max_first_10000'] = x[:10000].max()\n X_tr.loc[segment, 'max_last_10000'] = x[-10000:].max()\n \n X_tr.loc[segment, 'max_to_min'] = x.max() / np.abs(x.min())\n X_tr.loc[segment, 'max_to_min_diff'] = x.max() - np.abs(x.min())\n X_tr.loc[segment, 'count_big'] = len(x[np.abs(x) > 500])\n X_tr.loc[segment, 'sum'] = x.sum()\n \n X_tr.loc[segment, 'mean_change_rate_first_50000'] = calc_change_rate(x[:50000])\n X_tr.loc[segment, 'mean_change_rate_last_50000'] = calc_change_rate(x[-50000:])\n X_tr.loc[segment, 'mean_change_rate_first_10000'] = calc_change_rate(x[:10000])\n X_tr.loc[segment, 'mean_change_rate_last_10000'] = calc_change_rate(x[-10000:])\n \n X_tr.loc[segment, 'q95'] = np.quantile(x, 0.95)\n X_tr.loc[segment, 'q99'] = np.quantile(x, 0.99)\n X_tr.loc[segment, 'q05'] = np.quantile(x, 0.05)\n X_tr.loc[segment, 'q01'] = np.quantile(x, 0.01)\n \n X_tr.loc[segment, 'abs_q95'] = np.quantile(np.abs(x), 0.95)\n X_tr.loc[segment, 'abs_q99'] = np.quantile(np.abs(x), 0.99)\n X_tr.loc[segment, 'abs_q05'] = np.quantile(np.abs(x), 0.05)\n X_tr.loc[segment, 'abs_q01'] = np.quantile(np.abs(x), 0.01)\n \n X_tr.loc[segment, 'trend'] = add_trend_feature(x)\n X_tr.loc[segment, 'abs_trend'] = add_trend_feature(x, abs_values=True)\n X_tr.loc[segment, 'abs_mean'] = np.abs(x).mean()\n X_tr.loc[segment, 'abs_std'] = np.abs(x).std()\n \n X_tr.loc[segment, 'mad'] = x.mad()\n X_tr.loc[segment, 'kurt'] = x.kurtosis()\n X_tr.loc[segment, 'skew'] = x.skew()\n X_tr.loc[segment, 'med'] = x.median()\n \n X_tr.loc[segment, 'Hilbert_mean'] = np.abs(hilbert(x)).mean()\n X_tr.loc[segment, 'Hann_window_mean'] = (convolve(x, hann(150), mode='same') / sum(hann(150))).mean()\n X_tr.loc[segment, 'classic_sta_lta1_mean'] = classic_sta_lta(x, 500, 10000).mean()\n X_tr.loc[segment, 'classic_sta_lta2_mean'] = classic_sta_lta(x, 5000, 100000).mean()\n X_tr.loc[segment, 'classic_sta_lta3_mean'] = classic_sta_lta(x, 3333, 6666).mean()\n X_tr.loc[segment, 'classic_sta_lta4_mean'] = classic_sta_lta(x, 10000, 25000).mean()\n X_tr.loc[segment, 'classic_sta_lta5_mean'] = classic_sta_lta(x, 50, 1000).mean()\n X_tr.loc[segment, 'classic_sta_lta6_mean'] = classic_sta_lta(x, 100, 5000).mean()\n X_tr.loc[segment, 'classic_sta_lta7_mean'] = classic_sta_lta(x, 333, 666).mean()\n X_tr.loc[segment, 'classic_sta_lta8_mean'] = classic_sta_lta(x, 4000, 10000).mean()\n X_tr.loc[segment, 'Moving_average_700_mean'] = x.rolling(window=700).mean().mean(skipna=True)\n ewma = pd.Series.ewm\n X_tr.loc[segment, 'exp_Moving_average_300_mean'] = (ewma(x, span=300).mean()).mean(skipna=True)\n X_tr.loc[segment, 'exp_Moving_average_3000_mean'] = ewma(x, span=3000).mean().mean(skipna=True)\n X_tr.loc[segment, 'exp_Moving_average_30000_mean'] = ewma(x, span=30000).mean().mean(skipna=True)\n no_of_std = 3\n X_tr.loc[segment, 'MA_700MA_std_mean'] = x.rolling(window=700).std().mean()\n X_tr.loc[segment,'MA_700MA_BB_high_mean'] = (X_tr.loc[segment, 'Moving_average_700_mean'] + no_of_std * X_tr.loc[segment, 'MA_700MA_std_mean']).mean()\n X_tr.loc[segment,'MA_700MA_BB_low_mean'] = (X_tr.loc[segment, 'Moving_average_700_mean'] - no_of_std * X_tr.loc[segment, 'MA_700MA_std_mean']).mean()\n X_tr.loc[segment, 'MA_400MA_std_mean'] = x.rolling(window=400).std().mean()\n X_tr.loc[segment,'MA_400MA_BB_high_mean'] = (X_tr.loc[segment, 'Moving_average_700_mean'] + no_of_std * X_tr.loc[segment, 'MA_400MA_std_mean']).mean()\n X_tr.loc[segment,'MA_400MA_BB_low_mean'] = (X_tr.loc[segment, 'Moving_average_700_mean'] - no_of_std * X_tr.loc[segment, 'MA_400MA_std_mean']).mean()\n X_tr.loc[segment, 'MA_1000MA_std_mean'] = x.rolling(window=1000).std().mean()\n X_tr.drop('Moving_average_700_mean', axis=1, inplace=True)\n \n X_tr.loc[segment, 'iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n X_tr.loc[segment, 'q999'] = np.quantile(x,0.999)\n X_tr.loc[segment, 'q001'] = np.quantile(x,0.001)\n X_tr.loc[segment, 'ave10'] = stats.trim_mean(x, 0.1)\n\n for windows in [10, 100, 1000]:\n x_roll_std = x.rolling(windows).std().dropna().values\n x_roll_mean = x.rolling(windows).mean().dropna().values\n \n X_tr.loc[segment, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n X_tr.loc[segment, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n X_tr.loc[segment, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n X_tr.loc[segment, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n X_tr.loc[segment, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n X_tr.loc[segment, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n X_tr.loc[segment, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n X_tr.loc[segment, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n X_tr.loc[segment, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n X_tr.loc[segment, 'av_change_rate_roll_std_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n X_tr.loc[segment, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n \n X_tr.loc[segment, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n X_tr.loc[segment, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n X_tr.loc[segment, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n X_tr.loc[segment, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n X_tr.loc[segment, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n X_tr.loc[segment, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n X_tr.loc[segment, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n X_tr.loc[segment, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n X_tr.loc[segment, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n X_tr.loc[segment, 'av_change_rate_roll_mean_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n X_tr.loc[segment, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()\n\n# Fill missing and infinite values\nmeans_dict = {}\nfor col in X_tr.columns:\n if X_tr[col].isnull().any():\n mean_value = X_tr.loc[X_tr[col] != -np.inf, col].mean()\n X_tr.loc[X_tr[col] == -np.inf, col] = mean_value\n X_tr[col] = X_tr[col].fillna(mean_value)\n means_dict[col] = mean_value\n \nprint(\"Original shape: {}, final shape: {}\".format(train.shape, X_tr.shape))\nX_tr.head(3)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d2422ab1bce53e44b2c6e430c19439ddf23ea4ab"},"cell_type":"markdown","source":"Extract features for the test set:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"7bc7ef67b0ae5ba4812efa2d7c41f14ae3c85914"},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame(columns=X_tr.columns, dtype=np.float64, index=submission.index)\n\nfor i, seg_id in enumerate(X_test.index):\n seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n \n x = pd.Series(seg['acoustic_data'].values)\n X_test.loc[seg_id, 'mean'] = x.mean()\n X_test.loc[seg_id, 'std'] = x.std()\n X_test.loc[seg_id, 'max'] = x.max()\n X_test.loc[seg_id, 'min'] = x.min()\n \n X_test.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(x))\n X_test.loc[seg_id, 'mean_change_rate'] = calc_change_rate(x)\n X_test.loc[seg_id, 'abs_max'] = np.abs(x).max()\n X_test.loc[seg_id, 'std_first_50000'] = x[:50000].std()\n X_test.loc[seg_id, 'std_last_50000'] = x[-50000:].std()\n X_test.loc[seg_id, 'std_first_10000'] = x[:10000].std()\n X_test.loc[seg_id, 'std_last_10000'] = x[-10000:].std()\n \n X_test.loc[seg_id, 'avg_first_50000'] = x[:50000].mean()\n X_test.loc[seg_id, 'avg_last_50000'] = x[-50000:].mean()\n X_test.loc[seg_id, 'avg_first_10000'] = x[:10000].mean()\n X_test.loc[seg_id, 'avg_last_10000'] = x[-10000:].mean()\n \n X_test.loc[seg_id, 'min_first_50000'] = x[:50000].min()\n X_test.loc[seg_id, 'min_last_50000'] = x[-50000:].min()\n X_test.loc[seg_id, 'min_first_10000'] = x[:10000].min()\n X_test.loc[seg_id, 'min_last_10000'] = x[-10000:].min()\n \n X_test.loc[seg_id, 'max_first_50000'] = x[:50000].max()\n X_test.loc[seg_id, 'max_last_50000'] = x[-50000:].max()\n X_test.loc[seg_id, 'max_first_10000'] = x[:10000].max()\n X_test.loc[seg_id, 'max_last_10000'] = x[-10000:].max()\n \n X_test.loc[seg_id, 'max_to_min'] = x.max() / np.abs(x.min())\n X_test.loc[seg_id, 'max_to_min_diff'] = x.max() - np.abs(x.min())\n X_test.loc[seg_id, 'count_big'] = len(x[np.abs(x) > 500])\n X_test.loc[seg_id, 'sum'] = x.sum()\n \n X_test.loc[seg_id, 'mean_change_rate_first_50000'] = calc_change_rate(x[:50000])\n X_test.loc[seg_id, 'mean_change_rate_last_50000'] = calc_change_rate(x[-50000:])\n X_test.loc[seg_id, 'mean_change_rate_first_10000'] = calc_change_rate(x[:10000])\n X_test.loc[seg_id, 'mean_change_rate_last_10000'] = calc_change_rate(x[-10000:])\n \n X_test.loc[seg_id, 'q95'] = np.quantile(x,0.95)\n X_test.loc[seg_id, 'q99'] = np.quantile(x,0.99)\n X_test.loc[seg_id, 'q05'] = np.quantile(x,0.05)\n X_test.loc[seg_id, 'q01'] = np.quantile(x,0.01)\n \n X_test.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(x), 0.95)\n X_test.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(x), 0.99)\n X_test.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(x), 0.05)\n X_test.loc[seg_id, 'abs_q01'] = np.quantile(np.abs(x), 0.01)\n \n X_test.loc[seg_id, 'trend'] = add_trend_feature(x)\n X_test.loc[seg_id, 'abs_trend'] = add_trend_feature(x, abs_values=True)\n X_test.loc[seg_id, 'abs_mean'] = np.abs(x).mean()\n X_test.loc[seg_id, 'abs_std'] = np.abs(x).std()\n \n X_test.loc[seg_id, 'mad'] = x.mad()\n X_test.loc[seg_id, 'kurt'] = x.kurtosis()\n X_test.loc[seg_id, 'skew'] = x.skew()\n X_test.loc[seg_id, 'med'] = x.median()\n \n X_test.loc[seg_id, 'Hilbert_mean'] = np.abs(hilbert(x)).mean()\n X_test.loc[seg_id, 'Hann_window_mean'] = (convolve(x, hann(150), mode='same') / sum(hann(150))).mean()\n X_test.loc[seg_id, 'classic_sta_lta1_mean'] = classic_sta_lta(x, 500, 10000).mean()\n X_test.loc[seg_id, 'classic_sta_lta2_mean'] = classic_sta_lta(x, 5000, 100000).mean()\n X_test.loc[seg_id, 'classic_sta_lta3_mean'] = classic_sta_lta(x, 3333, 6666).mean()\n X_test.loc[seg_id, 'classic_sta_lta4_mean'] = classic_sta_lta(x, 10000, 25000).mean()\n X_test.loc[seg_id, 'classic_sta_lta5_mean'] = classic_sta_lta(x, 50, 1000).mean()\n X_test.loc[seg_id, 'classic_sta_lta6_mean'] = classic_sta_lta(x, 100, 5000).mean()\n X_test.loc[seg_id, 'classic_sta_lta7_mean'] = classic_sta_lta(x, 333, 666).mean()\n X_test.loc[seg_id, 'classic_sta_lta8_mean'] = classic_sta_lta(x, 4000, 10000).mean()\n X_test.loc[seg_id, 'Moving_average_700_mean'] = x.rolling(window=700).mean().mean(skipna=True)\n ewma = pd.Series.ewm\n X_test.loc[seg_id, 'exp_Moving_average_300_mean'] = (ewma(x, span=300).mean()).mean(skipna=True)\n X_test.loc[seg_id, 'exp_Moving_average_3000_mean'] = ewma(x, span=3000).mean().mean(skipna=True)\n X_test.loc[seg_id, 'exp_Moving_average_30000_mean'] = ewma(x, span=30000).mean().mean(skipna=True)\n no_of_std = 3\n X_test.loc[seg_id, 'MA_700MA_std_mean'] = x.rolling(window=700).std().mean()\n X_test.loc[seg_id,'MA_700MA_BB_high_mean'] = (X_test.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X_test.loc[seg_id, 'MA_700MA_std_mean']).mean()\n X_test.loc[seg_id,'MA_700MA_BB_low_mean'] = (X_test.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X_test.loc[seg_id, 'MA_700MA_std_mean']).mean()\n X_test.loc[seg_id, 'MA_400MA_std_mean'] = x.rolling(window=400).std().mean()\n X_test.loc[seg_id,'MA_400MA_BB_high_mean'] = (X_test.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X_test.loc[seg_id, 'MA_400MA_std_mean']).mean()\n X_test.loc[seg_id,'MA_400MA_BB_low_mean'] = (X_test.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X_test.loc[seg_id, 'MA_400MA_std_mean']).mean()\n X_test.loc[seg_id, 'MA_1000MA_std_mean'] = x.rolling(window=1000).std().mean()\n X_test.drop('Moving_average_700_mean', axis=1, inplace=True)\n \n X_test.loc[seg_id, 'iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n X_test.loc[seg_id, 'q999'] = np.quantile(x,0.999)\n X_test.loc[seg_id, 'q001'] = np.quantile(x,0.001)\n X_test.loc[seg_id, 'ave10'] = stats.trim_mean(x, 0.1)\n \n for windows in [10, 100, 1000]:\n x_roll_std = x.rolling(windows).std().dropna().values\n x_roll_mean = x.rolling(windows).mean().dropna().values\n \n X_test.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n X_test.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n X_test.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n X_test.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n X_test.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n X_test.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n X_test.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n X_test.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n X_test.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n X_test.loc[seg_id, 'av_change_rate_roll_std_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n X_test.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n \n X_test.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n X_test.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n X_test.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n X_test.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n X_test.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n X_test.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n X_test.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n X_test.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n X_test.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n X_test.loc[seg_id, 'av_change_rate_roll_mean_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n X_test.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()\n\nfor col in X_test.columns:\n if X_test[col].isnull().any():\n X_test.loc[X_test[col] == -np.inf, col] = means_dict[col]\n X_test[col] = X_test[col].fillna(means_dict[col])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"649591c22783ab1ca5c70a9d0e93dd9ac52b1b82"},"cell_type":"markdown","source":"

Baseline

\n\nWe will be using a gradient boosting decision tree model (lightgbm) as our main estimator in this notebook, which is implemented with the following class. Note that we are using a 10-fold validation scheme with the mean absolute error (competition metric). Our benchmark will be the validation score using all the 137 features."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"3957da7d73fe80afef261da08a22284b68763e09"},"cell_type":"code","source":"class BoostingModel(object):\n def __init__(self, params, verbose=False, early_stopping=50, num_folds=10):\n self.hyperparams = params\n self.stop = early_stopping\n self.num_folds = num_folds\n self.verbose = verbose\n self.seed = RANDOM_SEED\n self.hyperparams['random_seed'] = self.seed\n self.hyperparams['objective'] = 'regression_l1'\n \n def train(self, X, y):\n \"\"\"Train num_folds estimators (boosters) using KFold.\n \n Arguments:\n X: pandas dataframe with features\n y: series or dataframe with target values for regression\n \"\"\"\n self.estimators = []\n self.feat_importance = pd.DataFrame({'feature': X.columns})\n self.folds = KFold(self.num_folds, shuffle=True, random_state=self.seed)\n for i, (train_index, valid_index) in enumerate(self.folds.split(X, y)):\n estimator = lgb.LGBMRegressor(**self.hyperparams)\n estimator.fit(X.iloc[train_index], y.iloc[train_index],\n early_stopping_rounds=self.stop, verbose=self.verbose,\n eval_set=[(X.iloc[train_index], y.iloc[train_index]),\n (X.iloc[valid_index], y.iloc[valid_index])])\n \n fn = \"fold_\" + str(i+1)\n self.feat_importance[fn] = estimator.booster_.feature_importance(importance_type='gain')\n self.estimators.append(estimator)\n self.feat_importance['fold_mean'] = self.feat_importance.mean(axis=1)\n self.feat_importance['fold_std'] = self.feat_importance.std(axis=1)\n self.feat_importance.sort_values(by='fold_mean', ascending=False, inplace=True)\n \n def predict(self, X):\n \"\"\"Makes predictions for test set with each estimator and returns the average.\"\"\"\n predictions = np.zeros(X.shape[0])\n for estimator in self.estimators:\n predictions += estimator.predict(X) / len(self.estimators)\n return predictions\n \n def predict_oof(self, X):\n \"\"\"Returns the predictions for validation data in each fold.\"\"\"\n oof_predictions = np.zeros(X.shape[0])\n for i, (train_index, valid_index) in enumerate(self.folds.split(X)):\n oof_predictions[valid_index] = self.estimators[i].predict(X.iloc[valid_index])\n return oof_predictions\n \n def validation_score(self, X, y):\n oof_predictions = self.predict_oof(X)\n return mean_absolute_error(y, oof_predictions)\n\nparams = {\n 'learning_rate': 0.005,\n 'num_leaves': 8,\n 'max_depth': 8,\n 'feature_fraction': 0.8,\n 'subsample': 0.9,\n 'lambda_l1': 0,\n 'lambda_l2': 0.4,\n 'min_data_in_leaf': 40,\n 'min_gain_to_split': 0.001,\n 'boosting': 'gbdt',\n 'verbosity': -1,\n 'n_estimators': 20000\n}\nmodel = BoostingModel(params)\nmodel.train(X_tr, y_tr)\nbenchmark = model.validation_score(X_tr, y_tr)\nprint(\"Baseline score with all features (mae): {:.6f}\".format(benchmark))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5ba80d1e4e7f89239f8443d695017ebb88962589"},"cell_type":"markdown","source":"

1. Mutual Information

\n\nMutual information is a quantity that measures a relationship between two random variables. In particular, it measures how much information is communicated, on average, in one random variable about another. It's considered a filter method since we don't need any particular learning model to select features. It can be used with SelectKBest and a given number of features:"},{"metadata":{"trusted":true,"_uuid":"abd9b024d3ff0b41e2aa55f7a4d4b72f6c03212f","_kg_hide-input":false},"cell_type":"code","source":"model_ = BoostingModel(params)\nscores = []\nfor k in range(6, 60, 2):\n values = SelectKBest(mutual_info_regression, k=k).fit_transform(X_tr, y_tr)\n X = pd.DataFrame(values, columns=[i for i in range(values.shape[1])])\n model_.train(X, y_tr)\n scores.append(model_.validation_score(X, y_tr))\n\nlineplot(list(range(6, 60, 2)), scores, 'FS with Mutual information')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"252e6dd2edaa948c379befabfc8b2c850bea1c42"},"cell_type":"markdown","source":"

2. Feature importance

\n\nAnother simple approach is to select the top *n* features according to the feature importance in boosting trees or random forests."},{"metadata":{"trusted":true,"_uuid":"903e85d7ffaff20b0681162e27890965124140b9"},"cell_type":"code","source":"scores = []\nfor n in range(6, 60, 2):\n features = model.feat_importance[:n].feature\n model_.train(X_tr[features], y_tr)\n scores.append(model_.validation_score(X_tr[features], y_tr))\n \nlineplot(list(range(6, 60, 2)), scores)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"96fd6d4ef0fdfc461fcc461b6da38ce800d31e0b"},"cell_type":"markdown","source":"The x-axis represents the number of features according to their feature importance (e.g. at 10 the model is being trained with the 10 most important features)."},{"metadata":{"_uuid":"fef41be47786622cd030c9ff337d2610b136ac72"},"cell_type":"markdown","source":"

3. Permutation Importance

\n\nAs explained in this [article](https://explained.ai/rf-importance/index.html) [4], the standard feature importance in tree models is a good indicator for feature selection, but it can be biased. Features with different ranges of values or high cardinality might have their importance inflated. Permutation importance should solve these problems and works as following:\n\n* Record a baseline accuracy by passing a validation set or the out-of-bag samples through the model.\n* Permute the column values of a single predictor feature and then pass all test samples back through the model and recompute the accuracy or R2.\n* The importance of that feature is the difference between the baseline and the drop in overall accuracy or R2 caused by permuting the column.\n\nNote that this strategy does not require retraining the model after permuting each column; we just have to re-run the perturbed test samples through the already-trained model.\nThe permutation mechanism is much more computationally expensive than the mean decrease in impurity mechanism (standard feature importance), but the results are more reliable. It's also very simple to implement:"},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"d564f44d39a2902e0edb231948306a13e659f94e"},"cell_type":"code","source":"diff = {}\nfor col in X_tr.columns:\n save = X_tr[col].values.copy()\n X_tr[col] = np.random.permutation(X_tr[col])\n diff[col] = model.validation_score(X_tr, y_tr) - benchmark\n X_tr[col] = save","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a1946d0916cc882a4a3940a52633874f0145fbbd"},"cell_type":"markdown","source":"Now it's possible to rank the features according to their baseline score difference and train the model with different numbers of features:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"c5da3e1d398831d6f85720c58feab041bf36ef8e"},"cell_type":"code","source":"scores = []\nfor num_feats in range(6, 60, 2):\n features = sorted(diff, key=diff.get, reverse=True)[:num_feats]\n model_.train(X_tr[features], y_tr)\n scores.append(model_.validation_score(X_tr[features], y_tr))\n \nlineplot(list(range(6, 60, 2)), scores)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a133196eadbb434c8c722855dd674997f6fedcf6"},"cell_type":"markdown","source":"We have some improvement, so let's try a submission:"},{"metadata":{"trusted":true,"_uuid":"e502edffaf0a26f8ac1ae352dbc362faf913cb75"},"cell_type":"code","source":"num_feats = scores.index(min(scores))\nfeatures = sorted(diff, key=diff.get, reverse=True)[:num_feats]\nmodel_.train(X_tr[features], y_tr)\nsubmission['time_to_failure'] = model_.predict(X_test[features])\nsubmission.to_csv('submission_pi.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"efb946c43302b557b738bb6a93ff0e7516b8c6ac"},"cell_type":"markdown","source":"

4. Recursive Feature Elimination

\n\nGiven an estimator that assigns weights to features (e.g., the coefficients of a linear model), the goal of recursive feature elimination is to select features by recursively considering smaller and smaller sets of features. First, the estimator is trained on the initial set of features and the importance of each feature is obtained either through a coefficient or a feature importances attribute. Then, the least important features are pruned from current set of features. That procedure is recursively repeated on the pruned set until the desired number of features to select is eventually reached. This can also be implemented with the [RFE class](https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.RFE.html) in Scikit-learn."},{"metadata":{"trusted":true,"_uuid":"13c8a9dd489b72755bfc6549e3c2f66a47715ca6","_kg_hide-input":false},"cell_type":"code","source":"scores = []\nfeatures = []\nimportance = model.feat_importance.feature.values\nX = X_tr.copy()\nfor i in range(X_tr.shape[1] - 1):\n X.drop(importance[-1], axis=1, inplace=True)\n model_.train(X, y_tr)\n scores.append(model_.validation_score(X, y_tr))\n features.append(list(X.columns))\n importance = model_.feat_importance.feature.values\n\nlineplot(list(range(X_tr.shape[1], 1, -1)), scores)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1d3b6bf3f9b95274ee393c196aa9c872110b4df4"},"cell_type":"markdown","source":"Another submission:"},{"metadata":{"trusted":true,"_uuid":"0dd4619f8283417b4aece7e5f5e9d5d951f89b93"},"cell_type":"code","source":"features = features[scores.index(min(scores))]\nmodel_.train(X_tr[features], y_tr)\nsubmission['time_to_failure'] = model_.predict(X_test[features])\nsubmission.to_csv('submission_rfe.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bda00b0f40006bb78ed440b9a18271ec12198b50"},"cell_type":"markdown","source":"

5. Boruta

\n\nBoruta is a library for feature selection in R, but it was later implemented in python. It works with the concept of shadow features, which are similar to the permutation importance.\n\n* Firstly, it adds randomness to the given data set by creating shuffled copies of all features (which are called shadow features). \n* Then, it trains a estimator on the extended dataset and applies a feature importance measure to evaluate the importance of each feature where higher means more important.\n* At every iteration, it checks whether a real feature has a higher importance than the best of its shadow features (i.e. whether the feature has a higher Z-score than the maximum Z-score of its shadow features) and constantly removes features which are deemed highly unimportant.\n* Finally, the algorithm stops either when all features get confirmed or rejected or it reaches a specified limit of random forest runs."},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"8fd40868eb262e2471f2b7c9a76629ed2bdcc303"},"cell_type":"code","source":"estimator = RandomForestRegressor(n_estimators=500, n_jobs=4, max_depth=4)\nboruta_selector = BorutaPy(estimator, n_estimators='auto', verbose=0)\nboruta_selector.fit(X_tr.values, y_tr.values.flatten())\n\nboruta_df = pd.DataFrame({'feature': X_tr.columns, 'rank': boruta_selector.ranking_})\nfeats = boruta_df[boruta_df['rank'] == 1].feature\nmodel_.train(X_tr[feats], y_tr)\nprint(\"Boruta validation score: {:.6f}\".format(model_.validation_score(X_tr[feats], y_tr)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6a2df1b2fcf867e55dc9db501d13932b56dea749"},"cell_type":"markdown","source":"

6. Null importances

\n\nThe key idea here is not to shuffle or drop a feature, but the target variable. We start by fitting the model over several runs on a shuffled version of the target and saving the feature importance. This gives us the null importance distribution, which will be compared to the original importance (without target permutation). For a complete description check the [original notebook](https://www.kaggle.com/ogrellier/feature-selection-with-null-importances) from olivier or white paper [5]."},{"metadata":{"trusted":true,"_uuid":"4e39fe54de4532889c4e44b1b3ee531135b12346"},"cell_type":"code","source":"actual_importance = model.feat_importance[['feature', 'fold_mean']]\nnull_importance = pd.DataFrame()\nfor i in range(200):\n # Get current run importances\n shuffled_target = y_tr.copy().sample(frac=1.0)\n model_.train(X_tr, shuffled_target)\n importance = model_.feat_importance[['feature', 'fold_mean']]\n importance['run'] = i + 1\n # Concat the latest importances with the old ones\n null_importance = pd.concat([null_importance, importance])\n \nnull_importance.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9f2cee04e7f8dd84721ffd39f8530b075e286fc5","_kg_hide-input":true},"cell_type":"code","source":"feature_scores = []\nfor _f in actual_importance['feature'].unique():\n f_null_imps_gain = null_importance.loc[null_importance['feature'] == _f, 'fold_mean'].values\n f_act_imps_gain = actual_importance.loc[actual_importance['feature'] == _f, 'fold_mean'].mean()\n gain_score = np.log(1e-10 + f_act_imps_gain / (1 + np.percentile(f_null_imps_gain, 75))) # Avoid didvide by zero\n feature_scores.append((_f, gain_score))\n\nscores_df = pd.DataFrame(feature_scores, columns=['feature', 'gain_score'])\nscores_df.sort_values('gain_score', ascending=False, inplace=True)\n\nfor num_feats in range(6, 60, 2):\n features = scores_df.loc[:num_feats, 'feature']\n model_.train(X_tr[features], y_tr)\n scores.append(model_.validation_score(X_tr[features], y_tr))\n \nlineplot(list(range(6, 60, 2)), scores)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6514f16b3fac973777d0af59475b0c36b0082f0c"},"cell_type":"markdown","source":"

7. Lasso

\n\nLasso is a linear model penalized with the L1 norm, which has sparse solutions: many of their estimated coefficients are zero. It can be used along with SelectFromModel to select the non-zero coefficients and reduce the number of features in a dataset."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"991bac46fc506efa16bc82bb38209fe43584f295"},"cell_type":"code","source":"scores = []\nfor num_feats in range(6, 60, 2):\n values = SelectFromModel(LassoCV(cv=5), max_features=num_feats,\n threshold=-np.inf).fit_transform(X_tr, y_tr)\n X = pd.DataFrame(values, columns=[i for i in range(num_feats)])\n model_.train(X, y_tr)\n scores.append(model_.validation_score(X, y_tr))\n \nlineplot(list(range(6, 60, 2)), scores, \"Feature selection with Lasso\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0df474b85c518b74739a62fda1641313c9acabec"},"cell_type":"markdown","source":"

References

\n\n[1] [An Introduction to Variable and Feature Selection](http://www.jmlr.org/papers/volume3/guyon03a/guyon03a.pdf), JMLR, 2003.\n\n[2] Computational Methods of Feature Selection, Huan Liu, Hiroshi Motoda, 2017.\n\n[3] [Introduction to Feature Selection Methods](https://www.analyticsvidhya.com/blog/2016/12/introduction-to-feature-selection-methods-with-an-example-or-how-to-select-the-right-variables/), Analytics Vidhya, 2016.\n\n[4] [Beware Default Random Forest Importances](https://explained.ai/rf-importance/index.html), 2018.\n\n[5] [Permutation importance: a corrected feature importance measure](https://academic.oup.com/bioinformatics/article/26/10/1340/193348), 2010."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}