<jupyter_start><jupyter_text>Lambda School Data Science*Unit 2, Sprint 2, Module 4*---  Classification Metrics- get and interpret the **confusion matrix** for classification models- use classification metrics: **precision, recall**- understand the relationships between precision, recall, **thresholds, and predicted probabilities**, to help **make decisions and allocate budgets**- Get **ROC AUC** (Receiver Operating Characteristic, Area Under the Curve)  SetupRun the code cell below. You can work locally (follow the [local setup instructions](https://lambdaschool.github.io/ds/unit2/local/)) or on Colab.Libraries- category_encoders- ipywidgets- matplotlib- numpy- pandas- scikit-learn- seaborn<jupyter_code>%%capture
import sys

# If you're on Colab:
if 'google.colab' in sys.modules:
    DATA_PATH = 'https://raw.githubusercontent.com/LambdaSchool/DS-Unit-2-Kaggle-Challenge/master/data/'
    !pip install category_encoders==2.*

# If you're working locally:
else:
    DATA_PATH = '../data/'<jupyter_output><empty_output><jupyter_text>Get and interpret the confusion matrix for classification models  Overview First, load the Tanzania Waterpumps data and fit a model. (This code isn't new, we've seen it all before.)<jupyter_code>%matplotlib inline
# import category_encoders as ce
from category_encoders import OneHotEncoder, OrdinalEncoder
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.impute import SimpleImputer
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, plot_confusion_matrix, classification_report
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler


def wrangle(X):
    """Wrangles train, validate, and test sets in the same way"""
    X = X.copy()
    # drop repeat rows(observations)
    X.drop_duplicates(inplace=True)
    # drop repeat cols (features)
#     X.drop(columns=['quantity_group', 'extraction_type_group'], inplace=True)
    # drop constant features
#     X.drop(columns=['recorded_by'], inplace=True)
    # feature engineering
#     X['pump_age'] = X['date_recorded'].dt.year - X['construction_year']
#     X['date_recorded'] = X['date_recorded'].dt.year
    # drop datetime col
#     X.drop(columns='date_recorded', inplace=True)
    # drop high-cardinality
    drop_cols = [col for col in X.select_dtypes('object').columns if X[col].nunique() > 100]
    X.drop(columns=drop_cols, inplace=True)
    # drop features with lots of NaN values
    X.dropna(axis=1, thresh=len(X)*0.8, inplace=True)
    
    df['needs_repair'] = df['status_group'].apply(lambda x: 0 if x == 'functional' else 1)
    df.drop(columns='status_group', inplace=True)

#     # Convert date_recorded to datetime
#     X['date_recorded'] = pd.to_datetime(X['date_recorded'], infer_datetime_format=True)
    
#     # Extract components from date_recorded, then drop the original column
#     X['year_recorded'] = X['date_recorded'].dt.year
#     X['month_recorded'] = X['date_recorded'].dt.month
#     X['day_recorded'] = X['date_recorded'].dt.day
#     X = X.drop(columns='date_recorded')
    
#     # Engineer feature: how many years from construction_year to date_recorded
#     X['years'] = X['year_recorded'] - X['construction_year']    
    
#     # Drop recorded_by (never varies) and id (always varies, random)
#     unusable_variance = ['recorded_by', 'id']
#     X = X.drop(columns=unusable_variance)
    
#     # Drop duplicate columns
#     duplicate_columns = ['quantity_group']
#     X = X.drop(columns=duplicate_columns)
    
#     # About 3% of the time, latitude has small values near zero,
#     # outside Tanzania, so we'll treat these like null values
#     X['latitude'] = X['latitude'].replace(-2e-08, np.nan)
    
#     # When columns have zeros and shouldn't, they are like null values
#     cols_with_zeros = ['construction_year', 'longitude', 'latitude', 'gps_height', 'population']
#     for col in cols_with_zeros:
#         X[col] = X[col].replace(0, np.nan)
        
    return X


# Merge train_features.csv & train_labels.csv
df = pd.merge(pd.read_csv(DATA_PATH+'waterpumps/train_features.csv', parse_dates=['date_recorded'], na_values=[0, -2.000000e-08]), 
                 pd.read_csv(DATA_PATH+'waterpumps/train_labels.csv')).set_index('id')
df = wrangle(train)
# Read test_features.csv & sample_submission.csv
X_test = pd.read_csv(DATA_PATH+'waterpumps/test_features.csv', parse_dates=['date_recorded'], na_values=[0, -2.000000e-08], index_col='id')
X_test = wrangle(test)
# sample_submission = pd.read_csv(DATA_PATH+'waterpumps/sample_submission.csv')

# Split train into train & val. Make val the same size as test.
# target = 'status_group'
# train, val = train_test_split(train, test_size=len(test),  
#                               stratify=train[target], random_state=42)

# # Arrange data into X features matrix and y target vector
# X_train = train.drop(columns=target)
# y_train = train[target]
# X_val = val.drop(columns=target)
# y_val = val[target]
# X_test = test

# # Make pipeline!
# pipeline = make_pipeline(
#     FunctionTransformer(wrangle, validate=False), 
#     ce.OrdinalEncoder(), 
#     SimpleImputer(strategy='mean'), 
#     RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1)
# )

# # Fit on train, score on val
# pipeline.fit(X_train, y_train)
# y_pred = pipeline.predict(X_val)
# print('Validation Accuracy', accuracy_score(y_val, y_pred))<jupyter_output><empty_output><jupyter_text>EDA * How can we transform our target so that this is a binary classification problem?<jupyter_code># print(df.info())
# print(df.shape)
X_test.head()
df.head()
# df['status_group'].value_counts()
# np.where(df['status_group'] == 'functional', 0, 1)
# df['needs_repair'] = df['status_group'].apply(lambda x: 0 if x == 'functional' else 1)<jupyter_output><empty_output><jupyter_text>Split Data<jupyter_code># split tv / fm
target = 'needs_repair'
y = df[target]
X = df.drop(columns=target)

# train-val split
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42)

# sanity check
assert len(X_train) + len(X_val) == len(X)<jupyter_output><empty_output><jupyter_text>Baseline<jupyter_code>print('Baseline Accuracy:', y_train.value_counts(normalize=True).max())<jupyter_output>Baseline Accuracy: 0.5437530959919019<jupyter_text>Build Model* `OrdinalEncoder`* `SimpleImputer`* `RandomForestClassifier`<jupyter_code>model = make_pipeline(
    OrdinalEncoder(),
    SimpleImputer(),
    RandomForestClassifier(
        n_estimators=50,
        n_jobs=-1,
        random_state=42
    )
)
model.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>Check Metrics**Accuracy**<jupyter_code>print('Training Accuracy:', model.score(X_train, y_train))
print('Validation Accuracy:', model.score(X_val, y_val))<jupyter_output>Training Accuracy: 0.9982554758674161
Validation Accuracy: 0.8095279117849759<jupyter_text>**Confusion Matrix**<jupyter_code>plot_confusion_matrix(
    model,
    X_val, # use validation data only
    y_val,
    values_format='.0f',
    display_labels=['no repair needed', 'needs repair']
)<jupyter_output>C:\Users\jeffkang\AppData\Local\Programs\Python\Python39\lib\site-packages\sklearn\utils\deprecation.py:87: FutureWarning: Function plot_confusion_matrix is deprecated; Function `plot_confusion_matrix` is deprecated in 1.0 and will be removed in 1.2. Use one of the class methods: ConfusionMatrixDisplay.from_predictions or ConfusionMatrixDisplay.from_estimator.
  warnings.warn(msg, category=FutureWarning)<jupyter_text>**Recall:** Of those pumps that actually needed repair, what proportion did you correctly predict as needing repair?<jupyter_code>print('Recall Score:', 4042 / (4042 + 1262))<jupyter_output>Recall Score: 0.7620663650075414<jupyter_text>**Precision:** Of all the pumps that you predicted as needing repair, what proportion actually needed repair?<jupyter_code>print('Precision Score:', 4042 / (4042 + 949))<jupyter_output>Precision Score: 0.8098577439390904<jupyter_text>**Classification Report**<jupyter_code>print(classification_report(y_val, model.predict(X_val), target_names=['no repair needed', 'needs repair']))<jupyter_output>precision    recall  f1-score   support

no repair needed       0.81      0.85      0.83      6304
    needs repair       0.81      0.76      0.79      5304

        accuracy                           0.81     11608
       macro avg       0.81      0.81      0.81     11608
    weighted avg       0.81      0.81      0.81     11608<jupyter_text>Case StudyLet's say that it costs the Tanzanian goverment $100 to inspect a water pump, and there is only funding for 2000 pumps inspections<jupyter_code>n_inspections = 2000<jupyter_output><empty_output><jupyter_text>Scenario 1: Choose pumps randomly<jupyter_code>repair_prob = y_val.value_counts(normalize=True).min()
print('Inspections conducted:', n_inspections)
print('Pumps repaired:', n_inspections * (1-repair_prob) * 100)<jupyter_output>Inspections conducted: 2000
Pumps repaired: 108614.74844934528<jupyter_text>Scenario 2: Using our model "out of the box"<jupyter_code>data = {'y_val': y_val, 'y_pred': model.predict(X_val)}

results = pd.DataFrame(data)

mask = results['y_pred'] == 1
sample = results[mask].sample(n_inspections)

sample.head()

print('Inspections conducted:', n_inspections)
print('Pumps repaired:', sample['y_val'].sum())
print('Funds wasted:', (n_inspections - sample.y_val.sum()) * 100)<jupyter_output>Inspections conducted: 2000
Pumps repaired: 1642
Funds wasted: 35800<jupyter_text>Scenario 3: We emphasize **precision** in our model, and select only select pumps that our model is very certain (>0.9) need repair.<jupyter_code>data = {'y_val': y_val, 'y_pred_proba': model.predict_proba(X_val)[:,-1]}
results = pd.DataFrame(data)
display(results.head())

threshold = 0.9

mask = results['y_pred_proba'] > threshold
sample = results[mask].sample(n_inspections)
print('Inspections conducted:',n_inspections)
print('Pumps repaired:',sample['y_val'].sum())
print('Funds wasted:',(n_inspections - sample.y_val.sum()) * 100)<jupyter_output><empty_output><jupyter_text>Interlude: Beware of leakageIf you leave 'status_group' in your feature matrix, you'll have **leakage**  Follow AlongScikit-learn added a [**`plot_confusion_matrix`**](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.plot_confusion_matrix.html) function in version 0.22!  How many correct predictions were made?  How many total predictions were made?  What was the classification accuracy?  Use classification metrics: precision, recall  Overview[Scikit-Learn User Guide — Classification Report](https://scikit-learn.org/stable/modules/model_evaluation.htmlclassification-report)  Wikipedia, [Precision and recall](https://en.wikipedia.org/wiki/Precision_and_recall)> Both precision and recall are based on an understanding and measure of relevance.> Suppose a computer program for recognizing dogs in photographs identifies 8 dogs in a picture containing 12 dogs and some cats. Of the 8 identified as dogs, 5 actually are dogs (true positives), while the rest are cats (false positives). The program's precision is 5/8 while its recall is 5/12.> High precision means that an algorithm returned substantially more relevant results than irrelevant ones, while high recall means that an algorithm returned most of the relevant results.  Follow Along  [We can get precision & recall from the confusion matrix](https://en.wikipedia.org/wiki/Precision_and_recallDefinition_(classification_context))  How many correct predictions of "non functional"?  How many total predictions of "non functional"?  What's the precision for "non functional"?  How many actual "non functional" waterpumps?  What's the recall for "non functional"?  Understand the relationships between precision, recall, thresholds, and predicted probabilities, to help make decisions and allocate budgets  Overview  Imagine this scenario...Suppose there are over 14,000 waterpumps that you _do_ have some information about, but you _don't_ know whether they are currently functional, or functional but need repair, or non-functional.<jupyter_code>len(test)<jupyter_output><empty_output><jupyter_text>**You have the time and resources to go to just 2,000 waterpumps for proactive maintenance.** You want to predict, which 2,000 are most likely non-functional or in need of repair, to help you triage and prioritize your waterpump inspections.You have historical inspection data for over 59,000 other waterpumps, which you'll use to fit your predictive model.<jupyter_code>len(train) + len(val)<jupyter_output><empty_output><jupyter_text>Based on this historical data, if you randomly chose waterpumps to inspect, then about 46% of the waterpumps would need repairs, and 54% would not need repairs.<jupyter_code>y_train.value_counts(normalize=True)
random_inspections = 2000
print(f'With {random_inspections} random inspections, we expect to repair {0.46*random_inspections} waterpumps')<jupyter_output><empty_output><jupyter_text>**Can you do better than random at prioritizing inspections?** In this scenario, we should define our target differently. We want to identify which waterpumps are non-functional _or_ are functional but needs repair:<jupyter_code>y_train = y_train != 'functional'
y_val = y_val != 'functional'
y_train.value_counts(normalize=True)<jupyter_output><empty_output><jupyter_text>We already made our validation set the same size as our test set.<jupyter_code>len(val) == len(test)<jupyter_output><empty_output><jupyter_text>We can refit our model, using the redefined target.Then make predictions for the validation set.<jupyter_code>pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_val)<jupyter_output><empty_output><jupyter_text>Follow Along  Look at the confusion matrix:  How many total predictions of "True" ("non functional" or "functional needs repair") ?  We don't have "budget" to take action on all these predictions- But we can get predicted probabilities, to rank the predictions. - Then change the threshold, to change the number of positive predictions, based on our budget.  Get predicted probabilities and plot the distribution  Change the threshold  Or, get exactly 2,000 positive predictions Identify the 2,000 waterpumps in the validation set with highest predicted probabilities. Most of these top 2,000 waterpumps will be relevant recommendations, meaning `y_val==True`, meaning the waterpump is non-functional or needs repairs.Some of these top 2,000 waterpumps will be irrelevant recommendations, meaning `y_val==False`, meaning the waterpump is functional and does not need repairs.Let's look at a random sample of 50 out of these top 2,000: So how many of our recommendations were relevant? ... What's the precision for this subset of 2,000 predictions?  In this scenario ... Accuracy _isn't_ the best metric!Instead, change the threshold, to change the number of positive predictions, based on the budget. (You have the time and resources to go to just 2,000 waterpumps for proactive maintenance.)Then, evaluate with the precision for "non functional"/"functional needs repair".This is conceptually like **Precision@K**, where k=2,000.Read more here: [Recall and Precision at k for Recommender Systems: Detailed Explanation with examples](https://medium.com/@m_n_malaeb/recall-and-precision-at-k-for-recommender-systems-618483226c54)> Precision at k is the proportion of recommended items in the top-k set that are relevant> Mathematically precision@k is defined as: `Precision@k = ( of recommended items @k that are relevant) / ( of recommended items @k)`> In the context of recommendation systems we are most likely interested in recommending top-N items to the user. So it makes more sense to compute precision and recall metrics in the first N items instead of all the items. Thus the notion of precision and recall at k where k is a user definable integer that is set by the user to match the top-N recommendations objective.We asked, can you do better than random at prioritizing inspections?If we had randomly chosen waterpumps to inspect, we estimate that only 920 waterpumps would be repaired after 2,000 maintenance visits. (46%)But using our predictive model, in the validation set, we succesfully identified over 1,900 waterpumps in need of repair!So we will use this predictive model with the dataset of over 14,000 waterpumps that we _do_ have some information about, but we _don't_ know whether they are currently functional, or functional but need repair, or non-functional.We will predict which 2,000 are most likely non-functional or in need of repair.We estimate that approximately 1,900 waterpumps will be repaired after these 2,000 maintenance visits.So we're confident that our predictive model will help triage and prioritize waterpump inspections.  But ...This metric (~1,900 waterpumps repaired after 2,000 maintenance visits) is specific for _one_ classification problem and _one_ possible trade-off.Can we get an evaluation metric that is generic for _all_ classification problems and _all_ possible trade-offs?Yes — the most common such metric is **ROC AUC.**  Get ROC AUC (Receiver Operating Characteristic, Area Under the Curve)[Wikipedia explains,](https://en.wikipedia.org/wiki/Receiver_operating_characteristic) "A receiver operating characteristic curve, or ROC curve, is a graphical plot that illustrates the diagnostic ability of a binary classifier system as its discrimination threshold is varied. **The ROC curve is created by plotting the true positive rate (TPR) against the false positive rate (FPR) at various threshold settings.**"ROC AUC is the area under the ROC curve. [It can be interpreted](https://stats.stackexchange.com/questions/132777/what-does-auc-stand-for-and-what-is-it) as "the expectation that a uniformly drawn random positive is ranked before a uniformly drawn random negative." ROC AUC measures **how well a classifier ranks predicted probabilities.** So, when you get your classifier’s ROC AUC score, you need to **use predicted probabilities, not discrete predictions.**ROC AUC ranges **from 0 to 1.** Higher is better. A naive majority class **baseline** will have an ROC AUC score of **0.5.**  Scikit-Learn docs- [User Guide: Receiver operating characteristic (ROC)](https://scikit-learn.org/stable/modules/model_evaluation.htmlreceiver-operating-characteristic-roc)- [sklearn.metrics.roc_curve](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_curve.html)- [sklearn.metrics.roc_auc_score](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_auc_score.html) More links- [ROC curves and Area Under the Curve explained](https://www.dataschool.io/roc-curves-and-auc-explained/)- [The philosophical argument for using ROC curves](https://lukeoakdenrayner.wordpress.com/2018/01/07/the-philosophical-argument-for-using-roc-curves/)<jupyter_code># "The ROC curve is created by plotting the true positive rate (TPR) 
# against the false positive rate (FPR) 
# at various threshold settings."

# Use scikit-learn to calculate TPR & FPR at various thresholds
from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(y_val, y_pred_proba)
# See the results in a table
pd.DataFrame({
    'False Positive Rate': fpr, 
    'True Positive Rate': tpr, 
    'Threshold': thresholds
})
# See the results on a plot. 
# This is the "Receiver Operating Characteristic" curve
plt.scatter(fpr, tpr)
plt.title('ROC curve')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate');
# Use scikit-learn to calculate the area under the curve.
from sklearn.metrics import roc_auc_score
roc_auc_score(y_val, y_pred_proba)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09738562091503268, 'fraction_numerical': 0.022473604826546003, 'mean_word_length': 5.350893997445722, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 90, 'https://': 14, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '536334', 'license': 'MIT', 'n_tokens_mistral': 6387, 'n_tokens_neox': 5836, 'n_words': 1936, 'path': 'module4-classification-metrics/LS_DS_224.ipynb', 'repo_name': 'jeffyjkang/DS-Unit-2-Kaggle-Challenge', 'size': '19890'}