{"cells":[{"metadata":{},"cell_type":"markdown","source":"## General information\n\nIn this kernel I focus on various ML approaches:\n- adversarial validation for selecting stable features;\n- comparing several validation schemes;\n- trying several different models to find which ones are better for this competition;\n- some approaches to feature selection and model interpretation.\n\n\n\n*Work in progress*"},{"metadata":{},"cell_type":"markdown","source":"Importing libraries"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install -U vega_datasets notebook vega","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom tqdm import tqdm_notebook\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.svm import NuSVR, SVR\nfrom sklearn.metrics import mean_absolute_error\npd.options.display.precision = 15\n\nimport lightgbm as lgb\nimport xgboost as xgb\nimport time\nimport datetime\nfrom catboost import CatBoostRegressor\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold, GroupKFold, GridSearchCV, train_test_split\nfrom sklearn import metrics\nfrom sklearn import linear_model\nimport gc\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nimport eli5\nimport shap\nfrom IPython.display import HTML\nimport json\nimport altair as alt\n\nimport networkx as nx\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nalt.renderers.enable('notebook')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Functions used in this kernel\nThey are in the hidden cell below."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import os\nimport time\nimport datetime\nimport json\nimport gc\nfrom numba import jit\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm import tqdm_notebook\n\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom catboost import CatBoostRegressor, CatBoostClassifier\nfrom sklearn import metrics\n\nfrom itertools import product\n\nimport altair as alt\nfrom altair.vega import v5\nfrom IPython.display import HTML\n\n# using ideas from this kernel: https://www.kaggle.com/notslush/altair-visualization-2018-stackoverflow-survey\ndef prepare_altair():\n \"\"\"\n Helper function to prepare altair for working.\n \"\"\"\n\n vega_url = 'https://cdn.jsdelivr.net/npm/vega@' + v5.SCHEMA_VERSION\n vega_lib_url = 'https://cdn.jsdelivr.net/npm/vega-lib'\n vega_lite_url = 'https://cdn.jsdelivr.net/npm/vega-lite@' + alt.SCHEMA_VERSION\n vega_embed_url = 'https://cdn.jsdelivr.net/npm/vega-embed@3'\n noext = \"?noext\"\n \n paths = {\n 'vega': vega_url + noext,\n 'vega-lib': vega_lib_url + noext,\n 'vega-lite': vega_lite_url + noext,\n 'vega-embed': vega_embed_url + noext\n }\n \n workaround = f\"\"\" requirejs.config({{\n baseUrl: 'https://cdn.jsdelivr.net/npm/',\n paths: {paths}\n }});\n \"\"\"\n \n return workaround\n \n\ndef add_autoincrement(render_func):\n # Keep track of unique IDs\n cache = {}\n def wrapped(chart, id=\"vega-chart\", autoincrement=True):\n if autoincrement:\n if id in cache:\n counter = 1 + cache[id]\n cache[id] = counter\n else:\n cache[id] = 0\n actual_id = id if cache[id] == 0 else id + '-' + str(cache[id])\n else:\n if id not in cache:\n cache[id] = 0\n actual_id = id\n return render_func(chart, id=actual_id)\n # Cache will stay outside and \n return wrapped\n \n\n@add_autoincrement\ndef render(chart, id=\"vega-chart\"):\n \"\"\"\n Helper function to plot altair visualizations.\n \"\"\"\n chart_str = \"\"\"\n \n \"\"\"\n return HTML(\n chart_str.format(\n id=id,\n chart=json.dumps(chart) if isinstance(chart, dict) else chart.to_json(indent=None)\n )\n )\n \n\ndef reduce_mem_usage(df, verbose=True):\n numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n start_mem = df.memory_usage().sum() / 1024**2\n for col in df.columns:\n col_type = df[col].dtypes\n if col_type in numerics:\n c_min = df[col].min()\n c_max = df[col].max()\n if str(col_type)[:3] == 'int':\n if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n df[col] = df[col].astype(np.int8)\n elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n df[col] = df[col].astype(np.int16)\n elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n df[col] = df[col].astype(np.int32)\n elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n df[col] = df[col].astype(np.int64)\n else:\n c_prec = df[col].apply(lambda x: np.finfo(x).precision).max()\n if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max and c_prec == np.finfo(np.float16).precision:\n df[col] = df[col].astype(np.float16)\n elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max and c_prec == np.finfo(np.float32).precision:\n df[col] = df[col].astype(np.float32)\n else:\n df[col] = df[col].astype(np.float64)\n end_mem = df.memory_usage().sum() / 1024**2\n if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n return df\n \n\n@jit\ndef fast_auc(y_true, y_prob):\n \"\"\"\n fast roc_auc computation: https://www.kaggle.com/c/microsoft-malware-prediction/discussion/76013\n \"\"\"\n y_true = np.asarray(y_true)\n y_true = y_true[np.argsort(y_prob)]\n nfalse = 0\n auc = 0\n n = len(y_true)\n for i in range(n):\n y_i = y_true[i]\n nfalse += (1 - y_i)\n auc += y_i * nfalse\n auc /= (nfalse * (n - nfalse))\n return auc\n\n\ndef eval_auc(y_true, y_pred):\n \"\"\"\n Fast auc eval function for lgb.\n \"\"\"\n return 'auc', fast_auc(y_true, y_pred), True\n\n\ndef group_mean_log_mae(y_true, y_pred, types, floor=1e-9):\n \"\"\"\n Fast metric computation for this competition: https://www.kaggle.com/c/champs-scalar-coupling\n Code is from this kernel: https://www.kaggle.com/uberkinder/efficient-metric\n \"\"\"\n maes = (y_true-y_pred).abs().groupby(types).mean()\n return np.log(maes.map(lambda x: max(x, floor))).mean()\n \n\ndef train_model_regression(X, X_test, y, params, folds=None, model_type='lgb', eval_metric='mae', columns=None, plot_feature_importance=False, model=None,\n verbose=10000, early_stopping_rounds=200, n_estimators=50000, splits=None, n_folds=3):\n \"\"\"\n A function to train a variety of regression models.\n Returns dictionary with oof predictions, test predictions, scores and, if necessary, feature importances.\n \n :params: X - training data, can be pd.DataFrame or np.ndarray (after normalizing)\n :params: X_test - test data, can be pd.DataFrame or np.ndarray (after normalizing)\n :params: y - target\n :params: folds - folds to split data\n :params: model_type - type of model to use\n :params: eval_metric - metric to use\n :params: columns - columns to use. If None - use all columns\n :params: plot_feature_importance - whether to plot feature importance of LGB\n :params: model - sklearn model, works only for \"sklearn\" model type\n \n \"\"\"\n columns = X.columns if columns is None else columns\n X_test = X_test[columns]\n splits = folds.split(X) if splits is None else splits\n n_splits = folds.n_splits if splits is None else n_folds\n \n # to set up scoring parameters\n metrics_dict = {'mae': {'lgb_metric_name': 'mae',\n 'catboost_metric_name': 'MAE',\n 'sklearn_scoring_function': metrics.mean_absolute_error},\n 'group_mae': {'lgb_metric_name': 'mae',\n 'catboost_metric_name': 'MAE',\n 'scoring_function': group_mean_log_mae},\n 'mse': {'lgb_metric_name': 'mse',\n 'catboost_metric_name': 'MSE',\n 'sklearn_scoring_function': metrics.mean_squared_error}\n }\n\n \n result_dict = {}\n \n # out-of-fold predictions on train data\n oof = np.zeros(len(X))\n \n # averaged predictions on train data\n prediction = np.zeros(len(X_test))\n \n # list of scores on folds\n scores = []\n feature_importance = pd.DataFrame()\n \n # split and train on folds\n for fold_n, (train_index, valid_index) in enumerate(splits):\n if verbose:\n print(f'Fold {fold_n + 1} started at {time.ctime()}')\n if type(X) == np.ndarray:\n X_train, X_valid = X[columns][train_index], X[columns][valid_index]\n y_train, y_valid = y[train_index], y[valid_index]\n else:\n X_train, X_valid = X[columns].iloc[train_index], X[columns].iloc[valid_index]\n y_train, y_valid = y.iloc[train_index], y.iloc[valid_index]\n \n if model_type == 'lgb':\n model = lgb.LGBMRegressor(**params, n_estimators = n_estimators, n_jobs = -1)\n model.fit(X_train, y_train, \n eval_set=[(X_train, y_train), (X_valid, y_valid)], eval_metric=metrics_dict[eval_metric]['lgb_metric_name'],\n verbose=verbose, early_stopping_rounds=early_stopping_rounds)\n \n y_pred_valid = model.predict(X_valid)\n y_pred = model.predict(X_test, num_iteration=model.best_iteration_)\n \n if model_type == 'xgb':\n train_data = xgb.DMatrix(data=X_train, label=y_train, feature_names=X.columns)\n valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X.columns)\n\n watchlist = [(train_data, 'train'), (valid_data, 'valid_data')]\n model = xgb.train(dtrain=train_data, num_boost_round=20000, evals=watchlist, early_stopping_rounds=200, verbose_eval=verbose, params=params)\n y_pred_valid = model.predict(xgb.DMatrix(X_valid, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n \n if model_type == 'sklearn':\n model = model\n model.fit(X_train, y_train)\n \n y_pred_valid = model.predict(X_valid).reshape(-1,)\n score = metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid)\n print(f'Fold {fold_n}. {eval_metric}: {score:.4f}.')\n print('')\n \n y_pred = model.predict(X_test).reshape(-1,)\n \n if model_type == 'cat':\n model = CatBoostRegressor(iterations=20000, eval_metric=metrics_dict[eval_metric]['catboost_metric_name'], **params,\n loss_function=metrics_dict[eval_metric]['catboost_metric_name'])\n model.fit(X_train, y_train, eval_set=(X_valid, y_valid), cat_features=[], use_best_model=True, verbose=False)\n\n y_pred_valid = model.predict(X_valid)\n y_pred = model.predict(X_test)\n \n oof[valid_index] = y_pred_valid.reshape(-1,)\n if eval_metric != 'group_mae':\n scores.append(metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid))\n else:\n scores.append(metrics_dict[eval_metric]['scoring_function'](y_valid, y_pred_valid, X_valid['type']))\n\n prediction += y_pred \n \n if model_type == 'lgb' and plot_feature_importance:\n # feature importance\n fold_importance = pd.DataFrame()\n fold_importance[\"feature\"] = columns\n fold_importance[\"importance\"] = model.feature_importances_\n fold_importance[\"fold\"] = fold_n + 1\n feature_importance = pd.concat([feature_importance, fold_importance], axis=0)\n\n prediction /= n_splits\n print('CV mean score: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores), np.std(scores)))\n \n result_dict['oof'] = oof\n result_dict['prediction'] = prediction\n result_dict['scores'] = scores\n \n if model_type == 'lgb':\n if plot_feature_importance:\n feature_importance[\"importance\"] /= n_splits\n cols = feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n by=\"importance\", ascending=False)[:50].index\n\n best_features = feature_importance.loc[feature_importance.feature.isin(cols)]\n\n plt.figure(figsize=(16, 12));\n sns.barplot(x=\"importance\", y=\"feature\", data=best_features.sort_values(by=\"importance\", ascending=False));\n plt.title('LGB Features (avg over folds)');\n \n result_dict['feature_importance'] = feature_importance\n \n return result_dict\n \n\n\ndef train_model_classification(X, X_test, y, params, folds, model_type='lgb', eval_metric='auc', columns=None, plot_feature_importance=False, model=None,\n verbose=10000, early_stopping_rounds=200, n_estimators=50000, splits=None, n_folds=3):\n \"\"\"\n A function to train a variety of regression models.\n Returns dictionary with oof predictions, test predictions, scores and, if necessary, feature importances.\n \n :params: X - training data, can be pd.DataFrame or np.ndarray (after normalizing)\n :params: X_test - test data, can be pd.DataFrame or np.ndarray (after normalizing)\n :params: y - target\n :params: folds - folds to split data\n :params: model_type - type of model to use\n :params: eval_metric - metric to use\n :params: columns - columns to use. If None - use all columns\n :params: plot_feature_importance - whether to plot feature importance of LGB\n :params: model - sklearn model, works only for \"sklearn\" model type\n \n \"\"\"\n columns = X.columns if columns is None else columns\n n_splits = folds.n_splits if splits is None else n_folds\n X_test = X_test[columns]\n \n # to set up scoring parameters\n metrics_dict = {'auc': {'lgb_metric_name': eval_auc,\n 'catboost_metric_name': 'AUC',\n 'sklearn_scoring_function': metrics.roc_auc_score},\n }\n \n result_dict = {}\n \n # out-of-fold predictions on train data\n oof = np.zeros((len(X), len(set(y.values))))\n \n # averaged predictions on train data\n prediction = np.zeros((len(X_test), oof.shape[1]))\n \n # list of scores on folds\n scores = []\n feature_importance = pd.DataFrame()\n \n # split and train on folds\n for fold_n, (train_index, valid_index) in enumerate(folds.split(X)):\n print(f'Fold {fold_n + 1} started at {time.ctime()}')\n if type(X) == np.ndarray:\n X_train, X_valid = X[columns][train_index], X[columns][valid_index]\n y_train, y_valid = y[train_index], y[valid_index]\n else:\n X_train, X_valid = X[columns].iloc[train_index], X[columns].iloc[valid_index]\n y_train, y_valid = y.iloc[train_index], y.iloc[valid_index]\n \n if model_type == 'lgb':\n model = lgb.LGBMClassifier(**params, n_estimators=n_estimators, n_jobs = -1)\n model.fit(X_train, y_train, \n eval_set=[(X_train, y_train), (X_valid, y_valid)], eval_metric=metrics_dict[eval_metric]['lgb_metric_name'],\n verbose=verbose, early_stopping_rounds=early_stopping_rounds)\n \n y_pred_valid = model.predict_proba(X_valid)\n y_pred = model.predict_proba(X_test, num_iteration=model.best_iteration_)\n \n if model_type == 'xgb':\n train_data = xgb.DMatrix(data=X_train, label=y_train, feature_names=X.columns)\n valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X.columns)\n\n watchlist = [(train_data, 'train'), (valid_data, 'valid_data')]\n model = xgb.train(dtrain=train_data, num_boost_round=n_estimators, evals=watchlist, early_stopping_rounds=early_stopping_rounds, verbose_eval=verbose, params=params)\n y_pred_valid = model.predict(xgb.DMatrix(X_valid, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n \n if model_type == 'sklearn':\n model = model\n model.fit(X_train, y_train)\n \n y_pred_valid = model.predict(X_valid).reshape(-1,)\n score = metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid)\n print(f'Fold {fold_n}. {eval_metric}: {score:.4f}.')\n print('')\n \n y_pred = model.predict_proba(X_test)\n \n if model_type == 'cat':\n model = CatBoostClassifier(iterations=n_estimators, eval_metric=metrics_dict[eval_metric]['catboost_metric_name'], **params,\n loss_function=metrics_dict[eval_metric]['catboost_metric_name'])\n model.fit(X_train, y_train, eval_set=(X_valid, y_valid), cat_features=[], use_best_model=True, verbose=False)\n\n y_pred_valid = model.predict(X_valid)\n y_pred = model.predict(X_test)\n \n oof[valid_index] = y_pred_valid\n scores.append(metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid[:, 1]))\n\n prediction += y_pred \n \n if model_type == 'lgb' and plot_feature_importance:\n # feature importance\n fold_importance = pd.DataFrame()\n fold_importance[\"feature\"] = columns\n fold_importance[\"importance\"] = model.feature_importances_\n fold_importance[\"fold\"] = fold_n + 1\n feature_importance = pd.concat([feature_importance, fold_importance], axis=0)\n\n prediction /= n_splits\n \n print('CV mean score: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores), np.std(scores)))\n \n result_dict['oof'] = oof\n result_dict['prediction'] = prediction\n result_dict['scores'] = scores\n \n if model_type == 'lgb':\n if plot_feature_importance:\n feature_importance[\"importance\"] /= n_splits\n cols = feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n by=\"importance\", ascending=False)[:50].index\n\n best_features = feature_importance.loc[feature_importance.feature.isin(cols)]\n\n plt.figure(figsize=(16, 12));\n sns.barplot(x=\"importance\", y=\"feature\", data=best_features.sort_values(by=\"importance\", ascending=False));\n plt.title('LGB Features (avg over folds)');\n \n result_dict['feature_importance'] = feature_importance\n \n return result_dict\n\n# setting up altair\nworkaround = prepare_altair()\nHTML(\"\".join((\n \"\",\n)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data loading and short overview"},{"metadata":{"trusted":true},"cell_type":"code","source":"file_folder = '../input/champs-scalar-coupling' if 'champs-scalar-coupling' in os.listdir('../input/') else '../input'\ntrain = pd.read_csv(f'{file_folder}/train.csv')\ntest = pd.read_csv(f'{file_folder}/test.csv')\nsub = pd.read_csv(f'{file_folder}/sample_submission.csv')\nstructures = pd.read_csv(f'{file_folder}/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I have already done EDA in my previous [kernel](https://www.kaggle.com/artgor/molecular-properties-eda-and-models), so no need to repeat it.\n\nStill I want to show two plots again to explain what I'll do further in my kernel."},{"metadata":{"trusted":true},"cell_type":"code","source":"type_count = train['type'].value_counts().reset_index().rename(columns={'type': 'count', 'index': 'type'})\nchart = alt.Chart(type_count).mark_bar().encode(\n x=alt.X(\"type:N\", axis=alt.Axis(title='type')),\n y=alt.Y('count:Q', axis=alt.Axis(title='Count')),\n tooltip=['type', 'count']\n).properties(title=\"Counts of type\", width=350).interactive()\nrender(chart)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.violinplot(x='type', y='scalar_coupling_constant', data=train);\nplt.title('Violinplot of scalar_coupling_constant by type');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we already know there are 8 types of coupling and they have very different values of scalar coupling constant. While it is possible to build one single model on this dataset, training 8 separate models seems to be more reasonable. And it could be worth working on different types separately.\n\nI'm going to generate my usual features on the whole dataset, but after this I'll work only with three types: `1JHN`, `2JNH` and `3JHN`. They are the smallest and will allow to do things faster. Let's see what insights will we able to get from this data."},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Feature generation\n\nIn the hidden cell below I generate features as in my kernel: https://www.kaggle.com/artgor/brute-force-feature-engineering\nThe only difference is that I don't apply label encoding for `type` variable."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"def map_atom_info(df, atom_idx):\n df = pd.merge(df, structures, how = 'left',\n left_on = ['molecule_name', f'atom_index_{atom_idx}'],\n right_on = ['molecule_name', 'atom_index'])\n \n df = df.drop('atom_index', axis=1)\n df = df.rename(columns={'atom': f'atom_{atom_idx}',\n 'x': f'x_{atom_idx}',\n 'y': f'y_{atom_idx}',\n 'z': f'z_{atom_idx}'})\n return df\n\ntrain = map_atom_info(train, 0)\ntrain = map_atom_info(train, 1)\n\ntest = map_atom_info(test, 0)\ntest = map_atom_info(test, 1)\n\ntrain_p_0 = train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = test[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = test[['x_1', 'y_1', 'z_1']].values\n\ntrain['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ntest['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)\ntrain['dist_x'] = (train['x_0'] - train['x_1']) ** 2\ntest['dist_x'] = (test['x_0'] - test['x_1']) ** 2\ntrain['dist_y'] = (train['y_0'] - train['y_1']) ** 2\ntest['dist_y'] = (test['y_0'] - test['y_1']) ** 2\ntrain['dist_z'] = (train['z_0'] - train['z_1']) ** 2\ntest['dist_z'] = (test['z_0'] - test['z_1']) ** 2\n\ntrain['type_0'] = train['type'].apply(lambda x: x[0])\ntest['type_0'] = test['type'].apply(lambda x: x[0])\n\ntrain = reduce_mem_usage(train)\ntest = reduce_mem_usage(test)\n\ndef create_features(df):\n df['molecule_couples'] = df.groupby('molecule_name')['id'].transform('count')\n df['molecule_dist_mean'] = df.groupby('molecule_name')['dist'].transform('mean')\n df['molecule_dist_min'] = df.groupby('molecule_name')['dist'].transform('min')\n df['molecule_dist_max'] = df.groupby('molecule_name')['dist'].transform('max')\n df['atom_0_couples_count'] = df.groupby(['molecule_name', 'atom_index_0'])['id'].transform('count')\n df['atom_1_couples_count'] = df.groupby(['molecule_name', 'atom_index_1'])['id'].transform('count')\n \n df[f'molecule_atom_index_0_x_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['x_1'].transform('std')\n df[f'molecule_atom_index_0_y_1_mean'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('mean')\n df[f'molecule_atom_index_0_y_1_mean_diff'] = df[f'molecule_atom_index_0_y_1_mean'] - df['y_1']\n df[f'molecule_atom_index_0_y_1_mean_div'] = df[f'molecule_atom_index_0_y_1_mean'] / df['y_1']\n df[f'molecule_atom_index_0_y_1_max'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('max')\n df[f'molecule_atom_index_0_y_1_max_diff'] = df[f'molecule_atom_index_0_y_1_max'] - df['y_1']\n df[f'molecule_atom_index_0_y_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('std')\n df[f'molecule_atom_index_0_z_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['z_1'].transform('std')\n df[f'molecule_atom_index_0_dist_mean'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('mean')\n df[f'molecule_atom_index_0_dist_mean_diff'] = df[f'molecule_atom_index_0_dist_mean'] - df['dist']\n df[f'molecule_atom_index_0_dist_mean_div'] = df[f'molecule_atom_index_0_dist_mean'] / df['dist']\n df[f'molecule_atom_index_0_dist_max'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('max')\n df[f'molecule_atom_index_0_dist_max_diff'] = df[f'molecule_atom_index_0_dist_max'] - df['dist']\n df[f'molecule_atom_index_0_dist_max_div'] = df[f'molecule_atom_index_0_dist_max'] / df['dist']\n df[f'molecule_atom_index_0_dist_min'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('min')\n df[f'molecule_atom_index_0_dist_min_diff'] = df[f'molecule_atom_index_0_dist_min'] - df['dist']\n df[f'molecule_atom_index_0_dist_min_div'] = df[f'molecule_atom_index_0_dist_min'] / df['dist']\n df[f'molecule_atom_index_0_dist_std'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('std')\n df[f'molecule_atom_index_0_dist_std_diff'] = df[f'molecule_atom_index_0_dist_std'] - df['dist']\n df[f'molecule_atom_index_0_dist_std_div'] = df[f'molecule_atom_index_0_dist_std'] / df['dist']\n df[f'molecule_atom_index_1_dist_mean'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('mean')\n df[f'molecule_atom_index_1_dist_mean_diff'] = df[f'molecule_atom_index_1_dist_mean'] - df['dist']\n df[f'molecule_atom_index_1_dist_mean_div'] = df[f'molecule_atom_index_1_dist_mean'] / df['dist']\n df[f'molecule_atom_index_1_dist_max'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('max')\n df[f'molecule_atom_index_1_dist_max_diff'] = df[f'molecule_atom_index_1_dist_max'] - df['dist']\n df[f'molecule_atom_index_1_dist_max_div'] = df[f'molecule_atom_index_1_dist_max'] / df['dist']\n df[f'molecule_atom_index_1_dist_min'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('min')\n df[f'molecule_atom_index_1_dist_min_diff'] = df[f'molecule_atom_index_1_dist_min'] - df['dist']\n df[f'molecule_atom_index_1_dist_min_div'] = df[f'molecule_atom_index_1_dist_min'] / df['dist']\n df[f'molecule_atom_index_1_dist_std'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('std')\n df[f'molecule_atom_index_1_dist_std_diff'] = df[f'molecule_atom_index_1_dist_std'] - df['dist']\n df[f'molecule_atom_index_1_dist_std_div'] = df[f'molecule_atom_index_1_dist_std'] / df['dist']\n df[f'molecule_atom_1_dist_mean'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('mean')\n df[f'molecule_atom_1_dist_min'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('min')\n df[f'molecule_atom_1_dist_min_diff'] = df[f'molecule_atom_1_dist_min'] - df['dist']\n df[f'molecule_atom_1_dist_min_div'] = df[f'molecule_atom_1_dist_min'] / df['dist']\n df[f'molecule_atom_1_dist_std'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('std')\n df[f'molecule_atom_1_dist_std_diff'] = df[f'molecule_atom_1_dist_std'] - df['dist']\n df[f'molecule_type_0_dist_std'] = df.groupby(['molecule_name', 'type_0'])['dist'].transform('std')\n df[f'molecule_type_0_dist_std_diff'] = df[f'molecule_type_0_dist_std'] - df['dist']\n df[f'molecule_type_dist_mean'] = df.groupby(['molecule_name', 'type'])['dist'].transform('mean')\n df[f'molecule_type_dist_mean_diff'] = df[f'molecule_type_dist_mean'] - df['dist']\n df[f'molecule_type_dist_mean_div'] = df[f'molecule_type_dist_mean'] / df['dist']\n df[f'molecule_type_dist_max'] = df.groupby(['molecule_name', 'type'])['dist'].transform('max')\n df[f'molecule_type_dist_min'] = df.groupby(['molecule_name', 'type'])['dist'].transform('min')\n df[f'molecule_type_dist_std'] = df.groupby(['molecule_name', 'type'])['dist'].transform('std')\n df[f'molecule_type_dist_std_diff'] = df[f'molecule_type_dist_std'] - df['dist']\n\n df = reduce_mem_usage(df)\n return df\n\ntrain = create_features(train)\ntest = create_features(test)\n\ngood_columns = [\n'molecule_atom_index_0_dist_min',\n'molecule_atom_index_0_dist_max',\n'molecule_atom_index_1_dist_min',\n'molecule_atom_index_0_dist_mean',\n'molecule_atom_index_0_dist_std',\n'dist',\n'molecule_atom_index_1_dist_std',\n'molecule_atom_index_1_dist_max',\n'molecule_atom_index_1_dist_mean',\n'molecule_atom_index_0_dist_max_diff',\n'molecule_atom_index_0_dist_max_div',\n'molecule_atom_index_0_dist_std_diff',\n'molecule_atom_index_0_dist_std_div',\n'atom_0_couples_count',\n'molecule_atom_index_0_dist_min_div',\n'molecule_atom_index_1_dist_std_diff',\n'molecule_atom_index_0_dist_mean_div',\n'atom_1_couples_count',\n'molecule_atom_index_0_dist_mean_diff',\n'molecule_couples',\n'atom_index_1',\n'molecule_dist_mean',\n'molecule_atom_index_1_dist_max_diff',\n'molecule_atom_index_0_y_1_std',\n'molecule_atom_index_1_dist_mean_diff',\n'molecule_atom_index_1_dist_std_div',\n'molecule_atom_index_1_dist_mean_div',\n'molecule_atom_index_1_dist_min_diff',\n'molecule_atom_index_1_dist_min_div',\n'molecule_atom_index_1_dist_max_div',\n'molecule_atom_index_0_z_1_std',\n'y_0',\n'molecule_type_dist_std_diff',\n'molecule_atom_1_dist_min_diff',\n'molecule_atom_index_0_x_1_std',\n'molecule_dist_min',\n'molecule_atom_index_0_dist_min_diff',\n'molecule_atom_index_0_y_1_mean_diff',\n'molecule_type_dist_min',\n'molecule_atom_1_dist_min_div',\n'atom_index_0',\n'molecule_dist_max',\n'molecule_atom_1_dist_std_diff',\n'molecule_type_dist_max',\n'molecule_atom_index_0_y_1_max_diff',\n'molecule_type_0_dist_std_diff',\n'molecule_type_dist_mean_diff',\n'molecule_atom_1_dist_mean',\n'molecule_atom_index_0_y_1_mean_div',\n'molecule_type_dist_mean_div']\n\nfor f in ['atom_1', 'type_0']:\n if f in good_columns:\n lbl = LabelEncoder()\n lbl.fit(list(train[f].values) + list(test[f].values))\n train[f] = lbl.transform(list(train[f].values))\n test[f] = lbl.transform(list(test[f].values))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train[good_columns].copy()\ny = train['scalar_coupling_constant']\nX_test = test[good_columns].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_1JHN = train.loc[train['type'] == '1JHN', good_columns]\ny_1JHN = train.loc[train['type'] == '1JHN', 'scalar_coupling_constant']\nX_test_1JHN = test.loc[test['type'] == '1JHN', good_columns]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Adversarial validation\n\nAdversarial validation is a commonly used approach for selecting stable features. The idea is to combine train and test data and build a binary classifier. If there are features which work well to separate train and test data, they will likely be useless or harmful.\n\nI limit n_estimators and some other parameters to make training faster."},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'num_leaves': 128,\n 'min_data_in_leaf': 79,\n 'objective': 'binary',\n 'max_depth': 5,\n 'learning_rate': 0.3,\n \"boosting\": \"gbdt\",\n \"bagging_freq\": 1,\n \"bagging_fraction\": 0.8,\n \"bagging_seed\": 11,\n \"metric\": 'auc',\n \"verbosity\": -1,\n 'reg_alpha': 0.1,\n 'reg_lambda': 0.3,\n 'feature_fraction': 1.0\n }\n\nn_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)\n\nfeature_importance = {}\nfor t in ['1JHN', '2JHN', '3JHN']:\n print(f'Type: {t}. {time.ctime()}')\n X = train.loc[train['type'] == t, good_columns]\n y = train.loc[train['type'] == t, 'scalar_coupling_constant']\n X_test = test.loc[test['type'] == t, good_columns]\n \n features = X.columns\n X['target'] = 0\n X_test['target'] = 1\n\n train_test = pd.concat([X, X_test], axis=0)\n target = train_test['target']\n \n result_dict_lgb = train_model_classification(train_test, X_test, target, params, folds, model_type='lgb',\n columns=features, plot_feature_importance=True, model=None, verbose=500,\n early_stopping_rounds=200, n_estimators=500)\n plt.show();\n \n feature_importance[t] = result_dict_lgb['feature_importance']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that auc is ~0.7 (even though we didn't train model for a long time), so it would be worth dropping top features."},{"metadata":{"trusted":true},"cell_type":"code","source":"bad_advers_columns = {}\nfor t in feature_importance.keys():\n cols = feature_importance[t][[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n by=\"importance\", ascending=False).index[:7]\n bad_advers_columns[t] = cols","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Validation\n\nCreating a stable validation is one of serious challenges in each competition. Let's compare several cross-validation schemes and see how they work. I'll try the following schemes:\n\n- simple KFold with shuffling and without it;\n- strafitied KFold by by atom_index_0;\n- group KFold by molecules names and by atom_index_0;\n\nFor each scheme I'll train models and then compare their scores on cross-validation to see which one is more stable. Another way is to make a holdout set and compare model quality on cross-validation and on this holdout set. Or even better - compare it to the leaderboard."},{"metadata":{"trusted":true},"cell_type":"code","source":"val_folds = {'kfold_shuffle': KFold(n_splits=n_fold, shuffle=True, random_state=11),\n 'kfold_no_shuffle': KFold(n_splits=n_fold, shuffle=False, random_state=11),\n 'stratified_atom_index_0': StratifiedKFold(n_splits=n_fold, shuffle=False, random_state=11),\n 'group_molecules': GroupKFold(n_splits=n_fold),\n 'group_atom_index_0': GroupKFold(n_splits=n_fold)}\n\nparams = {'num_leaves': 128,\n 'min_child_samples': 79,\n 'objective': 'regression',\n 'max_depth': 9,\n 'learning_rate': 0.2,\n \"boosting_type\": \"gbdt\",\n \"subsample_freq\": 1,\n \"subsample\": 0.9,\n \"bagging_seed\": 11,\n \"metric\": 'mae',\n \"verbosity\": -1,\n 'reg_alpha': 0.1,\n 'reg_lambda': 0.3,\n 'colsample_bytree': 1.0\n }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"val_scores = {}\nfor t in ['1JHN', '2JHN', '3JHN']:\n val_scores[t] = {}\n print(f'Type: {t}. {time.ctime()}')\n \n X = train.loc[train['type'] == t, [col for col in good_columns if col not in bad_advers_columns[t]]]\n y = train.loc[train['type'] == t, 'scalar_coupling_constant']\n X_test = test.loc[test['type'] == t, [col for col in good_columns if col not in bad_advers_columns[t]]]\n\n for k, v in val_folds.items():\n print(f'Validation type: {k}')\n if 'kfold' in k:\n splits = v.split(X)\n elif 'stratified' in k:\n if 'molecules' in k:\n splits = v.split(X, train.loc[train['type'] == t, 'molecule_name'])\n else:\n splits = v.split(X, train.loc[train['type'] == t, 'atom_index_0'])\n elif 'group' in k:\n if 'molecules' in k:\n splits = v.split(X, y, train.loc[train['type'] == t, 'molecule_name'])\n else:\n splits = v.split(X, y, train.loc[train['type'] == t, 'atom_index_0'])\n\n result_dict_lgb = train_model_regression(X=X, X_test=X_test, y=y, params=params, model_type='lgb', eval_metric='mae', plot_feature_importance=False,\n verbose=False, early_stopping_rounds=100, n_estimators=500, splits=splits, n_folds=n_fold)\n\n val_scores[t][k] = result_dict_lgb['scores']\n print()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (18, 6))\nfor i, k in enumerate(val_scores.keys()):\n val_scores_df = pd.DataFrame(val_scores[k]).T.unstack().reset_index().rename(columns={'level_1': 'validation_scheme', 0: 'scores'})\n plt.subplot(1, 3, i + 1);\n sns.boxplot(data=val_scores_df, x='validation_scheme', y='scores');\n plt.xticks(rotation=45);\n plt.title(k);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is quite interesting!\n\n- stratified splits have a wide range of values as well as group kfold on `atom_index_0`. This feature has some correlation with the target, but it seems that splitting by it isn't really useful;\n- strangely kfold without shuffling works quite bad. I suppose this is because we need to have diverse information about different molecules in each fold;\n- kfold with shuffling and group kfold by molecules work well!\n\nSo we can should try group kfold by molecules. But for now let's still use simple kfold"},{"metadata":{},"cell_type":"markdown","source":"## Comparing several models"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_1JHN = train.loc[train['type'] == '1JHN', [col for col in good_columns if col not in bad_advers_columns['1JHN']]].fillna(-999)\ny_1JHN = train.loc[train['type'] == '1JHN', 'scalar_coupling_constant']\nX_test_1JHN = test.loc[test['type'] == '1JHN', [col for col in good_columns if col not in bad_advers_columns['1JHN']]].fillna(-999)\n\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.ensemble import ExtraTreesRegressor, AdaBoostRegressor, RandomForestRegressor\nfrom sklearn.kernel_ridge import KernelRidge\nfrom sklearn.neighbors import NearestNeighbors","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\netr = ExtraTreesRegressor()\n\nparameter_grid = {'n_estimators': [100, 300],\n 'max_depth': [15, 50]\n }\n\ngrid_search = GridSearchCV(etr, param_grid=parameter_grid, cv=folds, scoring='neg_mean_absolute_error', n_jobs=-1)\ngrid_search.fit(X_1JHN, y_1JHN, groups=train.loc[train['type'] == '1JHN', 'molecule_name'])\nprint('Best score: {}'.format(grid_search.best_score_))\nprint('Best parameters: {}'.format(grid_search.best_params_))\netr = ExtraTreesRegressor(**grid_search.best_params_)\nresult_dict_etr = train_model_regression(X=X_1JHN, X_test=X_test_1JHN, y=y_1JHN, params=params, model_type='sklearn', model=etr, eval_metric='mae', plot_feature_importance=False,\n verbose=False, early_stopping_rounds=100, n_estimators=500, folds=folds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nada = AdaBoostRegressor()\n\nparameter_grid = {'n_estimators': [50, 200]}\n\ngrid_search = GridSearchCV(ada, param_grid=parameter_grid, cv=folds, scoring='neg_mean_absolute_error', n_jobs=-1)\ngrid_search.fit(X_1JHN, y_1JHN, groups=train.loc[train['type'] == '1JHN', 'molecule_name'])\nprint('Best score: {}'.format(grid_search.best_score_))\nprint('Best parameters: {}'.format(grid_search.best_params_))\nada = AdaBoostRegressor(**grid_search.best_params_)\nresult_dict_ada = train_model_regression(X=X_1JHN, X_test=X_test_1JHN, y=y_1JHN, params=params, model_type='sklearn', model=ada, eval_metric='mae', plot_feature_importance=False,\n verbose=False, early_stopping_rounds=100, n_estimators=500, folds=folds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nrfr = RandomForestRegressor()\n\nparameter_grid = {'n_estimators': [100, 300],\n 'max_depth': [15, 50]\n }\n\ngrid_search = GridSearchCV(rfr, param_grid=parameter_grid, cv=folds, scoring='neg_mean_absolute_error', n_jobs=-1)\ngrid_search.fit(X_1JHN, y_1JHN, groups=train.loc[train['type'] == '1JHN', 'molecule_name'])\nprint('Best score: {}'.format(grid_search.best_score_))\nprint('Best parameters: {}'.format(grid_search.best_params_))\nrfr = RandomForestRegressor(**grid_search.best_params_)\nresult_dict_rfr = train_model_regression(X=X_1JHN, X_test=X_test_1JHN, y=y_1JHN, params=params, folds=folds, model_type='sklearn', model=rfr, eval_metric='mae', plot_feature_importance=False,\n verbose=10, early_stopping_rounds=100, n_estimators=500)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_1JHN = train.loc[train['type'] == '1JHN', [col for col in good_columns if col not in bad_advers_columns['1JHN']]].fillna(0)\ny_1JHN = train.loc[train['type'] == '1JHN', 'scalar_coupling_constant']\nX_test_1JHN = test.loc[test['type'] == '1JHN', [col for col in good_columns if col not in bad_advers_columns['1JHN']]].fillna(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nridge = linear_model.Ridge(normalize=True)\n\nparameter_grid = {'alpha': [0.01, 0.1, 1.0]}\n\ngrid_search = GridSearchCV(ridge, param_grid=parameter_grid, cv=folds, scoring='neg_mean_absolute_error', n_jobs=-1)\ngrid_search.fit(X_1JHN, y_1JHN, groups=train.loc[train['type'] == '1JHN', 'molecule_name'])\nprint('Best score: {}'.format(grid_search.best_score_))\nprint('Best parameters: {}'.format(grid_search.best_params_))\nridge = linear_model.Ridge(**grid_search.best_params_, normalize=True)\nresult_dict_ridge = train_model_regression(X=X_1JHN, X_test=X_test_1JHN, y=y_1JHN, params=params, folds=folds, model_type='sklearn', model=ridge, eval_metric='mae', plot_feature_importance=False,\n verbose=10, early_stopping_rounds=100, n_estimators=500)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nlasso = linear_model.Lasso(normalize=True)\n\nparameter_grid = {'alpha': [0.01, 0.1, 1.0]}\n\ngrid_search = GridSearchCV(lasso, param_grid=parameter_grid, cv=folds, scoring='neg_mean_absolute_error', n_jobs=-1)\ngrid_search.fit(X_1JHN, y_1JHN, groups=train.loc[train['type'] == '1JHN', 'molecule_name'])\nprint('Best score: {}'.format(grid_search.best_score_))\nprint('Best parameters: {}'.format(grid_search.best_params_))\nlasso = linear_model.Lasso(**grid_search.best_params_, normalize=True)\nresult_dict_lasso = train_model_regression(X=X_1JHN, X_test=X_test_1JHN, y=y_1JHN, params=params, model_type='sklearn', folds=folds, model=lasso, eval_metric='mae', plot_feature_importance=False,\n verbose=False, early_stopping_rounds=100, n_estimators=500)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(12, 8));\nscores_df = pd.DataFrame({'RandomForestRegressor': result_dict_rfr['scores']})\nscores_df['ExtraTreesRegressor'] = result_dict_etr['scores']\nscores_df['AdaBoostRegressor'] = result_dict_ada['scores']\n# scores_df['KNN'] = result_dict_knn['scores']\nscores_df['Ridge'] = result_dict_ridge['scores']\nscores_df['Lasso'] = result_dict_lasso['scores']\n\nsns.boxplot(data=scores_df);\nplt.xticks(rotation=45);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Main model training"},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = [set(i) for i in bad_advers_columns.values()]\nbad_cols = list(cols[0].intersection(cols[1]).intersection(cols[2]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_short = pd.DataFrame({'ind': list(train.index), 'type': train['type'].values, 'oof': [0] * len(train), 'target': train['scalar_coupling_constant'].values})\nX_short_test = pd.DataFrame({'ind': list(test.index), 'type': test['type'].values, 'prediction': [0] * len(test)})\nfor t in train['type'].unique():\n print(f'Training of type {t}')\n X_t = train.loc[train['type'] == t, [col for col in good_columns if col not in bad_cols]]\n X_test_t = test.loc[test['type'] == t, [col for col in good_columns if col not in bad_cols]]\n y_t = X_short.loc[X_short['type'] == t, 'target']\n splits = folds.split(X_t, y_t, train.loc[train['type'] == t, 'molecule_name'])\n result_dict_lgb = train_model_regression(X=X_t, X_test=X_test_t, y=y_t, params=params, folds=folds, model_type='lgb', eval_metric='mae', plot_feature_importance=True,\n verbose=1000, early_stopping_rounds=200, n_estimators=2000, splits=splits, n_folds=n_fold)\n plt.show();\n X_short.loc[X_short['type'] == t, 'oof'] = result_dict_lgb['oof']\n X_short_test.loc[X_short_test['type'] == t, 'prediction'] = result_dict_lgb['prediction']\n \n \nsub['scalar_coupling_constant'] = X_short_test['prediction']\nsub.to_csv('submission_t.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model interpretation"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_1JHN = train.loc[train['type'] == '1JHN', [col for col in good_columns if col not in bad_advers_columns['1JHN']]].fillna(-999)\ny_1JHN = train.loc[train['type'] == '1JHN', 'scalar_coupling_constant']\nX_test_1JHN = test.loc[test['type'] == '1JHN', [col for col in good_columns if col not in bad_advers_columns['1JHN']]].fillna(-999)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_valid, y_train, y_valid = train_test_split(X_1JHN, y_1JHN, test_size=0.1)\n\nparams = {'num_leaves': 30,\n 'min_data_in_leaf': 20,\n 'objective': 'regression',\n 'max_depth': 6,\n 'learning_rate': 0.01,\n \"boosting\": \"gbdt\",\n \"feature_fraction\": 0.9,\n \"bagging_freq\": 1,\n \"bagging_fraction\": 0.9,\n \"bagging_seed\": 11,\n \"metric\": 'mse',\n \"lambda_l1\": 0.2,\n \"verbosity\": -1}\nmodel1 = lgb.LGBMRegressor(**params, n_estimators = 20000, nthread = 4, n_jobs = -1)\nmodel1.fit(X_train, y_train, \n eval_set=[(X_train, y_train), (X_valid, y_valid)], eval_metric='rmse',\n verbose=1000, early_stopping_rounds=200)\n\neli5.show_weights(model1, feature_filter=lambda x: x != '')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"explainer = shap.TreeExplainer(model1, X_train)\nshap_values = explainer.shap_values(X_train)\n\nshap.summary_plot(shap_values, X_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"top_cols = X_train.columns[np.argsort(shap_values.std(0))[::-1]][:5]\nfor col in top_cols:\n shap.dependence_plot(col, shap_values, X_train)","execution_count":null,"outputs":[]}],"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}