{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# EDA which makes sense ⭐️⭐️⭐️⭐️⭐️\n\nThis notebook shows what kind of insight an EDA should give and how this insight is applied in creating a model and for correctly cross-validating it.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom matplotlib import pyplot as plt\nimport seaborn as sns\nimport warnings\nfrom colorama import Fore, Back, Style\nimport scipy.stats\n\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.preprocessing import OneHotEncoder, OrdinalEncoder, StandardScaler\nfrom sklearn.experimental import enable_iterative_imputer\nfrom sklearn.impute import SimpleImputer, IterativeImputer, KNNImputer\nfrom sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.linear_model import LogisticRegression, LinearRegression\nfrom sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis, LinearDiscriminantAnalysis\nfrom sklearn.metrics import roc_auc_score, roc_curve\n\nfrom lightgbm import LGBMClassifier, early_stopping, log_evaluation\n\nnp.set_printoptions(linewidth=150)\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:29:36.629423Z","iopub.execute_input":"2022-08-06T18:29:36.630118Z","iopub.status.idle":"2022-08-06T18:29:37.787665Z","shell.execute_reply.started":"2022-08-06T18:29:36.629949Z","shell.execute_reply":"2022-08-06T18:29:37.786158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading the data\n\nThe train file contains 26570 rows and the test file 20775 rows. \n\n**Insight:** The public leaderboard is based on a quarter of the test file, i.e. only 5000 rows. We have to make sure that we don't overfit to the public leaderboard. We should trust our cross-validation (which is based on 26570 rows) much more than the public leaderboard.","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv('../input/tabular-playground-series-aug-2022/train.csv',\n                    index_col='id')\ntest = pd.read_csv('../input/tabular-playground-series-aug-2022/test.csv',\n                    index_col='id')\ndisplay(train)\ndisplay(test)\nboth = pd.concat([train[test.columns], test])\n","metadata":{"execution":{"iopub.status.busy":"2022-08-06T18:29:37.790235Z","iopub.execute_input":"2022-08-06T18:29:37.790726Z","iopub.status.idle":"2022-08-06T18:29:38.047434Z","shell.execute_reply.started":"2022-08-06T18:29:37.790680Z","shell.execute_reply":"2022-08-06T18:29:38.046495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The target column\n\nOf the 26570 products tested, 79 % are good and 21 % fail.\n\n**Insight:** This looks like a slightly imbalanced classification. It might be good to stratify the train-test splits.","metadata":{}},{"cell_type":"code","source":"print(train.failure.value_counts() / len(train))","metadata":{"execution":{"iopub.status.busy":"2022-08-06T18:29:38.048851Z","iopub.execute_input":"2022-08-06T18:29:38.049117Z","iopub.status.idle":"2022-08-06T18:29:38.055715Z","shell.execute_reply.started":"2022-08-06T18:29:38.049091Z","shell.execute_reply":"2022-08-06T18:29:38.054873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The float columns\n\nThe data have 16 float columns. All columns have missing values (up to 10 % of the column).\n\n**Insight:** We'll need to impute the missing values (unless we use a classifier which deals with missing values automatically). The simplest solution is filling the missing values with the column's mean, but this simple solution won't win the competition. A more sophisticated solution might use the imputers from [sklearn.impute](https://scikit-learn.org/stable/modules/classes.html#module-sklearn.impute) or even some customized imputation scheme.\n\n","metadata":{}},{"cell_type":"code","source":"float_cols = [f for f in train.columns if train[f].dtype == float]\npd.concat([train[float_cols].isna().sum().rename('missing values in train'),\n           test[float_cols].isna().sum().rename('missing values in test')],\n          axis=1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:29:38.057626Z","iopub.execute_input":"2022-08-06T18:29:38.057847Z","iopub.status.idle":"2022-08-06T18:29:38.077466Z","shell.execute_reply.started":"2022-08-06T18:29:38.057824Z","shell.execute_reply":"2022-08-06T18:29:38.076545Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"More than half the rows have at least one missing value:","metadata":{}},{"cell_type":"code","source":"print(f\"{both[float_cols].isna().any(axis=1).sum() / len(both):.0%}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-06T18:29:38.078477Z","iopub.execute_input":"2022-08-06T18:29:38.078781Z","iopub.status.idle":"2022-08-06T18:29:38.094550Z","shell.execute_reply.started":"2022-08-06T18:29:38.078756Z","shell.execute_reply":"2022-08-06T18:29:38.093752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Insight:** With that many missing values, good value imputation will be the key for winning the competition.","metadata":{}},{"cell_type":"markdown","source":"Let's look at the distribution of the float features. We plot the train and test histograms in the same diagram - train is blue, test is orange.\n\nWe see that the first feature, loading, has a skewed distribution, perhaps log-normal; the other features are normally distributed. The first eight features have the same distribution in train and test; from measurement_10 onwards the distributions differ slightly.\n\nIn the same diagrams, we show the failure probabilities with magenta dots. The top left diagram shows that higher loadings imply a higher failure probability. The bottom right diagram shows that measurement_17 is positively correlated to the target as well. All other float features seem to be uncorrelated to the failure probability (the failure probability is 21 % independent of the feature value).\n\n**Insight:**\n- Maybe we can apply a log-transformation to the loading to make the distribution more symmetric.\n- As the measurements are uncorrelated to the failure probability, they are useless for a linear classifier. We need more complex classifiers which can deal with feature interactions (decision trees, neural networks, ...) - or good feature engineering.","metadata":{}},{"cell_type":"code","source":"_, axs = plt.subplots(4, 4, figsize=(12,12))\nfor f, ax in zip(float_cols, axs.ravel()):\n    mi = min(train[f].min(), test[f].min())\n    ma = max(train[f].max(), test[f].max())\n    bins = np.linspace(mi, ma, 50)\n    ax.hist(train[f], bins=bins, alpha=0.5, density=True, label='train')\n    ax.hist(test[f], bins=bins, alpha=0.5, density=True, label='test')\n    ax.set_xlabel(f)\n    if ax == axs[0, 0]: ax.legend(loc='lower left')\n        \n    ax2 = ax.twinx()\n    total, _ = np.histogram(train[f], bins=bins)\n    failures, _ = np.histogram(train[f][train.failure == 1], bins=bins)\n    with warnings.catch_warnings(): # ignore divide by zero for empty bins\n        warnings.filterwarnings('ignore', category=RuntimeWarning)\n        ax2.scatter((bins[1:] + bins[:-1]) / 2, failures / total,\n                    color='m', s=10, label='failure probability')\n    ax2.set_ylim(0, 0.5)\n    ax2.tick_params(axis='y', colors='m')\n    if ax == axs[0, 0]: ax2.legend(loc='upper right')\nplt.tight_layout(w_pad=1)\nplt.suptitle('Train and test distributions of the continuous features', fontsize=20, y=1.02)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:29:38.095542Z","iopub.execute_input":"2022-08-06T18:29:38.095776Z","iopub.status.idle":"2022-08-06T18:29:44.167662Z","shell.execute_reply.started":"2022-08-06T18:29:38.095752Z","shell.execute_reply":"2022-08-06T18:29:44.166919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Did it occur to you that a missing measurement might be an early indicator of product failure? Maybe the cause of a product failure triggers the failure of a measurement device. How can we test this idea? Obviously we have to calculate the conditional product failure rate given the measurement is missing E\\[product fails | measurement is missing\\] and compare it to the unconditional product failure rate, which is 0.212608.\n\nOf course, the deviations will be small and we should test them for significance against a null hypothesis which says that the failure count is binomially distributed with p = 0.212608. For simplicity, we can approximate the binomial distribution by a normal distribution, calculate the [z-score](https://en.wikipedia.org/wiki/Standard_score) and the [p-value](https://en.wikipedia.org/wiki/P-value):\n","metadata":{}},{"cell_type":"code","source":"# Start by plotting the bell curve\nplt.figure(figsize=(12, 4))\nz_ticks = np.linspace(-3.5, 3.5, 61)\npdf = scipy.stats.norm.pdf(z_ticks)\nplt.plot(z_ticks, pdf)\n\n# Calculate the conditional failure rate for every missing feature\n# Print the values and plot them\nprint('feature           fail   miss   failure rate       z    p-value')\nfor f in train.columns:\n    if train[f].isna().sum() > 0:\n        total = train[f].isna().sum()\n        fail = train[train[f].isna()].failure.sum()\n        z = (fail / total - 0.212608) / (np.sqrt(0.212608 * (1-0.212608)) / np.sqrt(total))\n        plt.scatter([z], [scipy.stats.norm.pdf(z)], c='r' if abs(z) > 2 else 'g', s=100)\n        print(f\"{f:15} : {fail:4} / {total:4} = {fail/total:.3f}          {z:5.2f}      {2*scipy.stats.norm.cdf(-abs(z)):.3f}\")\n        if abs(z) > 1: plt.annotate(f\"{f}: {fail / total:.3f}\",\n                                    (z, scipy.stats.norm.pdf(z)),\n                                    xytext=(0,10), \n                                    textcoords='offset points', ha='left' if z > 0 else 'right',\n                                    color='r' if abs(z) > 2 else 'g')\n            \n# Annotage the center (z=0)\nplt.vlines([0], 0, 0.05, color='g')\nplt.annotate(f\"z_score = 0\\naverage failure rate: {0.212608:.3f}\",\n                                    (0, 0.05),\n                                    xytext=(0,10), \n                                    textcoords='offset points', ha='center',\n                                    color='g')\nplt.title('Failure rate when feature is missing')\nplt.yticks([])\nplt.xlabel('z_score')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-06T18:34:53.477705Z","iopub.execute_input":"2022-08-06T18:34:53.478083Z","iopub.status.idle":"2022-08-06T18:34:53.931674Z","shell.execute_reply.started":"2022-08-06T18:34:53.478038Z","shell.execute_reply":"2022-08-06T18:34:53.930791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- When measurement_3 is missing, the failure rate is 0.160 (much lower than average).\n- When measurement_5 is missing, the failure rate is 0.254 (much higher than average).\n\nWith abs(z) > 2.5 and pvalue < 2 %, the conditional failure rates of missing measurement_3 and missing measurement_5 deviate significantly from the average failure rate and we can use the features m_3_missing and m_5_missing in our models:\n\n```\nX['m_3_missing'] = X.measurement_3.isna()\nX['m_5_missing'] = X.measurement_5.isna()\n```","metadata":{}},{"cell_type":"markdown","source":"# The integer columns\n\nThere are five integer columns. We are lucky: the columns are complete, without missing values. Two of the columns are called *attribute* and three *measurement*.","metadata":{}},{"cell_type":"code","source":"int_cols = [f for f in train.columns if train[f].dtype == int and f != 'failure']\npd.concat([train[int_cols].isna().sum().rename('missing values in train'),\n           test[int_cols].isna().sum().rename('missing values in test')],\n          axis=1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:29:44.479707Z","iopub.execute_input":"2022-08-06T18:29:44.480023Z","iopub.status.idle":"2022-08-06T18:29:44.495029Z","shell.execute_reply.started":"2022-08-06T18:29:44.479994Z","shell.execute_reply":"2022-08-06T18:29:44.494235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's plot the distributions of the five integer features for train and test:","metadata":{}},{"cell_type":"code","source":"_, axs = plt.subplots(2, 3, figsize=(12, 8))\nfor f, ax in zip(int_cols, axs.ravel()):\n    temp1 = train.failure.groupby(train[f]).agg(['mean', 'size'])\n    ax.bar(temp1.index, temp1['size'] / len(train), alpha=0.5, label='train')\n    temp2 = test[f].value_counts()\n    ax.bar(temp2.index, temp2 / len(test), alpha=0.5, label='test')\n    ax.set_xlabel(f)\n    ax.set_ylabel('frequency')\n\n    ax2 = ax.twinx()\n    ax2.scatter(temp1.index, temp1['mean'],\n                color='m', label='failure probability')\n    ax2.set_ylim(0, 0.5)\n    ax2.tick_params(axis='y', colors='m')\n    if ax == axs[0, 0]: ax2.legend(loc='upper right')\n\naxs[0, 0].legend()\naxs[1, 2].axis('off')\nplt.tight_layout(w_pad=1)\nplt.suptitle('Train and test distributions of the integer features', fontsize=20, y=1.02)\nplt.show()\ndel temp1, temp2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:29:44.496141Z","iopub.execute_input":"2022-08-06T18:29:44.496971Z","iopub.status.idle":"2022-08-06T18:29:46.107856Z","shell.execute_reply.started":"2022-08-06T18:29:44.496941Z","shell.execute_reply":"2022-08-06T18:29:46.106891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We see that attribute_2 has two values (5 and 8) which occur only in the training data, and another value (7) occurs only in the test data. attribute_3 is similar.\nMeasurement_2 has positive correlation with the target, at least for values > 10.\n\n**Insight**:\n- The two attributes could be categorical features. Perhaps we should one-hot encode them and make sure that our classifier can deal with the values which occur only in test.\n- As measurement_2 is correlated to the target only for values above 10, linear classifiers will profit if we clip all values below 10.\n- Again, most features are uncorrelated to the failure probability. If a linear classifier wins this competition, it will win only with very good feature engineering.\n","metadata":{}},{"cell_type":"markdown","source":"# The string columns\n\nThree columns contain strings. They have no missing values.","metadata":{}},{"cell_type":"code","source":"string_cols = [f for f in train.columns if train[f].dtype == object]\npd.concat([train[string_cols].isna().sum().rename('missing values in train'),\n           test[string_cols].isna().sum().rename('missing values in test')],\n          axis=1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:29:46.110878Z","iopub.execute_input":"2022-08-06T18:29:46.111235Z","iopub.status.idle":"2022-08-06T18:29:46.134897Z","shell.execute_reply.started":"2022-08-06T18:29:46.111204Z","shell.execute_reply":"2022-08-06T18:29:46.133664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_, axs = plt.subplots(1, 3, figsize=(12, 4))\nfor f, ax in zip(string_cols, axs.ravel()):\n    temp1 = train[f].value_counts(dropna=False, normalize=True)\n    temp2 = test[f].value_counts(dropna=False, normalize=True)\n    values = sorted(set(temp1.index).union(temp2.index))\n    temp1 = temp1.reindex(values)\n    temp2 = temp2.reindex(values)\n    ax.bar(range(len(values)), temp1, alpha=0.5, label='train')\n    ax.bar(range(len(values)), temp2, alpha=0.5, label='test')\n    ax.set_xlabel(f)\n    ax.set_ylabel('frequency')\n    ax.set_xticks(range(len(values)), values)\n    \n    temp1 = train.failure.groupby(train[f]).agg(['mean', 'size'])\n    temp1 = temp1.reindex(values)\n    ax2 = ax.twinx()\n    ax2.scatter(range(len(values)), temp1['mean'],\n                color='m', label='failure probability')\n    ax2.tick_params(axis='y', colors='m')\n    ax2.set_ylim(0, 0.5)\n    if ax == axs[0]: ax2.legend(loc='lower right')\n\naxs[0].legend()\nplt.suptitle('Train and test distributions of the string features', fontsize=20, y=0.96)\nplt.tight_layout(w_pad=1)\nplt.show()\ndel temp1, temp2   \n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:29:46.136735Z","iopub.execute_input":"2022-08-06T18:29:46.137243Z","iopub.status.idle":"2022-08-06T18:29:46.813790Z","shell.execute_reply.started":"2022-08-06T18:29:46.137206Z","shell.execute_reply":"2022-08-06T18:29:46.812724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Insight:**\n- The product codes of train and test are disjoint: A through E are the training products, F through I are the test products. We want to create a classifier which predicts correct probabilities for previously unseen products. To validate such a classifier, we have to simulate this situation by splitting the data so that the validation set contains other products than the training set. The correct method here is a five-fold cross-validation where every fold uses four products for training and the fifth product for validation ([GroupKFold](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GroupKFold.html)). Forget what I wrote about imbalanced classes and stratified folds above - getting the product codes right is much more important!\n- Feature engineering: We can use the product codes for feature engineering by adding aggregate statistics of the measurements, grouped by product code, as new features.\n- Attribute_0 and attribute_1 are categorical features which should be one-hot encoded.","metadata":{}},{"cell_type":"markdown","source":"# Product codes and attributes\n\n@takanashihumbert has [discovered a dependency among the attributes](https://www.kaggle.com/code/takanashihumbert/interesting-patterns-found-in-category-features): All four attributes are completely determined by the product code. If we wanted to store the data in a relational database, we'd normalize it by storing the product attributes in a separate table.","metadata":{}},{"cell_type":"code","source":"both[string_cols + ['attribute_2', 'attribute_3']].drop_duplicates().set_index('product_code')","metadata":{"execution":{"iopub.status.busy":"2022-08-06T18:29:46.815622Z","iopub.execute_input":"2022-08-06T18:29:46.816017Z","iopub.status.idle":"2022-08-06T18:29:46.845830Z","shell.execute_reply.started":"2022-08-06T18:29:46.815985Z","shell.execute_reply":"2022-08-06T18:29:46.844873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Baseline model with cross-validation\n\nIn the following, we'll see how to correctly implement preprocessing and cross-validation for a simple model. You can play around with the hyperparameters - I haven't really optimized them yet.\n\nThe code contains a choice of three classifiers: RandomForestClassifier, ExtraTreesClassifier and LogisticRegression. The saved version of the notebook uses LogisticRegression and prints the weights which the classifier optimizes.\n\nWhen I did this experiment for the first time, I saw that most regression coefficients were near zero, which means that the corresponding features are noise. To verify this thought, I configured the classifier for a l1 penalty. The l1 penalty tries to select as few features as possible for the model.\n\nThe pipeline consists of the following steps:\n1. Split train and validation ([GroupKFold](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GroupKFold.html) on product_code)\n2. One-hot encode attribute_0 and attribute_1\n3. Impute the missing values ([KNNImputer](https://scikit-learn.org/stable/modules/generated/sklearn.impute.KNNImputer.html))\n4. Clip measurement_2\n5. Scale the data (StandardScaler)\n6. Logistic regression with l1 penalty\n7. Evaluate feature importances (regression coefficients)\n","metadata":{}},{"cell_type":"code","source":"auc_list = []\ntest_pred_list = []\nimportance_list = []\nkf = GroupKFold(n_splits=5) # must be 5 because of the 5 product codes\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(train, train.failure, train.product_code)):\n    X_tr = train.iloc[idx_tr][test.columns]\n    X_va = train.iloc[idx_va][test.columns]\n    X_te = test.copy()\n    y_tr = train.iloc[idx_tr].failure\n    y_va = train.iloc[idx_va].failure\n\n    # We one-hot encode attribute_0 and attribute_1\n    ohe_attributes = ['attribute_0', 'attribute_1']\n    ohe_output = ['ohe0_7', 'ohe1_6', 'ohe1_8']\n    ohe = OneHotEncoder(categories=[['material_5', 'material_7'],\n                                    ['material_5', 'material_6', 'material_8']],\n                        drop='first', sparse=False, handle_unknown='ignore')\n    ohe.fit(X_tr[ohe_attributes])\n    for df in [X_tr, X_va, X_te]:\n        with warnings.catch_warnings(): # ignore \"Found unknown categories\"\n            warnings.filterwarnings('ignore', category=UserWarning)\n            df[ohe_output] = ohe.transform(df[ohe_attributes])\n        df.drop(columns=ohe_attributes, inplace=True)\n\n    # We add the indicators for missing values\n    for df in [X_tr, X_va, X_te]:\n        df['m_3_missing'] = df.measurement_3.isna()\n        df['m_5_missing'] = df.measurement_5.isna()\n\n    # We fill the missing values\n    features = [f for f in X_tr.columns if f == 'loading' or f.startswith('measurement')]\n    imputer = KNNImputer(n_neighbors=3)\n    imputer.fit(X_tr[features])\n    for df in [X_tr, X_va, X_te]:\n        df[features] = imputer.transform(df[features])\n\n    # The EDA diagram of measurement 2 shows that the feature is correlated\n    # to the target only for values above 10. For this reason, we clip\n    # all values below 11.\n    for df in [X_tr, X_va, X_te]:\n        df['measurement_2'] = df['measurement_2'].clip(11, None)\n\n    # We fit a model\n    features2 = [f for f in X_tr.columns if f != 'product_code']\n    #model = RandomForestClassifier(n_estimators=200, max_depth=8, min_samples_leaf=100, n_jobs=-1, random_state=1)\n    #model = ExtraTreesClassifier(n_estimators=100, max_depth=8, min_samples_leaf=100, max_features=5, n_jobs=-1, random_state=1)\n    model = make_pipeline(StandardScaler(), \n                          LogisticRegression(penalty='l1', C=0.01,\n                                             solver='liblinear', random_state=1))\n    model.fit(X_tr[features2], y_tr)\n    importance_list.append(model.named_steps['logisticregression'].coef_.ravel())\n\n    # We validate the model\n    y_va_pred = model.predict_proba(X_va[features2])[:,1]\n    score = roc_auc_score(y_va, y_va_pred)\n    print(f\"Fold {fold}: auc = {score:.5f}\")\n    auc_list.append(score)\n\n    test_pred_list.append(model.predict_proba(X_te[features2])[:,1])\n\n# Show overall score\nprint(f\"{Fore.GREEN}{Style.BRIGHT}Average auc = {sum(auc_list) / len(auc_list):.5f}{Style.RESET_ALL}\")\n\n# Show feature importances\nimportance_df = pd.DataFrame(np.array(importance_list).T, index=features2)\nimportance_df['mean'] = importance_df.mean(axis=1).abs()\nimportance_df['feature'] = features2\nimportance_df = importance_df.sort_values('mean', ascending=False).reset_index().head(10)\nplt.figure(figsize=(14, 4))\nplt.barh(importance_df.index, importance_df['mean'], color='lightgreen')\nplt.gca().invert_yaxis()\nplt.yticks(ticks=importance_df.index, labels=importance_df['feature'])\nplt.title('LogisticRegression feature importances')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-06T18:29:46.847861Z","iopub.execute_input":"2022-08-06T18:29:46.848466Z","iopub.status.idle":"2022-08-06T18:32:01.485671Z","shell.execute_reply.started":"2022-08-06T18:29:46.848419Z","shell.execute_reply":"2022-08-06T18:32:01.484622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Manual feature selection\n\nIf you look at the logistic regression coefficients above (the green bar chart), you'll see that the first coefficient is 0.27 or 0.28 and all other coefficients are near zero. Loading, the first feature, is the most important feature. A few manual tests show that selecting the five features with the highest logistic regression weights ('loading', 'attribute_3', 'measurement_2', 'measurement_4', 'measurement_17') increases the cv score to 0.59256.\n\nAt this point I suggest that you scroll back to the diagrams which show the correlation between the features and the target (the magenta scatterplots). Can you see why these five features get selected?\n\nIf you scroll back even more (up to the missing value count), you'll see that measurement_17 (second highest importance) is the feature with the most missing values. ","metadata":{}},{"cell_type":"code","source":"# Same code as in the cell above except for the assignment to features2\nauc_list = []\ntest_pred_list = []\nkf = GroupKFold(n_splits=5) # must be 5 because of the 5 product codes\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(train, train.failure, train.product_code)):\n    X_tr = train.iloc[idx_tr][test.columns]\n    X_va = train.iloc[idx_va][test.columns]\n    X_te = test.copy()\n    y_tr = train.iloc[idx_tr].failure\n    y_va = train.iloc[idx_va].failure\n    \n    # We one-hot encode attribute_0 and attribute_1\n    ohe_attributes = ['attribute_0', 'attribute_1']\n    ohe_output = ['ohe0_7', 'ohe1_6', 'ohe1_8']\n    ohe = OneHotEncoder(categories=[['material_5', 'material_7'],\n                                    ['material_5', 'material_6', 'material_8']],\n                        drop='first', sparse=False, handle_unknown='ignore')\n    ohe.fit(X_tr[ohe_attributes])\n    for df in [X_tr, X_va, X_te]:\n        with warnings.catch_warnings(): # ignore \"Found unknown categories\"\n            warnings.filterwarnings('ignore', category=UserWarning)\n            df[ohe_output] = ohe.transform(df[ohe_attributes])\n        df.drop(columns=ohe_attributes, inplace=True)\n\n    # We add the indicators for missing values\n    for df in [X_tr, X_va, X_te]:\n        df['m_3_missing'] = df.measurement_3.isna()\n        df['m_5_missing'] = df.measurement_5.isna()\n\n    # We fill the missing values\n    features = [f for f in X_tr.columns if f == 'loading' or f.startswith('measurement')]\n    imputer = KNNImputer(n_neighbors=3)\n    imputer.fit(X_tr[features])\n    for df in [X_tr, X_va, X_te]:\n        df[features] = imputer.transform(df[features])\n                \n    # The EDA diagram of measurement 2 shows that the feature is correlated\n    # to the target only for values above 10. For this reason, we clip\n    # all values below 11.\n    for df in [X_tr, X_va, X_te]:\n        df['measurement_2'] = df['measurement_2'].clip(11, None)\n    \n    # We fit a model using only the most important features\n    features2 = ['loading', 'attribute_3', 'measurement_2', 'measurement_4', 'measurement_17', 'm_3_missing', 'm_5_missing']\n    #model = RandomForestClassifier(n_estimators=200, max_depth=3, min_samples_leaf=100, n_jobs=-1, random_state=1)\n    #model = ExtraTreesClassifier(n_estimators=1000, max_depth=9, min_samples_leaf=100, n_jobs=-1, random_state=1)\n    model = make_pipeline(StandardScaler(), LogisticRegression())\n    model.fit(X_tr[features2], y_tr)\n    with np.printoptions(linewidth=150, precision=2, suppress=True):\n        print(model.named_steps['logisticregression'].coef_)\n    \n    # We validate the model\n    y_va_pred = model.predict_proba(X_va[features2])[:,1]\n    score = roc_auc_score(y_va, y_va_pred)\n    print(f\"Fold {fold}: auc = {score:.5f}\")\n    auc_list.append(score)\n\n    test_pred_list.append(model.predict_proba(X_te[features2])[:,1])\nprint(f\"{Fore.GREEN}{Style.BRIGHT}Average auc = {sum(auc_list) / len(auc_list):.5f}{Style.RESET_ALL}\")\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:32:01.486886Z","iopub.execute_input":"2022-08-06T18:32:01.487235Z","iopub.status.idle":"2022-08-06T18:34:15.271641Z","shell.execute_reply.started":"2022-08-06T18:32:01.487202Z","shell.execute_reply":"2022-08-06T18:34:15.268137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If you want to learn more about feature importance in linear models, [read here](https://scikit-learn.org/stable/auto_examples/inspection/plot_linear_model_coefficient_interpretation.html).","metadata":{}},{"cell_type":"markdown","source":"# ROC curve\n\nTo get a feeling for the auc metric, we plot the roc curve. The light red area under the curve corresponds to the AUC score.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(5, 5))\nfpr, tpr, _ = roc_curve(y_va, y_va_pred)\nplt.plot(fpr, tpr, color='#c00000', lw=3, label=f\"(auc (fold 4) = {roc_auc_score(y_va, y_va_pred):.5f})\") # curve\nplt.fill_between(fpr, tpr, color='#ffc0c0') # area under the curve\nplt.plot([0, 1], [0, 1], color=\"navy\", lw=1, linestyle=\"--\") # diagonal\nplt.gca().set_aspect('equal')\nplt.xlim([0.0, 1.0])\nplt.ylim([0.0, 1.0])\nplt.xlabel(\"False Positive Rate\")\nplt.ylabel(\"True Positive Rate\")\nplt.title(\"Receiver operating characteristic\")\nplt.legend(loc=\"lower right\")\nplt.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:34:15.273442Z","iopub.execute_input":"2022-08-06T18:34:15.273856Z","iopub.status.idle":"2022-08-06T18:34:15.442207Z","shell.execute_reply.started":"2022-08-06T18:34:15.273817Z","shell.execute_reply":"2022-08-06T18:34:15.441124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission\n\nWe finally create a submission file. I won't submit it myself because I trust the cv scores more than the public leaderboard.","metadata":{}},{"cell_type":"code","source":"submission = pd.DataFrame({'id': test.index,\n                           'failure': sum(test_pred_list)/len(test_pred_list)})\nsubmission.to_csv('submission_four_features.csv', index=False)\nsubmission","metadata":{"execution":{"iopub.status.busy":"2022-08-06T18:34:15.443865Z","iopub.execute_input":"2022-08-06T18:34:15.444249Z","iopub.status.idle":"2022-08-06T18:34:15.496057Z","shell.execute_reply.started":"2022-08-06T18:34:15.444215Z","shell.execute_reply":"2022-08-06T18:34:15.494935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To get a feeling for the predicted probabilities and as a sanity check, we plot histograms of the validation and test predictions. If the two histograms deviate from each other, this could indicate a bug in the model.","metadata":{}},{"cell_type":"code","source":"plt.hist(y_va_pred, bins=np.linspace(0, 0.6, 50), density=True, label='validation')\nplt.hist(submission.failure, bins=np.linspace(0, 0.6, 50), density=True,\n         color='orange', rwidth=0.4, label='test')\nplt.xlabel('y_pred')\nplt.ylabel('density')\nplt.title('Distribution of the predicted probabilities')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-06T18:34:15.497980Z","iopub.execute_input":"2022-08-06T18:34:15.500388Z","iopub.status.idle":"2022-08-06T18:34:15.767900Z","shell.execute_reply.started":"2022-08-06T18:34:15.500358Z","shell.execute_reply":"2022-08-06T18:34:15.767100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature correlations per product\n\nLet's look at the correlations of the continuous features for every product. Green matrix entries are positive correlations, red ones are negative. In most browsers, you can right-click on the matrices to zoom in.\n\nWe see:\n1. Most correlations are zero.\n2. loading isn't correlated with any measurement.\n3. The correlations matrices of the products differ from each other.\n4. Measurement_17 (the row with many green entries) seems to depend on measurement_3 to measurement_9.\n5. There are some correlations among measurement_10 to measurement_16.\n6. There are some correlations among the integer features measurement_0 to measurement_2.\n\n**Insight:** People like to show correlation matrices in EDAs, but they are often useless. In this competition, however, the correlations may become important: they will help impute missing values well, in particular for measurement_17.","metadata":{}},{"cell_type":"code","source":"_, axs = plt.subplots(3, 3, figsize=(18, 18))\nfor product, ax in zip(np.unique(both.product_code), axs.ravel()):\n    corr = both[float_cols + ['measurement_0', 'measurement_1', 'measurement_2']][both.product_code == product].corr()\n    mask = np.triu(np.ones_like(corr, dtype=bool))\n    sns.heatmap(corr*10, mask=mask, linewidth=0.0, fmt='.0f', \n                annot=True, annot_kws={'size': 8}, \n                cmap='PiYG', center=0, ax=ax, cbar=False)\n    ax.set_title(product)\nplt.tight_layout(w_pad=0.5)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:34:15.768933Z","iopub.execute_input":"2022-08-06T18:34:15.769205Z","iopub.status.idle":"2022-08-06T18:34:27.956086Z","shell.execute_reply.started":"2022-08-06T18:34:15.769181Z","shell.execute_reply":"2022-08-06T18:34:27.954999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Per-product distributions of the features\n\nOne can ask how the feature distributions depend on the product codes. The diagram shows that the first eight features are independent of the product code. Measurement_10 through measurement_17 and measurement_0 through measurement_2 depend on the product.\n\n**Insight:** When we impute the missing values for these features, the imputed value should depend on the product code. SimpleImputer and KNNImputer do not take this dependence into account.","metadata":{}},{"cell_type":"code","source":"_, axs = plt.subplots(4, 4, figsize=(12,12))\nfor f, ax in zip(float_cols, axs.ravel()):\n    mi = both[f].min()\n    ma = both[f].max()\n    bins = np.linspace(mi, ma, 40)\n    for product in np.unique(both.product_code):\n        h, edges = np.histogram(both[f][both.product_code == product], bins=bins)\n        ax.plot((edges[:-1] + edges[1:]) / 2, h, label=product, lw=3)\n    ax.set_xlabel(f)\n    if ax == axs[0, 0]: ax.legend(loc='upper right')\nplt.tight_layout(w_pad=1)\nplt.suptitle('Distributions of the continuous features for every product', fontsize=20, y=1.02)\nplt.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:34:27.957770Z","iopub.execute_input":"2022-08-06T18:34:27.958925Z","iopub.status.idle":"2022-08-06T18:34:30.444044Z","shell.execute_reply.started":"2022-08-06T18:34:27.958877Z","shell.execute_reply":"2022-08-06T18:34:30.443042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prop_cycle = plt.rcParams['axes.prop_cycle']\n\nfig, axs = plt.subplots(1, 3, figsize=(12, 4))\naxs = axs.ravel()\nfor(ax, f) in zip(axs, ['measurement_0', 'measurement_1', 'measurement_2']):\n    for i, product in enumerate(np.unique(both.product_code)):\n        uv, uc = np.unique(both[f][both.product_code == product], return_counts=True)\n        ax.plot(uv, uc, alpha=1, color=prop_cycle.by_key()['color'][i % 10],\n                lw=3, label=product)\n    ax.set_xlabel(f)\n    if ax == axs[0]: ax.legend()\nplt.suptitle('Distributions of the integer features for every product', y=1.01, fontsize=20)\nplt.tight_layout(w_pad=1)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-06T18:34:30.445186Z","iopub.execute_input":"2022-08-06T18:34:30.445539Z","iopub.status.idle":"2022-08-06T18:34:31.022922Z","shell.execute_reply.started":"2022-08-06T18:34:30.445508Z","shell.execute_reply":"2022-08-06T18:34:31.021794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}