<jupyter_start><jupyter_text>Combined Data This notebook was loaded with:```bashPYSPARK_DRIVER_PYTHON=jupyter PYSPARK_DRIVER_PYTHON_OPTS=notebook ./dse/bin/dse pyspark --num-executors 5 --driver-memory 8g --executor-memory 8g```At this point, we've got several sets of data processed and cleaned. We also have discovered several fields we can use for joining:- license_id- longitude, latitudeLongitude and latitude are great candidates for joining crime, sanitation, weather, and inspections. The problem is that it's not reasonable to expect them to fall on exactly the same coordinate.Suppose we divided the city up into a grid and determined the coordinates for the center of each cell. Then, we could determine which sanitation complaints and crimes were committed in the cell, and connect that to inspections.<jupyter_code>%pylab inline 
from pyspark.sql.functions import lit, col, coalesce, when
from pyspark.ml.classification import RandomForestClassifier, DecisionTreeClassifier, LogisticRegression, GBTClassifier
from pyspark.ml.feature import VectorIndexer, VectorAssembler, StringIndexer
from pyspark.ml import Pipeline
import pyspark

from pyspark.mllib.evaluation import MulticlassMetrics, BinaryClassificationMetrics

import seaborn as sn
import matplotlib.pyplot as plt
import pandas as pd<jupyter_output><empty_output><jupyter_text>Cross Validation Now that we're ready to build a model, let's think about our approach. How do we know if the model will be any good? How do we know it will even be any better than just randomly guessing?No matter which model we choose, we'll need to be able to compare them. There are many measures, but one simple one in classification problems is accuracy. Accuracy is simply the proportion of how predictions were right vs wrong. The approach we'll opt to use is a simple test-train split. That means we use some of our data to train the model, then use the other data we set aside to test its accuracy. We don't want to check accuracy against the training data, becsause if we overfit (meaning we are trained to only be good at pridicting this set of data, and not future data), we'll never know.We'll start by taking our full set and splitting it 75/25<jupyter_code>df_full = sqlContext.read.format('parquet').load('full_set')
df_full.dtypes
df_full = df_full.drop("city_grid").drop("y")
df_train, df_test = df_full.randomSplit([.75, .25], 24)
df_train.count()
df_test.count()
df_train.columns<jupyter_output><empty_output><jupyter_text>Baseline Model As we mentioned above, we still have to ask... sure that number is great... is it any better than just guessing? So, as a baseline model, let's just say we will predict that everyone fails.<jupyter_code>accuracy = df_test.select("y_fail").withColumn("y_pred", lit(1)).filter(col("y_fail") == col("y_pred")).count() / float(19170) * 100
print "Baseline Model Accuracy: %.2f%%" % accuracy<jupyter_output>Baseline Model Accuracy: 20.47%<jupyter_text>If we say that everyone fails, we'll only be right 20.47% of the time.  Preparing the  Model We've got a little work to do to get spark.ml to be able to process our data. First, we'll index the label (our actual pass/fail) that we're training on.<jupyter_code>labelIndexer = StringIndexer(inputCol="y_fail", outputCol="label")<jupyter_output><empty_output><jupyter_text>Next, we'll record our columns. We'll use one for the indexing the features. Remember when we one-hot encoded a few columns? We'll need another list of columns so that we can decode them back (the second array). We could broadcast these, but we'll keep it simple since this training is a one-time thing.<jupyter_code>#broadcast this?
feature_cols = ['canvass',
 'complaint',
 'cumulative_failures',
 'cumulative_inspections',
 'days_since_last_inspection',
 'ever_failed',
 'inspection_type',
 'license_related',
 'liquor',
 'prev_fail',
 'proportion_past_failures',
 'recent_inspection',
 'reinspection',
 'special_event',
 'task_force',
 'monthVec',
 'weekdayVec',
 'riskVec',
 'conditional_approval',
 'wardVec',
 'police_districtVec',
 'crime_count_theft',
 'crime_count_burglary',
 'crime_count_other',
 'sanitation_violation_count',
 'TMAX',
 'TMIN',
 'TMAX_3',
 'TMAX_5']
full_feature_set = [
   'canvass',
 'complaint',
 'cumulative_failures',
 'cumulative_inspections',
 'days_since_last_inspection',
 'ever_failed',
 'inspection_type',
 'license_related',
 'liquor',
 'prev_fail',
 'proportion_past_failures',
 'recent_inspection',
 'reinspection',
 'special_event',
 'task_force',
 'month_jan','month_feb','month_march','month_april','month_may','month_june','month_july','month_august','month_sept','month_oct','month_nov','month_dec',
 'weekday_mon','weekday_tues','weekday_weds','weekday_thurs','weekday_fri','weekday_sat',
 'risk_high','risk_med','risk_low',
 'conditional_approval',
 'ward_1',
 'ward_2',
 'ward_3',
 'ward_4',
 'ward_5',
 'ward_6',
 'ward_7',
 'ward_8',
 'ward_9',
 'ward_10',
 'ward_11',
 'ward_12',
 'ward_13',
 'ward_14',
 'ward_15',
 'ward_16',
 'ward_17',
 'ward_18',
 'ward_19',
 'ward_20',
 'ward_21',
 'ward_22',
 'ward_23',
 'ward_24',
 'ward_25',
 'ward_26',
 'ward_27',
 'ward_28',
 'ward_29',
 'ward_30',
 'ward_31',
 'ward_32',
 'ward_33',
 'ward_34',
 'ward_35',
 'ward_36',
 'ward_37',
 'ward_38',
 'ward_39',
 'ward_40',
 'ward_41',
 'ward_42',
 'ward_43',
 'ward_44',
 'ward_45',
 'ward_46',
 'ward_47',
 'ward_48',
 'ward_49', 'ward_50','district_1',
 'district_2',
 'district_3',
 'district_4',
 'district_5',
 'district_6',
 'district_7',
 'district_8',
 'district_9',
 'district_10',
 'district_11',
 'district_12',
 'district_13',
 'district_14',
 'district_15',
 'district_16',
 'district_17',
 'district_18',
 'district_19',
 'district_20',
 'district_21',
 'district_22',
 'district_23',
 'district_24', 'district_25',
 'crime_count_theft',
 'crime_count_burglary',
 'crime_count_other',
 'sanitation_violation_count',
 'TMAX',
 'TMIN',
 'TMAX_3',
 'TMAX_5' 
]<jupyter_output><empty_output><jupyter_text>Now, we assemble our feature columns into a single Sparse Vector (https://spark.apache.org/docs/1.1.0/mllib-data-types.html)<jupyter_code>featureAssembler = VectorAssembler(inputCols=feature_cols, outputCol='features')<jupyter_output><empty_output><jupyter_text>A pipeline is used to compose our functions. Remember that Spark doesn't execute these until an action is taken... that will be `fit` below<jupyter_code>pipeline = Pipeline(stages=[featureAssembler, labelIndexer])
data = pipeline.fit(df_full).transform(df_full)
data.cache()<jupyter_output><empty_output><jupyter_text>Now, we'll split the data again, using the same random seed, which should give us the same random selections as above.<jupyter_code>train, test = data.randomSplit([0.75,0.25], seed=24) # need to ensure same split for each time (and with the baseline)<jupyter_output><empty_output><jupyter_text>Normally, we really want the labeled classes for our data to be evenly split. 50/50 would be ideal. Let's see what we have.<jupyter_code>train.filter(col("y_fail")==1).count()/float(train.count())
train.filter(col("y_fail")==0).count()/float(train.count())<jupyter_output><empty_output><jupyter_text>So, we've got many more passes than fails. That's reasonable or else the city would have a crisis on its hands. We have a few options. One way is to randomly add duplicates from the failed group until we have an even proportion (called Bagging... see https://en.wikipedia.org/wiki/Bootstrap_aggregating). Another approach is to weight them we'll try that below.<jupyter_code>def reportAccuracy(metrics):
    print metrics.accuracy
    matrix = metrics.confusionMatrix().toArray()
    print matrix
    plt.figure(figsize = (3,3))
    sn.heatmap(matrix, annot=True, fmt="g")
    plt.xlabel("True")
    plt.ylabel("Predicted")
def reportResults(rdd, model):
    fig = plt.figure(figsize=(12, 4))
    ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
    ax2= plt.subplot2grid((1, 5), (0, 3), colspan=3)
    
    metrics = MulticlassMetrics(rdd)

    print metrics.accuracy
    matrix = metrics.confusionMatrix().toArray()
    print matrix
    
    sn.heatmap(matrix, annot=True, fmt="g", ax=ax1)
    ax1.set_xlabel("True")
    ax1.set_ylabel("Predicted")

    # Include only the 10 most important predictors
    df_importance = pd.DataFrame(model.featureImportances.toArray())
    df_importance.columns = ["importance"]
    df_importance["names"] = full_feature_set
    
    df_importance.sort_values("importance", ascending=False, inplace=True)
    df_importance = df_importance.head(10)
    ax2.barh(range(1,11),df_importance["importance"],align='center', color = "darkblue", alpha = 0.6)
    ax2.set_yticklabels(df_importance["names"])
    ax2.set_yticks(range(1,11))
    ax2.set_xlabel('Relative Importance in the Model')
    ax2.set_ylabel('Features')
    ax2.set_title('Relative importance of Each Feature')<jupyter_output><empty_output><jupyter_text>Decision Tree A decision tree is a classification method that builds a decision tree for predicting the label. (https://en.wikipedia.org/wiki/Decision_tree_learning).<jupyter_code>dt = DecisionTreeClassifier(labelCol="label", featuresCol="features")
model = dt.fit(train)
predictions = model.transform(test)
reportResults(predictions.select("prediction", "label").rdd, model)<jupyter_output>0.913145539906
[[ 14808.    438.]
 [  1227.   2697.]]<jupyter_text>As we can see, it does a reasonable job. For the first split, it used `proportion_past_failures`, then `days_since_last_inspection` and so on. You can see why the intuition that beyond 5 levels deep, we really aren't adding much value. That makes sense. If you fail a lot, then you'll probably fail again... that's probably the gut feeling inspectors get. But, why not split on prev_fail first? If the data looked slightly differently, we may have.The matrix on the left is called a Confusion Matrix. It shows us how often we got each classification right and wrong. For example, the lower right is out "false positive" quadrant. It's the number of times we said it was a fail (1) but was actually a pass (0).  Random Forest Classification Random Forest (https://en.wikipedia.org/wiki/Random_forest) addresses that concern by creating many trees, and combining their results in a sort of "voting" process. That makes it less susceptible to overfitting.<jupyter_code>rf = RandomForestClassifier(labelCol="label", featuresCol="features", numTrees=500)
model = rf.fit(train)
predictions = model.transform(test)
reportResults(predictions.select("prediction", "label").rdd, model)<jupyter_output>0.865414710485
[[  1.52450000e+04   1.00000000e+00]
 [  2.57900000e+03   1.34500000e+03]]<jupyter_text>Logistic Regression Logistic Regression doesn't use trees, but rather uses "priors" to determine a probability of being in one of the labeled classes (https://en.wikipedia.org/wiki/Logistic_regression). That might be useful for us if we want to sort them by likelihood of being in a class.<jupyter_code>lr = LogisticRegression(maxIter=10, regParam=0.3, elasticNetParam=0.8)
model = lr.fit(train)
predictions = model.transform(test)
metrics = MulticlassMetrics(predictions.select("prediction", "label").rdd)
reportAccuracy(metrics)<jupyter_output>0.795305164319
[[ 15246.      0.]
 [  3924.      0.]]<jupyter_text>Unfortunately, it performs terribly. It never gets fails right (those are the ones we care about most). That's because our data is imbalanced, and since we have more passes, it's going to lean toward passing. Let's see if we can weight fails more and passes less to get a better prediction.<jupyter_code>count = predictions.count()
l = predictions.select("label").toPandas()
count / (2 * np.bincount(l["label"].astype("int64")))
weighted_train = train.withColumn("logistic_weights", when(col("label")==0.0, 0.25).otherwise(0.75))
lr = LogisticRegression(maxIter=10, regParam=0.3, elasticNetParam=0.8, weightCol="logistic_weights")
model = lr.fit(weighted_train)
predictions = model.transform(test)
metrics = MulticlassMetrics(predictions.select("prediction", "label").rdd)
reportAccuracy(metrics)<jupyter_output>0.875169535733
[[ 15065.    181.]
 [  2212.   1712.]]<jupyter_text>Much better!  Gradient Boosted Trees I'll leave you to read about GBT (https://en.wikipedia.org/wiki/Gradient_boostingGradient_tree_boosting), but it's no accident that I left it for last. It performed best.<jupyter_code>gbt = GBTClassifier(labelCol="label", featuresCol="features", maxIter=100)
gbt_model = gbt.fit(train)
predictions = gbt_model.transform(test)
reportResults(predictions.select("prediction", "label").rdd, gbt_model)<jupyter_output>0.914919144497
[[ 14816.    430.]
 [  1201.   2723.]]<jupyter_text>Choosing our model Purely using accuracy in our case is problematic. We don't want a perfectly accurate model (in fact if it is, you should be worried). What we're trying to accomplish is to determine an inspection order that would reduce the liklihood of people getting sick. That means we want to catch as many failures as soon as we can. That requires a balance. We're not all that worried about over-predicting false negatives, unless they spam our inspection list. GBT did the best overall, but let's take a look at how it got there. It predicted 3924 fails, but missed 430 fails. Weighted logistic regression did just as well at predicting those 3924 fails, but it only missed 181 false negatives. So, even though it had a lower total accuracy, it had a more impactful one for our case... and as a bonus, we can use the prior probability as a sort order for inspections. We can also use proportion of past failures, which we would have to have used in the tree approach.<jupyter_code>model.featureImportances
gbt_model.save("final_model")
model.save("final_model_lr")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11127241053930464, 'fraction_numerical': 0.03179211729694418, 'mean_word_length': 5.753921568627451, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 59, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '522332', 'license': 'CNRI-Python', 'n_tokens_mistral': 4721, 'n_tokens_neox': 4320, 'n_words': 1386, 'path': 'Step8-Modeling.ipynb', 'repo_name': 'angelok1/chicago-inspections-dse', 'size': '13777'}