{"cells":[{"metadata":{"papermill":{"duration":0.005774,"end_time":"2020-09-08T02:02:36.044631","exception":false,"start_time":"2020-09-08T02:02:36.038857","status":"completed"},"tags":[]},"cell_type":"markdown","source":"The objective of this notebook is to look for patterns across the 206 scored targets. It might be possible to find patterns that trace back to features and give some incremental reduction in logloss. \n\nUpdate: There is a pattern that shows improvement in logloss for one of the scored targets.I haven't tried it on the full target set. It might be useful if you are ensembling high-scoring deep learning models with a tabular-type model."},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:36.059426Z","iopub.status.busy":"2020-09-08T02:02:36.058498Z","iopub.status.idle":"2020-09-08T02:02:36.06142Z","shell.execute_reply":"2020-09-08T02:02:36.060922Z"},"papermill":{"duration":0.012054,"end_time":"2020-09-08T02:02:36.061539","exception":false,"start_time":"2020-09-08T02:02:36.049485","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"import pandas as pd","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.004214,"end_time":"2020-09-08T02:02:36.070518","exception":false,"start_time":"2020-09-08T02:02:36.066304","status":"completed"},"tags":[]},"cell_type":"markdown","source":"# Part 1. Scored targets\n\n\n## The control group pattern\n\nIn this [notebook](https://www.kaggle.com/artgor/lish-moa-baseline-approach), 2x Grandmaster @artfor mentions that observations in the control group show no response in the train set targets. Let's confirm that claim."},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:36.09015Z","iopub.status.busy":"2020-09-08T02:02:36.084015Z","iopub.status.idle":"2020-09-08T02:02:40.064748Z","shell.execute_reply":"2020-09-08T02:02:40.06402Z"},"papermill":{"duration":3.987769,"end_time":"2020-09-08T02:02:40.064877","exception":false,"start_time":"2020-09-08T02:02:36.077108","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"targets =  pd.read_csv(\"../input/lish-moa/train_targets_scored.csv\",\n                 index_col=['sig_id'])\ntargets_0_idx = targets[targets.sum(axis=1)==0].index\n\ntrain_features =  pd.read_csv(\"../input/lish-moa/train_features.csv\",\n                 index_col=['sig_id'])\ncontrol_idx = train_features.query('cp_type==\"ctl_vehicle\"').index\n\ndiffs = len(set(control_idx) - set(targets_0_idx))\n","execution_count":null,"outputs":[]},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:40.083377Z","iopub.status.busy":"2020-09-08T02:02:40.08248Z","iopub.status.idle":"2020-09-08T02:02:40.753031Z","shell.execute_reply":"2020-09-08T02:02:40.753529Z"},"papermill":{"duration":0.683655,"end_time":"2020-09-08T02:02:40.753683","exception":false,"start_time":"2020-09-08T02:02:40.070028","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"test_features =  pd.read_csv(\"../input/lish-moa/test_features.csv\",\n                 index_col=['sig_id'])\ntest_control_idx = test_features.query('cp_type==\"ctl_vehicle\"').index\n\nctrl_pct = len(set(test_control_idx))/len(test_features)\n\nprint(diffs, ctrl_pct)","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.004161,"end_time":"2020-09-08T02:02:40.762426","exception":false,"start_time":"2020-09-08T02:02:40.758265","status":"completed"},"tags":[]},"cell_type":"markdown","source":"Yes, it is the case that all ids in the train set's control group show no response. So there is our first observable pattern in the targets that translates to the feature set: All ids in the control group - about 9% of the test set - will likely have all 0 targets."},{"metadata":{"papermill":{"duration":0.004022,"end_time":"2020-09-08T02:02:40.770966","exception":false,"start_time":"2020-09-08T02:02:40.766944","status":"completed"},"tags":[]},"cell_type":"markdown","source":"## Correlation in the treated cases\n\nOne way to find aggregated patterns in targets is of course to look at correlations. Here is a correlation matrix heatmap that shows points for any two targets with above 70% correlation (positive or negative). There are 4 such pairs with two of those pairs above 90%. (The heatmap shows two permutations per pair, one above and one below the diagonal.) "},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:40.786195Z","iopub.status.busy":"2020-09-08T02:02:40.785499Z","iopub.status.idle":"2020-09-08T02:02:44.33917Z","shell.execute_reply":"2020-09-08T02:02:44.33963Z"},"papermill":{"duration":3.564548,"end_time":"2020-09-08T02:02:44.339774","exception":false,"start_time":"2020-09-08T02:02:40.775226","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"import numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\n\ntargets_treated = targets.drop(targets_0_idx)\ncorr_mtx = abs(targets_treated.corr())\n\ncorr_map = corr_mtx[corr_mtx>=.7]\nplt.figure(figsize=(12,8))\nsns.heatmap(corr_map, cmap=\"viridis\")","execution_count":null,"outputs":[]},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:44.357792Z","iopub.status.busy":"2020-09-08T02:02:44.357109Z","iopub.status.idle":"2020-09-08T02:02:44.370035Z","shell.execute_reply":"2020-09-08T02:02:44.369517Z"},"papermill":{"duration":0.024181,"end_time":"2020-09-08T02:02:44.37015","exception":false,"start_time":"2020-09-08T02:02:44.345969","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"pairs = (corr_mtx.where(np.triu(np.ones(corr_mtx.shape), k=1).astype(np.bool))\n                     .stack()\n                     .sort_values(ascending=False))\npairs[:4]","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.005631,"end_time":"2020-09-08T02:02:44.382065","exception":false,"start_time":"2020-09-08T02:02:44.376434","status":"completed"},"tags":[]},"cell_type":"markdown","source":"Correlation drops right off and is very weak for the remaining 185 pairs. The histogram below shows how most test cases show 0 or 1 response across all targets. "},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:44.399269Z","iopub.status.busy":"2020-09-08T02:02:44.398673Z","iopub.status.idle":"2020-09-08T02:02:44.644371Z","shell.execute_reply":"2020-09-08T02:02:44.643771Z"},"papermill":{"duration":0.256561,"end_time":"2020-09-08T02:02:44.644479","exception":false,"start_time":"2020-09-08T02:02:44.387918","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"targets.sum(axis=1).hist()","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.005891,"end_time":"2020-09-08T02:02:44.656361","exception":false,"start_time":"2020-09-08T02:02:44.65047","status":"completed"},"tags":[]},"cell_type":"markdown","source":"Looking again at the 4 pairs above, we see the 3 pairs below the top pair are all related. Here's a distance graph of those along with 16 more pairs. The networkx package is magic in this application!"},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:44.675685Z","iopub.status.busy":"2020-09-08T02:02:44.674876Z","iopub.status.idle":"2020-09-08T02:02:45.399879Z","shell.execute_reply":"2020-09-08T02:02:45.400352Z"},"papermill":{"duration":0.737903,"end_time":"2020-09-08T02:02:45.400491","exception":false,"start_time":"2020-09-08T02:02:44.662588","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"import networkx as nx\n\npairs_df = (1-pairs).reset_index()\nG = nx.from_pandas_edgelist(pairs_df[:20], source='level_0', target='level_1', edge_attr=0)\n\ngraph_opts = dict(arrows=False,\n                  node_size=5,\n                  width=2,\n                  alpha=0.8,\n                  font_size=12,\n                  font_color='darkblue',\n                  edge_color='darkgray'\n                 )\n\nfig= plt.figure(figsize=(12,10))\nnx.draw_spring(G, with_labels=True, **graph_opts)","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.007924,"end_time":"2020-09-08T02:02:45.416009","exception":false,"start_time":"2020-09-08T02:02:45.408085","status":"completed"},"tags":[]},"cell_type":"markdown","source":"## Another pattern\n\nThe next step is to look at individual samples and see which sets look similar. We could go group by group through the above network, or we can look for similarities among cases from a different angle. Remember, we want to trace these patterns in the target set back to the feature set. \n\nHere I use the [missingno](https://github.com/ResidentMario/missingno) package by Notebooks Master @residentmario to help find patterns. There are a few interesting things at which to look (check that proper grammar!). Notice the dense bands at the top leftish part of the chart."},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:45.437431Z","iopub.status.busy":"2020-09-08T02:02:45.436746Z","iopub.status.idle":"2020-09-08T02:02:46.336808Z","shell.execute_reply":"2020-09-08T02:02:46.337234Z"},"papermill":{"duration":0.913888,"end_time":"2020-09-08T02:02:46.337405","exception":false,"start_time":"2020-09-08T02:02:45.423517","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"import missingno as msno\n\ncols_sorted = targets_treated.sum().sort_values(ascending=False).index\ntargets_visible = targets_treated[cols_sorted].replace(0, np.nan)\n                                                # you can also use pd.NA with pandas v1+\n\nmsno.matrix(targets_visible.iloc[:, :50].sample(n=1000), sort='descending', color=(1,0,0))","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.008193,"end_time":"2020-09-08T02:02:46.354417","exception":false,"start_time":"2020-09-08T02:02:46.346224","status":"completed"},"tags":[]},"cell_type":"markdown","source":"We see this picture upon drilling down."},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:46.375312Z","iopub.status.busy":"2020-09-08T02:02:46.374657Z","iopub.status.idle":"2020-09-08T02:02:46.658156Z","shell.execute_reply":"2020-09-08T02:02:46.657606Z"},"papermill":{"duration":0.295122,"end_time":"2020-09-08T02:02:46.658274","exception":false,"start_time":"2020-09-08T02:02:46.363152","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"top_cols, idx = np.unique(pairs_df.iloc[:4, :2].values.flatten(), return_index=True)\n                                                    # use df.to_numpy() for 1.0+\nmsno.matrix(targets_visible.loc[targets_treated[top_cols].any(axis=1), \n                top_cols[idx]].sort_values(['flt3_inhibitor', 'pdgfr_inhibitor', 'kit_inhibitor']), \n                sort=None, color=(1,0,0))\n","execution_count":null,"outputs":[]},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:46.683711Z","iopub.status.busy":"2020-09-08T02:02:46.682557Z","iopub.status.idle":"2020-09-08T02:02:46.68658Z","shell.execute_reply":"2020-09-08T02:02:46.686992Z"},"papermill":{"duration":0.019684,"end_time":"2020-09-08T02:02:46.687136","exception":false,"start_time":"2020-09-08T02:02:46.667452","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"targets.proteasome_inhibitor.sum()/len(targets)","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.008816,"end_time":"2020-09-08T02:02:46.705534","exception":false,"start_time":"2020-09-08T02:02:46.696718","status":"completed"},"tags":[]},"cell_type":"markdown","source":"Well this is rather interesting. One way to characterize it is that if a case shows a positive response for the proteasome inhibitor, it will almost always show a positive response for the nfkb_inhibitor as well. Furthermore, the case is very likely to show no repsonse to the other 3 inhibitors in the chart. It's 3% of the populaton, which is potentially nice."},{"metadata":{"papermill":{"duration":0.008966,"end_time":"2020-09-08T02:02:46.723641","exception":false,"start_time":"2020-09-08T02:02:46.714675","status":"completed"},"tags":[]},"cell_type":"markdown","source":"## More patterns?\n\nThe network graph above showed several relationships in the targets. Missingno has a dendrogram tool that does something similar and unlike the matrix, doesn't depend on sort order. Here is the dendrogram for the cases receiving real treatment.\n\nShorter bars indicate tighter relationships."},{"metadata":{"execution":{"iopub.execute_input":"2020-09-08T02:02:46.751878Z","iopub.status.busy":"2020-09-08T02:02:46.751253Z","iopub.status.idle":"2020-09-08T02:02:49.875665Z","shell.execute_reply":"2020-09-08T02:02:49.876141Z"},"papermill":{"duration":3.142763,"end_time":"2020-09-08T02:02:49.87628","exception":false,"start_time":"2020-09-08T02:02:46.733517","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"msno.dendrogram(targets_visible)","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.013713,"end_time":"2020-09-08T02:02:49.903984","exception":false,"start_time":"2020-09-08T02:02:49.890271","status":"completed"},"tags":[]},"cell_type":"markdown","source":"# Part 2. Features\n\n## Feature importance\n\nDetecting patterns is uesful in the competition only if we can trace it back to patterns in the features. And then, it depends on whether the model we use found those patterns already. In my experience machine learning models usually find the patterns before humans do - that's why we're all here. Sometimes though we can find incremental improvements that clue us in to new feature combinations.\n\nThe relationship of the 5 inhibitors seems a good place to start. I'll use a RandomForestClassifier on the proteasome inhibitor target as a simple baseline for the model and for finding important features. The imbalanced data makes results very sensitive to data splits and random seeds, so we want something that can generalize. I'm using the ELI5 library and permutation importance based on several iterations of a cross-validated model. "},{"metadata":{"trusted":true},"cell_type":"code","source":"from multiprocessing import cpu_count\nimport eli5\nfrom eli5.sklearn import PermutationImportance\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.model_selection import cross_val_score\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.metrics import log_loss, make_scorer\n\nnproc = cpu_count()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train_features\ny = targets.proteasome_inhibitor\n\nX[['cp_type', 'cp_dose']] = X[['cp_type', 'cp_dose']].astype('category').apply(lambda x: x.cat.codes)\n\nrf_model = RandomForestClassifier(n_estimators=100, random_state=10, verbose=True, n_jobs=nproc)\nscorer = make_scorer(log_loss)\nskf = StratifiedKFold(n_splits=5, random_state=24)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"%%time\n\nperm = PermutationImportance(rf_model, scoring=scorer, cv=skf)\nperm.fit(X, y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weights = eli5.explain_weights_dfs(perm, feature_names=X.columns.tolist())\nweights_df = weights['feature_importances']\nweights_df[:15]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature engineering\n\nNow we can add some of the important features in as combinations of each other. I'll use pairwise combinations here. This technique usually has more of a benefit for linear models than for tree-based models. I didn't notice any competitive linear models in the notebooks though and am checking it out this way."},{"metadata":{"trusted":true},"cell_type":"code","source":"from itertools import combinations\n\ntrain_features_with = train_features.copy()\nimportant = weights_df.loc[:5, 'feature'].tolist()\nfor pair in combinations(important, 2):\n    col = \"_\".join(pair)\n    train_features_with[col] = train_features_with[pair[0]] * train_features_with[pair[1]]\n\ntrain_features_with[:5]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_with = train_features_with\nrf_model_with = RandomForestClassifier(n_estimators=100, random_state=10, verbose=True, n_jobs=nproc)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Part 3. Models and scores\n\n## Comparison\nHere are loglosses for the baseline and augmented set. Again, there are several iterations with the same 5-fold CV splits within an iteration."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nlosses = np.zeros((2, 3, 5))\nfor i in range(3):\n    skf = StratifiedKFold(n_splits=5, random_state=i*8)\n    scorer = make_scorer(log_loss, greater_is_better=False, needs_proba=True)\n    losses[0,i] = cross_val_score(rf_model, X, y, scoring=scorer, cv=skf, n_jobs=nproc)\n    losses[1,i] = cross_val_score(rf_model_with, X_with, y, scoring=scorer, cv=skf, n_jobs=nproc)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f\"baseline: {np.mean(losses[0])} mean, {np.std(losses[0])} std dev \"\n      f\"with features: {np.mean(losses[1])} mean, {np.std(losses[1])} std dev\"\n      )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"markdown","source":"## Closing thoughts\n\nThere is about an 8% improvement in the mean when using the features, which is promising. However, the standard deviation runs at ~10%. Even a bump on the public leaderboard may not translate to gains on the private LB.\n\nAlso there's the matter of the other 205 targets! I hope at least the notebook gives you some ideas for feature engineering based on patterns in the targets. \n\nGood luck for the duration of the contest!"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}