{"cells":[{"metadata":{"_uuid":"f1c226f0755e6517f4c6bc7cf3f20ccf25fdc94b"},"cell_type":"markdown","source":"# Missingness Data != Missing Any Information"},{"metadata":{"trusted":true,"_uuid":"8564283fd813631b45f67d86cfc77287a7803004"},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom matplotlib import style\nsns.set_style('whitegrid')\nplt.rcParams['figure.figsize'] = (12, 8)\n%matplotlib inline\nimport re\nfrom itertools import product\nimport itertools\n\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.metrics import matthews_corrcoef, roc_curve, auc, roc_auc_score, classification_report, confusion_matrix\nfrom sklearn.model_selection import train_test_split\n\nfrom scipy import stats","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ffe88ee88b6dd1d2cc15c793530d3fa33fe5d8fa"},"cell_type":"markdown","source":"In this notebook, we will explore the relationship between the missingness of data fields and their values."},{"metadata":{"trusted":true,"_uuid":"2e7a8834c8baccc271e2d3d54c436ccc6732dcb2"},"cell_type":"code","source":"dat = pd.read_pickle(\"../input/gstore-revenue-data-preprocessing/train.pkl\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c022b63a00ea12cd8d5935b69ed629d5703a37a6"},"cell_type":"markdown","source":"If we look at the counts of NA values in each attribute below, we can see that for many of the attributes, the data is missing for a significant portion of the rows. For many of these attributes, \"missing\" does not exactly mean \"we know nothing about it\". In fact, in many cases, \"missing\" can be a source of useful insight just as valid values do. For instance, if we failed to record the number of pages viewed, does this mean the user did not view anything?\nMissingness can also reveal which attributes were collected together and could potentially dependent on each other."},{"metadata":{"trusted":true,"_uuid":"2a65069334e3152e41c8ce63623086e694c755ca"},"cell_type":"code","source":"dat.apply(lambda x: np.sum(pd.isna(x)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"263af40b4e6e7a30ac812dabcb9555107382e26c"},"cell_type":"markdown","source":"### Independence Test of Column Values vs Column Missingness"},{"metadata":{"_uuid":"fb8d41527c6780aca4fb4105a39d5a18d4c4475d"},"cell_type":"markdown","source":"As a first step, we would like to know if some of the columns are dependent on the missingness of other columns. For categorical columns, we can apply the chi-square test to determine if they are independent from column missingnesses. If a significant number of columns appear to be dependent on column missingness, there is a good reason to include these missingnesses in models based on them.\n\nTo do so, we will first have to find all the categorical columns as well as columns with missing values, and create a new dataframe including both data columns and missingness indicator columns."},{"metadata":{"trusted":false,"_uuid":"feddea57021c56d1702694f66742a44893eda8cd"},"cell_type":"code","source":"cat_columns = [c for c in dat.columns if str(dat[c].dtype) == 'category']","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"2253cb874f1e87ca0727f5a854c33fbac957e30f"},"cell_type":"code","source":"missing_count = dat.apply(lambda x: np.sum(pd.isna(x)))\ncol_w_missing = list(missing_count[missing_count > 0].index)\ncol_w_missing","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"c062c0fd4c309a17d1d876f264d18f01557c117a"},"cell_type":"code","source":"missing = dat.copy()\nfor col in col_w_missing:\n    missing['miss_' + col] = pd.isnull(dat[col])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2e11624a140dea479bbf4012e74eeb676ae87e98"},"cell_type":"markdown","source":"Since in our preprocessing step, we converted missing revenue values to 0, here we add the revenue missingness column back."},{"metadata":{"trusted":false,"_uuid":"9a0a3957294cb6711628b43cd21638772e77b291"},"cell_type":"code","source":"zero_revenue = missing['totals.transactionRevenue'] == 0\nmissing['miss_totals.transactionRevenue'] = zero_revenue\ncol_w_missing.append('totals.transactionRevenue')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6f1920d76ae31d00990009a5512b807cbfb64631"},"cell_type":"markdown","source":"Now we can perform a pairwise chi2 independence test for categorical columns vs column missingness:"},{"metadata":{"trusted":false,"_uuid":"0abfe29bc61a192f578578bfc7bffcf4ceea7423"},"cell_type":"code","source":"ind_miss_p = np.full((len(cat_columns), len(col_w_missing)), np.nan)\nfor i, j in product(\n        range(len(cat_columns)), range(len(col_w_missing))):\n    chi2, p, dof, ex = stats.chi2_contingency(\n        missing.groupby([cat_columns[i], 'miss_' + col_w_missing[j]\n                         ]).size().unstack().fillna(0).astype(np.int))\n    ind_miss_p[i, j] = p\n    \nmiss_ind_test_output = pd.DataFrame(\n    ind_miss_p,\n    index=cat_columns,\n    columns=['miss_' + c for c in col_w_missing])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"2c4e5edd580757e3e1abed6bcf97fe5c063a67a9"},"cell_type":"code","source":"f, ax = plt.subplots(figsize=(12, 9))\nsns.heatmap(data=miss_ind_test_output, ax=ax, linewidths=0.01)\nax.set_title(\"p-values of chi2 independence test of categorical values vs missingness\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2ab6dd1ac242a7f93b79b27b29284faeb40695eb"},"cell_type":"markdown","source":"As we see above, the value of many of the columns appear to be dependent on many other columns' missingness (e.g. chi2 statistic large or p-value sufficiently small), we know that when considering the whole dataset, the values are not missing at random. There is potential information to be extracted from the missing values or the relationship between existing values and missing values. A missing value might indicate a specific state of the user or session that has an effect on the existing values, or an existing value might give away clues on what a missing value should have been if it were not missing."},{"metadata":{"_uuid":"da53b539c63290e31476a606583fe15e531dd48c"},"cell_type":"markdown","source":"### Independence Test of Column Missingness"},{"metadata":{"_uuid":"d56d6b0b689c70a43ad05022ad08b69e38d21d04"},"cell_type":"markdown","source":"Now that we know some of the columns are dependent on column missingness, what about the relationship between the missingness of different columns? Here we perform the same chi2 test, except with only cloumn missingness and between themselves:"},{"metadata":{"trusted":false,"_uuid":"8bb5dae9f6afbf6ae8fff55b3c035042e01f3532"},"cell_type":"code","source":"ind_miss2miss_p = np.full((len(col_w_missing), len(col_w_missing)), 0.)\nfor i, j in product(range(len(col_w_missing)), range(len(col_w_missing))):\n    if i < j:\n        chi2, p, dof, ex = stats.chi2_contingency(\n            missing.groupby([\n                'miss_' + col_w_missing[i], 'miss_' + col_w_missing[j]\n            ]).size().unstack().fillna(0).astype(np.int))\n        ind_miss2miss_p[i, j] = p\n        ind_miss2miss_p[j, i] = ind_miss2miss_p[i, j]\n    elif i == j:\n        ind_miss2miss_p[i, j] = 0\n\nmiss2miss_p_output = pd.DataFrame(\n    ind_miss2miss_p,\n    index=['miss_' + c for c in col_w_missing],\n    columns=['miss_' + c for c in col_w_missing])\n\ng = sns.clustermap(\n    data=miss2miss_p_output, figsize=(12, 12), linewidths=0.01)\ng.ax_col_dendrogram.set_title(\"pairwise p-value of column missingness independence test\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9bb439605ff0e971cf598f708bf713e8f039ce82"},"cell_type":"markdown","source":"As we see here, there are two major clusters in the pairwise p-value heatmap. Remember that larger p-value (brighter colour) indicates that the pair is more likely to be independent. The upper left corner has four columns that are dependent on each other but mostly independent from other columns **(device.browser, trafficSource.source, totals.pageviews, trafficSource.medium)**, and the lower right corner has a large number of columns that are all dependent on each other.\n\nThe analysis above tells us whether there are relationships between the missingnesses, but not how they are related to each other. Apart from independence, we would also like to know if the missingness of different columns are \"in sync\" which each other, are good predictors of each other, or at least offer much information about each other.\n\nWe first analyse if some of the missingnesses are \"in sync\", e.g. tend to happen together. It is roughly the same as asking if one missingness is a good predictor of another. Here we will be using pairwise Matthews correlation coefficient, a common measure for binary classification evaluation. A coefficient of +1 represents a perfect prediction, 0 no better than random prediction and −1 indicates total disagreement between prediction and observation. We also perform heatmap clustering to identify clusters of columns that are closely related to each other."},{"metadata":{"trusted":false,"_uuid":"7b9e2df0982c3866fcf022b42f8c055c7e6b18d9"},"cell_type":"code","source":"ind_miss2miss_mcc = np.full((len(col_w_missing), len(col_w_missing)), 0.)\nfor i, j in product(range(len(col_w_missing)), range(len(col_w_missing))):\n    if i < j:\n        ind_miss2miss_mcc[i, j] = matthews_corrcoef(\n            missing['miss_' + col_w_missing[i]],\n            missing['miss_' + col_w_missing[j]])\n        ind_miss2miss_mcc[j, i] = ind_miss2miss_mcc[i, j]\n    elif i == j:\n        ind_miss2miss_mcc[i, j] = 1\n\nmiss2miss_mcc_output = pd.DataFrame(\n    ind_miss2miss_mcc,\n    index=['miss_' + c for c in col_w_missing],\n    columns=['miss_' + c for c in col_w_missing])\nmiss2miss_mcc_output.index.name = 'predicted'\nmiss2miss_mcc_output.columns.name = 'input'\n\ng = sns.clustermap(\n    data=miss2miss_mcc_output, figsize=(12, 12), linewidths=0.01)\ng.ax_col_dendrogram.set_title(\"pairwise MCC score of column missingness\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e47a99d4664ece8f81871eedd0083da297971fa6"},"cell_type":"markdown","source":"Here we can clearly see four clusters of cloumns and eight other relatively isolated columns. \n\nIn the top left we see seven columns related to the ad contents. Among them, the **trafficSource.adwardsClickInfo** attributes are more closely related with each other than the others (they are always missing at the same time). \n\nThen we see a separation between the **geoNetwork** attributes, with **(country, continent, SubContinent)** in one cluster, always appearing togther and with **(metro, city, region)** in another cluster, mostly appearing together. This is a good indicator that there might be two separate sources of data for these attributes, and potential (actually proven to exist) conflicts between the two clusters can be explained that way. We also notice that **networkDomain**, despite not in any of these two clusters, appear to have higher score with the first clutser than the second, indicating that it is less likely related to the second cluster. \n\nThen we find that **medium** and **source** are related in missingness.\n\nThis is all useful, but sometimes we are not too concerned about whether one column's missingness is a good predictor of another. We just want to know if one column can tell us some information about another, even if it is very noisy information. We need some other measures that are more about information gain or \"doing better than random\", such as AUC, entropy, etc."},{"metadata":{"trusted":false,"_uuid":"7cc0b630da72bb5c8ad1020986d8e2bdb58d4cca"},"cell_type":"code","source":"ind_miss2miss_auc = np.full((len(col_w_missing), len(col_w_missing)), 0.)\nfor i, j in product(range(len(col_w_missing)), range(len(col_w_missing))):\n        score1 = roc_auc_score(missing['miss_' + col_w_missing[i]],\n                                          missing['miss_' + col_w_missing[j]])\n        score2 = roc_auc_score(missing['miss_' + col_w_missing[i]],\n                                          ~missing['miss_' + col_w_missing[j]])\n        ind_miss2miss_auc[i, j] = max(score1, score2)\n        \nmiss2miss_auc_output = pd.DataFrame(\n    ind_miss2miss_auc,\n    index=['miss_' + c for c in col_w_missing],\n    columns=['miss_' + c for c in col_w_missing])\nmiss2miss_auc_output.index.name = 'predicted'\nmiss2miss_auc_output.columns.name = 'input'\n\ng = sns.clustermap(data=miss2miss_auc_output, figsize=(12, 12), linewidths=0.01)\ng.ax_col_dendrogram.set_title(\"pairwise AUC score of column missingness\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d35ec9fa0efffd292413592efb5159ca5b0151f5"},"cell_type":"markdown","source":"This pairwise AUC heatmap tells us whether the missingness of attributes at the bottom gives us useful information about attributes on the right. Some observations are expected, such as the two **geoNetwork** clusters still being present here, as well as the presense of the ad contents cluster. However, there are something unexpected as well, such as **totals.transactionRevenue**, the attribute we are most interested about, actually leaks some information about it in several other columns! The missingness of **totals.bounces** appear to tell us a great deal about whether revenue exist, but the missingnesses of **geoNetwork.metro** and **totals.newVisits** also reveals a little. Let us plot the ROC graph of predicting **transcationRevenue** missingness with these columns:"},{"metadata":{"trusted":false,"_uuid":"6a7d8dc1b31553c9354159c6a53e783b47f99778"},"cell_type":"code","source":"cur_dict = dict()\ncols = [c for c in col_w_missing if c != 'totals.transactionRevenue']\nfor c in cols:\n    fpr_p, tpr_p, _ = roc_curve(~missing['miss_totals.transactionRevenue'],\n                                missing['miss_' + c])\n    fpr_n, tpr_n, _ = roc_curve(~missing['miss_totals.transactionRevenue'],\n                                ~missing['miss_' + c])\n    auc_p, auc_n = auc(fpr_p, tpr_p), auc(fpr_n, tpr_n)\n    if auc_p >= 0.55:\n        cur_dict[c] = [fpr_p, tpr_p, auc_p]\n    elif auc_n >= 0.55:\n        cur_dict[c] = [fpr_n, tpr_n, auc_n]\n\nplt.figure(figsize=(12, 12))\nlw = 2\nfor c, v in cur_dict.items():\n    plt.plot(v[0], v[1], lw=lw, label=\"{0}  AUC={1}\".format(c, v[2]))\nplt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\nplt.xlabel('False Positive Rate')\nplt.ylabel('True Positive Rate')\nplt.legend(loc=\"lower right\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f70408d7172d0a14dbd63984ef3fa2d4c5d02c81"},"cell_type":"markdown","source":"Using **totals.bounces** missingness to catch cases with positive revenue is surprisingly good, reaching 100% success with only about 45% false positive. As we see below, there are no cases where **bounces** and **transactionRevenue** are both present, so it appears that **transactionRevenue** can only be positive if **bounces** is missing. This alone will not make a good predictor though, as the misclassification of cases where revenue = 0 will be exceedingly high. Nevertheless, this is some information we can use with almost certainty, and is definitely better than no information."},{"metadata":{"trusted":false,"_uuid":"031acd9e6ea642aa253da510a460b7d8cd096395"},"cell_type":"code","source":"missing.groupby(['miss_totals.transactionRevenue', 'miss_totals.bounces']).size().unstack().fillna(0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9ce23f48003086bfaf18d4523b7166cf84e72348"},"cell_type":"markdown","source":"What if we try to use *all* the column missingness to predict the missingness of revenue? Here we go, using a random forest:"},{"metadata":{"trusted":false,"_uuid":"cd7dbca447974a738fe591f37d343c34167378f5"},"cell_type":"code","source":"X = missing.loc[:, [\n    c for c in missing.columns if re.match(r'miss_', c) is not None\n    and c != 'miss_totals.transactionRevenue'\n]]\ny = ~missing['miss_totals.transactionRevenue']","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"8c8242109efb58161487d0c8ad2211c03f16fd54"},"cell_type":"code","source":"train_X, test_X, train_y, test_y = train_test_split(X, y, test_size=0.3, random_state=7777)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"e4c9405a0ce0345f789ef34387e61b4e0c27c258"},"cell_type":"code","source":"clf = RandomForestClassifier(n_estimators=100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"375d34767e25104c07202a91e1c71f2da76c2232"},"cell_type":"code","source":"clf.fit(train_X, train_y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"bed7f47c85c9da5cac466dc87103f27e14c2e60f"},"cell_type":"code","source":"preds = clf.predict(test_X)\nprobs = clf.predict_proba(test_X)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5af34bb0f670a7ef637014b7631ef0418938f600"},"cell_type":"markdown","source":"As we see, due to the extreme imbalance between revenue = 0 and revenue > 0 cases, the classifier does not learn to classify the positive class (revenue > 0), unfortunately."},{"metadata":{"trusted":false,"_uuid":"b7d6b7e7cad4c0ff71fbd65ba3c66f507d84174f"},"cell_type":"code","source":"print(classification_report(test_y, preds))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c6c0d3c6e94b92dcd716db8a3c754fc486245117"},"cell_type":"markdown","source":"However, the classifier does learn to find likely suspects of the positive class, if we loosen the threshold a little (read: a lot). Yes, we end up with very poor precision, but remember that this is with just column missingness without touching the actual data, and we already found a way to exclude many rows that cannot be in the positive class."},{"metadata":{"trusted":false,"_uuid":"e3da6df5db29b7a02ccaa81f6d2e82371b796188"},"cell_type":"code","source":"print(classification_report(test_y, probs[:, 1] > 0.01))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"17d8be428a07dabae00cb72268cb5ce578b9122c"},"cell_type":"code","source":"pd.DataFrame(\n    confusion_matrix(test_y, probs[:, 1] > 0.01),\n    columns=['pred_miss', 'pred_exist'],\n    index=['miss', 'exist'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c3a3c221e8684fe54f429ccd02f1613f792eb5ee"},"cell_type":"markdown","source":"With the ROC curve below, we can see that with the combined might of all column missingness, we can do better than just using the most informative column **bounces**."},{"metadata":{"trusted":false,"_uuid":"9df82b3ab18e271200a098927f72813282c1aaa2"},"cell_type":"code","source":"fpr, tpr, _ = roc_curve(test_y, probs[:, 1])\nauc_score = auc(fpr, tpr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"16656ad35ae61471201ac7c34d2450adacb9a242"},"cell_type":"code","source":"plt.figure(figsize=(12, 12))\nlw = 2\nc = 'totals.bounces'\nv = cur_dict[c]\nplt.plot(v[0], v[1], lw=lw, label=\"{0}  AUC={1}\".format(c, v[2]))\nplt.plot(fpr, tpr, lw=lw, label=\"{0}  AUC={1}\".format('RF classifier', auc_score))\nplt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\nplt.xlabel('False Positive Rate')\nplt.ylabel('True Positive Rate')\nplt.legend(loc=\"lower right\")\nplt.title('ROC of RF Classifier Based on Missingness')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"e4b4efbabebb1a5e70bdba4c3a05a5fb537cd9e7"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"hide_input":false,"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.0"}},"nbformat":4,"nbformat_minor":1}