a-complete-introduction-and-walkthrough.ipynb

{"cells":[{"metadata":{"_uuid":"cbb767a6ce7a389511a40311cd4f01f1b09ad270"},"cell_type":"markdown","source":"# Costa Rican Household Poverty Level Prediction\n\n\nWelcome to another Kaggle challenge! The objective of the Costa Rican Household Poverty Level Prediction contest is to develop a machine learning model that can predict the poverty level of households using both individual and household characteristics. This \"data science for good\" project offers the opportunity to put our skills towards a task more beneficial to society than getting people to click on ads!\n\nIn this notebook, we will walk through a complete machine learning solution: first, get introduced to the problem, then perform a thorough Exploratory Data Analysis of the dataset, work on feature engineering, try out multiple machine learning models, select a model, work to optimize the model, and finally, inspect the outputs of the model and draw conclusions. __While this notebook may not get us to the top of the leaderboard, it is meant to be used as a teaching tool to give you a solid foundation to build on for future machine learning projects. Kaggle projects can teach us a lot about machine learning, but several of the strategies used to get to the very top of the leaderboard are not best practices, so here we'll stick to building a very good - although not quite first place - machine learning solution. While Kaggle projects are competitions, I think they are best described as \"a machine learning education\" disguised as a contest!\"\n\nIf you are looking to follow-up on this work, I have additional work including a [kernel on using Automated Feature Engineering](https://www.kaggle.com/willkoehrsen/featuretools-for-good) with [Featuretools](https://docs.featuretools.com/#minute-quick-start) for this problem (with slightly higher leaderboard score). (If you enjoy my writing style and explanations, I write for [Towards Data Science](http://medium.com/@williamkoehrsen/))\n\n## Problem and Data Explanation\n\nThe data for this competition is provided in two files: `train.csv` and `test.csv`. The training set has 9557 rows and 143 columns while the testing set has 23856 rows and 142 columns. Each row represents __one individual__ and each column is a __feature, either unique to the individual, or for the household of the individual__. The training set has one additional column, `Target`, which represents the poverty level on a 1-4 scale and is the label for the competition. A value of 1 is the most extreme poverty. \n\nThis is a __supervised multi-class classification machine learning problem__:\n\n* __Supervised__: provided with the labels for the training data\n* __Multi-class classification__: Labels are discrete values with 4 classes\n\n### Objective\n\nThe objective is to predict poverty on a __household level__. We are given data on the individual level with each individual having unique features but also information about their household. In order to create a dataset for the task, we'll have to perform some _aggregations of the individual data_ for each household. Moreover, we have to make a prediction for every individual in the test set, but _\"ONLY the heads of household are used in scoring\"_ which means we want to predict poverty on a household basis. \n\n__Important note: while all members of a household should have the same label in the training data, there are errors where individuals in the same household have different labels. In these cases, we are told to use the label for the head of each household, which can be identified by the rows where `parentesco1 == 1.0`.__ We will cover how to correct this in the notebook (for more info take a look at the [competition main discussion](https://www.kaggle.com/c/costa-rican-household-poverty-prediction/discussion/61403)).\n\nThe `Target` values represent poverty levels as follows:\n\n    1 = extreme poverty \n    2 = moderate poverty \n    3 = vulnerable households \n    4 = non vulnerable households\n\nThe explanations for all 143 columns can be found in the [competition documentation](https://www.kaggle.com/c/costa-rican-household-poverty-prediction/data), but a few to note are below:\n\n* __Id__: a unique identifier for each individual, this should not be a feature that we use! \n* __idhogar__: a unique identifier for each household. This variable is not a feature, but will be used to group individuals by household as all individuals in a household will have the same identifier.\n* __parentesco1__: indicates if this person is the head of the household.\n* __Target__: the label, which should be equal for all members in a household\n\nWhen we make a model, we'll train on a household basis with the label for each household _the poverty level of the head of household_. The raw data contains a mix of both household and individual characteristics and for the individual data, we will have to find a way to aggregate this for each household. Some of the individuals belong to a household with _no head of household_ which means that unfortunately we can't use this data for training. These issues with the data are completely typical of __real-world__ data and hence this problem is great preparation for the datasets you'll encounter in a data science job! \n\n### Metric\n\nUltimately we want to build a machine learning model that can predict the integer poverty level of a household. Our predictions will be assessed by the __Macro F1 Score.__ You may be familiar with the [standard F1 score](http://blog.exsilio.com/all/accuracy-precision-recall-f1-score-interpretation-of-performance-measures/) for binary classification problems which is the harmonic mean of precision and recall:\n\n$$F_1 = \\frac{2}{\\tfrac{1}{\\mathrm{recall}} + \\tfrac{1}{\\mathrm{precision}}} = 2 \\cdot \\frac{\\mathrm{precision} \\cdot \\mathrm{recall}}{\\mathrm{precision} + \\mathrm{recall}}$$\n\nFor mutli-class problems, we have to average the F1 scores for each class. The macro F1 score averages the F1 score for each class _without taking into account label imbalances_. \n\n$$\\text{Macro F1} = \\frac{\\text{F1 Class 1} + \\text{F1 Class 2} + \\text{F1 Class 3} + \\text{F1 Class 4}}{4}$$\n\nIn other words, the number of occurrences of each label does not figure into the calculation when using macro (while it does when using the \"weighted\" score). (For more information on the differences, look at the [Scikit-Learn Documention for F1 Score](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.f1_score.html#sklearn.metrics.f1_score) or this [Stack Exchange question and answers](https://datascience.stackexchange.com/q/15989/42908). If we want to assess our performance, we can use the code:\n\n```\nfrom sklearn.metrics import f1_score\nf1_score(y_true, y_predicted, average = 'macro`)\n```\n\nFor this problem, the labels are imbalanced, which makes it a little strange to use macro averaging for the evaluation metric, but that's a decision made by the organizers and not something we can change! In your own work, you want to be aware of label imbalances and choose a metric accordingly. \n\n## Roadmap\n\nThe end objective is a machine learning model that can predict the poverty level of a household. However, before we get carried away with modeling, it's important to understand the problem and data. Also, we want to evaluate numerous models before choosing one as the \"best\" and after building a model, we want to investigate the predictions. Our roadmap is therefore as follows:\n\n1. Understand the problem (we're almost there already)\n2. Exploratory Data Analysis\n3. Feature engineering to create a dataset for machine learning\n4. Compare several baseline machine learning models\n5. Try more complex machine learning models\n6. Optimize the selected model\n7. Investigate model predictions in context of problem\n6. Draw conclusions and lay out next steps \n\nThe steps laid out above are iterative meaning that while we will go through them one at a time, we might go back to an earlier step and revisit some of our decisions. In general, data science is a non-linear pracice where we are constantly evaluating our past decisions and making improvements. In particular, feature engineering, modeling, and optimization are steps that we often repeat because we never know if we got them right the first time!\n\n# Getting Started\n\nWe have a pretty good grasp of the problem, so we'll move into the Exploratory Data Analysis (EDA) and feature engineering. For the EDA we'll examine any interesting anomalies, trends, correlations, or patterns that can be used for feature engineering and for modeling. We'll make sure to investigate our data both quantitatively (with statistics) and visually (with figures). \n\nOnce we have a good grasp of the data and any potentially useful relationships, we can do some feature engineering (the most important part of the machine learning pipeline) and establish a baseline model. This won't get us to the top of the leaderboard, but it will provide a strong foundation to build on! \n\nWith all that info in mind (don't worry if you haven't got all the details), let's get started! "},{"metadata":{"_uuid":"df5bf2d5d559636504db46b783c45e77d9bb8776"},"cell_type":"markdown","source":"#### Imports\n\nWe'll use a familiar stack of data science libraries: `Pandas`, `numpy`, `matplotlib`, `seaborn`, and eventually `sklearn` for modeling. "},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"# Data manipulation\nimport pandas as pd\nimport numpy as np\n\n# Visualization\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Set a few plotting defaults\n%matplotlib inline\nplt.style.use('fivethirtyeight')\nplt.rcParams['font.size'] = 18\nplt.rcParams['patch.edgecolor'] = 'k'","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9c541a166077ce34eadcde7a702b4e1f0d3bbba6"},"cell_type":"markdown","source":"### Read in Data and Look at Summary Information"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"pd.options.display.max_columns = 150\n\n# Read in data\ntrain = pd.read_csv('../input/train.csv')\ntest = pd.read_csv('../input/test.csv')\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8ae7f4f9dc01ad847283246327c7db467ba76498"},"cell_type":"markdown","source":"That gives us a look at all of the columns which don't appear to be in any order. To get a quick overview of the data we use `df.info()`."},{"metadata":{"_uuid":"b180722a8e9a1688463f0340e4830e37a5586c4a","trusted":true},"cell_type":"code","source":"train.info()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"35a4c874b82d853208f32d0d1dbaed422045f914"},"cell_type":"markdown","source":"This tells us there are 130 integer columns, 8 float (numeric) columns, and 5 object columns. The integer columns probably represent Boolean variables (that take on either 0 or 1) or [ordinal variables](https://www.ma.utexas.edu/users/mks/statmistakes/ordinal.html) with discrete ordered values. The object columns might pose an issue because they cannot be fed directly into a machine learning model.\n\nLet's glance at the test data which has many more rows (individuals) than the train. It does have one fewer column because there's no Target!"},{"metadata":{"_uuid":"61687f25325218b77d4d75804044c7083da85eb7","trusted":true},"cell_type":"code","source":"test.info()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3c0d0640f3867d2c290ac58a55639fefb08875c0"},"cell_type":"markdown","source":"#### Integer Columns\n\nLet's look at the distribution of unique values in the integer columns. For each column, we'll count the number of unique values and show the result in a bar plot."},{"metadata":{"_uuid":"7db394756abedd215885bc7d6b9dc2c767964b35","trusted":true},"cell_type":"code","source":"train.select_dtypes(np.int64).nunique().value_counts().sort_index().plot.bar(color = 'blue', \n                                                                             figsize = (8, 6),\n                                                                            edgecolor = 'k', linewidth = 2);\nplt.xlabel('Number of Unique Values'); plt.ylabel('Count');\nplt.title('Count of Unique Values in Integer Columns');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1a9db9a6caa3b77b6a6f0184644c1e11b2ed53d5"},"cell_type":"markdown","source":"The columns with only 2 unique values represent Booleans (0 or 1). In a lot of cases, this boolean information is already on a household level. For example, the `refrig` column says whether or not the household has a refrigerator. When it comes time to make features from the Boolean columns that are on the household level, we will _not need to aggregate_ these. However, the Boolean columns that are on the individual level will need to be aggregated. \n\n#### Float Columns\n\nAnother column type is floats which represent continuous variables. We can make a quick distribution plot to show the distribution of all float columns. We'll use an [`OrderedDict`](https://pymotw.com/2/collections/ordereddict.html) to map the poverty levels to colors because this keeps the keys and values in the same order as we specify (unlike a regular Python dictionary).\n\nThe following graphs shows the distributions of the `float` columns colored by the value of the `Target`. With these plots, we can see if there is a significant difference in the variable distribution depending on the household poverty level."},{"metadata":{"_uuid":"c22a1876cf1d8a717b0baab8564acb4a7f531b51","scrolled":false,"trusted":true},"cell_type":"code","source":"from collections import OrderedDict\n\nplt.figure(figsize = (20, 16))\nplt.style.use('fivethirtyeight')\n\n# Color mapping\ncolors = OrderedDict({1: 'red', 2: 'orange', 3: 'blue', 4: 'green'})\npoverty_mapping = OrderedDict({1: 'extreme', 2: 'moderate', 3: 'vulnerable', 4: 'non vulnerable'})\n\n# Iterate through the float columns\nfor i, col in enumerate(train.select_dtypes('float')):\n    ax = plt.subplot(4, 2, i + 1)\n    # Iterate through the poverty levels\n    for poverty_level, color in colors.items():\n        # Plot each poverty level as a separate line\n        sns.kdeplot(train.loc[train['Target'] == poverty_level, col].dropna(), \n                    ax = ax, color = color, label = poverty_mapping[poverty_level])\n        \n    plt.title(f'{col.capitalize()} Distribution'); plt.xlabel(f'{col}'); plt.ylabel('Density')\n\nplt.subplots_adjust(top = 2)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"79ec245f541336c564fdce3d013b0b988627b6e5"},"cell_type":"markdown","source":"Later on we'll calculate correlations between the variables and the `Target` to gauge the relationships between the features, but these plots can already give us a sense of which variables may be most \"relevant\" to a model. For example, the `meaneduc`, representing the average education of the adults in the household appears to be related to the poverty level: __a higher average adult education leads to higher values of the target which are less severe levels of poverty__. The theme of the importance of education is one we will come back to again and again in this notebook! "},{"metadata":{"_uuid":"355fd9746ec30af6661cf97a1243790730a5df99"},"cell_type":"markdown","source":"#### Object Columns\n\nThe last column type is `object` which we can view as follows."},{"metadata":{"_uuid":"9b110b76084ba14045c4c69ce909148d91ce0ecc","trusted":true},"cell_type":"code","source":"train.select_dtypes('object').head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d1e660b7e06b81b61bec81a82be4a5854115d09a"},"cell_type":"markdown","source":"The `Id` and `idhogar` object types make sense because these are identifying variables. However, the other columns seem to be a mix of strings and numbers which we'll need to address before doing any machine learning. According to the documentation for these columns:\n\n* `dependency`: Dependency rate, calculated = (number of members of the household younger than 19 or older than 64)/(number of member of household between 19 and 64)\n* `edjefe`: years of education of male head of household, based on the interaction of escolari (years of education), head of household and gender, yes=1 and no=0\n* `edjefa`: years of education of female head of household, based on the interaction of escolari (years of education), head of household and gender, yes=1 and no=0\n\nThese explanations clear up the issue. For these three variables, __\"yes\" = 1__ and __\"no\" = 0__. We can correct the variables using a mapping and convert to floats. "},{"metadata":{"_uuid":"acf4a663dd5d6713dad025ef6674b547cc02690f","trusted":true},"cell_type":"code","source":"mapping = {\"yes\": 1, \"no\": 0}\n\n# Apply same operation to both train and test\nfor df in [train, test]:\n    # Fill in the values with the correct mapping\n    df['dependency'] = df['dependency'].replace(mapping).astype(np.float64)\n    df['edjefa'] = df['edjefa'].replace(mapping).astype(np.float64)\n    df['edjefe'] = df['edjefe'].replace(mapping).astype(np.float64)\n\ntrain[['dependency', 'edjefa', 'edjefe']].describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"55af5988853fd538f674b07f93d5655c133ed744","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (16, 12))\n\n# Iterate through the float columns\nfor i, col in enumerate(['dependency', 'edjefa', 'edjefe']):\n    ax = plt.subplot(3, 1, i + 1)\n    # Iterate through the poverty levels\n    for poverty_level, color in colors.items():\n        # Plot each poverty level as a separate line\n        sns.kdeplot(train.loc[train['Target'] == poverty_level, col].dropna(), \n                    ax = ax, color = color, label = poverty_mapping[poverty_level])\n        \n    plt.title(f'{col.capitalize()} Distribution'); plt.xlabel(f'{col}'); plt.ylabel('Density')\n\nplt.subplots_adjust(top = 2)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e2fb233b0d5876e71ba33e6b5a8d29e2429fab94"},"cell_type":"markdown","source":"These variables are now correctly represented as numbers and can be fed into a machine learning model. \n\nTo make operations like that above a little easier, we'll join together the training and testing dataframes. This is important once we start feature engineering because we want to apply the same operations to both dataframes so we end up with the same features. Later we can separate out the sets based on the `Target`."},{"metadata":{"_uuid":"59a54206932e6414745226a101097b69b7845107","trusted":true},"cell_type":"code","source":"# Add null Target column to test\ntest['Target'] = np.nan\ndata = train.append(test, ignore_index = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"423ee29cdbb0975eb9d8f6f729e1595f401acbc3"},"cell_type":"markdown","source":"## Exploring Label Distribution\n\nNext, we can get an idea of how imbalanced the problem is by looking at the distribution of labels. There are four possible integer levels, indicating four different levels of poverty. To look at the correct labels, we'll subset only to the columns where `parentesco1 == 1` because this is the head of household, the correct label for each household.\n\nThe bar plot below shows the distribution of training labels (since there are no testing labels)."},{"metadata":{"_uuid":"6ed7b957df6b758a5bf586c682289edcb1a78f66","trusted":true},"cell_type":"code","source":"# Heads of household\nheads = data.loc[data['parentesco1'] == 1].copy()\n\n# Labels for training\ntrain_labels = data.loc[(data['Target'].notnull()) & (data['parentesco1'] == 1), ['Target', 'idhogar']]\n\n# Value counts of target\nlabel_counts = train_labels['Target'].value_counts().sort_index()\n\n# Bar plot of occurrences of each label\nlabel_counts.plot.bar(figsize = (8, 6), \n                      color = colors.values(),\n                      edgecolor = 'k', linewidth = 2)\n\n# Formatting\nplt.xlabel('Poverty Level'); plt.ylabel('Count'); \nplt.xticks([x - 1 for x in poverty_mapping.keys()], \n           list(poverty_mapping.values()), rotation = 60)\nplt.title('Poverty Level Breakdown');\n\nlabel_counts","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2876c882d7be2bc7ab58bb0dd0dc09c09bcc65a5"},"cell_type":"markdown","source":"We are dealing with an imbalanced class problem (which makes it intriguing why the contest organizers choose the _macro_ F1 score as the metric instead of _weighted_ F1!). There are many more households that classify as _non vulnerable_ than in any other category. The _extreme_ poverty class is the smallest (I guess this should make us optimistic!).\n\nOne problem with imbalanced classification problems is that the machine learning model can have a difficult time predicting the minority classes because it sees far less examples. Think about this in human terms: if we are classifiying poverty and we see far more cases of no poverty than extreme poverty, it will make it more difficult for us to identify the high poverty households because of less exposure. One potential method to address class imbalanceds is through oversampling  (which is covered in more advanced notebooks)."},{"metadata":{"_uuid":"4b4e7cebd091b34469e78c123eac55410dfb8e77"},"cell_type":"markdown","source":"## Addressing Wrong Labels\n\nAs with any realistic dataset, the Costa Rican Poverty data has some issues. Typically, 80% of a data science project will be spent cleaning data and fixing anomalies/errors. These can be either human entry errors, measurement errors, or sometimes just extreme values that are correct but stand out. For this problem, some of the labels are not correct because _individuals in the same household have a different poverty level_. We're not told why this may be the case, but we are told to use the head of household as the true label. \n\nThat information makes our job much easier, but in a real-world problem, we would have to figure out the reason _Why_ the labels are wrong and how to address the issue on our own. This section fixes the issue with the labels although it is not strictly necessary: I kept it in the notebook just to show how we may deal with this issue.\n\n### Identify Errors\n\nFirst we need to find the errors before we can correct them. To find the households with different labels for  family members, we can group the data by the household and then check if there is only one unique value of the `Target`."},{"metadata":{"_uuid":"81b7d45bfe8fd0feb09c307c5f2c9f00a6f8bc73","trusted":true},"cell_type":"code","source":"# Groupby the household and figure out the number of unique values\nall_equal = train.groupby('idhogar')['Target'].apply(lambda x: x.nunique() == 1)\n\n# Households where targets are not all equal\nnot_equal = all_equal[all_equal != True]\nprint('There are {} households where the family members do not all have the same target.'.format(len(not_equal)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25af32cb8437379e65c2100e4b1d2f26ef5f0410"},"cell_type":"markdown","source":"Let's look at one example."},{"metadata":{"_uuid":"c3c17e4f9a62990e08d9704d51917b0ace16d54f","trusted":true},"cell_type":"code","source":"train[train['idhogar'] == not_equal.index[0]][['idhogar', 'parentesco1', 'Target']]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3c934f81190e614dd6910c75828bc65d71e08233"},"cell_type":"markdown","source":"The organizers tell us that the correct label is that for the head of household, where `parentesco1 == 1`. For this household, the correct label is __3__ for all members. We can correct this (as shown later) by reassigning all the individuals in this household the correct poverty level. In the real-world, you might have to make the tough decision of how to address the problem by yourself (or with the help of your team)."},{"metadata":{"_uuid":"5f3bfc8ada58df621ab0d9bb789a161e9db834b9"},"cell_type":"markdown","source":"### Families without Heads of Household\n\nWe can correct all the label discrepancies by assigning the individuals in the same household the label of the head of household. But wait, you may ask: \"What if there are households without a head of household? And what if the members of those households have differing values of the label?\" \n\nWell, since you asked, let's investigate exactly that question! "},{"metadata":{"_uuid":"32c355c2537617ddb4806d4b2b1300817d5b7d99","trusted":true},"cell_type":"code","source":"households_leader = train.groupby('idhogar')['parentesco1'].sum()\n\n# Find households without a head\nhouseholds_no_head = train.loc[train['idhogar'].isin(households_leader[households_leader == 0].index), :]\n\nprint('There are {} households without a head.'.format(households_no_head['idhogar'].nunique()))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"317ab2763756312d66fcbe92d8d320878e560c1b","trusted":true},"cell_type":"code","source":"# Find households without a head and where labels are different\nhouseholds_no_head_equal = households_no_head.groupby('idhogar')['Target'].apply(lambda x: x.nunique() == 1)\nprint('{} Households with no head have different labels.'.format(sum(households_no_head_equal == False)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"88929859b4ba55bed14c3ec95f0bc7e60aa5c5c7"},"cell_type":"markdown","source":"Well that's a relief! This means that we don't have to worry about a household both where there is no head __AND__ the members have different values of the label! For this problem, according to the organizers, __if a household does not have a head, then there is no true label. Therefore, we actually won't use any of the households without a head for training__ Nonetheless, it's still a good exercise to go through this process of investigating the data! \n\n### Correct Errors\n\nNow we can correct labels for the households that do have a head __AND__ the members have different poverty levels."},{"metadata":{"_uuid":"8b3a4200ee6425d54d6321503fdd14bc7c8f2756","trusted":true},"cell_type":"code","source":"# Iterate through each household\nfor household in not_equal.index:\n    # Find the correct label (for the head of household)\n    true_target = int(train[(train['idhogar'] == household) & (train['parentesco1'] == 1.0)]['Target'])\n    \n    # Set the correct label for all members in the household\n    train.loc[train['idhogar'] == household, 'Target'] = true_target\n    \n    \n# Groupby the household and figure out the number of unique values\nall_equal = train.groupby('idhogar')['Target'].apply(lambda x: x.nunique() == 1)\n\n# Households where targets are not all equal\nnot_equal = all_equal[all_equal != True]\nprint('There are {} households where the family members do not all have the same target.'.format(len(not_equal)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d1eb1771bc629d6f02452dc4e3edf669bffd1b3d"},"cell_type":"markdown","source":"Since we are only going to use the heads of household for the labels, __this step is not completely necessary but it shows a workflow for correcting data errors like you may encounter in real life__. Don't consider it extra work, just practice for your career! "},{"metadata":{"_uuid":"5c309e0206872b450409ec6bc78e74f49bcaf206"},"cell_type":"markdown","source":"## Missing Variables\n\nOne of the most important steps of exploratory data analysis is finding missing values in the data and determining how to handle them. Missing values have to be filled in before we use a machine learning model and we need to think of the best strategy for filling them in based on the feature: this is where we'll have to start digging into the data definitions.\n\nFirst we can look at the percentage of missing values in each column."},{"metadata":{"_uuid":"47dd23a8222bd1d89512662dcb299b2ee6b6b04e","trusted":true},"cell_type":"code","source":"# Number of missing in each column\nmissing = pd.DataFrame(data.isnull().sum()).rename(columns = {0: 'total'})\n\n# Create a percentage missing\nmissing['percent'] = missing['total'] / len(data)\n\nmissing.sort_values('percent', ascending = False).head(10).drop('Target')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7134602d0d3cd3045fb97eaaa724592b2d23553a"},"cell_type":"markdown","source":"We don't have to worry about the `Target` becuase we made that `NaN` for the test data. However, we do need to address the other 3 columns with a high percentage of missing values.\n\n__v18q1__: Number of tablets\n\nLet's start with `v18q1` which indicates the number of tablets owned by a family. We can look at the value counts of this variable. Since this is a household variable, it only makes sense to look at it on a household level, so we'll only select the rows for the head of household.\n\n#### Function to Plot Value Counts\n\nSince we might want to plot value counts for different columns, we can write a simple function that will do it for us!"},{"metadata":{"_uuid":"1eadbc5c6266518b2f0a047449966662ef8ea538","trusted":true},"cell_type":"code","source":"def plot_value_counts(df, col, heads_only = False):\n    \"\"\"Plot value counts of a column, optionally with only the heads of a household\"\"\"\n    # Select heads of household\n    if heads_only:\n        df = df.loc[df['parentesco1'] == 1].copy()\n        \n    plt.figure(figsize = (8, 6))\n    df[col].value_counts().sort_index().plot.bar(color = 'blue',\n                                                 edgecolor = 'k',\n                                                 linewidth = 2)\n    plt.xlabel(f'{col}'); plt.title(f'{col} Value Counts'); plt.ylabel('Count')\n    plt.show();","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7d5df5fccf64cd3ccecbb1c1bf2c2133a77c31e5","trusted":true},"cell_type":"code","source":"plot_value_counts(heads, 'v18q1')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cd53530fcdb192d6f7b6512c33ca134cf7984bed"},"cell_type":"markdown","source":"It looks like the most common number of tablets to own is 1 if we go only by the data that is present. However, we also need to think about the data that is missing. In this case, it could be that families with a `nan` in this category just do not own a tablet! If we look at the data definitions, we see that `v18q` indicates whether or not a family owns a tablet. We should investigate this column combined with the number of tablets to see if our hypothesis holds.\n\nWe can `groupby` the value of `v18q` (which is 1 for owns a tablet and 0 for does not) and then calculate the number of null values for `v18q1`. This will tell us if the null values represent that the family does not own a tablet."},{"metadata":{"_uuid":"19c44193187c0db95649a60447b6e0fe43703795","trusted":true},"cell_type":"code","source":"heads.groupby('v18q')['v18q1'].apply(lambda x: x.isnull().sum())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"adff3153e886d7adc1ced34c157a747fd0c4e9a1"},"cell_type":"markdown","source":"Well, that solves the issue! Every family that has `nan` for `v18q1` does not own a tablet. Therefore, we can fill in this missing value with zero."},{"metadata":{"_uuid":"20291bd2766b006bca3255034c4d1432e8024b0d","trusted":true},"cell_type":"code","source":"data['v18q1'] = data['v18q1'].fillna(0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"42662b00e496efdd18cf3fc9cf497604cf873896"},"cell_type":"markdown","source":"__v2a1__: Monthly rent payment\n\nThe next missing column is `v2a1` which represents the montly rent payment. "},{"metadata":{"_uuid":"a425dc1a313ce8b75f64d4aa034513b84694c6ba"},"cell_type":"markdown","source":"In addition to looking at the missing values of the monthly rent payment, it will be interesting to also look at the distribution of `tipovivi_`, the columns showing the ownership/renting status of the home. For this plot, we show the ownership status of those homes with a `nan` for the monthyl rent payment."},{"metadata":{"_uuid":"136b2e5c08622f91d56a275d278e9a4731a7396f","trusted":true},"cell_type":"code","source":"# Variables indicating home ownership\nown_variables = [x for x in data if x.startswith('tipo')]\n\n\n# Plot of the home ownership variables for home missing rent payments\ndata.loc[data['v2a1'].isnull(), own_variables].sum().plot.bar(figsize = (10, 8),\n                                                                        color = 'green',\n                                                              edgecolor = 'k', linewidth = 2);\nplt.xticks([0, 1, 2, 3, 4],\n           ['Owns and Paid Off', 'Owns and Paying', 'Rented', 'Precarious', 'Other'],\n          rotation = 60)\nplt.title('Home Ownership Status for Households Missing Rent Payments', size = 18);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d072b3df8b6486cb90c1e9832094c7e13680a856"},"cell_type":"markdown","source":"The meaning of the home ownership variables is below:\n\n    tipovivi1, =1 own and fully paid house\n    tipovivi2, \"=1 own,  paying in installments\"\n    tipovivi3, =1 rented\n    tipovivi4, =1 precarious\n    tipovivi5, \"=1 other(assigned,  borrowed)\"\n    \nWe've solved the issue! Well, mostly: the households that do not have a monthly rent payment generally own their own home. In a few other situations, we are not sure of the reason for the missing information. \n\nFor the houses that are owned and have a missing monthly rent payment, we can set the value of the rent payment to zero. For the other homes, we can leave the missing values to be imputed but we'll add a flag (Boolean) column indicating that these households had missing values."},{"metadata":{"_uuid":"6ecfd2a5e15efa4ffe5436f43e84dbb5ea929314","trusted":true},"cell_type":"code","source":"# Fill in households that own the house with 0 rent payment\ndata.loc[(data['tipovivi1'] == 1), 'v2a1'] = 0\n\n# Create missing rent payment column\ndata['v2a1-missing'] = data['v2a1'].isnull()\n\ndata['v2a1-missing'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2e147ff8825964c53c634b05ad6653ebf8390539"},"cell_type":"markdown","source":"__rez_esc__: years behind in school\n\nThe last column with a high percentage of missing values is `rez_esc` indicating years behind in school. For the families with a null value, is possible that they have no children currently in school. Let's test this out by finding the ages of those who have a missing value in this column and the ages of those who do not have a missing value."},{"metadata":{"_uuid":"d9576b2de07fbab360435f26d3d869f4f3b8a269","trusted":true},"cell_type":"code","source":"data.loc[data['rez_esc'].notnull()]['age'].describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b7f67c2d888db26770ff3c0e108d003b34137d4a"},"cell_type":"markdown","source":"What this tells us is that the oldest age with a missing value is 17. For anyone older than this, maybe we can assume that they are simply not in school. Let's look at the ages of those who have a missing value. "},{"metadata":{"_uuid":"b8b91108b4f7029b0642dd68052cc4b6d5906316","trusted":true},"cell_type":"code","source":"data.loc[data['rez_esc'].isnull()]['age'].describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fdc641191009f43a7d3407f46c41dc89bb0fc551"},"cell_type":"markdown","source":"If we read through some of the [discussions for this competition](https://www.kaggle.com/c/costa-rican-household-poverty-prediction/discussion/61403), we learn that this variable is only defined for individuals between 7 and 19. Anyone younger or older than this range presumably has no years behind and therefore the value should be set to 0. \nFor this variable, if the individual is over 19 and they have a missing value, or if they are younger than 7 and have a missing value we can set it to zero. For anyone else, we'll leave the value to be imputed and add a boolean flag."},{"metadata":{"_uuid":"791b94a06ade9896fd5a7d6dc38152f831642a6a","trusted":true},"cell_type":"code","source":"# If individual is over 19 or younger than 7 and missing years behind, set it to 0\ndata.loc[((data['age'] > 19) | (data['age'] < 7)) & (data['rez_esc'].isnull()), 'rez_esc'] = 0\n\n# Add a flag for those between 7 and 19 with a missing value\ndata['rez_esc-missing'] = data['rez_esc'].isnull()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bce88b0ebef5a8d3eeb3ca4f978a0d4a69b1e185"},"cell_type":"markdown","source":"There is also one outlier in the `rez_esc` column. Again, if we read through the competition discussions, we learn that the maximum value for this variable is 5. Therefore, any values above 5 should be set to 5."},{"metadata":{"_uuid":"18a972cc48d9e4818777afb7b06eab7e1f97855a","trusted":true},"cell_type":"code","source":"data.loc[data['rez_esc'] > 5, 'rez_esc'] = 5","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b3163f5a385c3f7d8fb7e3d38ea93526bf73e788"},"cell_type":"markdown","source":"## Plot Two Categorical Variables\n\nTo show how two categorical variables interact with one another, there are a number of plotting options: scatterplots, faceted bar plots, boxplots, etc. I wasn't satisfied with any of these choices so I wrote the function below, which essentially is a scatterplot of two categoricals where the size of the points represent the percentage of a given y-value represented by each x-value. "},{"metadata":{"_uuid":"bc5a09a7c9d504394785aa144a14547c31d3ab2f","trusted":true},"cell_type":"code","source":"def plot_categoricals(x, y, data, annotate = True):\n    \"\"\"Plot counts of two categoricals.\n    Size is raw count for each grouping.\n    Percentages are for a given value of y.\"\"\"\n    \n    # Raw counts \n    raw_counts = pd.DataFrame(data.groupby(y)[x].value_counts(normalize = False))\n    raw_counts = raw_counts.rename(columns = {x: 'raw_count'})\n    \n    # Calculate counts for each group of x and y\n    counts = pd.DataFrame(data.groupby(y)[x].value_counts(normalize = True))\n    \n    # Rename the column and reset the index\n    counts = counts.rename(columns = {x: 'normalized_count'}).reset_index()\n    counts['percent'] = 100 * counts['normalized_count']\n    \n    # Add the raw count\n    counts['raw_count'] = list(raw_counts['raw_count'])\n    \n    plt.figure(figsize = (14, 10))\n    # Scatter plot sized by percent\n    plt.scatter(counts[x], counts[y], edgecolor = 'k', color = 'lightgreen',\n                s = 100 * np.sqrt(counts['raw_count']), marker = 'o',\n                alpha = 0.6, linewidth = 1.5)\n    \n    if annotate:\n        # Annotate the plot with text\n        for i, row in counts.iterrows():\n            # Put text with appropriate offsets\n            plt.annotate(xy = (row[x] - (1 / counts[x].nunique()), \n                               row[y] - (0.15 / counts[y].nunique())),\n                         color = 'navy',\n                         s = f\"{round(row['percent'], 1)}%\")\n        \n    # Set tick marks\n    plt.yticks(counts[y].unique())\n    plt.xticks(counts[x].unique())\n    \n    # Transform min and max to evenly space in square root domain\n    sqr_min = int(np.sqrt(raw_counts['raw_count'].min()))\n    sqr_max = int(np.sqrt(raw_counts['raw_count'].max()))\n    \n    # 5 sizes for legend\n    msizes = list(range(sqr_min, sqr_max,\n                        int(( sqr_max - sqr_min) / 5)))\n    markers = []\n    \n    # Markers for legend\n    for size in msizes:\n        markers.append(plt.scatter([], [], s = 100 * size, \n                                   label = f'{int(round(np.square(size) / 100) * 100)}', \n                                   color = 'lightgreen',\n                                   alpha = 0.6, edgecolor = 'k', linewidth = 1.5))\n        \n    # Legend and formatting\n    plt.legend(handles = markers, title = 'Counts',\n               labelspacing = 3, handletextpad = 2,\n               fontsize = 16,\n               loc = (1.10, 0.19))\n    \n    plt.annotate(f'* Size represents raw count while % is for a given y value.',\n                 xy = (0, 1), xycoords = 'figure points', size = 10)\n    \n    # Adjust axes limits\n    plt.xlim((counts[x].min() - (6 / counts[x].nunique()), \n              counts[x].max() + (6 / counts[x].nunique())))\n    plt.ylim((counts[y].min() - (4 / counts[y].nunique()), \n              counts[y].max() + (4 / counts[y].nunique())))\n    plt.grid(None)\n    plt.xlabel(f\"{x}\"); plt.ylabel(f\"{y}\"); plt.title(f\"{y} vs {x}\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a4dd0f5cb0567694952aeb24995cdf5d0f0a0cd3","trusted":true},"cell_type":"code","source":"plot_categoricals('rez_esc', 'Target', data);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f25ffc2c2a6c8a0c37217dacdccee3f3f7f5313f"},"cell_type":"markdown","source":"The size of the markers represents the raw count. To read the plot, choose a given y-value and then read across the row. For example, with a poverty level of 1, 93% of individuals have no years behind with a total count of around 800 individuals and about 0.4% of individuals are 5 years behind with about 50 total individuals in this category. This plot attempts to show both the overall counts and the within category proportion; it's not perfect , but I gave it a shot!"},{"metadata":{"_uuid":"0bebbe6d0673a374c7d7c0dcb029c4dd92687ba5","trusted":true},"cell_type":"code","source":"plot_categoricals('escolari', 'Target', data, annotate = False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6c3444e8f8baea05b29fe2490ee09bd8b255512f"},"cell_type":"markdown","source":"The remaining missing values in each column will be filled in, a process known as `Imputation`. There are several types of imputation commonly used, and one of the simplest and most effective methods is to fill in the missing values with the `median` of the column."},{"metadata":{"_uuid":"4929266f920c5fc1bdc9b5ae14fbd54454a99cb3"},"cell_type":"markdown","source":"As a final step with the missing values, we can plot the distribution of target for the case where either of these values are missing. "},{"metadata":{"_uuid":"aab147c98444a34aee0259bcd505d878db00c86f","trusted":true},"cell_type":"code","source":"plot_value_counts(data[(data['rez_esc-missing'] == 1)], \n                  'Target')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5ac0324614abd5ee24f5f32277ca08b9bcbd72e6"},"cell_type":"markdown","source":"The distribution here seems to match that for all the data at large."},{"metadata":{"_uuid":"d6e74fff4f86c84e961debf2cce2cf54c2c7e87b","trusted":true},"cell_type":"code","source":"plot_value_counts(data[(data['v2a1-missing'] == 1)], \n                  'Target')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a1b4f1d3869a4a216aaf1dee9c949d8b317c11f1"},"cell_type":"markdown","source":"This looks like it could be an indicator of more poverty given the higher prevalence of 2: moderate poverty. \n\n__This represents an important point__: sometimes the missing information is just as important as the information you are given. "},{"metadata":{"_uuid":"989a4051c0166d6d910ff526fc30f50a3b58f969"},"cell_type":"markdown","source":"# Feature Engineering\n\nThere is plenty more exploratory data analysis we can do, but first we should work on consolidating our data at a household level. We already have some of the information for each household, but for training, we will need _all_ of the information summarized for each household. This means grouping the individuals in a house (`groupby`) and performing an aggregation (`agg`) of the individual variables. \n\nIn another notebook, I show how we can use automated feature engineering to do this, and __automated feature engineering__ should be a standard part of the machine learning workflow. Right now, we'll stick to doing this by hand, but definitely take a look at __automated feature engineering in Featuretools__."},{"metadata":{"_uuid":"81fa4da4a8733fdfdd756e0a946ffb8a799fc51c"},"cell_type":"markdown","source":"## Column Definitions\n\nSometimes in data science we have to get our hands dirty digging through the data or do tedious tasks that take a lot of time. This is that part of the analysis: we have to define the columns that are at an individual level and at a household level using the [data decsriptions](https://www.kaggle.com/c/costa-rican-household-poverty-prediction/data). There is simply no other way to identify which variables at are the household level than to go through the variables themselves in the data description. Except, I've already done this for you, so all you have to do is copy and paste!\n\nWe'll define different variables because we need to treat some of them in a different manner. Once we have the variables defined on each level, we can work to start aggregating them as needed.\n\nThe process is as follows\n\n1. Break variables into household level and invididual level\n2. Find suitable aggregations for the individual level data\n    * Ordinal variables can use statistical aggregations\n    * Boolean variables can also be aggregated but with fewer stats\n3. Join the individual aggregations to the household level data\n\n### Define Variable Categories\n\nThere are several different categories of variables:\n\n1. Individual Variables: these are characteristics of each individual rather than the household\n    * Boolean: Yes or No (0 or 1)\n    * Ordered Discrete: Integers with an ordering\n2. Household variables\n    * Boolean: Yes or No\n    * Ordered Discrete: Integers with an ordering\n    * Continuous numeric\n3. Squared Variables: derived from squaring variables in the data\n4. Id variables: identifies the data and should not be used as features\n\nBelow we manually define the variables in each category. This is a little tedious, but also necessary. "},{"metadata":{"_uuid":"9e5a3b03ca064e4dc055e785730d1a43af037277","trusted":true},"cell_type":"code","source":"id_ = ['Id', 'idhogar', 'Target']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dc6f4e6c5f343e917d65a1d8e6c1043236bcdce1","trusted":true},"cell_type":"code","source":"ind_bool = ['v18q', 'dis', 'male', 'female', 'estadocivil1', 'estadocivil2', 'estadocivil3', \n            'estadocivil4', 'estadocivil5', 'estadocivil6', 'estadocivil7', \n            'parentesco1', 'parentesco2',  'parentesco3', 'parentesco4', 'parentesco5', \n            'parentesco6', 'parentesco7', 'parentesco8',  'parentesco9', 'parentesco10', \n            'parentesco11', 'parentesco12', 'instlevel1', 'instlevel2', 'instlevel3', \n            'instlevel4', 'instlevel5', 'instlevel6', 'instlevel7', 'instlevel8', \n            'instlevel9', 'mobilephone', 'rez_esc-missing']\n\nind_ordered = ['rez_esc', 'escolari', 'age']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b382ddbea1cd6ff25b71f13b9263a6565ebadd57","trusted":true},"cell_type":"code","source":"hh_bool = ['hacdor', 'hacapo', 'v14a', 'refrig', 'paredblolad', 'paredzocalo', \n           'paredpreb','pisocemento', 'pareddes', 'paredmad',\n           'paredzinc', 'paredfibras', 'paredother', 'pisomoscer', 'pisoother', \n           'pisonatur', 'pisonotiene', 'pisomadera',\n           'techozinc', 'techoentrepiso', 'techocane', 'techootro', 'cielorazo', \n           'abastaguadentro', 'abastaguafuera', 'abastaguano',\n            'public', 'planpri', 'noelec', 'coopele', 'sanitario1', \n           'sanitario2', 'sanitario3', 'sanitario5',   'sanitario6',\n           'energcocinar1', 'energcocinar2', 'energcocinar3', 'energcocinar4', \n           'elimbasu1', 'elimbasu2', 'elimbasu3', 'elimbasu4', \n           'elimbasu5', 'elimbasu6', 'epared1', 'epared2', 'epared3',\n           'etecho1', 'etecho2', 'etecho3', 'eviv1', 'eviv2', 'eviv3', \n           'tipovivi1', 'tipovivi2', 'tipovivi3', 'tipovivi4', 'tipovivi5', \n           'computer', 'television', 'lugar1', 'lugar2', 'lugar3',\n           'lugar4', 'lugar5', 'lugar6', 'area1', 'area2', 'v2a1-missing']\n\nhh_ordered = [ 'rooms', 'r4h1', 'r4h2', 'r4h3', 'r4m1','r4m2','r4m3', 'r4t1',  'r4t2', \n              'r4t3', 'v18q1', 'tamhog','tamviv','hhsize','hogar_nin',\n              'hogar_adul','hogar_mayor','hogar_total',  'bedrooms', 'qmobilephone']\n\nhh_cont = ['v2a1', 'dependency', 'edjefe', 'edjefa', 'meaneduc', 'overcrowding']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2771c641c24bb3ed21eb655320a4ad85066a2cd6","trusted":true},"cell_type":"code","source":"sqr_ = ['SQBescolari', 'SQBage', 'SQBhogar_total', 'SQBedjefe', \n        'SQBhogar_nin', 'SQBovercrowding', 'SQBdependency', 'SQBmeaned', 'agesq']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0f6effa461740a3436b7af6c3e866014fe7742f6"},"cell_type":"markdown","source":"Let's make sure we covered all of the variables and didn't repeat any. "},{"metadata":{"_uuid":"8aa08338c60ae323e4e4d6654f1c4380befd5c76","trusted":true},"cell_type":"code","source":"x = ind_bool + ind_ordered + id_ + hh_bool + hh_ordered + hh_cont + sqr_\n\nfrom collections import Counter\n\nprint('There are no repeats: ', np.all(np.array(list(Counter(x).values())) == 1))\nprint('We covered every variable: ', len(x) == data.shape[1])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"39d839a711f0783285ff0ca570366f045eb995f1"},"cell_type":"markdown","source":"#### Squared Variables\n\nFirst, the easiest step: we'll remove all of the squared variables. Sometimes variables are squared or transformed as part of feature engineering because it can help linear models learn relationships that are non-linear. However, since we will be using more complex models, these squared features are redundant. They are highly correlated with the non-squared version, and hence can actually hurt our model by adding irrelevant information and also slowing down training.\n\nFor an example, let's take a look at `SQBage` vs `age`."},{"metadata":{"_uuid":"5c89fc990cf416705e577fda9c086900f2c3abef","trusted":true},"cell_type":"code","source":"sns.lmplot('age', 'SQBage', data = data, fit_reg=False);\nplt.title('Squared Age versus Age');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"489d3c75283d701f3ea74bc14db815f212b546bb"},"cell_type":"markdown","source":"These variables are highly correlated, and we don't need to keep both in our data."},{"metadata":{"_uuid":"b258b0875757bc29016c020164b21021b63a7bcc","trusted":true},"cell_type":"code","source":"# Remove squared variables\ndata = data.drop(columns = sqr_)\ndata.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"68b5104b73966a357e6c278fccb562f52ba29863"},"cell_type":"markdown","source":"## Id Variables\n\nThese are pretty simple: they will be kept as is in the data since we need them for identification.\n\n## Household Level Variables\n\nFirst let's subset to the heads of household and then to the household level variables."},{"metadata":{"_uuid":"f3c177da83980642540c763666f6859d59d506aa","trusted":true},"cell_type":"code","source":"heads = data.loc[data['parentesco1'] == 1, :]\nheads = heads[id_ + hh_bool + hh_cont + hh_ordered]\nheads.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"730b9df91ea83a5fce8097745fbc66e768c5d321"},"cell_type":"markdown","source":"For most of the household level variables, we can simply keep them as is: since we want to make predictions for each household, we use these variables as features. However, we can also remove some redundant variables and also add in some more features derived from existing data. "},{"metadata":{"_uuid":"ff1b80ef402231ab5191a5f65e2c47ae7bcaf696"},"cell_type":"markdown","source":"### Redundant Household Variables\n\nLet's take a look at the correlations between all of the household variables. If there are any that are too highly correlated, then we might want to remove one of the pair of highly correlated variables.\n\nThe following code identifies any variables with a greater than 0.95 absolute magnitude correlation."},{"metadata":{"_uuid":"ff9ed5e6364dadd624f5717fbc2b3ee37195e6f3","trusted":true},"cell_type":"code","source":"# Create correlation matrix\ncorr_matrix = heads.corr()\n\n# Select upper triangle of correlation matrix\nupper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))\n\n# Find index of feature columns with correlation greater than 0.95\nto_drop = [column for column in upper.columns if any(abs(upper[column]) > 0.95)]\n\nto_drop","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1ae294fb76f5bcaaaf61de9cc0abffefaf881cdf"},"cell_type":"markdown","source":"These show one out of each pair of correlated variables. To find the other pair, we can subset the `corr_matrix`. "},{"metadata":{"_uuid":"a4c8970c265d7eff827abac3356ddd1d12b5572a","trusted":true},"cell_type":"code","source":"corr_matrix.loc[corr_matrix['tamhog'].abs() > 0.9, corr_matrix['tamhog'].abs() > 0.9]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ae006c90765293d8a09f4fea910420cf77dde8db"},"cell_type":"code","source":"sns.heatmap(corr_matrix.loc[corr_matrix['tamhog'].abs() > 0.9, corr_matrix['tamhog'].abs() > 0.9],\n            annot=True, cmap = plt.cm.autumn_r, fmt='.3f');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"76b4979c18c990a853e907de0d5e68982eaabdc5"},"cell_type":"markdown","source":"There are several variables here having to do with the size of the house:\n\n* r4t3, Total persons in the household\n* tamhog, size of the household\n* tamviv, number of persons living in the household\n* hhsize, household size\n* hogar_total, # of total individuals in the household\n\nThese variables are all highly correlated with one another. In fact, `hhsize` has a perfect correlation with `tamhog` and `hogar_total`. We will remove these two variables because the information is redundant. We can also remove `r4t3` because it has a near perfect correlation with `hhsize`.\n\n`tamviv` is not necessarily the same as `hhsize` because there might be family members that are not living in the household. Let's visualize this difference in a scatterplot."},{"metadata":{"_uuid":"b3517e3d0f11f886d82f2fcc7b07c3b03e455c8d","trusted":true},"cell_type":"code","source":"heads = heads.drop(columns = ['tamhog', 'hogar_total', 'r4t3'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cd03278e9ae4244036f60a223f9fb7f16da18b02","trusted":true},"cell_type":"code","source":"sns.lmplot('tamviv', 'hhsize', data, fit_reg=False, size = 8);\nplt.title('Household size vs number of persons living in the household');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"907c00eb54bf710b7edf1a2729e6afd8bf4c276c"},"cell_type":"markdown","source":"We see for a number of cases, there are more people living in the household than there are in the family. This gives us a good idea for a new feature: __the difference between these two measurements!__\n\nLet's make this new feature."},{"metadata":{"_uuid":"2cbebfb6c3274e08bb76feb9bd311f79f1f4b2fd","trusted":true},"cell_type":"code","source":"heads['hhsize-diff'] = heads['tamviv'] - heads['hhsize']\nplot_categoricals('hhsize-diff', 'Target', heads)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a46229591b686ffbc89d3c3e6308b4505bb2ce85"},"cell_type":"markdown","source":"Even though most households do not have a difference, there are a few that have more people living in the household than are members of the household.\n\nLet's move on to the other redundant variables. First we can look at `coopele`"},{"metadata":{"_uuid":"27ce4ffc97c6a745a85545cc2bd9fb4a5f0d4071","trusted":true},"cell_type":"code","source":"corr_matrix.loc[corr_matrix['coopele'].abs() > 0.9, corr_matrix['coopele'].abs() > 0.9]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"562839a941bac3265a5d0d4ce9eccc12e416bf0e"},"cell_type":"markdown","source":"These variables indicate where the electricity in the home is coming from. There are four options, and the families that don't have one of these two options either have no electricity (`noelec`) or get it from a private plant (`planpri`). \n\n#### Creating Ordinal Variable\n\nI'm going to compress these four variables into one by creating an ordinal variable. I'm going to choose the mapping myself, based on the data decriptions: \n\n    0: No electricity\n    1: Electricity from cooperative\n    2: Electricity from CNFL, ICA, ESPH/JASEC\n    3: Electricity from private plant\n\nAn ordered variable has an inherent ordering, and for this we choose our own based on the domain knowledge. After we create this new ordered variable, we can drop the four others. There are several households that do not have a variable here, so we will use a `nan` (which will be filled in during imputation) and add a Boolean column indicating there was no measure for this variable."},{"metadata":{"_uuid":"a9bd357a09eeccc0f0c479651470331f47af2648","trusted":true},"cell_type":"code","source":"elec = []\n\n# Assign values\nfor i, row in heads.iterrows():\n    if row['noelec'] == 1:\n        elec.append(0)\n    elif row['coopele'] == 1:\n        elec.append(1)\n    elif row['public'] == 1:\n        elec.append(2)\n    elif row['planpri'] == 1:\n        elec.append(3)\n    else:\n        elec.append(np.nan)\n        \n# Record the new variable and missing flag\nheads['elec'] = elec\nheads['elec-missing'] = heads['elec'].isnull()\n\n# Remove the electricity columns\n# heads = heads.drop(columns = ['noelec', 'coopele', 'public', 'planpri'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5d97954ff9c3420f9b44a8865e05fc473a17a1a5","trusted":true},"cell_type":"code","source":"plot_categoricals('elec', 'Target', heads)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0276d599b977349915d0ac88e02d25b3cffd5357"},"cell_type":"markdown","source":"We can see that for every value of the Target, the most common source of electricity is from one of the listed providers. "},{"metadata":{"_uuid":"5a2f14086e364572f84a05d21daad869cabe5378"},"cell_type":"markdown","source":"The final redundant column is `area2`. This means the house is in a rural zone, but it's redundant because we have a column indicating if the house is in a urban zone. Therefore, we can drop this column."},{"metadata":{"_uuid":"8417d12ca211bd0aab75ba94d0215161000ccf81","trusted":true},"cell_type":"code","source":"heads = heads.drop(columns = 'area2')\n\nheads.groupby('area1')['Target'].value_counts(normalize = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4400e7ac30ef507308ea9c7a24581eff0ad570ef"},"cell_type":"markdown","source":"It seems like households in an urban area (value of 1) are more likely to have lower poverty levels than households in a rural area (value of 0)."},{"metadata":{"_uuid":"0e7aa155b0c310c6ec4ac25c88c90af5e007e247"},"cell_type":"markdown","source":"### Creating Ordinal Variables\n\nFor the walls, roof, and floor of the house, there are three columns each: the first indicating 'bad', the second 'regular', and the third 'good'. We could leave the variables as booleans, but to me it makes more sense to turn them into ordinal variables because there is an inherent order: `bad < regular < good`. To do this, we can simply find whichever column is non-zero for each household using `np.argmax`. \n\nOnce we have created the ordinal variable, we are able to drop the original variables."},{"metadata":{"_uuid":"8c776cec90c7f4829eda309d7c925e6c751008ab","trusted":true},"cell_type":"code","source":"# Wall ordinal variable\nheads['walls'] = np.argmax(np.array(heads[['epared1', 'epared2', 'epared3']]),\n                           axis = 1)\n\n# heads = heads.drop(columns = ['epared1', 'epared2', 'epared3'])\nplot_categoricals('walls', 'Target', heads)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5df42c066b57cac634891f08570f33811961c42a","trusted":true},"cell_type":"code","source":"# Roof ordinal variable\nheads['roof'] = np.argmax(np.array(heads[['etecho1', 'etecho2', 'etecho3']]),\n                           axis = 1)\nheads = heads.drop(columns = ['etecho1', 'etecho2', 'etecho3'])\n\n# Floor ordinal variable\nheads['floor'] = np.argmax(np.array(heads[['eviv1', 'eviv2', 'eviv3']]),\n                           axis = 1)\n# heads = heads.drop(columns = ['eviv1', 'eviv2', 'eviv3'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"747de840f4fe3374f8ef95637b9948a1be5a737b"},"cell_type":"markdown","source":"## Feature Construction\n\nIn addition to mapping variables to ordinal features, we can also create entirely new features from the existing data, known as feature construction. For example, we can add up the previous three features we just created to get an overall measure of the quality of the house's structure."},{"metadata":{"_uuid":"f698c5f8bb5dcff40f1711939b9cdefd060cf8b0","trusted":true},"cell_type":"code","source":"# Create new feature\nheads['walls+roof+floor'] = heads['walls'] + heads['roof'] + heads['floor']\n\nplot_categoricals('walls+roof+floor', 'Target', heads, annotate=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ce2862b49ece2347ee94cce38c8c6ba316dfeb1d"},"cell_type":"markdown","source":"This new feature may be useful because it seems like a Target of 4 (the lowest poverty level) tends to have higher values of the 'house quality' variable. We can also look at this in a table to get the fine-grained details."},{"metadata":{"_uuid":"a0e0570396accb529420c25eb42224363daedd54","trusted":true},"cell_type":"code","source":"counts = pd.DataFrame(heads.groupby(['walls+roof+floor'])['Target'].value_counts(normalize = True)).rename(columns = {'Target': 'Normalized Count'}).reset_index()\ncounts.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"038c8d7ad491a74016bd47cfdbb8ea568c318515"},"cell_type":"markdown","source":"The next variable will be a `warning` about the quality of the house. It will be a negative value, with -1 point each for no toilet, electricity, floor, water service, and ceiling."},{"metadata":{"_uuid":"1c69b16c9508d5c453fb3582a0731bee60116938","trusted":true},"cell_type":"code","source":"# No toilet, no electricity, no floor, no water service, no ceiling\nheads['warning'] = 1 * (heads['sanitario1'] + \n                         (heads['elec'] == 0) + \n                         heads['pisonotiene'] + \n                         heads['abastaguano'] + \n                         (heads['cielorazo'] == 0))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"40acec7571f26575f2001c08675ce88206ffb014"},"cell_type":"markdown","source":"We can keep using our `plot_categoricals` function to visualize these relationships, but `seaborn` also has a number of plotting options that can work with categoricals. One is the `violinplot` which shows the distribution of a variable on the y axis with the width of each plot showing the number of observations in that category. "},{"metadata":{"_uuid":"44e3169ba43c0a59828965e1a091abbacea15f72","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (10, 6))\nsns.violinplot(x = 'warning', y = 'Target', data = heads);\nplt.title('Target vs Warning Variable');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c8861533da1cdf69b1a2f0d08808b5c4cffd99ee","trusted":true},"cell_type":"code","source":"plot_categoricals('warning', 'Target', data = heads)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"15d1d11a7c70bbf2cedba08a6a8af28e8a699d88"},"cell_type":"markdown","source":"The violinplot is not great here because it smooths out the categorical variable with the effect that it looks as if the Target can take on lesser and greater values than in reality. Nonetheless, we can see a high concentration of households that have no warning signs and have the lowest level of poverty. It looks as if this may be a useful feature, but we can't know for sure until we get to modeling!"},{"metadata":{"_uuid":"9884e666ffb55adacc01931675170f39ab0886cf"},"cell_type":"markdown","source":"The final household feature we can make for now is a `bonus` where a family gets a point for having a refrigerator, computer, tablet, or television."},{"metadata":{"_uuid":"edf1cf9043391a9af16c45e19e84bbc91dce052a","trusted":true},"cell_type":"code","source":"# Owns a refrigerator, computer, tablet, and television\nheads['bonus'] = 1 * (heads['refrig'] + \n                      heads['computer'] + \n                      (heads['v18q1'] > 0) + \n                      heads['television'])\n\nsns.violinplot('bonus', 'Target', data = heads,\n                figsize = (10, 6));\nplt.title('Target vs Bonus Variable');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"97bfbf33dddaed18073ba92bda62305cc38a770f"},"cell_type":"markdown","source":"## Per Capita Features\n\nAdditional features we can make calculate the number of certain measurements for each person in the household."},{"metadata":{"_uuid":"cc7c411527e8fd9955e8ab83bcec4c78322c6355","trusted":true},"cell_type":"code","source":"heads['phones-per-capita'] = heads['qmobilephone'] / heads['tamviv']\nheads['tablets-per-capita'] = heads['v18q1'] / heads['tamviv']\nheads['rooms-per-capita'] = heads['rooms'] / heads['tamviv']\nheads['rent-per-capita'] = heads['v2a1'] / heads['tamviv']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"615275fff580b041d45a9d358468f3fd8900169d"},"cell_type":"markdown","source":"## Exploring Household Variables\n\nAfter going to all the trouble of getting our features in order, now we can take a look at them in relation to the Target. We've already done a little of this, but now we can try to quantify relationships.\n\n### Measuring Relationships\n\nThere are many ways for measuring relationships between two variables. Here we will examine two of these:\n\n1. The Pearson Correlation: from -1 to 1 measuring the linear relationship between two variables\n2. The Spearman Correlation: from -1 to 1 measuring the monotonic relationship between two variables\n\nThe Spearman correlation is 1 if as one variable increases, the other does as well, even if the relationship is not linear. On the other hand, the Pearson correlation can only be one if the increase is exactly linear. These are best illustrated by example."},{"metadata":{"_uuid":"37b0c3e16ea4e9055bca656ce590cc532a5c3d92","trusted":true},"cell_type":"code","source":"from scipy.stats import spearmanr","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6827b8e42f0462baed0d74409386e1ec255b0f46","trusted":true},"cell_type":"code","source":"def plot_corrs(x, y):\n    \"\"\"Plot data and show the spearman and pearson correlation.\"\"\"\n    \n    # Calculate correlations\n    spr = spearmanr(x, y).correlation\n    pcr = np.corrcoef(x, y)[0, 1]\n    \n    # Scatter plot\n    data = pd.DataFrame({'x': x, 'y': y})\n    plt.figure( figsize = (6, 4))\n    sns.regplot('x', 'y', data = data, fit_reg = False);\n    plt.title(f'Spearman: {round(spr, 2)}; Pearson: {round(pcr, 2)}');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"eb3d0ab8df4bdfab62e817d55437ab0857690784","trusted":true},"cell_type":"code","source":"x = np.array(range(100))\ny = x ** 2\n\nplot_corrs(x, y)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dbe1eac29065813b705c37448f01e7ac660ff784"},"cell_type":"markdown","source":"The Spearman correlation is often considered to be better for ordinal variables such as the Target or the years of education. Most relationshisp in the real world aren't linear, and although the Pearson correlation can be an approximation of how related two variables are, it's inexact and not the best method of comparison."},{"metadata":{"_uuid":"9ff7af6f2f0c34d67db3e352df99066de8343091","trusted":true},"cell_type":"code","source":"x = np.array([1, 1, 1, 2, 3, 3, 4, 4, 4, 5, 5, 6, 7, 8, 8, 9, 9, 9])\ny = np.array([1, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 3, 3, 2, 4, 2, 2, 4])\n\nplot_corrs(x, y)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"391b0659f036f8a6509640eaf5c3a4ef6906566d"},"cell_type":"markdown","source":"In most cases, the values are very similar. "},{"metadata":{"_uuid":"d1a3e7c95f4d24c1253110c0965e1cc8fd43480e","trusted":true},"cell_type":"code","source":"x = np.array(range(-19, 20))\ny = 2 * np.sin(x)\n\nplot_corrs(x, y)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e1074fba48fed2034ce1f176cf2b6652ef661348"},"cell_type":"markdown","source":"First, we'll calculate the Pearson correlation of every variable with the Target."},{"metadata":{"_uuid":"41762a06dff718cdc452c4181df92f352e24b9df","trusted":true},"cell_type":"code","source":"# Use only training data\ntrain_heads = heads.loc[heads['Target'].notnull(), :].copy()\n\npcorrs = pd.DataFrame(train_heads.corr()['Target'].sort_values()).rename(columns = {'Target': 'pcorr'}).reset_index()\npcorrs = pcorrs.rename(columns = {'index': 'feature'})\n\nprint('Most negatively correlated variables:')\nprint(pcorrs.head())\n\nprint('\\nMost positively correlated variables:')\nprint(pcorrs.dropna().tail())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b54fc898fde336625149bda09c5c6bfbfb82b033"},"cell_type":"markdown","source":"For the negative correlations, as we increase the variable, the Target decreases indicating the poverty severity increases. Therefore, as the `warning` increases, the poverty level also increases which makes sense because this was meant to show potential bad signs about a house. The `hogar_nin` is the number of children 0 - 19 in the family which also makes sense: younger children can be financial source of stress on a family leading to higher levels of poverty. Or, families with lower socioeconomic status have more children in the hopes that one of them will be able to succeed. Whatever the explanation, there is a [real link between family size and poverty](https://www.adb.org/sites/default/files/publication/157217/adbi-rp68.pdf)\n\nOn the other hand, for the positive correlations, a higher value means a higher value of Target indicating the poverty severity decreases. The most highly correlated household level variable is `meaneduc`, the average education level of the adults in the household. This relationship between education and poverty intuitively makes sense: [greater levels of education generally correlate with lower levels of poverty](https://www.childfund.org/poverty-and-education/). We don't necessarily know which causes which, but we do know these tend to move in the same direction.\n\nThe general guidelines for correlation values are below, but these will change depending on who you ask ([source](http://www.statstutor.ac.uk/resources/uploaded/pearsons.pdf) for these):\n\n*  .00-.19 “very weak”\n*  .20-.39 “weak”\n*  .40-.59 “moderate”\n*  .60-.79 “strong”\n*  .80-1.0 “very strong”\n\nWhat these correlations show is that there are some weak relationships that hopefully our model will be able to use to learn a mapping from the features to the Target.\n\nNow we can move on to the Spearman correlation. "},{"metadata":{"_uuid":"57ff58385cf395ad96537cf775bfa80a61a45979","trusted":true},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore', category = RuntimeWarning)\n\nfeats = []\nscorr = []\npvalues = []\n\n# Iterate through each column\nfor c in heads:\n    # Only valid for numbers\n    if heads[c].dtype != 'object':\n        feats.append(c)\n        \n        # Calculate spearman correlation\n        scorr.append(spearmanr(train_heads[c], train_heads['Target']).correlation)\n        pvalues.append(spearmanr(train_heads[c], train_heads['Target']).pvalue)\n\nscorrs = pd.DataFrame({'feature': feats, 'scorr': scorr, 'pvalue': pvalues}).sort_values('scorr')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f804082f69658b46ad0c9a134f84122856dd37df"},"cell_type":"markdown","source":"The Spearman correlation coefficient calculation also comes with a `pvalue` indicating the significance level of the relationship. Any `pvalue` less than 0.05 is genearally regarded as significant, although since we are doing multiple comparisons, we want to divide the p-value by the number of comparisons, a process known as the Bonferroni correction. "},{"metadata":{"_uuid":"2aa018089c7e349611c3a0ddb7247d40de9a4351","trusted":true},"cell_type":"code","source":"print('Most negative Spearman correlations:')\nprint(scorrs.head())\nprint('\\nMost positive Spearman correlations:')\nprint(scorrs.dropna().tail())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"65312c90b8be134c59ad0c1b182305f0cc6738b2"},"cell_type":"markdown","source":"For the most part, the two methods of calculating correlations are in agreement. Just out of curiousity, we can look for the values that are furthest apart."},{"metadata":{"_uuid":"9cf2669563f444e88edb40a22e77ac3c0debef63","trusted":true},"cell_type":"code","source":"corrs = pcorrs.merge(scorrs, on = 'feature')\ncorrs['diff'] = corrs['pcorr'] - corrs['scorr']\n\ncorrs.sort_values('diff').head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e8246b539b238758630826413d7a191ddf30488b","trusted":true},"cell_type":"code","source":"corrs.sort_values('diff').dropna().tail()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5a5a5f7ce58dfbd7652c97a877a027a053b66e26"},"cell_type":"markdown","source":"The largest discrepancy in the correlations is `dependency`. We can make a scatterplot of the `Target` versus the `dependency` to visualize the relationship. We'll add a little jitter to the plot because these are both discrete variables."},{"metadata":{"_uuid":"8361ec53ddf8642b79a3e206bb9fee47514792ae","trusted":true},"cell_type":"code","source":"sns.lmplot('dependency', 'Target', fit_reg = True, data = train_heads, x_jitter=0.05, y_jitter=0.05);\nplt.title('Target vs Dependency');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1575020ee3771a475db3243ececb639e5f7ef8e7"},"cell_type":"markdown","source":"It's hard to see the relationship, but it's slightly negative: as the `dependency` increases, the value of the `Target` decreases. This makes sense: the `dependency` is the number of dependent individuals divided by the number of non-dependents. As we increase this value, the poverty severty tends to increase: having more dependent family members (who usually are non-working) leads to higher levels of poverty because they must be supported by the non-dependent family members. "},{"metadata":{"trusted":true,"_uuid":"c45feb5e458547023a97560328141a5e41a22ea7"},"cell_type":"code","source":"sns.lmplot('rooms-per-capita', 'Target', fit_reg = True, data = train_heads, x_jitter=0.05, y_jitter=0.05);\nplt.title('Target vs Rooms Per Capita');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"46e99d738f30a21006a91209721f03b085fa1fa2"},"cell_type":"markdown","source":"#### Correlation Heatmap \n\nOne of my favorite plots is the correlation heatmap because it shows a ton of info in one image. For the heatmap, we'll pick 7 variables and show the correlations between themselves and with the target. "},{"metadata":{"_uuid":"b3535cc3f355808f220d78862360e7aede64d92c","trusted":true},"cell_type":"code","source":"variables = ['Target', 'dependency', 'warning', 'walls+roof+floor', 'meaneduc',\n             'floor', 'r4m1', 'overcrowding']\n\n# Calculate the correlations\ncorr_mat = train_heads[variables].corr().round(2)\n\n# Draw a correlation heatmap\nplt.rcParams['font.size'] = 18\nplt.figure(figsize = (12, 12))\nsns.heatmap(corr_mat, vmin = -0.5, vmax = 0.8, center = 0, \n            cmap = plt.cm.RdYlGn_r, annot = True);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"05b815c1126418ac5821be1053e602d36c02e0d8"},"cell_type":"markdown","source":"This plot shows us that there are a number of variables that have a weak correlation with the `Target`. There are also high correlations between some variables (such as `floor` and `walls+roof+floor`) which could pose an issue because of collinearity. \n\n### Features Plot\n\nFor the final exploration of the household level data, we can make a plot of some of the most correlated variables with the Target. This shows scatterplots on the upper triangle, kernel density estimate (kde) plots on the diagonal, and 2D KDE plots on the lower triangle."},{"metadata":{"_uuid":"012c5afa076a04df5d833c725e2fdcd097b783a5","trusted":true},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')\n\n# Copy the data for plotting\nplot_data = train_heads[['Target', 'dependency', 'walls+roof+floor',\n                         'meaneduc', 'overcrowding']]\n\n# Create the pairgrid object\ngrid = sns.PairGrid(data = plot_data, size = 4, diag_sharey=False,\n                    hue = 'Target', hue_order = [4, 3, 2, 1], \n                    vars = [x for x in list(plot_data.columns) if x != 'Target'])\n\n# Upper is a scatter plot\ngrid.map_upper(plt.scatter, alpha = 0.8, s = 20)\n\n# Diagonal is a histogram\ngrid.map_diag(sns.kdeplot)\n\n# Bottom is density plot\ngrid.map_lower(sns.kdeplot, cmap = plt.cm.OrRd_r);\ngrid = grid.add_legend()\nplt.suptitle('Feature Plots Colored By Target', size = 32, y = 1.05);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"058c36fdb6108c2a84045d4592ad8e53002a3b14"},"cell_type":"markdown","source":"We'll leave the feature engineering of the household variables for now. Later, we can come back to this step if we are not pleased with the model performance."},{"metadata":{"_uuid":"cb4e91d2ce09ff60e61ec1f3be16932ed7d0acd5","trusted":true},"cell_type":"code","source":"household_feats = list(heads.columns)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"de174018faf038f680d405baa68a2aba93ab1b77"},"cell_type":"markdown","source":"# Individual Level Variables\n\nThere are two types of individual level variables: Boolean (1 or 0 for True or False) and ordinal (discrete values with a meaningful ordering). "},{"metadata":{"_uuid":"88e7cc6cd8acab811702215a2dc2b08b8e203dd2","trusted":true},"cell_type":"code","source":"ind = data[id_ + ind_bool + ind_ordered]\nind.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"931587721d828d6a534d41f34f9c6a561a454cf4"},"cell_type":"markdown","source":"## Redundant Individual Variables\n\nWe can do the same process we did with the household level variables to identify any redundant individual variables. We'll focus on any variables that have an absolute magnitude of the correlation coefficient greater than 0.95."},{"metadata":{"_uuid":"7b8c7059c959d69372304663e1c4a8d839d27eac","trusted":true},"cell_type":"code","source":"# Create correlation matrix\ncorr_matrix = ind.corr()\n\n# Select upper triangle of correlation matrix\nupper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))\n\n# Find index of feature columns with correlation greater than 0.95\nto_drop = [column for column in upper.columns if any(abs(upper[column]) > 0.95)]\n\nto_drop","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b5feba5e8b3843a9e255497497184b5e8ab4754f"},"cell_type":"markdown","source":"This is simply the opposite of male! We can remove the male flag."},{"metadata":{"_uuid":"e34832557ff343d282c11259eb1e1e4c477b0158","trusted":true},"cell_type":"code","source":"ind = ind.drop(columns = 'male')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f6d5fcca07c194b82ba1ea89b8290bbd71ef632b"},"cell_type":"markdown","source":"### Creating Ordinal Variables\n\nMuch as we did with the household level data, we can map existing columns to an ordinal variable. Here we will focus on the `instlevel_` variables which indicate the amount of education an individual has from `instlevel1`: no level of education to `instlevel9`: postgraduate education. \n\nTo create the ordinal variable, for each individual, we will simply find which column is non-zero. The education has an inherent ordering (higher is better) so this conversion to an ordinal variable makes sense in the problem context."},{"metadata":{"_uuid":"870c1271305145593a6bee0fec5d8dd64141f384","trusted":true},"cell_type":"code","source":"ind[[c for c in ind if c.startswith('instl')]].head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f0cc456f9ec85144e32702a98f2c82871d72d268","trusted":true},"cell_type":"code","source":"ind['inst'] = np.argmax(np.array(ind[[c for c in ind if c.startswith('instl')]]), axis = 1)\n\nplot_categoricals('inst', 'Target', ind, annotate = False);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2005b023394839c61300ad2f131fd55e7065b36e"},"cell_type":"markdown","source":"Higher levels of education seem to correspond to less extreme levels of poverty. We do need to keep in mind this is on an individual level though and we eventually will have to aggregate this data at the household level. "},{"metadata":{"_uuid":"36dea7591a21564e8762547d098ae7c4ca2010b0","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (10, 8))\nsns.violinplot(x = 'Target', y = 'inst', data = ind);\nplt.title('Education Distribution by Target');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bf9faa15211503eba64468154970d620b8c0126a","trusted":true},"cell_type":"code","source":"# Drop the education columns\n# ind = ind.drop(columns = [c for c in ind if c.startswith('instlevel')])\nind.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"132986030dc283109318780a2f6330a75b72e3aa"},"cell_type":"markdown","source":"### Feature Construction\n\nWe can make a few features using the existing data. For example, we can divide the years of schooling by the age."},{"metadata":{"_uuid":"f673326c23b7d9b0806924f05af22c0601a07798","trusted":true},"cell_type":"code","source":"ind['escolari/age'] = ind['escolari'] / ind['age']\n\nplt.figure(figsize = (10, 8))\nsns.violinplot('Target', 'escolari/age', data = ind);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"35dddba8210e8fe298f21a61812a8492ea9697d4"},"cell_type":"markdown","source":"We can also take our new variable, `inst`, and divide this by the age. The final variable we'll name `tech`: this represents the combination of tablet and mobile phones."},{"metadata":{"_uuid":"eee1804b83ecf9ac83fb6aef014fe6ac5d3f93f8","trusted":true},"cell_type":"code","source":"ind['inst/age'] = ind['inst'] / ind['age']\nind['tech'] = ind['v18q'] + ind['mobilephone']\nind['tech'].describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fbaef320687392da25e028e2021e5d5be574bb7e"},"cell_type":"markdown","source":"## Feature Engineering through Aggregations\n\nIn order to incorporate the individual data into the household data, we need to aggregate it for each household. The simplest way to do this is to `groupby` the family id `idhogar` and then `agg` the data. For the aggregations for ordered or continuous variables, we can use six, five of which are built in to pandas, and one of which we define ourselves `range_`. The boolean aggregations can be the same, but this will create many redundant columns which we will then need to drop. For this case, we'll use the same aggregations and then go back and drop the redundant columns."},{"metadata":{"_uuid":"f59c9612532cfc67c2551268370b662d03fe84ed","trusted":true},"cell_type":"code","source":"# Define custom function\nrange_ = lambda x: x.max() - x.min()\nrange_.__name__ = 'range_'\n\n# Group and aggregate\nind_agg = ind.drop(columns = 'Target').groupby('idhogar').agg(['min', 'max', 'sum', 'count', 'std', range_])\nind_agg.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1eb6430fcba5d962ade9aa7e866bb6f16892057e"},"cell_type":"markdown","source":"With just that one line, we go from 30 features to 180. Next we can rename the columns to make it easier to keep track."},{"metadata":{"_uuid":"b444e21422941ca1c26b0b3404c3e2a16d5af991","trusted":true},"cell_type":"code","source":"# Rename the columns\nnew_col = []\nfor c in ind_agg.columns.levels[0]:\n    for stat in ind_agg.columns.levels[1]:\n        new_col.append(f'{c}-{stat}')\n        \nind_agg.columns = new_col\nind_agg.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f196a265f84ec22177da1822d0d732ff6f886d13","trusted":true},"cell_type":"code","source":"ind_agg.iloc[:, [0, 1, 2, 3, 6, 7, 8, 9]].head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"aa8ae6e58944998218d3ddae84390577a9979104"},"cell_type":"markdown","source":"### Feature Selection \n\nAs a first round of feature selection, we can remove one out of every pair of variables with a correlation greater than 0.95. "},{"metadata":{"_uuid":"6dd33798e4f99f82ac5c1ec659906bedeb12febf","trusted":true},"cell_type":"code","source":"# Create correlation matrix\ncorr_matrix = ind_agg.corr()\n\n# Select upper triangle of correlation matrix\nupper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))\n\n# Find index of feature columns with correlation greater than 0.95\nto_drop = [column for column in upper.columns if any(abs(upper[column]) > 0.95)]\n\nprint(f'There are {len(to_drop)} correlated columns to remove.')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d1f7bc5fc7635bfbee47507666fa844fbeef8c50"},"cell_type":"markdown","source":"We'll drop the columns and then merge with the `heads` data to create a final dataframe."},{"metadata":{"_uuid":"bccbcdbad813f7c1596fd8209de3df46e18db199","trusted":true},"cell_type":"code","source":"ind_agg = ind_agg.drop(columns = to_drop)\nind_feats = list(ind_agg.columns)\n\n# Merge on the household id\nfinal = heads.merge(ind_agg, on = 'idhogar', how = 'left')\n\nprint('Final features shape: ', final.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3e0b6e903c0780560d5e3b4283dcca8be3a79cf3","trusted":true},"cell_type":"code","source":"final.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8b31ac05762ca5c3a29bbe9739e47d688fd5fc35"},"cell_type":"markdown","source":"### Final Data Exploration\n\nWe'll do a little bit of exploration. "},{"metadata":{"_uuid":"32accf5f15275477da10863997bac77b1bb6113e","trusted":true},"cell_type":"code","source":"corrs = final.corr()['Target']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2e75e94d93a61cefe0f4c9fa6ac26cb285bafc0e","trusted":true},"cell_type":"code","source":"corrs.sort_values().head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e6c8a457da110b83db2fe47cf5c631c9d3bcce58","trusted":true},"cell_type":"code","source":"corrs.sort_values().dropna().tail()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5f3a895a60ad479fc2665b7ce1e88ddabed1f4c6"},"cell_type":"markdown","source":"We can see some of the variables that we made are highly correlated with the Target. Whether these variables are actually _useful_ will be determined in the modeling stage. "},{"metadata":{"_uuid":"7b2ed3ed02f7921fab8a9b2b57d98bdaf6f3ff03","trusted":true},"cell_type":"code","source":"plot_categoricals('escolari-max', 'Target', final, annotate=False);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d1c19def77a277533c2f663564eb99c1a29b314d","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (10, 6))\nsns.violinplot(x = 'Target', y = 'escolari-max', data = final);\nplt.title('Max Schooling by Target');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2adb7da9565c8e85f2899c5a8992d2ce13617a35","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (10, 6))\nsns.boxplot(x = 'Target', y = 'escolari-max', data = final);\nplt.title('Max Schooling by Target');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"83d009bf8a4f8f4cf18085524e4a5fdebc043286","scrolled":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (10, 6))\nsns.boxplot(x = 'Target', y = 'meaneduc', data = final);\nplt.xticks([0, 1, 2, 3], poverty_mapping.values())\nplt.title('Average Schooling by Target');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c3c0a518e219d7e2bafe0dba04c0fae4c0f4c7be","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (10, 6))\nsns.boxplot(x = 'Target', y = 'overcrowding', data = final);\nplt.xticks([0, 1, 2, 3], poverty_mapping.values())\nplt.title('Overcrowding by Target');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"18187d9d6e401a27362372a3ec361a40bc157d25"},"cell_type":"markdown","source":"One other feature that might be useful is the gender of the head of household. Since we aggregated the data, we'll have to go back to the individual level data and find the gender for the head of household."},{"metadata":{"_uuid":"bc12c2723bb2867b7890bba79fe0be1887ba3a4f","trusted":true},"cell_type":"code","source":"head_gender = ind.loc[ind['parentesco1'] == 1, ['idhogar', 'female']]\nfinal = final.merge(head_gender, on = 'idhogar', how = 'left').rename(columns = {'female': 'female-head'})","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1ba7d757545ead1b25ab27bd2bf149ac9819e598","trusted":true},"cell_type":"code","source":"final.groupby('female-head')['Target'].value_counts(normalize=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b251e03cfb73b70ecb1563150ed6ce5a5318cdb3"},"cell_type":"markdown","source":"It looks like households where the head is female are slightly more likely to have a severe level of poverty."},{"metadata":{"_uuid":"03b2c39589c3488a9f1241aed8760718fb9f3176","trusted":true},"cell_type":"code","source":"sns.violinplot(x = 'female-head', y = 'Target', data = final);\nplt.title('Target by Female Head of Household');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"601b30f76124044c2949bb3e95d58067b3be2e00"},"cell_type":"markdown","source":"We can also look at the difference in average education by whether or not the family has a female head of household."},{"metadata":{"_uuid":"e81f71c5576a83dc59b43d53fa79f1b422705e4b","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 8))\nsns.boxplot(x = 'Target', y = 'meaneduc', hue = 'female-head', data = final);\nplt.title('Average Education by Target and Female Head of Household', size = 16);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"796eba22d6d1ba7cc68ba0ab300613934bbb44bf"},"cell_type":"markdown","source":"It looks like at every value of the `Target`, households with female heads have higher levels of education. Yet, we saw that overall, households with female heads are more likely to have severe poverty. "},{"metadata":{"_uuid":"4691aa23a0407a63dec6da789c6da64c9c0186ba","trusted":true},"cell_type":"code","source":"final.groupby('female-head')['meaneduc'].agg(['mean', 'count'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d0d75824293f57247e833e2c1386ff807d56585d"},"cell_type":"markdown","source":"Overall, the average education of households with female heads is slightly higher than those with male heads. I'm not too sure what to make of this, but it seems right to me."},{"metadata":{"_uuid":"7d16a374bc37c67c427a2dac99d5be12fed49ee3"},"cell_type":"markdown","source":"# Machine Learning Modeling\n\nOnce feature engineering/construction is done, we can get started with the machine learning! All of our data (both training and testing) is aggregated for each household and so can be directly used in a model. To first show the process of modeling, we'll use the capable Random Forest Classifier in Scikit-Learn. This probably won't get us to the top of the leaderboard, but it will allow us to establish a baseline. Later we'll try several other models including the powerful Gradient Boosting Machine.\n\nTo assess our model, we'll use 10-fold cross validation on the training data. This will essentially _train and test the model 10 times_ using different splits of the training data. 10-fold cross validation is an effective method for estimating the performance of a model on the test set. We want to look at the average performance in cross validation as well as the standard deviation to see how much scores change between the folds. We use the `F1 Macro` measure to evaluate performance."},{"metadata":{"_uuid":"a6f26316452144bfeff36474c36c32d1334668ba","trusted":true},"cell_type":"code","source":"from sklearn.ensemble import RandomForestClassifier\nfrom sklearn.metrics import f1_score, make_scorer\nfrom sklearn.model_selection import cross_val_score\nfrom sklearn.preprocessing import Imputer\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.pipeline import Pipeline\n\n# Custom scorer for cross validation\nscorer = make_scorer(f1_score, greater_is_better=True, average = 'macro')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8037345f25c39fdaa811adfd38f65714678359a2","trusted":true},"cell_type":"code","source":"# Labels for training\ntrain_labels = np.array(list(final[final['Target'].notnull()]['Target'].astype(np.uint8)))\n\n# Extract the training data\ntrain_set = final[final['Target'].notnull()].drop(columns = ['Id', 'idhogar', 'Target'])\ntest_set = final[final['Target'].isnull()].drop(columns = ['Id', 'idhogar', 'Target'])\n\n# Submission base which is used for making submissions to the competition\nsubmission_base = test[['Id', 'idhogar']].copy()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"44fc64eb3f3be49735f6198209b0971e2bc9e072"},"cell_type":"markdown","source":"Because we are going to be comparing different models, we want to scale the features (limit the range of each column to between 0 and 1). For many ensemble models this is not necessary, but when we use models that depend on a distance metric, such as KNearest Neighbors or the Support Vector Machine, feature scaling is an absolute necessity. When comparing different models, it's always safest to scale the features. We also impute the missing values with the median of the feature.\n\nFor imputing missing values and scaling the features in one step, we can make a pipeline. This will be fit on the training data and used to transform the training and testing data."},{"metadata":{"_uuid":"e4459fe1d50077ba09995fb01664ade706905d10","trusted":true},"cell_type":"code","source":"features = list(train_set.columns)\n\npipeline = Pipeline([('imputer', Imputer(strategy = 'median')), \n                      ('scaler', MinMaxScaler())])\n\n# Fit and transform training data\ntrain_set = pipeline.fit_transform(train_set)\ntest_set = pipeline.transform(test_set)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d288a06b9b05d4868c3771e1e35c7658eb678e3f"},"cell_type":"markdown","source":"The data has no missing values and is scaled between zero and one. This means it can be directly used in any Scikit-Learn model. "},{"metadata":{"_uuid":"b1b0fde920b9cb7caa2653e417b9cfb7cb3db075","trusted":true},"cell_type":"code","source":"model = RandomForestClassifier(n_estimators=100, random_state=10, \n                               n_jobs = -1)\n# 10 fold cross validation\ncv_score = cross_val_score(model, train_set, train_labels, cv = 10, scoring = scorer)\n\nprint(f'10 Fold Cross Validation F1 Score = {round(cv_score.mean(), 4)} with std = {round(cv_score.std(), 4)}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b68367d25404464a1e62f49eab871f4fd97d75e8"},"cell_type":"markdown","source":"That score is not great, but it will serve as a baseline and leaves us plenty of room to improve! "},{"metadata":{"_uuid":"d93a45c46832d0b7804a112925d19bb5739e586b"},"cell_type":"markdown","source":"## Feature Importances\n\nWith a tree-based model, we can look at the feature importances which show a relative ranking of the usefulness of features in the model. These represent the sum of the reduction in impurity at nodes that used the variable for splitting, but we don't have to pay much attention to the absolute value. Instead we'll focus on relative scores.\n\nIf we want to view the feature importances, we'll have to train a model on the whole training set. Cross validation does not return the feature importances. "},{"metadata":{"_uuid":"76e8463384bd4b81f6d6ba8288b0684571bac380","trusted":true},"cell_type":"code","source":"model.fit(train_set, train_labels)\n\n# Feature importances into a dataframe\nfeature_importances = pd.DataFrame({'feature': features, 'importance': model.feature_importances_})\nfeature_importances.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a474c17b1394ef6c008c6b9133fd40b4d8f04e98"},"cell_type":"markdown","source":"Below is a short function we'll use to plot the feature importances. I use this function a lot and often copy and paste it between scripts. I hope the documentation makes sense! "},{"metadata":{"_uuid":"d9912b9bdd0204aea423ec6e610f7788726ed038","trusted":true},"cell_type":"code","source":"def plot_feature_importances(df, n = 10, threshold = None):\n    \"\"\"Plots n most important features. Also plots the cumulative importance if\n    threshold is specified and prints the number of features needed to reach threshold cumulative importance.\n    Intended for use with any tree-based feature importances. \n    \n    Args:\n        df (dataframe): Dataframe of feature importances. Columns must be \"feature\" and \"importance\".\n    \n        n (int): Number of most important features to plot. Default is 15.\n    \n        threshold (float): Threshold for cumulative importance plot. If not provided, no plot is made. Default is None.\n        \n    Returns:\n        df (dataframe): Dataframe ordered by feature importances with a normalized column (sums to 1) \n                        and a cumulative importance column\n    \n    Note:\n    \n        * Normalization in this case means sums to 1. \n        * Cumulative importance is calculated by summing features from most to least important\n        * A threshold of 0.9 will show the most important features needed to reach 90% of cumulative importance\n    \n    \"\"\"\n    plt.style.use('fivethirtyeight')\n    \n    # Sort features with most important at the head\n    df = df.sort_values('importance', ascending = False).reset_index(drop = True)\n    \n    # Normalize the feature importances to add up to one and calculate cumulative importance\n    df['importance_normalized'] = df['importance'] / df['importance'].sum()\n    df['cumulative_importance'] = np.cumsum(df['importance_normalized'])\n    \n    plt.rcParams['font.size'] = 12\n    \n    # Bar plot of n most important features\n    df.loc[:n, :].plot.barh(y = 'importance_normalized', \n                            x = 'feature', color = 'darkgreen', \n                            edgecolor = 'k', figsize = (12, 8),\n                            legend = False, linewidth = 2)\n\n    plt.xlabel('Normalized Importance', size = 18); plt.ylabel(''); \n    plt.title(f'{n} Most Important Features', size = 18)\n    plt.gca().invert_yaxis()\n    \n    \n    if threshold:\n        # Cumulative importance plot\n        plt.figure(figsize = (8, 6))\n        plt.plot(list(range(len(df))), df['cumulative_importance'], 'b-')\n        plt.xlabel('Number of Features', size = 16); plt.ylabel('Cumulative Importance', size = 16); \n        plt.title('Cumulative Feature Importance', size = 18);\n        \n        # Number of features needed for threshold cumulative importance\n        # This is the index (will need to add 1 for the actual number)\n        importance_index = np.min(np.where(df['cumulative_importance'] > threshold))\n        \n        # Add vertical line to plot\n        plt.vlines(importance_index + 1, ymin = 0, ymax = 1.05, linestyles = '--', colors = 'red')\n        plt.show();\n        \n        print('{} features required for {:.0f}% of cumulative importance.'.format(importance_index + 1, \n                                                                                  100 * threshold))\n    \n    return df\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c0d04a9495e7be357599a1c91ea196908efd0624","trusted":true},"cell_type":"code","source":"norm_fi = plot_feature_importances(feature_importances, threshold=0.95)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5403f4fdc5a5bd8192a2ec3a09f062d8573f3c1b"},"cell_type":"markdown","source":"__Education reigns supreme!__ The most important variable is the average amount of education in the household, followed by the maximum education of anyone in the household.  I have a suspicion these variables are highly correlated (collinear) which means we may want to remove one of them from the data. The other most important features are a combination of variables we created and variables that were already present in the data. \n\nIt's interesting that we only need 106 of the ~180 features to account for 90% of the importance. This tells us that we may be able to remove some of the features. However, feature importances don't tell us which direction of the feature is important (for example, we can't use these to tell whether more or less education leads to more severe poverty) they only tell us which features the model considered relevant. "},{"metadata":{"_uuid":"540771805485f245fa5a0a7d18f5fc35b7169e2a","trusted":true},"cell_type":"code","source":"def kde_target(df, variable):\n    \"\"\"Plots the distribution of `variable` in `df` colored by the `Target` column\"\"\"\n    \n    colors = {1: 'red', 2: 'orange', 3: 'blue', 4: 'green'}\n\n    plt.figure(figsize = (12, 8))\n    \n    df = df[df['Target'].notnull()]\n    \n    for level in df['Target'].unique():\n        subset = df[df['Target'] == level].copy()\n        sns.kdeplot(subset[variable].dropna(), \n                    label = f'Poverty Level: {level}', \n                    color = colors[int(subset['Target'].unique())])\n\n    plt.xlabel(variable); plt.ylabel('Density');\n    plt.title('{} Distribution'.format(variable.capitalize()));","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3db2aab481d2c7a641b722c87ef9ff9eded0d102","trusted":true},"cell_type":"code","source":"kde_target(final, 'meaneduc')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0b74cce70fba4268ca56a5202f879309045c9eb0","trusted":true},"cell_type":"code","source":"kde_target(final, 'escolari/age-range_')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fb3574b0f4d56a904729cc7c800e5619e3e89e8f"},"cell_type":"markdown","source":"# Model Selection\n\nNow that we have a good set of features, it's time to get into the modeling. We already tried one basic model, the Random Forest Classifier which delivered a best macro F1 of 0.35. However, in machine learning, there is no way to know ahead of time which model will work best for a given dataset. The following plot shows that __there are some problems where even Gaussian Naive Bayes will outperform a gradient boosting machine__. This is from [an excellent paper by Randal Olson that discusses many points of machine learning](https://psb.stanford.edu/psb-online/proceedings/psb18/olson.pdf)\n\n![algorithm_comparison](https://raw.githubusercontent.com/WillKoehrsen/Machine-Learning-Projects/master/algorithm_comparison.png)\n\nWhat this plot tells us is that we have to try out a number of different models to see which is optimal. Most people eventually settle on the __gradient boosting machine__ and we will try that out, but for now we'll take a look at some of the other options.  There are literally dozens (maybe hundreds) of multi-class machine learning models if we look at the [Scikit-Learn documentation](http://scikit-learn.org/stable/modules/multiclass.html). We don't have to try them all, but we should sample from the options.\n\nWhat we want to do is write a function that can evaluate a model. This will be pretty simple since we already wrote most of the code. In addition to the Random Forest Classifier, we'll try eight other Scikit-Learn models. Luckily, this dataset is relatively small and we can rapidly iterate through the models. We will make a dataframe to hold the results and the function will add a row to the dataframe for each model."},{"metadata":{"_uuid":"78a7710f39bab271c4320f42b9bbc7e3f9a4bd2f","trusted":true},"cell_type":"code","source":"# Model imports\nfrom sklearn.svm import LinearSVC\nfrom sklearn.naive_bayes import GaussianNB\nfrom sklearn.neural_network import MLPClassifier\nfrom sklearn.linear_model import LogisticRegressionCV, RidgeClassifierCV\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis\nfrom sklearn.neighbors import KNeighborsClassifier","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"665cd410e664a6092b756a888727f166fbb5a5ed","trusted":true},"cell_type":"code","source":"import warnings \nfrom sklearn.exceptions import ConvergenceWarning\n\n# Filter out warnings from models\nwarnings.filterwarnings('ignore', category = ConvergenceWarning)\nwarnings.filterwarnings('ignore', category = DeprecationWarning)\nwarnings.filterwarnings('ignore', category = UserWarning)\n\n# Dataframe to hold results\nmodel_results = pd.DataFrame(columns = ['model', 'cv_mean', 'cv_std'])\n\ndef cv_model(train, train_labels, model, name, model_results=None):\n    \"\"\"Perform 10 fold cross validation of a model\"\"\"\n    \n    cv_scores = cross_val_score(model, train, train_labels, cv = 10, scoring=scorer, n_jobs = -1)\n    print(f'10 Fold CV Score: {round(cv_scores.mean(), 5)} with std: {round(cv_scores.std(), 5)}')\n    \n    if model_results is not None:\n        model_results = model_results.append(pd.DataFrame({'model': name, \n                                                           'cv_mean': cv_scores.mean(), \n                                                            'cv_std': cv_scores.std()},\n                                                           index = [0]),\n                                             ignore_index = True)\n\n        return model_results","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ab07945c331294a68d600de58e4ed9c2739473bf","trusted":true},"cell_type":"code","source":"model_results = cv_model(train_set, train_labels, LinearSVC(), \n                         'LSVC', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2ea1ec18bda63d33c6e6482de5361d1728f41bbd"},"cell_type":"markdown","source":"That's one model to cross off the list (although we didn't perform hyperparameter tuning so the actual performance could possibly be improved)."},{"metadata":{"_uuid":"e5c13ea9fb13134e3b25d59bedf5b9f974e1517c","trusted":true},"cell_type":"code","source":"model_results = cv_model(train_set, train_labels, \n                         GaussianNB(), 'GNB', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a73eca1b3b7be54e3d6bae663169e67854f91a06"},"cell_type":"markdown","source":"That performance is very poor. I don't think we need to revisit the Gaussian Naive Bayes method (although there are problems on which it can outperform the Gradient Boosting Machine)."},{"metadata":{"_uuid":"8b1dca57fa07e406f1fe81fa320b3b97ab078d8e","trusted":true},"cell_type":"code","source":"model_results = cv_model(train_set, train_labels, \n                         MLPClassifier(hidden_layer_sizes=(32, 64, 128, 64, 32)),\n                         'MLP', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"84fb16fdda905c410ef60cef95016c85306fe81c"},"cell_type":"markdown","source":"The multi-layer perceptron (a deep neural network) has decent performance. This might be an option if we are able to hyperparameter tune the network. However, the limited amount of data could be an issue with a neural network as these generally require hundreds of thousands of examples to learn effectively."},{"metadata":{"_uuid":"2310f6d79426e16c9e3af6f8d4bfe03aa847d1b7","trusted":true},"cell_type":"code","source":"model_results = cv_model(train_set, train_labels, \n                          LinearDiscriminantAnalysis(), \n                          'LDA', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"44862a70f227c4a2ec0b73431f0fdeeeebf09f0a"},"cell_type":"markdown","source":"__If you run `LinearDiscriminantAnalysis` without filtering out the `UserWarning`s, you get many messages saying \"Variables are collinear.\"__ This might give us a hint that we want to remove some collinear features! We might want to try this model again after removing the collinear variables because the score is comparable to the random forest."},{"metadata":{"_uuid":"5c5c3d651e96cfe01af34f93e6faa501fa35a94f","trusted":true},"cell_type":"code","source":"model_results = cv_model(train_set, train_labels, \n                         RidgeClassifierCV(), 'RIDGE', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"94459bb16bb7ddf3f49ee1a547818dd52efbe764"},"cell_type":"markdown","source":"The linear model (with ridge regularization) does surprisingly well. This might indicate that a simple model can go a long way in this problem (although we'll probably end up using a more powerful method)."},{"metadata":{"_uuid":"b10b93b0caed61b6b54df8acfdda3e144c6499cf","trusted":true},"cell_type":"code","source":"for n in [5, 10, 20]:\n    print(f'\\nKNN with {n} neighbors\\n')\n    model_results = cv_model(train_set, train_labels, \n                             KNeighborsClassifier(n_neighbors = n),\n                             f'knn-{n}', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"49639627d4e1513ca07395b9a471295afc0a03c6"},"cell_type":"markdown","source":"As one more attempt, we'll consider the ExtraTreesClassifier, a variant on the random forest using ensembles of decision trees as well."},{"metadata":{"_uuid":"eeab91bfb4b953e0c7aa5112625b65268537c662","trusted":true},"cell_type":"code","source":"from sklearn.ensemble import ExtraTreesClassifier\n\nmodel_results = cv_model(train_set, train_labels, \n                         ExtraTreesClassifier(n_estimators = 100, random_state = 10),\n                         'EXT', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4f4a3dc32b62b3fcc625e54f8d92afd9f1de2e7d"},"cell_type":"markdown","source":"## Comparing Model Performance\n\nWith the modeling results in a dataframe, we can plot them to see which model does the best."},{"metadata":{"_uuid":"d8c2a9a2df39625e6b1066ca950e44b5fe68862e","trusted":true},"cell_type":"code","source":"model_results = cv_model(train_set, train_labels,\n                          RandomForestClassifier(100, random_state=10),\n                              'RF', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8cbd9d6729288337c5c3cace44af254ebd480fce","trusted":true},"cell_type":"code","source":"model_results.set_index('model', inplace = True)\nmodel_results['cv_mean'].plot.bar(color = 'orange', figsize = (8, 6),\n                                  yerr = list(model_results['cv_std']),\n                                  edgecolor = 'k', linewidth = 2)\nplt.title('Model F1 Score Results');\nplt.ylabel('Mean F1 Score (with error bar)');\nmodel_results.reset_index(inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"876ff3ad8202066da47d7aea5d1784a981c81edb"},"cell_type":"markdown","source":"The most likely candidate seems to be the Random Forest because it does best right out of the box.  While we didn't tune any of the hyperparameters so the comparison between models is not perfect, these results reflect those of many other Kaggle competitiors finding that tree-based ensemble methods (including the Gradient Boosting Machine) perform very well on structured datasets. Hyperparameter performance does improve the performance of machine learning models, but we don't have time to try all possible combinations of settings for all models. The graph below ([from the paper by Randal Olson](https://psb.stanford.edu/psb-online/proceedings/psb18/olson.pdf)) shows the effect of hyperparameter tuning versus the default values in Scikit-Learn.\n\n![hyperparameter_improvement](https://raw.githubusercontent.com/WillKoehrsen/Machine-Learning-Projects/master/hyperparameter_improvement.png)\n\nIn most cases the accuracy gain is less than 10% so the worst model is probably not suddenly going to become the best model through tuning. "},{"metadata":{"_uuid":"3e9b271483c048ab27e584577a94f0dc4a2e47e6"},"cell_type":"markdown","source":"For now we'll say the random forest does the best. Later we'll look at using the Gradient Boosting Machine, although not implemented in Scikit-Learn. Instead we'll be using the more powerful [LightGBM version](http://lightgbm.readthedocs.io/en/latest/). Now, let's turn to making a submission using the random forest."},{"metadata":{"_uuid":"de34b3435f0ff50a2c437799867de0689fc91c7c"},"cell_type":"markdown","source":"# Making a Submission\n\nIn order to make a submission, we need the test data. Fortunately, we have the test data formatted in exactly the same manner as the train data. \n\nThe format of a testing submission is shown below. Although we are making predictions for each household, we actually need one row per individual (identified by the `Id`) but only the prediction for the head of household is scored. \n\n```\nId,Target\nID_2f6873615,1\nID_1c78846d2,2\nID_e5442cf6a,3\nID_a8db26a79,4\nID_a62966799,4 \n```\n"},{"metadata":{"_uuid":"3e38e6c120980df2d11f1fc5a2377b674b066481"},"cell_type":"markdown","source":"The `submission_base` will have all the individuals in the test set since we have to have a \"prediction\" for each individual while the `test_ids` will only contain the `idhogar` from the heads of households. When predicting, we only predict for each household and then we merge the `predictions` dataframe with all of the individuals on the household id (`idhogar`). This will set the `Target` to the same value for everyone in a household. For the test households without a head of household, we can just set these predictions to 4 since they will not be scored. "},{"metadata":{"_uuid":"98fa2a883bf9365159dbbc51a15fb9c7ed1c1c10","trusted":true},"cell_type":"code","source":"test_ids = list(final.loc[final['Target'].isnull(), 'idhogar'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5c4072b5bb0210abddfd1652eb123f988048c081"},"cell_type":"markdown","source":"The function below takes in a model, a training set, the training labels, and a testing set and performs the following operations:\n\n* Trains the model on the training data using `fit`\n* Makes predictions on the test data using `predict`\n* Creates a `submission` dataframe that can be saved and uploaded to the competition"},{"metadata":{"_uuid":"444972fb20122f5a73f868789ecd170e2ffdced8","trusted":true},"cell_type":"code","source":"def submit(model, train, train_labels, test, test_ids):\n    \"\"\"Train and test a model on the dataset\"\"\"\n    \n    # Train on the data\n    model.fit(train, train_labels)\n    predictions = model.predict(test)\n    predictions = pd.DataFrame({'idhogar': test_ids,\n                               'Target': predictions})\n\n     # Make a submission dataframe\n    submission = submission_base.merge(predictions, \n                                       on = 'idhogar',\n                                       how = 'left').drop(columns = ['idhogar'])\n    \n    # Fill in households missing a head\n    submission['Target'] = submission['Target'].fillna(4).astype(np.int8)\n\n    return submission ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5107d84e364842331dd23e3112835646e22ce004"},"cell_type":"markdown","source":"Let's make a submission with the Random Forest."},{"metadata":{"_uuid":"a0eeba88ccae80a22fad6d62ce924f6b517183d6","trusted":true},"cell_type":"code","source":"rf_submission = submit(RandomForestClassifier(n_estimators = 100, \n                                              random_state=10, n_jobs = -1), \n                         train_set, train_labels, test_set, test_ids)\n\nrf_submission.to_csv('rf_submission.csv', index = False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"eee428be91dd21820654f91b9225384b3d7f8527"},"cell_type":"markdown","source":"These predictions score __0.370__ when submitted to the competition."},{"metadata":{"_uuid":"0ffae83f745053c9d244165c24c3ad714419e934"},"cell_type":"markdown","source":"# Feature Selection\n\nOne potential method for improving model performance is feature selection. This is the process where we try to keep only the most useful features for our model. \"Most useful\" can mean many different things, and there are numerous heuristics for selecting the most important features. For feature selection in this notebook, we'll first remove any columns with greater than 0.95 correlation (we already did some of this during feature engineering) and then we'll apply recursive feature elimination with the Scikit-Learn library. \n\nFirst up are the correlations. 0.95 is an arbitrary threshold - feel free to change the values and see how the performance changes! "},{"metadata":{"_uuid":"719a2dcad1caa87ee515e7abcd53a6e55488ce15","trusted":true},"cell_type":"code","source":"train_set = pd.DataFrame(train_set, columns = features)\n\n# Create correlation matrix\ncorr_matrix = train_set.corr()\n\n# Select upper triangle of correlation matrix\nupper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))\n\n# Find index of feature columns with correlation greater than 0.95\nto_drop = [column for column in upper.columns if any(abs(upper[column]) > 0.95)]\n\nto_drop","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f2dbecdc74055adee9ebc41b4c1e2acb19b80133","trusted":true},"cell_type":"code","source":"train_set = train_set.drop(columns = to_drop)\ntrain_set.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e160a64b76cd10632d0bce2f1df28d6a2f687656","trusted":true},"cell_type":"code","source":"test_set = pd.DataFrame(test_set, columns = features)\ntrain_set, test_set = train_set.align(test_set, axis = 1, join = 'inner')\nfeatures = list(train_set.columns)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b90d72e28f73d59887d0f192d25389f5c5876f7e"},"cell_type":"markdown","source":"### Recursive Feature Elimination with Random Forest\n\nThe `RFECV` in Sklearn stands for Recursive Feature Elimination with Cross Validation. The selector operates using a model with feature importances in an iterative manner. At each iteration, it removes either a fraction of features or a set number of features. The iterations continue until the cross validation score no longer improves. \n\nTo create the `selector` object, we pass in the the model, the number of features to remove at each iteration, the cross validation folds, our custom scorer, and any other parameters to guide the selection."},{"metadata":{"_uuid":"e43d0465fa522ba0980943c0a76ae7db0bb71439","trusted":true},"cell_type":"code","source":"from sklearn.feature_selection import RFECV\n\n# Create a model for feature selection\nestimator = RandomForestClassifier(random_state = 10, n_estimators = 100,  n_jobs = -1)\n\n# Create the object\nselector = RFECV(estimator, step = 1, cv = 3, scoring= scorer, n_jobs = -1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"89f75c8dfabf846495fe8b79bfed99924537d7e0"},"cell_type":"markdown","source":"Then we `fit` the selector on the training data as with any other sklearn model. This will continue the feature selection until the cross validation scores no longer improve. "},{"metadata":{"_uuid":"d43cca4b4fc02df588f538e3cd8c3e0e85cc3a4f","trusted":true},"cell_type":"code","source":"selector.fit(train_set, train_labels)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fd03ca2be245681dd58dfefc95e436ceb0a82957"},"cell_type":"markdown","source":"We can investigate the object to see the training scores for each iteration. The following code will plot the validation scores versus the number of features for the training."},{"metadata":{"_uuid":"0227baf12a496ff6c9f25e4d0f1444eb62df5fb4","trusted":true},"cell_type":"code","source":"plt.plot(selector.grid_scores_);\n\nplt.xlabel('Number of Features'); plt.ylabel('Macro F1 Score'); plt.title('Feature Selection Scores');\nselector.n_features_","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"268b17a9c705b6406b0fa72158942d32133ebdbe"},"cell_type":"markdown","source":"We can see that the score improves as we add features up until 96 features. According to the selector, this is the optimal number of features.\n\nThe rankings of each feature can be found by inspecting the trained object. These represent essentially the importance of features averaged over the iterations. Features can share the same ranking, and only features with a rank of 1 are retained."},{"metadata":{"_uuid":"60ac7a6b5b0f7a3056ca7a41b39a63454fc708ce","trusted":true},"cell_type":"code","source":"rankings = pd.DataFrame({'feature': list(train_set.columns), 'rank': list(selector.ranking_)}).sort_values('rank')\nrankings.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e180e020470d68f36fa5f127e004aec9404456fa"},"cell_type":"markdown","source":"Finally, we select the features and then evaluate in cross validation."},{"metadata":{"_uuid":"9d098b899a3aaf81f9f567e21285920ee12ee6e2","trusted":true},"cell_type":"code","source":"train_selected = selector.transform(train_set)\ntest_selected = selector.transform(test_set)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2f58a618719f3e2a3c1c50e74073b5303fb89399","trusted":true},"cell_type":"code","source":"# Convert back to dataframe\nselected_features = train_set.columns[np.where(selector.ranking_==1)]\ntrain_selected = pd.DataFrame(train_selected, columns = selected_features)\ntest_selected = pd.DataFrame(test_selected, columns = selected_features)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3fb2c4635b15a519dc83452fa77a991cf4aa23c9","trusted":true},"cell_type":"code","source":"model_results = cv_model(train_selected, train_labels, model, 'RF-SEL', model_results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e7da5fe95ad404363260a5e34f312be358b9aca2","trusted":true},"cell_type":"code","source":"model_results.set_index('model', inplace = True)\nmodel_results['cv_mean'].plot.bar(color = 'orange', figsize = (8, 6),\n                                  yerr = list(model_results['cv_std']),\n                                 edgecolor = 'k', linewidth = 2)\nplt.title('Model F1 Score Results');\nplt.ylabel('Mean F1 Score (with error bar)');\nmodel_results.reset_index(inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0d6511ebfcab536f17dcd23ec59a7757c696280a"},"cell_type":"markdown","source":"The model with selected features does slightly better in cross-validation. We can keep both sets of features for use in our next model, the Gradient Boosting Machine."},{"metadata":{"_uuid":"22bee5cc58b40acd2badce0579283671ef488ed6"},"cell_type":"markdown","source":"# Upgrading Our Model: Gradient Boosting Machine\n\nAfter using the Random Forest and getting decent scores, it's time to step up and use the gradient boosting machine. If you spend any time on Kaggle, you'll notice that the Gradient Boosting Machine (GBM) wins a high percentage of competitions where the data is structured (in tables) and the datasets are not that large (less than a million observations). \n\nI won't go too much into the details here, but instead will focus on the implementation. We'll use the GBM in LightGBM, although there are also options in Scikit-Learn, XGBOOST, and CatBoost. The first set of hyperparameters we'll use were based on those I've found have worked well for other problems.\n\n### Choosing Number of Estimators with Early Stopping\n\nTo choose the number of estimators (the number of decision trees in the ensemble, called `n_estimators` or `num_boost_rounds`), we'll use early stopping with 5-fold cross validation. This will keep adding estimators until the performance as measured by the Macro F1 Score has not increased for 100 training rounds. To use this metric, we'll have to define a custom metric."},{"metadata":{"_uuid":"956832879bc14aef10d7f676e1148173d4e16a97","trusted":true},"cell_type":"code","source":"def macro_f1_score(labels, predictions):\n    # Reshape the predictions as needed\n    predictions = predictions.reshape(len(np.unique(labels)), -1 ).argmax(axis = 0)\n    \n    metric_value = f1_score(labels, predictions, average = 'macro')\n    \n    # Return is name, value, is_higher_better\n    return 'macro_f1', metric_value, True","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bdfde0db12b0e51f6e31d404126ab8036201f132"},"cell_type":"markdown","source":"# Light Gradient Boosting Machine Implementation\n\nThe function below implements training the gradient boosting machine with Stratified Kfold cross validation and early stopping to prevent overfitting to the training data (although this can still occur). The function performs training with cross validation and records the predictions in probability for each fold. To see how this works, we can return the predictions from each fold and then we'll return a submission to upload to the competition.\n\nChoosing hyperparameters for the Gradient Boosting Machine can be tough and generally is done through model optimization. In this notebook, we'll use a set of hyperparameters that I've found work well on previous problems (although they will not necessarily translate to this competition). \n\nWe set the `n_estimators` to 10000 but we won't actually reach this number because we are using early stopping which will quit training estimators when the cross validation metric does not improve for `early_stopping_rounds`. There's a lot going on in this function, and read through it carefully to make sure you have it all! I've tried to make the comments and code straightforward. (The `display` is used to show custom information during training in combination with `%%capture` so we don't have to see all the LightGBM information during training)."},{"metadata":{"_uuid":"fac36cb103709ca76cf85e874bb2b1c7a378d9a1","trusted":true},"cell_type":"code","source":"from sklearn.model_selection import StratifiedKFold\nimport lightgbm as lgb\nfrom IPython.display import display\n\ndef model_gbm(features, labels, test_features, test_ids, \n              nfolds = 5, return_preds = False, hyp = None):\n    \"\"\"Model using the GBM and cross validation.\n       Trains with early stopping on each fold.\n       Hyperparameters probably need to be tuned.\"\"\"\n    \n    feature_names = list(features.columns)\n\n    # Option for user specified hyperparameters\n    if hyp is not None:\n        # Using early stopping so do not need number of esimators\n        if 'n_estimators' in hyp:\n            del hyp['n_estimators']\n        params = hyp\n    \n    else:\n        # Model hyperparameters\n        params = {'boosting_type': 'dart', \n                  'colsample_bytree': 0.88, \n                  'learning_rate': 0.028, \n                   'min_child_samples': 10, \n                   'num_leaves': 36, 'reg_alpha': 0.76, \n                   'reg_lambda': 0.43, \n                   'subsample_for_bin': 40000, \n                   'subsample': 0.54, \n                   'class_weight': 'balanced'}\n    \n    # Build the model\n    model = lgb.LGBMClassifier(**params, objective = 'multiclass', \n                               n_jobs = -1, n_estimators = 10000,\n                               random_state = 10)\n    \n    # Using stratified kfold cross validation\n    strkfold = StratifiedKFold(n_splits = nfolds, shuffle = True)\n    \n    # Hold all the predictions from each fold\n    predictions = pd.DataFrame()\n    importances = np.zeros(len(feature_names))\n    \n    # Convert to arrays for indexing\n    features = np.array(features)\n    test_features = np.array(test_features)\n    labels = np.array(labels).reshape((-1 ))\n    \n    valid_scores = []\n    \n    # Iterate through the folds\n    for i, (train_indices, valid_indices) in enumerate(strkfold.split(features, labels)):\n        \n        # Dataframe for fold predictions\n        fold_predictions = pd.DataFrame()\n        \n        # Training and validation data\n        X_train = features[train_indices]\n        X_valid = features[valid_indices]\n        y_train = labels[train_indices]\n        y_valid = labels[valid_indices]\n        \n        # Train with early stopping\n        model.fit(X_train, y_train, early_stopping_rounds = 100, \n                  eval_metric = macro_f1_score,\n                  eval_set = [(X_train, y_train), (X_valid, y_valid)],\n                  eval_names = ['train', 'valid'],\n                  verbose = 200)\n        \n        # Record the validation fold score\n        valid_scores.append(model.best_score_['valid']['macro_f1'])\n        \n        # Make predictions from the fold as probabilities\n        fold_probabilitites = model.predict_proba(test_features)\n        \n        # Record each prediction for each class as a separate column\n        for j in range(4):\n            fold_predictions[(j + 1)] = fold_probabilitites[:, j]\n            \n        # Add needed information for predictions \n        fold_predictions['idhogar'] = test_ids\n        fold_predictions['fold'] = (i+1)\n        \n        # Add the predictions as new rows to the existing predictions\n        predictions = predictions.append(fold_predictions)\n        \n        # Feature importances\n        importances += model.feature_importances_ / nfolds   \n        \n        # Display fold information\n        display(f'Fold {i + 1}, Validation Score: {round(valid_scores[i], 5)}, Estimators Trained: {model.best_iteration_}')\n\n    # Feature importances dataframe\n    feature_importances = pd.DataFrame({'feature': feature_names,\n                                        'importance': importances})\n    \n    valid_scores = np.array(valid_scores)\n    display(f'{nfolds} cross validation score: {round(valid_scores.mean(), 5)} with std: {round(valid_scores.std(), 5)}.')\n    \n    # If we want to examine predictions don't average over folds\n    if return_preds:\n        predictions['Target'] = predictions[[1, 2, 3, 4]].idxmax(axis = 1)\n        predictions['confidence'] = predictions[[1, 2, 3, 4]].max(axis = 1)\n        return predictions, feature_importances\n    \n    # Average the predictions over folds\n    predictions = predictions.groupby('idhogar', as_index = False).mean()\n    \n    # Find the class and associated probability\n    predictions['Target'] = predictions[[1, 2, 3, 4]].idxmax(axis = 1)\n    predictions['confidence'] = predictions[[1, 2, 3, 4]].max(axis = 1)\n    predictions = predictions.drop(columns = ['fold'])\n    \n    # Merge with the base to have one prediction for each individual\n    submission = submission_base.merge(predictions[['idhogar', 'Target']], on = 'idhogar', how = 'left').drop(columns = ['idhogar'])\n        \n    # Fill in the individuals that do not have a head of household with 4 since these will not be scored\n    submission['Target'] = submission['Target'].fillna(4).astype(np.int8)\n    \n    # return the submission and feature importances along with validation scores\n    return submission, feature_importances, valid_scores","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"31a574ac08cd9d10be79d86d6ef40fd1b06beb8e"},"cell_type":"markdown","source":"### Cross Validation with Early Stopping Notes\n\nCross validation with early stopping is one of the most effective methods for preventing overfitting on the training set because it prevents us from continuing to add model complexity once it is clear that validation scores are not improving. Repeating this process across multiple folds helps to reduce the bias that comes from using a single fold. Early stopping also lets us train the model much quicker. Overall, __early stopping with cross validation__ is the best method to select the number of estimators in the Gradient Boosting Machine and should be our default technique when we desig an implementation."},{"metadata":{"_uuid":"99c07a92300f43c7df45e7eed1e175b1adc74304","trusted":true},"cell_type":"code","source":"%%capture --no-display\npredictions, gbm_fi = model_gbm(train_set, train_labels, test_set, test_ids, return_preds=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e0ed4f9d0ce7c6c3a5ced03017702e706907b02b"},"cell_type":"markdown","source":"The power of the Gradient Boosting Machine can be seen here! The cross validation score blows away anything we've done previously. \n\nLet's take a look at the predictions to understand what is going on with the predictions in each fold."},{"metadata":{"_uuid":"583d73e53388439db5fdf62f750c045c8338b069","trusted":true},"cell_type":"code","source":"predictions.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cfaa253938c54c139df95ddba27ce260f76c829f"},"cell_type":"markdown","source":"For each fold, the `1, 2, 3, 4` columns represent the probability for each `Target`. The `Target` is the maximum of these with the `confidence` the probability. We have the predictions for all 5 folds, so we can plot the confidence in each `Target` for the different folds."},{"metadata":{"_uuid":"0691d5a568b072f7818e73e7924f959eab8167f2","trusted":true},"cell_type":"code","source":"plt.rcParams['font.size'] = 18\n\n# Kdeplot\ng = sns.FacetGrid(predictions, row = 'fold', hue = 'Target', size = 3, aspect = 4)\ng.map(sns.kdeplot, 'confidence');\ng.add_legend();\n\nplt.suptitle('Distribution of Confidence by Fold and Target', y = 1.05);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"191a6a91c7e030d7c2314b8444b7010391f0ed5e"},"cell_type":"markdown","source":"What we see here is that the confidence for each class if relatively low. It does appear that the model has greater confidence in `Target=4` predictions which makes sense because of the _class imbalance and the high prevalence of this label._ \n\nAnother way to look at the information is as a `violinplot`. This shows the same information, with the number of observations related to the width of the plot."},{"metadata":{"_uuid":"c87e224213269b6ff51c8321389810f33820cec3","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (24, 12))\nsns.violinplot(x = 'Target', y = 'confidence', hue = 'fold', data = predictions);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"90c5e31ebd374d163030bc331ca52f60fe076891"},"cell_type":"markdown","source":"Overall, these results show the issue with imbalanced class problems: our model cannot distinguish very well between the classes that are underrepresented. Later we'll look at predictions themselves and see where our model is \"confused\". For now, we can generate a submission file and submit it to the competition.\n\nWhen we actually make predictions for each household, we average the predictions from each of the folds. Therefore, we are essentially using multiple models since each one is trained on a slightly different fold of the data. The gradient boosting machine is already an ensemble machine learning model, and now we are using it almost as a meta-ensemble by averaging predictions from several gbms. \n\nThis process is shown in the code below."},{"metadata":{"_uuid":"1c2f0e0a3f821c665dc58dd04ee21446483f5954","trusted":true},"cell_type":"code","source":"# Average the predictions over folds\npredictions = predictions.groupby('idhogar', as_index = False).mean()\n\n# Find the class and associated probability\npredictions['Target'] = predictions[[1, 2, 3, 4]].idxmax(axis = 1)\npredictions['confidence'] = predictions[[1, 2, 3, 4]].max(axis = 1)\npredictions = predictions.drop(columns = ['fold'])\n\n# Plot the confidence by each target\nplt.figure(figsize = (10, 6))\nsns.boxplot(x = 'Target', y = 'confidence', data = predictions);\nplt.title('Confidence by Target');\n\nplt.figure(figsize = (10, 6))\nsns.violinplot(x = 'Target', y = 'confidence', data = predictions);\nplt.title('Confidence by Target');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"12baeff0594ef621a9770ff02eb62bcff1e5e979"},"cell_type":"markdown","source":"We can have the function instead return the actual submission file. This takes the average predictions across the five folds, in effectm combining 5 different models, each one trained on a slghtly different subset of the data."},{"metadata":{"_uuid":"a772baa2d12f38ee099e34cc8c65b073aac504fe","trusted":true},"cell_type":"code","source":"%%capture\nsubmission, gbm_fi, valid_scores = model_gbm(train_set, train_labels, \n                                             test_set, test_ids, return_preds=False)\n\nsubmission.to_csv('gbm_baseline.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ea54214d3a7c1a0a75e1260f44dcc2c916ed2a14","trusted":true},"cell_type":"code","source":"_ = plot_feature_importances(gbm_fi, threshold=0.95)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"864b27256383ca6809ef9149bbb9c8adf8f0c785"},"cell_type":"markdown","source":"The gbm seems to think the most important features are those derived from ages. The education variables also show up in the most important features. \n\n### Try Selected Features\n\nThe next step with the LightGBM is to try the features that were selected through recursive feature elimination."},{"metadata":{"_uuid":"f67a7ccddd80db1748340f64e243c263f44cda8d","trusted":true},"cell_type":"code","source":"%%capture --no-display\nsubmission, gbm_fi_selected, valid_scores_selected = model_gbm(train_selected, train_labels, \n                                                               test_selected, test_ids)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fd5389ce1ad7e010b4f79703d0b3eb9ba92e629b","trusted":true},"cell_type":"code","source":"model_results = model_results.append(pd.DataFrame({'model': [\"GBM\", \"GBM_SEL\"], \n                                                   'cv_mean': [valid_scores.mean(), valid_scores_selected.mean()],\n                                                   'cv_std':  [valid_scores.std(), valid_scores_selected.std()]}),\n                                                sort = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2ba52c5cc48753b0ff8926c83ab79afa4b9e6b86","trusted":true},"cell_type":"code","source":"model_results.set_index('model', inplace = True)\nmodel_results['cv_mean'].plot.bar(color = 'orange', figsize = (8, 6),\n                                  yerr = list(model_results['cv_std']),\n                                 edgecolor = 'k', linewidth = 2)\nplt.title('Model F1 Score Results');\nplt.ylabel('Mean F1 Score (with error bar)');\nmodel_results.reset_index(inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"176bf22d90e67d0bb100a47d2aa9fee3c1858d14"},"cell_type":"markdown","source":"The massive advantage of the gradient boosting machine is on display here. For a final step, let's try using 10-folds with both sets and add them to the plot."},{"metadata":{"_uuid":"e22bf3151fc05b4532b75c9596c2479c154e15ec","trusted":true},"cell_type":"code","source":"%%capture\nsubmission, gbm_fi, valid_scores = model_gbm(train_set, train_labels, test_set, test_ids, \n                                             nfolds=10, return_preds=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5b6c0a7da2493bc849b4b6112418ce39340a4a19"},"cell_type":"code","source":"submission.to_csv('gbm_10fold.csv', index = False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c169d4b51cc589d676feeb67b0518a4c4b5a40eb","trusted":true},"cell_type":"code","source":"%%capture\nsubmission, gbm_fi_selected, valid_scores_selected = model_gbm(train_selected, train_labels, test_selected, test_ids,\n                                                               nfolds=10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"eca7c58c81c8f5b640d0898c9ece6b05839b4609","trusted":true},"cell_type":"code","source":"submission.to_csv('gmb_10fold_selected.csv', index = False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4bce239a707127faf14a36a1e4e49b7097fa8ab4","trusted":true},"cell_type":"code","source":"model_results = model_results.append(pd.DataFrame({'model': [\"GBM_10Fold\", \"GBM_10Fold_SEL\"], \n                                                   'cv_mean': [valid_scores.mean(), valid_scores_selected.mean()],\n                                                   'cv_std':  [valid_scores.std(), valid_scores_selected.std()]}),\n                                    sort = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9a475f0430bd525bceb771e120d55cf9ae470225","trusted":true},"cell_type":"code","source":"model_results.set_index('model', inplace = True)\nmodel_results['cv_mean'].plot.bar(color = 'orange', figsize = (8, 6), \n                                  edgecolor = 'k', linewidth = 2,\n                                  yerr = list(model_results['cv_std']))\nplt.title('Model F1 Score Results');\nplt.ylabel('Mean F1 Score (with error bar)');\nmodel_results.reset_index(inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"384d2ca0f26ff075b4ae4b197714ebf5574003bc"},"cell_type":"markdown","source":"The best model seems to be the Gradient Boosting Machine trained with 10 folds on the selected features. This model has not yet been optimized, but we might be able to get a little more performance through optimization."},{"metadata":{"_uuid":"c6747458a39bd71603fc90b142cb1cef796b14ee","trusted":true},"cell_type":"code","source":"print(f\"There are {gbm_fi_selected[gbm_fi_selected['importance'] == 0].shape[0]} features with no importance.\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"96cbfb73dc9bdb8c9cc783e229647da1f1d0ac2a"},"cell_type":"markdown","source":"Well that should make us feel better! All of the features we are using have some importance to the Gradient Boosting Machine. It might be a good idea to go back and retry feature selection but with the GBM since that is the model we are using. \n\nThe next step to take is Model Optimization, the process of getting the most from a machine learning model. "},{"metadata":{"_uuid":"2bf31ba603c98390729788691d28b32b88faba4b"},"cell_type":"markdown","source":"# Model Optimization\n\nModel optimization is the process of extracting the best performance from a machine learning model by tuning the hyperparameters through cross-validation. This is necessary because the best model hyperparameters are different for every dataset. \n\nIn this notebook, we won't get too far into model tuning, but there are multiple options:\n\n1. Manual\n2. Grid Search\n3. Random Search\n4. Automated Optimization\n\nWe will stick to 4. because it is generally the most efficient method and can easily be implemented in a number of libraries, including Hyperopt, which uses a modified version of Bayesian Optimization with the Tree Parzen Estimator. "},{"metadata":{"_uuid":"517131f66e23cb747ac25785188020572a944bb0"},"cell_type":"markdown","source":"# Model Tuning with Hyperopt\n\nBayesian optimization requires 4 parts:\n\n1. Objective function: what we want to maximize (or minimize)\n2. Domain space: region over which to search \n3. Algorithm for choosing next hyperparameters: uses past results to suggest next values\n4. Results history: saves the past results\n\nI've written previously about using Hyperopt, so here we'll stick to the implementation. "},{"metadata":{"_uuid":"97a92863df46523c682fc234e0035ba6e73d4f5b","trusted":true},"cell_type":"code","source":"from hyperopt import hp, tpe, Trials, fmin, STATUS_OK\nfrom hyperopt.pyll.stochastic import sample","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dbc545643ccfac06754c27f2b24c089d3137559e","trusted":true},"cell_type":"code","source":"import csv\nimport ast\nfrom timeit import default_timer as timer","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9c6f77562e1222e3fc92994724de3e8c8d8002ab"},"cell_type":"markdown","source":"### 1. Objective Function\n\nThis takes in the model hyperparameters and returns the associated validation score. Hyperopt needs a score to minimize, so we return 1 - Macro F1 score. "},{"metadata":{"_uuid":"91346af05b05712bf425713a7c30b1a3af31440e","trusted":true},"cell_type":"code","source":"def objective(hyperparameters, nfolds=5):\n    \"\"\"Return validation score from hyperparameters for LightGBM\"\"\"\n    \n    # Keep track of evals\n    global ITERATION\n    ITERATION += 1\n    \n    # Retrieve the subsample\n    subsample = hyperparameters['boosting_type'].get('subsample', 1.0)\n    subsample_freq = hyperparameters['boosting_type'].get('subsample_freq', 0)\n    \n    boosting_type = hyperparameters['boosting_type']['boosting_type']\n    \n    if boosting_type == 'dart':\n        hyperparameters['drop_rate'] = hyperparameters['boosting_type']['drop_rate']\n    \n    # Subsample and subsample frequency to top level keys\n    hyperparameters['subsample'] = subsample\n    hyperparameters['subsample_freq'] = subsample_freq\n    hyperparameters['boosting_type'] = boosting_type\n    \n    # Whether or not to use limit maximum depth\n    if not hyperparameters['limit_max_depth']:\n        hyperparameters['max_depth'] = -1\n    \n    # Make sure parameters that need to be integers are integers\n    for parameter_name in ['max_depth', 'num_leaves', 'subsample_for_bin', \n                           'min_child_samples', 'subsample_freq']:\n        hyperparameters[parameter_name] = int(hyperparameters[parameter_name])\n\n    if 'n_estimators' in hyperparameters:\n        del hyperparameters['n_estimators']\n    \n    # Using stratified kfold cross validation\n    strkfold = StratifiedKFold(n_splits = nfolds, shuffle = True)\n    \n    # Convert to arrays for indexing\n    features = np.array(train_selected)\n    labels = np.array(train_labels).reshape((-1 ))\n    \n    valid_scores = []\n    best_estimators = []\n    run_times = []\n    \n    model = lgb.LGBMClassifier(**hyperparameters, class_weight = 'balanced',\n                               n_jobs=-1, metric = 'None',\n                               n_estimators=10000)\n    \n    # Iterate through the folds\n    for i, (train_indices, valid_indices) in enumerate(strkfold.split(features, labels)):\n        \n        # Training and validation data\n        X_train = features[train_indices]\n        X_valid = features[valid_indices]\n        y_train = labels[train_indices]\n        y_valid = labels[valid_indices]\n        \n        start = timer()\n        # Train with early stopping\n        model.fit(X_train, y_train, early_stopping_rounds = 100, \n                  eval_metric = macro_f1_score, \n                  eval_set = [(X_train, y_train), (X_valid, y_valid)],\n                  eval_names = ['train', 'valid'],\n                  verbose = 400)\n        end = timer()\n        # Record the validation fold score\n        valid_scores.append(model.best_score_['valid']['macro_f1'])\n        best_estimators.append(model.best_iteration_)\n        \n        run_times.append(end - start)\n    \n    score = np.mean(valid_scores)\n    score_std = np.std(valid_scores)\n    loss = 1 - score\n    \n    run_time = np.mean(run_times)\n    run_time_std = np.std(run_times)\n    \n    estimators = int(np.mean(best_estimators))\n    hyperparameters['n_estimators'] = estimators\n    \n    # Write to the csv file ('a' means append)\n    of_connection = open(OUT_FILE, 'a')\n    writer = csv.writer(of_connection)\n    writer.writerow([loss, hyperparameters, ITERATION, run_time, score, score_std])\n    of_connection.close()\n    \n    # Display progress\n    if ITERATION % PROGRESS == 0:\n        display(f'Iteration: {ITERATION}, Current Score: {round(score, 4)}.')\n    \n    return {'loss': loss, 'hyperparameters': hyperparameters, 'iteration': ITERATION,\n            'time': run_time, 'time_std': run_time_std, 'status': STATUS_OK, \n            'score': score, 'score_std': score_std}","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"353308cc0cf7c84d79ff2485cf0b588fda7bdfd8"},"cell_type":"markdown","source":"### 2. Search Space\n\nThe domain is the entire range of values over which we want to search. The only difficult part is the `subsample` ratio which must be set to 1.0 if the `boosting_type=\"goss\"`. Feel free to play around with the values here."},{"metadata":{"_uuid":"e8af56dde243886a20f729a44df1793e93d8bcc3","trusted":true},"cell_type":"code","source":"# Define the search space\nspace = {\n    'boosting_type': hp.choice('boosting_type', \n                              [{'boosting_type': 'gbdt', \n                                'subsample': hp.uniform('gdbt_subsample', 0.5, 1),\n                                'subsample_freq': hp.quniform('gbdt_subsample_freq', 1, 10, 1)}, \n                               {'boosting_type': 'dart', \n                                 'subsample': hp.uniform('dart_subsample', 0.5, 1),\n                                 'subsample_freq': hp.quniform('dart_subsample_freq', 1, 10, 1),\n                                 'drop_rate': hp.uniform('dart_drop_rate', 0.1, 0.5)},\n                                {'boosting_type': 'goss',\n                                 'subsample': 1.0,\n                                 'subsample_freq': 0}]),\n    'limit_max_depth': hp.choice('limit_max_depth', [True, False]),\n    'max_depth': hp.quniform('max_depth', 1, 40, 1),\n    'num_leaves': hp.quniform('num_leaves', 3, 50, 1),\n    'learning_rate': hp.loguniform('learning_rate', \n                                   np.log(0.025), \n                                   np.log(0.25)),\n    'subsample_for_bin': hp.quniform('subsample_for_bin', 2000, 100000, 2000),\n    'min_child_samples': hp.quniform('min_child_samples', 5, 80, 5),\n    'reg_alpha': hp.uniform('reg_alpha', 0.0, 1.0),\n    'reg_lambda': hp.uniform('reg_lambda', 0.0, 1.0),\n    'colsample_bytree': hp.uniform('colsample_by_tree', 0.5, 1.0)\n}","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"18cca2f9585960cddba12012ba5220fcb2cadbc5","trusted":true},"cell_type":"code","source":"sample(space)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"71ac1849fe8587b3df1e50ca8780da819a2be6bd"},"cell_type":"markdown","source":"### 3. Algorithm\n\nThe algorithm for choosing the next values is the Tree Parzen Estimator which uses Bayes rule for constructing a surrogate model of the objective function. Instead of maximizing the objective function, the algorithm maximizes the Expected Improvement (EI) of the surrogate model."},{"metadata":{"_uuid":"4e9a6b0bee169d1957b91423b94ac9376b5f4123","trusted":true},"cell_type":"code","source":"algo = tpe.suggest","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fc37b3f20efd8499eb449c7882b18abf96c7302a"},"cell_type":"markdown","source":"### 4. Results History\n\nWe'll use two different methods for recording results:\n\n1. Trials object: stores everything returned from the objective function\n2. Write to a csv file on every iteration\n\nI like using multiple methods for tracking progress because it means redundancy. One way may fail, but hopefully both will not! The csv file can be used to monitor the method while it is running and the Trials object can be saved and then reloaded to resume optimization."},{"metadata":{"_uuid":"100556e61a1034bd6ee4901cb5fc6c5b6c8c17a6","trusted":true},"cell_type":"code","source":"# Record results\ntrials = Trials()\n\n# Create a file and open a connection\nOUT_FILE = 'optimization.csv'\nof_connection = open(OUT_FILE, 'w')\nwriter = csv.writer(of_connection)\n\nMAX_EVALS = 100\nPROGRESS = 10\nN_FOLDS = 5\nITERATION = 0\n\n# Write column names\nheaders = ['loss', 'hyperparameters', 'iteration', 'runtime', 'score', 'std']\nwriter.writerow(headers)\nof_connection.close()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2aa7ee1bdf700fc885e9af247d4586696598d63a","trusted":true},"cell_type":"code","source":"%%capture --no-display\ndisplay(\"Running Optimization for {} Trials.\".format(MAX_EVALS))\n\n# Run optimization\nbest = fmin(fn = objective, space = space, algo = tpe.suggest, trials = trials,\n            max_evals = MAX_EVALS)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e5637749d4785ff30dff310bdf09caa26b34d24f"},"cell_type":"markdown","source":"To resume training, we can pass in the same trials object and increase the max number of iterations. For later use, the trials can be saved as json."},{"metadata":{"_uuid":"007d7501f06ad53a41e66e0c50d2fe98ca012648","trusted":true},"cell_type":"code","source":"import json\n\n# Save the trial results\nwith open('trials.json', 'w') as f:\n    f.write(json.dumps(str(trials)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8230379a9b7dcc0aa8b0fefc0de307f568877ada"},"cell_type":"markdown","source":"## Using Optimized Model\n\nOnce the optimization procedure has finished, we can use the best results for modeling."},{"metadata":{"_uuid":"32a508bfb26d4e5423662a0046d415928cb31ee6","trusted":true},"cell_type":"code","source":"results = pd.read_csv(OUT_FILE).sort_values('loss', ascending = True).reset_index()\nresults.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6886ffb4c403c32b44dba06b885f1b58119e005c","trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 6))\nsns.regplot('iteration', 'score', data = results);\nplt.title(\"Optimization Scores\");\nplt.xticks(list(range(1, results['iteration'].max() + 1, 3)));","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ce7673c2ec2cdba2da5c2512d5b4240447bb6647","trusted":true},"cell_type":"code","source":"best_hyp = ast.literal_eval(results.loc[0, 'hyperparameters'])\nbest_hyp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ed5e8d02c12ccd9bb5542da34a0c5030973e8106"},"cell_type":"code","source":"%%capture\nsubmission, gbm_fi, valid_scores = model_gbm(train_selected, train_labels, \n                                             test_selected, test_ids, \n                                             nfolds = 10, return_preds=False)\n\nmodel_results = model_results.append(pd.DataFrame({'model': [\"GBM_OPT_10Fold_SEL\"], \n                                                   'cv_mean': [valid_scores.mean()],\n                                                   'cv_std':  [valid_scores.std()]}),\n                                    sort = True).sort_values('cv_mean', ascending = False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"52ff99d3d8e0406a169a5d3f9d1a985bdd79d65b"},"cell_type":"code","source":"%%capture\nsubmission, gbm_fi, valid_scores = model_gbm(train_set, train_labels, \n                                             test_set, test_ids, \n                                             nfolds = 10, return_preds=False)\n\nmodel_results = model_results.append(pd.DataFrame({'model': [\"GBM_OPT_10Fold\"], \n                                                   'cv_mean': [valid_scores.mean()],\n                                                   'cv_std':  [valid_scores.std()]}),\n                                    sort = True).sort_values('cv_mean', ascending = False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2d2366e42ea09f8e4472509d02f34ad34c0aced6","trusted":true},"cell_type":"code","source":"model_results.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c92e262793c2963462e33e2f24884d17dda2a889","trusted":true},"cell_type":"code","source":"submission.to_csv('gbm_opt_10fold_selected.csv', index = False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"09639e1895de282d62dd8c17c4b9a0b0fd61069c"},"cell_type":"markdown","source":"At this point, to improve our performance, we can continue with the optimization, try more feature engineering, stack or ensemble additional models, or look at more experimental methods such as dimension reduction or oversampling. However, the scores we are getting are fairly good and I'm going to move on to a crucial part of machine learning: investigating the predictions to see where our model is wrong. \n\nWhile Kaggle is great for presenting realistic datasets, the methods needed to get to the very top of the leaderboard are not generally used in the real-world. Past some level of performance, accuracy takes a back seat to interpretability. People want to know why a model makes the predictions it does, much as they would want a human to be able to explain her decisions. We could work on squeezing some more performance from our model, but right now, our energy is better spent investigating our model. "},{"metadata":{"trusted":true,"_uuid":"7b5c1246ff4463d634f1e33b58c3c5908b22ee21"},"cell_type":"code","source":"_ = plot_feature_importances(gbm_fi)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"83ae9ce4fe74963bfbe39d60b46c49d64000c883"},"cell_type":"markdown","source":"# Investigate Predictions\n\nAs a first attempt at looking into our model, we can visualize the distribution of predicted labels on the test data. We would expect these to show the same distribution as on the training data. Since we are concerned with household predictions, we'll look at only the predictions for each house and compare with that in the training data.\n\nThe following histrograms are normalize meaning that they show the relative frequency instead of the absolute counts. This is necessary because the raw counts differ in the training and testing data."},{"metadata":{"_uuid":"7f7c53c96da69065d37cb4fb0399337eabf5b0b2","trusted":true},"cell_type":"code","source":"preds = submission_base.merge(submission, on = 'Id', how = 'left')\npreds = pd.DataFrame(preds.groupby('idhogar')['Target'].mean())\n\n# Plot the training labels distribution\nfig, axes = plt.subplots(1, 2, sharey = True, figsize = (12, 6))\nheads['Target'].sort_index().plot.hist(normed = True,\n                                       edgecolor = r'k',\n                                       linewidth = 2,\n                                       ax = axes[0])\n\naxes[0].set_xticks([1, 2, 3, 4]);\naxes[0].set_xticklabels(poverty_mapping.values(), rotation = 60)\naxes[0].set_title('Train Label Distribution')\n\n# Plot the predicted labels\npreds['Target'].sort_index().plot.hist(normed = True, \n                                       edgecolor = 'k',\n                                       linewidth = 2,\n                                       ax = axes[1])\naxes[1].set_xticks([1, 2, 3, 4]);\naxes[1].set_xticklabels(poverty_mapping.values(), rotation = 60)\nplt.subplots_adjust()\nplt.title('Predicted Label Distribution');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7f2246344f5f0515365d6bed470cc79b7f49f4e7","trusted":true},"cell_type":"code","source":"heads['Target'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"440d3eb0052f375f66e0b7135e23c4687558a6e1","trusted":true},"cell_type":"code","source":"preds['Target'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"84e1621be51031a05b3a08db2291aeb220457ce9"},"cell_type":"markdown","source":"The predicted distribution looks close to the training distribution although there are some differences. Depending on the run of the notebook, the results you see may change, but for this edition, the 4s are underrepresented in the predictions and the 3s are overrepresented.\n\nOne potentially method for dealing with imbalanced classification problems is [oversampling the minority class](http://contrib.scikit-learn.org/imbalanced-learn/stable/over_sampling.html), which is easy to do in Python using the imbalanced learn library. We won't explore that option here."},{"metadata":{"_uuid":"660a5bddb4983ed6e0a2f7c8ea42337ac5de49bb"},"cell_type":"markdown","source":"## Validation\n\nFor the test predictions, we can only compare the distribution with that found on the training data. If we want to compare predictions to actual answers, we'll have to split the training data into a separate validation set. We'll use 1000 examples for testing and then we can do operations like make the confusion matrix because we have the right answer."},{"metadata":{"_uuid":"ebf65f2c212515fab1b6458e398e36c17c9edf62","trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\n# Split into validation set\nX_train, X_valid, y_train, y_valid = train_test_split(train_selected,\n                                                      train_labels,\n                                                      test_size = 1000,\n                                                      random_state = 10)\n\n# Create model and train\nmodel = lgb.LGBMClassifier(**best_hyp, \n                           class_weight = 'balanced',\n                           random_state = 10)\nmodel.fit(X_train, y_train);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3b0cc94884fcef395ff18d1d5498f56c69b4cf81","trusted":true},"cell_type":"code","source":"# Make validation predictions\nvalid_preds = model.predict_proba(X_valid)\npreds_df = pd.DataFrame(valid_preds, columns = [1, 2, 3, 4])\n\n# Convert into predictions\npreds_df['prediction'] = preds_df[[1, 2, 3, 4]].idxmax(axis = 1)\npreds_df['confidence'] = preds_df[[1, 2, 3, 4]].max(axis = 1)\n\npreds_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"764f6f9aceb440ce7b40822e7d34819c96bb70c0","trusted":true},"cell_type":"code","source":"print('F1 score:', round(f1_score(y_valid, preds_df['prediction'], average = 'macro'), 5))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0cde5e69be941382f723d0ce620d5dd8616d88a6"},"cell_type":"markdown","source":"Let's use our predictions to plot a confusion matrix. This lets us see where the model is \"confused\" by showing the differences between predictions and true values. The code for this confusion matrix is adapted from: http://scikit-learn.org/stable/auto_examples/model_selection/plot_confusion_matrix.html"},{"metadata":{"_uuid":"67dd7d0504754e01ecafcd9f584a2aa55422feff","trusted":true},"cell_type":"code","source":"from sklearn.metrics import confusion_matrix\nimport itertools\n\ndef plot_confusion_matrix(cm, classes,\n                          normalize=False,\n                          title='Confusion matrix',\n                          cmap=plt.cm.Oranges):\n    \"\"\"\n    This function prints and plots the confusion matrix.\n    Normalization can be applied by setting `normalize=True`.\n    Source: http://scikit-learn.org/stable/auto_examples/model_selection/plot_confusion_matrix.html\n    \"\"\"\n    if normalize:\n        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n        print(\"Normalized confusion matrix\")\n    else:\n        print('Confusion matrix, without normalization')\n\n    print(cm)\n\n    plt.figure(figsize = (10, 10))\n    plt.imshow(cm, interpolation='nearest', cmap=cmap)\n    plt.title(title, size = 24)\n    plt.colorbar(aspect=4)\n    tick_marks = np.arange(len(classes))\n    plt.xticks(tick_marks, classes, rotation=45, size = 14)\n    plt.yticks(tick_marks, classes, size = 14)\n\n    fmt = '.2f' if normalize else 'd'\n    thresh = cm.max() / 2.\n    \n    # Labeling the plot\n    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):\n        plt.text(j, i, format(cm[i, j], fmt), fontsize = 20,\n                 horizontalalignment=\"center\",\n                 color=\"white\" if cm[i, j] > thresh else \"black\")\n        \n    plt.grid(None)\n    plt.tight_layout()\n    plt.ylabel('True label', size = 18)\n    plt.xlabel('Predicted label', size = 18)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"44eae33038e327a6893c287c97c49bad9158de3f","trusted":true},"cell_type":"code","source":"cm = confusion_matrix(y_valid, preds_df['prediction'])\n\nplot_confusion_matrix(cm, classes = ['Extreme', 'Moderate', 'Vulnerable', 'Non-Vulnerable'],\n                      title = 'Poverty Confusion Matrix')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4f18e9f20611fd09919aea8a120820ccf30635ff"},"cell_type":"markdown","source":"Here's how to read the confusion matrix: any of the values on the diagonal, the model got correct because the predicted value matches the true value. Anything not on the diagonal our model got wrong which we can assess by looking at the predicted value versus the actual value. \n\nFor example, our model correctly predicted 25 observations where the poverty was extreme. On the other hand, for 26 cases where the poverty was extreme, our model predicted it was moderate. For 95 cases where the poverty was Non-Vulnerable, our model predicted the poverty was Vulnerable. Overall, we see that our model is only very accurate at idenifying the non-vulnerable households.\n\nTo look at the percentage of each true label predicted in each class, we can normalize the confusion matrix for the true labels."},{"metadata":{"_uuid":"e917ee492b4cb0618f39d4368b160c2ce3cefc57","trusted":true},"cell_type":"code","source":"plot_confusion_matrix(cm, normalize = True,\n                      classes = ['Extreme', 'Moderate', 'Vulnerable', 'Non-Vulnerable'],\n                      title = 'Poverty Confusion Matrix')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b9906be9e3ce7549aa02e23c8b7f4142ad035a24"},"cell_type":"markdown","source":"Now we can see that our model really does not do that well for classes other than Non Vulnerable. It only correctly identifies 15% of the Vulnerable households, classifying more of them as moderate or non vulnerable. Overall, these results show that imbalanced classification problems with relatively few observations are very difficult. There are some methods we can take to try and counter this such as oversampling or training multiple models on different sections of the data, but at the end of the day, the most effective method may be to __gather more data__. "},{"metadata":{"_uuid":"86487bc964d262c6be1144282587d6205a0b1c47"},"cell_type":"markdown","source":"# Dimension Reduction\n\nAs a final exploration of the problem, we can apply a few different dimension reductions methods to the selected data set. These methods can be used for visualization or as a preprocessing method for machine learning. We'll look at four different methods:\n\n1. [PCA: Principal Components Analysis](https://en.wikipedia.org/wiki/Principal_component_analysis). Finds the dimensions of greatest variation in the data\n2. [ICA: Independent Components Analysis](https://en.wikipedia.org/wiki/Independent_component_analysis). Attempts to separate a mutltivariate signal into independent signals. \n3. [TSNE: T-distributed Stochastic Neighbor Embedding](http://www.jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf). Maps high-dimensional data to a low-dimensional manifold attempting to maintain the local structure within the data. It is a non-linear technique and generally only used for visualization.\n4. [UMAP: Uniform Manifold Approximation and Projection](https://arxiv.org/pdf/1802.03426.pdf): A relatively new technique that also maps data to a low-dimensional manifold but tries to preserve more global structure than TSNE. \n\nAll four of these methods are relatively simple to implement in Python. We'll map the selected features down to 3 dimensions for visualization and then also use PCA, ICA, and UMAP as features for modeling (TSNE has no `transform` method and hence cannot be used for preprocessing)."},{"metadata":{"_uuid":"1f2df04033cbf3d7bfcb265766fd2ec8085200fe","trusted":true},"cell_type":"code","source":"from umap import UMAP\nfrom sklearn.decomposition import PCA, FastICA\nfrom sklearn.manifold import TSNE\n\nn_components = 3\n\numap = UMAP(n_components=n_components)\npca = PCA(n_components=n_components)\nica = FastICA(n_components=n_components)\ntsne = TSNE(n_components=n_components)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"556fe32f610cac7bfcd9a73032520db7190c05dc","trusted":true},"cell_type":"code","source":"train_df = train_selected.copy()\ntest_df = test_selected.copy()\n\nfor method, name in zip([umap, pca, ica, tsne], \n                        ['umap', 'pca', 'ica', 'tsne']):\n    \n    # TSNE has no transform method\n    if name == 'tsne':\n        start = timer()\n        reduction = method.fit_transform(train_selected)\n        end = timer()\n    \n    else:\n        start = timer()\n        reduction = method.fit_transform(train_selected)\n        end = timer()\n        \n        test_reduction = method.transform(test_selected)\n    \n        # Add components to test data\n        test_df['%s_c1' % name] = test_reduction[:, 0]\n        test_df['%s_c2' % name] = test_reduction[:, 1]\n        test_df['%s_c3' % name] = test_reduction[:, 2]\n\n    # Add components to training data for visualization and modeling\n    train_df['%s_c1' % name] = reduction[:, 0]\n    train_df['%s_c2' % name] = reduction[:, 1]\n    train_df['%s_c3' % name] = reduction[:, 2]\n    \n    print(f'Method: {name} {round(end - start, 2)} seconds elapsed.')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c1b101fb8825b949f9287881023b02ae551dea7e","trusted":true},"cell_type":"code","source":"from mpl_toolkits.mplot3d import Axes3D\n\ndef discrete_cmap(N, base_cmap=None):\n    \"\"\"Create an N-bin discrete colormap from the specified input map\n    Source: https://gist.github.com/jakevdp/91077b0cae40f8f8244a\"\"\"\n\n    base = plt.cm.get_cmap(base_cmap)\n    color_list = base(np.linspace(0, 1, N))\n    cmap_name = base.name + str(N)\n    return base.from_list(cmap_name, color_list, N)\n\ncmap = discrete_cmap(4, base_cmap = plt.cm.RdYlBu)\n\ntrain_df['label'] = train_labels","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"682812f652af4190fad74d498fbd85cc1288fdd0","scrolled":false,"trusted":true},"cell_type":"code","source":"# Plot each method\nfor method, name in zip([umap, pca, ica, tsne], \n                        ['umap', 'pca', 'ica', 'tsne']):\n    \n    fig = plt.figure(figsize = (8, 8))\n    ax = fig.add_subplot(111, projection='3d')\n    \n    p = ax.scatter(train_df['%s_c1' % name], train_df['%s_c2'  % name], train_df['%s_c3'  % name], \n                   c = train_df['label'].astype(int), cmap = cmap)\n    \n    plt.title(f'{name.capitalize()}', size = 22)\n    fig.colorbar(p, aspect = 4, ticks = [1, 2, 3, 4])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f279a603d6fc6169b301828d95bd37cc05af7883"},"cell_type":"markdown","source":"It's difficult to see much clustering in these plots, which indicates that separating the poverty levels is hard given the available data. Machine learning models are only as effective as the data we give to them, and sometimes the data does not have enough information.\n\nAs a final step, we can train a model using the added reductions from PCA, ICA, and UMAP. We'll use the same Light GBM hyperparameters from the optimization."},{"metadata":{"_uuid":"fbc7f6b837225c6413d7020c2018dce15a450bf2","trusted":true},"cell_type":"code","source":"train_df, test_df = train_df.align(test_df, axis = 1, join = 'inner')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"82308142d2e46857fc624eb0de0d3cb687867295","trusted":true},"cell_type":"code","source":"%%capture\n\nsubmission, gbm_fi, valid_scores = model_gbm(train_df, train_labels, \n                                             test_df, test_ids, nfolds = 10,\n                                             hyp = best_hyp)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3fdd481a721f7389c5235fbf2668f37a1a4d56f2"},"cell_type":"code","source":"submission.to_csv('gbm_opt_10fold_dr.csv', index = False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6abb9055a095b319cf557a8fdde46a00f71a2838","trusted":true},"cell_type":"code","source":"model_results = model_results.append(pd.DataFrame({'model': [\"GBM_OPT_10Fold_DR\"], \n                                                   'cv_mean': [valid_scores.mean()],\n                                                   'cv_std':  [valid_scores.std()]}),\n                                    sort = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5426347cd681954a4075eede1f7edca91a8fca84","trusted":true},"cell_type":"code","source":"model_results = model_results.sort_values('cv_mean')\nmodel_results.set_index('model', inplace = True)\nmodel_results['cv_mean'].plot.bar(color = 'orange', figsize = (10, 8),\n                                  edgecolor = 'k', linewidth = 2,\n                                  yerr = list(model_results['cv_std']))\nplt.title('Model F1 Score Results');\nplt.ylabel('Mean F1 Score (with error bar)');\nmodel_results.reset_index(inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d136251bb48cd913c495a23019323d7ce2565d60"},"cell_type":"markdown","source":"The dimensionality reduction components slightly hurt the overall score of the model. It's possible that they only lead to overfitting on the training data. "},{"metadata":{"_uuid":"14e164da3c26cb2d784cf95c0a7334d5bb808042","trusted":true},"cell_type":"code","source":"_ = plot_feature_importances(gbm_fi)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1ce17ae3e979a16f78b2be6f9cf0f810876bac43"},"cell_type":"markdown","source":"Given the high importance of the dimension reduction components, it looks like they may be causing overfitting. The dimension reduction methods are not using any of the label information, which means they might not necessarily have any useful information for the model."},{"metadata":{"_uuid":"546633035184ff485bdccc277c72bcdef52aca7b"},"cell_type":"markdown","source":"# Visualize Single Decision Tree"},{"metadata":{"_uuid":"b7d236a044ea59dfde8d62ce0e032ff9e132690d"},"cell_type":"markdown","source":"As a final step, we can look at one decision tree in the random forest. First we'll limited the max depth for visability, and then we'll expand the tree all the way.\n\nThe first step is simply to train a random forest and extract one tree (we could also train a single decision tree)."},{"metadata":{"trusted":true,"_uuid":"e33a96fab47ee122cb70e010bc3c0b7555cf8018"},"cell_type":"code","source":"model = RandomForestClassifier(max_depth = 3, n_estimators=10)\nmodel.fit(train_selected, train_labels)\nestimator_limited = model.estimators_[5]\nestimator_limited","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"44ff393ec2f3eb0c1b9860606a9a81429033d708"},"cell_type":"markdown","source":"We take the trained tree and export it as a `.dot` file using `export_graphviz`."},{"metadata":{"trusted":true,"_uuid":"142b90d9228b181f623103168f2ad7870a38de62"},"cell_type":"code","source":"from sklearn.tree import export_graphviz\n\nexport_graphviz(estimator_limited, out_file='tree_limited.dot', feature_names = train_selected.columns,\n                class_names = ['extreme', 'moderate' , 'vulnerable', 'non-vulnerable'],\n                rounded = True, proportion = False, precision = 2, filled = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6efaf01303dfaec39f5eb0601e242043362b44e4"},"cell_type":"markdown","source":"Next we use a system command to convert dot file to png."},{"metadata":{"trusted":true,"_uuid":"19adcf8cee519b2b07b001b97058331a0430d59f"},"cell_type":"code","source":"!dot -Tpng tree_limited.dot -o tree_limited.png","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b64f6ea788a41db0e9206f14f7826cb6ddebf7ab"},"cell_type":"markdown","source":"Finally, we can view the tree in a Jupyter Notebok using `IPython.display` (this is the nicest way to render an image. It also happens to be quite simple."},{"metadata":{"trusted":true,"_uuid":"07ede9900f754d4128b1987b066196e284a58afa"},"cell_type":"code","source":"# Display in jupyter notebook\nfrom IPython.display import Image\nImage(filename = 'tree_limited.png')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"69ee01500eb060a3e85d8563fe3569500fe3c353"},"cell_type":"markdown","source":"#### Visualize Tree with No Maximum Depth\n\nWhen we don't limit the depth, the tree grows extremely deep! It's therefore usually helpful to do some limitation (unless you just like the look of the large tree)! "},{"metadata":{"trusted":true,"_uuid":"0e1e4f152cb6d76d2e0e782551b2de9e1c7519bf"},"cell_type":"code","source":"# No maximum depth\nmodel = RandomForestClassifier(max_depth = None, n_estimators=10)\nmodel.fit(train_selected, train_labels)\nestimator_nonlimited = model.estimators_[5]\n\nexport_graphviz(estimator_nonlimited, out_file='tree_nonlimited.dot', feature_names = train_selected.columns,\n                class_names = ['extreme', 'moderate' , 'vulnerable', 'non-vulnerable'],\n                rounded = True, proportion = False, precision = 2)\n\n!dot -Tpng tree_nonlimited.dot -o tree_nonlimited.png -Gdpi=600","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2f5a4613925b457ef4df71d47265ea12a2488795"},"cell_type":"code","source":"Image(filename = 'tree_nonlimited.png')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"680110aedf38c4954aba2ad91fcfdad4370869d3"},"cell_type":"markdown","source":"# Conclusions\n\nIn this notebook, we went through a step-by-step implementation of an entire data science solution to a real-world problem. Machine learning is really just a series of steps, each simple by themselves, with the overall result often extremely powerful. \n\nOur path was as follows:\n\n1. Understand the problem\n2. Exploratory Data Analysis\n    * Deal with data issues\n    * Fill in missing values\n3. Feature Engineering\n    * Aggregate data \n    * Feature selection in stages\n4. Model Selection\n    * Try many different models to see which one is most promising\n    * Feature selection can also come into play\n5. Model Optimization\n    * Choose the best performing model and tune\n6. Implementing best model\n7. Investigate predictions\n    * Identify model shortcomings\n8. Try new techniques\n    * Experiment and learn! \n  \nAs stated in the introduction, while these steps do have a general order, there is nothing preventing us from repeating them based on subsequent work. For example, we often go back to feature engineering/selection after modeling if we are not satisfied with our performance. After investigating predictions, we might return to the modeling stage and re-think our approach.\n\nIt's also important to keep in mind that machine learning is largely empirical: there are few established best practices so you need to constantly experiement to determine what works best. That's why we included the last step: we got to try out some new methods and even if they don't pay off for this project, they may be useful in the future. As data scientists we need to be constantly learning to do our work more effectively with the best tools possible.\n\nOur final model does well compared to others in the competition, but overall is not extremely accurate. There might be ways to improve performance, but overall, we might not have enough data to achieve exceptional metrics. That's an important point to remember: at the end of the day, the success or failure of a data science project will rest on the quality and quantity of available data. We have solved the problem of classifying household poverty in Porto Rica as well as we can with the available tools, and the next process if we were in charge would be figuring out what more data could be collected to yield a better solution.\n\n## Next Steps\n\nAlthough we tried many possible approaches, there are still some methods that we did not implement but might prove useful:\n\n* Further Hyperparameter tuning: we did not spend much time optimizing the model and there are other packages for optimization that we could try. \n* Further feature selection: we probably do not need to retain all of the features to get the same performance\n* [Oversampling the minority class](http://contrib.scikit-learn.org/imbalanced-learn/stable/over_sampling.html) / undersampling the majority class: an advanced techinque that might improve performance\n* Ensembling or stacking multiple models. We could train models on different sections of the data and then combine their predictions to try and better separate the classes. It might also be worthwhile to include models other than the gradient boosting machine. \n\nThe most valuable part of a Kaggle competition is never the rewards, but the knowledge gained along the way. These contests are great for learning how to solve realistic problems and I approach every new problem as an opportunity to improve as a data scientist. Thanks for reading!\n\nIf you're looking for more on this competition, I have a kernel using automated feature engineering for this problem: [Featuretools for Good](https://www.kaggle.com/willkoehrsen/featuretools-for-good). \n\nBest,\n\nWill \n"},{"metadata":{"_uuid":"c9d93c507f54c2c9e95d3cb82e5221290c089565","trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"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.6.4"}},"nbformat":4,"nbformat_minor":1}