![](\"machine_learning_xkcd.png\")
\n",
"\n",
"Image source: www.xkcd.com\n",
"![](\"ml_overview.png\")
\n",
"\n",
"Image source: www.wordstream.com\n",
"![](\"fruits.jpg\")
\n",
"\n",
"Image source: www.freepik.com"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"![](\"speed_dating.jpg\")
\n",
"\n",
"Image source: www.heartcenteredevents.com\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"data = pd.read_csv('speed_dating.csv', encoding='latin-1')\n",
"\n",
"print(data.shape)\n",
"print(data.columns)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Problem statement"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This dataset is the result of several speed dating events conducted by two professors from Columbia University (Fisman et al., 2006. Gender differences in mate selection: Evidence from a speed dating experiment).\n",
"\n",
"Each row contains a record of one speed dating encounter: the participant's ID, their partner's ID, and information about the participant. \n",
"\n", "The **features** include:\n", "\n", "- socioeconomic status/income\n", "- self-rated personality traits such as ambitiousness, fun-ness and sicnerity\n", "- self-rated attractiveness and importance of attractiveness in a partner\n", "- planned career paths (participants are students)\n", "- etc.\n", "\n", "The question we're trying to answer is:\n", "
\n", "**Given these features, can we predict whether two participants will be interested in one another romantically?**\n", "\n", "We have a binary target variable (=predictand): \n", "\n", "- **1** if both individuals would like to meet again\n", "- **0** if one or both decline\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# What percentage of pairs resulted in a match?\n", "data['match'].value_counts(normalize=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# There were several rounds (\"waves\") of experiments, not all of which used the \n", "# same scoring system. We remove the samples that are based on scores from 1 to 10\n", "data = data[~data['wave'].isin(range(6,12))]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Clean-up" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Income is given as e.g. 67,890.00 --> change this to 67890.\n", "data['income'] = pd.to_numeric(data['income'].str.replace(',', ''))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Remove all columns of type 'object' (mostly strings) - we won't \n", "# use those \n", "data = data.select_dtypes(exclude=['object'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dealing with missing values" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# How many missing values do we have?\n", "data.isnull().sum().sort_values(ascending=False)[:10]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# A crude approach to deal wih missing values would be to just drop all rows that \n", "# contain at least one missing value. If we want that, we can just do\n", "# data.dropna()\n", "\n", "# Or, we could keep only the rows with a certain minimum number of \n", "# valid entries:\n", "# data.dropna(thresh=100)\n", "\n", "# We will drop all rows with invalid entries of 'match' (our predictand) and 'iid'.\n", "data = data.dropna(subset=['iid', 'match'])\n", "\n", "# But since we want to keep as many samples as possible, we are going to replace \n", "# all other missing values with the median\n", "data = data.fillna(data.median())\n", "\n", "# More sophisticated approaches would fill the missing values with values \n", "# estimated by another predictive model - there's in example in the cell below.\n", "\n", "# In general: Which strategy \"makes most sense\" really depends on the dataset!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# from sklearn.neighbors import KNeighborsRegressor\n", "# from sklearn.metrics import r2_score\n", "\n", "# # Get null columns\n", "# null_cols = data.isnull().sum()\n", "# null_cols = list(null_cols[null_cols != 0].index)\n", "\n", "# for i in null_cols:\n", "# x = data.fillna(data.mean()).drop(['match', 'id'], axis=1)\n", "# y = x.pop(i)\n", " \n", "# xtrain, xtest, ytrain, ytest = train_test_split(x, y, test_size=0.2)\n", "# knn = KNeighborsRegressor()\n", "# knn.fit(xtrain, ytrain) \n", "# pred = knn.predict(xtest)\n", "# r2 = r2_score(ytest, pred)\n", " \n", "# # If we can reasonably predict these values, do so\n", "# if r2 > 0.2:\n", "# data['predicted'] = knn.predict(data.fillna(data.median()).drop([i, 'match', 'id'], axis=1))\n", "# data[i] = np.where(data[i].isnull(), data['predicted'], data[i])\n", "# del data['predicted']\n", " \n", "# # Otherwise, just take the median\n", "# else:\n", "# data[i] = data[i].fillna(data[i].median())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding more features " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can construct additional features that may have some predictive skill:\n", "\n", "For example, we can use the average rating of each participant (as determined by how often his or her partners wanted to see them again) as a measure of each person's desirability. \n", "\n", "Or, we can we can determine how \"choosy\" someone is -- were they interested in every person they met, or hardly anyone?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "desirability = data.groupby('iid').mean()['dec_o'].to_dict()\n", "data['desirability'] = data['iid'].map(desirability)\n", "\n", "choosiness = data.groupby('iid').mean()['dec'].to_dict()\n", "data['choosiness'] = data['iid'].map(choosiness)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We might also guess that not just absolute values of attributes matter, but also (or even more so) their differences in the two participants meeting each other." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def get_partner_data(pid, col):\n", " \"\"\"\n", " Looks up the person's partner.\n", " \"\"\"\n", " partner = data[data['iid'] == pid].head(1)[col].iloc[0]\n", " return partner\n", "\n", "# Income (take the log difference)\n", "data['partner_income'] = data['pid'].apply(get_partner_data, col='income')\n", "data['income_difference'] = np.log1p(np.abs(data.income - data.partner_income))\n", "\n", "# Age\n", "data['age_difference'] = data['age'] - data['age_o']" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Compare the each person's overall desirability to their partner's.\n", "data['partner_desirability'] = data['pid'].apply(get_partner_data, col='desirability')\n", "data['des_diff'] = data['desirability'] - data['partner_desirability']\n", "\n", "# Also for Choosiness\n", "data['partner_choosiness'] = data['pid'].apply(get_partner_data, col='choosiness')\n", "data['choose_diff'] = data['choosiness'] - data['partner_choosiness']" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# For the purpose of this tutorial, we will only use a\n", "# subset of the 195 features in order to minimize the \n", "# preprocessing work\n", "\n", "data = data[['order', 'int_corr', 'samerace', 'choosiness', 'choose_diff', 'age',\n", " 'income_difference', 'desirability', 'partner_desirability', \n", " 'age_difference', 'intel', 'go_out', 'match']]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create train-test split" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sklearn.model_selection import train_test_split\n", "\n", "X = data.drop('match', axis=1)\n", "y = data['match']\n", "\n", "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create classifier object and train the classifier (fit the estimator) using the training data\n", "#### A) Simple logistic regression" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sklearn.linear_model import LogisticRegression\n", "\n", "logreg = LogisticRegression(penalty='l2', tol=0.0001, C=1.0, fit_intercept=True, \n", " class_weight=None, random_state=0, solver='liblinear')\n", "\n", "logreg.fit(X_train, y_train)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### B) More fancy approach using logistic regression with cross validation and a pipeline" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sklearn.linear_model import LogisticRegressionCV # includes cross-validation\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.pipeline import make_pipeline\n", "\n", "C_values = [0.05, 0.1, 0.5, 1., 5.]\n", "n_folds = 10\n", "logregCV = LogisticRegressionCV(Cs=C_values, fit_intercept=True, \n", " cv=n_folds, dual=False, penalty='l2', \n", " scoring='roc_auc', solver='liblinear', \n", " tol=0.0001, max_iter=1000, class_weight=None, \n", " random_state=0)\n", "\n", "logregCV_pipe = make_pipeline(StandardScaler(), logregCV)\n", "logregCV_pipe.fit(X_train, y_train)\n", "print(logregCV_pipe.get_params().keys())" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "coef_B = logregCV_pipe.named_steps['logisticregressioncv'].coef_\n", "\n", "plt.figure(figsize=(10,8))\n", "coef_logregCV_pipe = pd.Series(coef_B.ravel(),\n", " index=X_train.columns) \n", "coef_logregCV_pipe.sort_values().plot(kind=\"barh\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Evaluate the performance of the model on the test data\n", "\n", "There are lots of different metrics available to evaluate the performance of your model: http://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from sklearn.metrics import roc_curve, classification_report, \\\n", " confusion_matrix, brier_score_loss" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "y_pred_A = logreg.predict(X_test) # the predicted classes (\"1\" for match, \"0\" for no match)\n", "y_prob_A = logreg.predict_proba(X_test)[:,1] # the predicted probabilities of class 1\n", "\n", "# Calculate performance metrics\n", "n_TN, n_FP, n_FN, n_TP = confusion_matrix(y_test, y_pred_A).ravel()\n", "\n", "# Matthews correlation coefficient\n", "MCC = (n_TP * n_TN - n_FP * n_FN) / np.sqrt((n_TP + n_FP) * (n_TP + n_FN) * (n_TN + n_FP) * (n_TN + n_FN))\n", "if np.isnan(MCC): \n", " MCC = 0\n", " \n", "# Brier score loss (the smaller, the better)\n", "BSL = brier_score_loss(y_test, y_prob_A)\n", " \n", "# Accuracy\n", "accuracy = float((n_TP + n_TN)) / (n_TP + n_FP + n_TN + n_FN)\n", "\n", "n_test = float(len(y_pred_A))\n", "print(\"TP: {0} ({1:2.1f}%)\".format(n_TP, n_TP / n_test * 100.))\n", "print(\"TN: {0} ({1:2.1f}%)\".format(n_TN, n_TN / n_test * 100.))\n", "print(\"FP: {0} ({1:2.1f}%)\".format(n_FP, n_FP / n_test * 100.))\n", "print(\"FN: {0} ({1:2.1f}%)\".format(n_FN, n_FN / n_test * 100.))\n", "print(\"\\nAccuracy: {0:2.2f}\".format(accuracy))\n", "print(\"Matthews correlation coefficient: {0:2.2f}\".format(MCC))\n", "print(\"Brier score loss: {0:2.4f}\".format(BSL))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "y_pred_B = logregCV_pipe.predict(X_test)\n", "y_prob_B = logregCV_pipe.predict_proba(X_test)[:,1]\n", "\n", "# Calculate performance metrics\n", "n_TN, n_FP, n_FN, n_TP = confusion_matrix(y_test, y_pred_B).ravel()\n", "\n", "# Matthews correlation coefficient\n", "MCC = (n_TP * n_TN - n_FP * n_FN) / np.sqrt((n_TP + n_FP) * (n_TP + n_FN) * (n_TN + n_FP) * (n_TN + n_FN))\n", "if np.isnan(MCC): \n", " MCC = 0\n", " \n", "BSL = brier_score_loss(y_test, y_prob_B)\n", "\n", "n_test = float(len(y_pred_B))\n", "print(\"TP: {0} ({1:2.1f}%)\".format(n_TP, n_TP / n_test * 100.))\n", "print(\"TN: {0} ({1:2.1f}%)\".format(n_TN, n_TN / n_test * 100.))\n", "print(\"FP: {0} ({1:2.1f}%)\".format(n_FP, n_FP / n_test * 100.))\n", "print(\"FN: {0} ({1:2.1f}%)\".format(n_FN, n_FN / n_test * 100.))\n", "print(\"\\nAccuracy: {0:2.2f}\".format(accuracy))\n", "print(\"Matthews correlation coefficient: {0:2.2f}\".format(MCC))\n", "print(\"Brier score loss: {0:2.4f}\".format(BSL))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "** In this example, using cross-validation to optimize the regularization strength only slightly increases the performance. However, applying cross-validation, scaling features, choosing different \"penalties\" in the regularization ('l1' vs. 'l2') etc. can greatly improve your model! **" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 3) From Love to Wine: A Regression Example\n", "\n", "Our goal is to predict wine quality." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "
![](\"wine_quality.jpg\")
\n",
"\n",
"Wine quality control in Adega Cooperativa de Borba winery, Alentejo, Portugal. Image by Manuel Ribeiro.\n",
"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Load dataset" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dataset_url = \"http://mlr.cs.umass.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv\"\n", "data = pd.read_csv(dataset_url, sep=';')\n", " \n", "print(data.head())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create train-test split" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "y = data['quality'] \n", "X = data.drop('quality', axis=1) \n", "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create classifier object and train the classifier (fit the estimator) using the training data\n", "\n", "#### A) Linear regression" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Sometimes (e.g. when feature values have very different ranges) it can be\n", "# useful to scale or normalize the training and test sets before using them\n", "# in a model, so let's do that here:\n", "scaler = StandardScaler()\n", "X_train_scaled = scaler.fit_transform(X_train)\n", "X_test_scaled = scaler.transform(X_test)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**!! A note about data leakage !!**\n", "\n", "The proper way to scale the data is *not* to scale the entire dataset with a single transform, since this will indirectly leak information into the training data about the whole dataset, including the test data. Instead, \n", "you should fit the scaler only to the training data, transform the training data, and then transform the test data with that same scaler.\n", "\n", "**\"Contaminated\" training data can lead to over-optimistic performance estimations on your testing data!**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sklearn.linear_model import LinearRegression\n", "from sklearn.metrics import explained_variance_score, r2_score\n", "\n", "linreg = LinearRegression(fit_intercept=True)\n", "linreg.fit(X_train, y_train)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = plt.figure(figsize=(10,8))\n", "coef_linreg = pd.Series(linreg.coef_.ravel(),\n", " index=X_train.columns) \n", "coef_linreg.sort_values().plot(kind=\"barh\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Evaluate the Linear regression model" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "y_pred = linreg.predict(X_test)\n", "expl_var = explained_variance_score(y_test, y_pred)\n", "r2 = r2_score(y_test, y_pred)\n", "print(\"*********\")\n", "print(\"Performance of linear regression:\")\n", "print(\"*********\")\n", "print(\"\\texplained variance: {0:2.3f}\".format(expl_var))\n", "print(\"\\tR2 score: {0:2.3f}\".format(r2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### B) Random Forest regression\n", "\n", "Random forest is an ensemble tool which takes a subset of observations and a subset of variables to build decision trees. It builds multiple such decision trees and combines them to get a more accurate and stable prediction.\n", "\n", "The model has a couple of levers, i.e., parameters we can tune to improve performance. Such parameters that define the \"higher-level\" structure of a model are called **hyperparameters**. \n", "\n", "E.g., in Random Forest regression, hyperparameters include the number of individual decision trees, the maximum number of features used in each of these trees, or their maximum depth. We can explore the hyperparameter space using a **grid search**, where we try different combinations of hyperparameters and then use the best combination for our final model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
![](\"random_forest.jpg\")
\n",
"\n",
"Image source: www.dataaspirant.com"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from sklearn.ensemble import RandomForestRegressor\n",
"from sklearn.model_selection import GridSearchCV\n",
"\n",
"\n",
"# Declare data preprocessing steps\n",
"pipeline = make_pipeline(StandardScaler(), \n",
" RandomForestRegressor())\n",
" \n",
"# Define hyperparameters to tune\n",
"hyperparameters = {'randomforestregressor__max_features' : ['auto', 'sqrt', 'log2'],\n",
" 'randomforestregressor__max_depth': [None, 5, 3, 1],\n",
" 'randomforestregressor__n_estimators': [50, 100]}\n",
" \n",
"# Tune model using cross-validation pipeline\n",
"forest_cv = GridSearchCV(pipeline, hyperparameters, cv=10)\n",
"forest_cv.fit(X_train, y_train)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Evaluate the Random Forest model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"y_pred = forest_cv.predict(X_test)\n",
"expl_var = explained_variance_score(y_test, y_pred)\n",
"r2 = r2_score(y_test, y_pred)\n",
"\n",
"print(\"*********\")\n",
"print(\"Performance of Random Forest regression:\")\n",
"print(\"*********\")\n",
"print(\"\\texplained variance: {0:2.3f}\".format(expl_var))\n",
"print(\"\\tR2 score: {0:2.3f}\".format(r2))"
]
}
],
"metadata": {
"anaconda-cloud": {},
"gist_info": {
"gist_id": null,
"gist_url": null
},
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}