<jupyter_start><jupyter_text>Modeling profitability (binary classification)  Imports<jupyter_code>import pickle
import pandas as pd 
import matplotlib.pyplot as plt 
import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import StratifiedKFold, GridSearchCV
from sklearn.metrics import confusion_matrix, classification_report, roc_curve, roc_auc_score, auc
import itertools
import os

pd.set_option('display.max_rows', 100)
pd.set_option('display.max_columns', 100)
pd.set_option('display.width', 100)
df = pd.read_csv('/Users/ryanrappa/Desktop/dsi/film-profit-prediction/csv_files/clean_data_v2.csv')<jupyter_output><empty_output><jupyter_text>Add'l features: Columns for avg budget and avg "adjusted budget" of competing films<jupyter_code>def get_competitor_avg_budget(df, id_col):
    '''
    Iterates over each film (each row)
    and generates new columns with the average
    revenue, profit, number of movies, and number of
    profitable movies of the casts/crews from competing movies
    (i.e. those released within +/- 2 weeks) have been in.
    --------
    INPUT
    df: dataframe with cast columns
    id_col: column name (str) with film id
    '''
    df['release_date'] = pd.to_datetime(df['release_date'])
    for i in range(len(df)):
        #1. getting a dataframe with competing films
        id_ = df.at[i, id_col]
        release_date = df.at[i, 'release_date']
        df_to_avg = df[(df['release_date'] >= (release_date - pd.to_timedelta(14, unit='d'))) & \
               (df['release_date'] <= (release_date + pd.to_timedelta(14, unit='d'))) & \
               (df['id'] != id_)]
        
        #2. creating new col for competing films' avg budget.
        df.at[i, 'compet_budget'] = df_to_avg['budget'].mean()
    
    return df
df = get_competitor_avg_budget(df, 'id')<jupyter_output><empty_output><jupyter_text>**"Adjusted budget"** is defined as the average of the cast + director's average box office revenue to date, minus the budget of the current film they're part of. In theory this should give us some information about whether the film's budget makes sense considering the past success of the cast and director. We already have this column pertaining to films themselves, now creating it for avg competitors' adjusted budget:<jupyter_code>df['compet_cast_dir_avg_rev'] = (df['compet_cast_rev'] + df['compet_dir_rev'])/2
df['compet_adj_budget'] = df['compet_cast_dir_avg_rev'] - df['compet_budget']<jupyter_output><empty_output><jupyter_text>Dropping non-numerical and unnecessary cols<jupyter_code>df.head()
cols_X = ['budget', 
          'runtime', 
          'releases', 
          'cast_rev', 
          'cast_prof', 
          'cast_films', 
          'cast_prof_films', 
          'dir_rev',
          'dir_prof', 
          'dir_films', 
          'dir_prof_films', 
          'writ_rev', 
          'writ_prof', 
          'writ_films', 
          'writ_prof_films', 
          'adj_budget',
          'cast_dir_avg_rev',
          'fall', 'spring', 'summer', 'winter',
          'Action', 'Adventure', 'Animation', 'Comedy', 'Crime', 'Drama', 'Family', 'Fantasy', 'History', 'Horror', 'Music', 'Mystery', 'None', 'Romance', 'Science Fiction', 'Thriller', 'War', 'Western',
          'compet_cast_rev', 'compet_cast_prof', 'compet_cast_films', 'compet_cast_prof_films', 'compet_dir_rev', 'compet_dir_prof', 'compet_dir_films', 'compet_dir_prof_films', 'compet_writ_rev', 'compet_writ_prof', 'compet_writ_films', 'compet_writ_prof_films', 
          'compet_budget',
          'compet_adj_budget',
          'compet_cast_dir_avg_rev']
col_y = ['made_money']
X = df.loc[:, cols_X].values
y = df.loc[:, col_y].values.ravel() 
#need to use ravel for y array to be correct shape for analysis, 
#otherwise np.mean(y_pred == y_test) does not work properly later
print(X.shape)
print(y.shape)<jupyter_output>(3600, 54)
(3600,)<jupyter_text>Functions to show confusion matrix  source: https://www.kaggle.com/dstuerzer/optimized-logistic-regression<jupyter_code>def plot_confusion_matrix(cm, classes,
                          normalize=False,
                          title='Confusion matrix',
                          cmap=plt.cm.Blues):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.
    """
    
    
    fig, ax = plt.subplots()
    
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
#     plt.title(title)  #my toggle
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=45)
    plt.yticks(tick_marks, classes)

    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        print("Normalized confusion matrix")
#    else:
#        print('Confusion matrix, without normalization')

#    print(cm)

    thresh = cm.max() / 2.
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        plt.text(j, i, cm[i, j],
                 horizontalalignment="center",
                 color="white" if cm[i, j] > thresh else "black")

    plt.tight_layout()
    plt.ylabel('True label', fontsize=14)
    ax.xaxis.tick_top()  #my addition
    plt.xlabel('Predicted label', fontsize=14)
    ax.xaxis.set_label_position('top')  #my addition
    plt.show()

def show_data(cm, print_res = 0):
    tp = cm[1,1]
    fn = cm[1,0]
    fp = cm[0,1]
    tn = cm[0,0]
    if print_res == 1:
        print('Accuracy =     {:.3f}'.format((tp+tn)/(tp+fp+tn+fn)))  #my addition
        print('Precision =     {:.3f}'.format(tp/(tp+fp)))
        print('Recall (TPR) =  {:.3f}'.format(tp/(tp+fn)))
        print('Fallout (FPR) = {:.3f}'.format(fp/(fp+tn)))
    return (tp+tn)/(tp+fp+tn+fn), tp/(tp+fp), tp/(tp+fn), fp/(fp+tn)<jupyter_output><empty_output><jupyter_text>Logistic regression  Data is imbalanced (has mostly profitable films), so using stratified train-test split<jupyter_code>#train-test split with scaled data - avoiding data leakage by scaling train and test sets separately
scaler = StandardScaler()
skf = StratifiedKFold(n_splits = 5, shuffle = True, random_state=313)
for train_index, test_index in skf.split(X, y):
    X_train_scaled, y_train = scaler.fit(X[train_index]).transform(X[train_index]), y[train_index]
    X_test_scaled, y_test = scaler.fit(X[test_index]).transform(X[test_index]), y[test_index]
    break<jupyter_output><empty_output><jupyter_text>Training the model<jupyter_code>lrn = LogisticRegression(penalty='l1')  #l1 penalty performs slightly better than l2 penalty

lrn.fit(X_train_scaled, y_train)
y_pred = lrn.predict(X_test_scaled)
cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)<jupyter_output>/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/logistic.py:433: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)<jupyter_text>Accuracy is the percentage of correct predictions (classifications) made by the model. Precision is the probability that a film that is classified as profitable actually made money. Recall (True Positive Rate) is the probability that a profitable film is classified as such. Fallout (False Positive Rate) is the probability that an unprofitable film is wrongly classified as profitable. **Really important to minimize fallout if we don't want to make unprofitable movies.**<jupyter_code>np.mean(y_pred == y_test) #accuracy checks out<jupyter_output><empty_output><jupyter_text>Parameter tuning - looking at ROC curves & confusion matrices  Functions to plot ROC. Source: https://www.kaggle.com/dstuerzer/optimized-logistic-regression<jupyter_code>def ROC(X, y, c, r):
#makes cross_validation for given parameters c,r. Returns FPR, TPR (averaged)
    dic_weight = {1:len(y)/(r*np.sum(y)), 0:len(y)/(len(y)-r*np.sum(y))} 
    #unfortunately this takes too long to run with l1 penalty, so looking at l2
    lrn = LogisticRegression(penalty = 'l2', C = c, class_weight = dic_weight,
                             solver='liblinear') #specifying default solver to avoid a zillion warnings printing
    
    N = 5      #how much k-fold
    N_iter = 3    #repeat how often (taking the mean)
    mean_tpr = 0.0
    mean_thresh = 0.0
    mean_fpr = np.linspace(0, 1, 50000)
    

    for it in range(N_iter):
        skf = StratifiedKFold(n_splits = N, shuffle = True)
        for train_index, test_index in skf.split(X, y):
            X_train, y_train = X[train_index], y[train_index]
            X_test, y_test = X[test_index], y[test_index]
         
            lrn.fit(X_train, y_train)
            y_prob = lrn.predict_proba(X_test)[:,lrn.classes_[1]]
            
            fpr, tpr, thresholds = roc_curve(y_test, y_prob)
            mean_tpr += np.interp(mean_fpr, fpr, tpr)
            mean_thresh += np.interp(mean_fpr, fpr, thresholds)
            mean_tpr[0] = 0.0

    mean_tpr /= (N*N_iter)
    mean_thresh /= (N*N_iter)
    mean_tpr[-1] = 1.0
    return mean_fpr, mean_tpr, roc_auc_score(y_test, y_prob), mean_thresh

def plot_roc(X,y, list_par_1, par_1 = 'C', par_2 = 1):

    f = plt.figure(figsize = (9,6));
    for p in list_par_1:
        if par_1 == 'C':
            c = p
            r = par_2
        else:
            r = p
            c = par_2
        list_FP, list_TP, AUC, mean_thresh = ROC(X, y, c, r)      
        plt.plot(list_FP, list_TP, label = 'C = {}, r = {}, TPR(3e-4) = {:.4f}, AUC = {:.4f}'.format(c,r,list_TP[10],AUC));
    plt.legend(title = 'values', loc='lower right')
    plt.xlim(0, 1)   #we are only interested in small values of FPR
    plt.ylim(0, 1)
    plt.xlabel('FPR')
    plt.ylabel('TPR')
    plt.title('ROC detail')
    desired_fpr = 0.2
    plt.axvline(desired_fpr, color='b', linestyle='dashed', linewidth=2) #line to show where desired FPR occurs 
    plt.show()
    plt.close()<jupyter_output><empty_output><jupyter_text>Casual/visual grid search 1: Look for optimal class_weight, 'r'<jupyter_code>plot_roc(X,y, [10, 1,.1,.01,.001], 'r', 1)
# looks like r = 1 (default) is best overall
# I think this makes sense because we did stratified train-test split earlier<jupyter_output><empty_output><jupyter_text>Casual/visual grid search 2: Look for optimal 'C' (inverse regularization strength)<jupyter_code>plot_roc(X,y, [0.001, 0.01, 0.1, 1, 10], 'C', 1)
# differences are so small it hardly matters
# running this repeatedly shows no pattern in terms of one C value standing out
# so we'll stick with the default C = 1.0<jupyter_output><empty_output><jupyter_text>Casual/visual grid search 3: checking different decision boundaries - 60, 65, 70, 80%<jupyter_code>y_pred = (lrn.predict_proba(X_test_scaled)[:,1] >= 0.6).astype(bool)

cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)
y_pred = (lrn.predict_proba(X_test_scaled)[:,1] >= 0.7).astype(bool)

cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)
y_pred = (lrn.predict_proba(X_test_scaled)[:,1] >= 0.8).astype(bool)

cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)
y_pred = (lrn.predict_proba(X_test_scaled)[:,1] >= 0.65).astype(bool)

cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)<jupyter_output><empty_output><jupyter_text>**Hard to say which decision boundary is optimal. We want to minimize FPR, but accuracy is getting trashed at higher thresholds. 65-70% boundary looks ok. Time to try other models: SVM, RF, GB...**  Support Vector Classification<jupyter_code>svc = SVC(kernel='sigmoid')  #tried all kernels, SVC just doesn't do well for this

svc.fit(X_train_scaled, y_train)
y_pred = svc.predict(X_test_scaled)
cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)<jupyter_output><empty_output><jupyter_text>Random Forest (RF)<jupyter_code># using scaled data (performs better):
rf = RandomForestClassifier()

rf.fit(X_train_scaled, y_train)
y_pred = rf.predict(X_test_scaled)
cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)<jupyter_output>/anaconda3/lib/python3.6/site-packages/sklearn/ensemble/forest.py:246: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.
  "10 in version 0.20 to 100 in 0.22.", FutureWarning)<jupyter_text>sklearn Gradient Boosting (GB)<jupyter_code># using scaled data (performs slightly better):

gb = GradientBoostingClassifier()

gb.fit(X_train_scaled, y_train)
y_pred = gb.predict(X_test_scaled)
cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)<jupyter_output><empty_output><jupyter_text>**GB does best. RF fallout only marginally lower.** GB feature importances:<jupyter_code>importances = gb.feature_importances_
indices = np.argsort(importances)[::-1]

# Print the feature ranking
print("Feature importance ranking:")

feat_imp_idxs = []
for f in range(X.shape[1]):
#     print("%d. feature %d (%f)" % (f + 1, indices[f], importances[indices[f]]))
    feat_imp_idxs.append(indices[f])
    if importances[indices[f]] < 0.01:
        break
        
count = 0
for i in feat_imp_idxs:
    count += 1
    print(count, "-", list(df.loc[:, cols_X])[i])

# Plot the feature importances of the forest
plt.figure(figsize=(12,3))
plt.title("Feature importances")
plt.bar(range(X.shape[1]), importances[indices],
       color="r", align="center")
plt.xticks(range(X.shape[1]), indices)
plt.xlim([-1, X.shape[1]])
plt.show()
df.loc[:, cols_X].head(1)<jupyter_output><empty_output><jupyter_text>checking different decision boundaries - 70, 75%<jupyter_code>y_pred = (gb.predict_proba(X_test_scaled)[:,1] >= 0.70).astype(bool)

cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)
y_pred = (gb.predict_proba(X_test_scaled)[:,1] >= 0.75).astype(bool)

cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)<jupyter_output><empty_output><jupyter_text>GB parameter tuning  Grid search steps run on EC2:<jupyter_code># parameters = {
#     "learning_rate": [0.1],  #decrease and graph this later
#     "min_samples_split": np.linspace(0.1, 0.5, 5),
#     "min_samples_leaf": np.linspace(0.1, 0.5, 5),
#     "max_depth":[3,5,8],
#     "max_features":["log2","sqrt"],
#     "criterion": ["friedman_mse",  "mae"],
#     "subsample":[0.6, 0.7, 0.8, 0.9, 1.0],
#     "n_estimators":[500]  #decrease and graph this later
#     }
# clf = GridSearchCV(GradientBoostingClassifier(), parameters, cv=10, n_jobs=-1)
## doing this step on EC2, 
## pickling to later inspect best parameters
## then will tune learning rate and n_estimators:

# clf.fit(X_train_scaled, y_train)

# pickle_out = open("gbc.pkl","wb")
# pickle.dump(clf, pickle_out)
# pickle_out.close()<jupyter_output><empty_output><jupyter_text>Checking out optimal parameters, then tuning learning rate and n_estimators as per: https://medium.com/all-things-ai/in-depth-parameter-tuning-for-gradient-boosting-3363992e9bae<jupyter_code>os.getcwd()
pickle_path = '/Users/ryanrappa/Desktop/dsi/film-profit-prediction/grid_search_files/gbc.pkl'
pickle_in = open(pickle_path, "rb")
gbc = pickle.load(pickle_in)<jupyter_output><empty_output><jupyter_text>best params:<jupyter_code>gbc.best_estimator_<jupyter_output><empty_output><jupyter_text>tuning learning rate:<jupyter_code>learning_rates = [1, 0.5, 0.25, 0.1, 0.05, 0.01]
train_results = []
test_results = []
for eta in learning_rates:
    model = GradientBoostingClassifier(learning_rate=eta, criterion='mae', init=None,
              loss='deviance', max_depth=8, max_features='log2',
              max_leaf_nodes=None, min_impurity_decrease=0.0,
              min_impurity_split=None, min_samples_leaf=0.1,
              min_samples_split=0.2, min_weight_fraction_leaf=0.0,
              n_estimators=500, n_iter_no_change=None, presort='auto',
              random_state=None, subsample=0.9, tol=0.0001,
              validation_fraction=0.1, verbose=0, warm_start=False)
    model.fit(X_train_scaled, y_train)
    train_pred = model.predict(X_train_scaled)
    false_positive_rate, true_positive_rate, thresholds = roc_curve(y_train, train_pred)
    roc_auc = auc(false_positive_rate, true_positive_rate)
    train_results.append(roc_auc)
    y_pred = model.predict(X_test_scaled)
    false_positive_rate, true_positive_rate, thresholds = roc_curve(y_test, y_pred)
    roc_auc = auc(false_positive_rate, true_positive_rate)
    test_results.append(roc_auc)
from matplotlib.legend_handler import HandlerLine2D
line1, = plt.plot(learning_rates, train_results, 'b', label='Train AUC')
line2, = plt.plot(learning_rates, test_results, 'r', label='Test AUC')
plt.legend(handler_map={line1: HandlerLine2D(numpoints=2)})
plt.ylabel('AUC score')
plt.xlabel('learning rate')
plt.show()

#looks like 0.1 is best<jupyter_output><empty_output><jupyter_text>tuning n_estimators:<jupyter_code>n_estimators = [1000, 500, 100, 50, 10]
train_results = []
test_results = []
for eta in n_estimators:
    model = GradientBoostingClassifier(learning_rate=0.1, criterion='mae', init=None,
              loss='deviance', max_depth=8, max_features='log2',
              max_leaf_nodes=None, min_impurity_decrease=0.0,
              min_impurity_split=None, min_samples_leaf=0.1,
              min_samples_split=0.2, min_weight_fraction_leaf=0.0,
              n_estimators=eta, n_iter_no_change=None, presort='auto',
              random_state=None, subsample=0.9, tol=0.0001,
              validation_fraction=0.1, verbose=0, warm_start=False)
    model.fit(X_train_scaled, y_train)
    train_pred = model.predict(X_train_scaled)
    false_positive_rate, true_positive_rate, thresholds = roc_curve(y_train, train_pred)
    roc_auc = auc(false_positive_rate, true_positive_rate)
    train_results.append(roc_auc)
    y_pred = model.predict(X_test_scaled)
    false_positive_rate, true_positive_rate, thresholds = roc_curve(y_test, y_pred)
    roc_auc = auc(false_positive_rate, true_positive_rate)
    test_results.append(roc_auc)
from matplotlib.legend_handler import HandlerLine2D
line1, = plt.plot(n_estimators, train_results, 'b', label='Train AUC')
line2, = plt.plot(n_estimators, test_results, 'r', label='Test AUC')
plt.legend(handler_map={line1: HandlerLine2D(numpoints=2)})
plt.ylabel('AUC score')
plt.xlabel('n_estimators')
plt.show()<jupyter_output><empty_output><jupyter_text>results w tuned params:<jupyter_code>final_gbc = gbc.best_estimator_
final_gbc.fit(X_train_scaled, y_train)


y_pred = (final_gbc.predict_proba(X_test_scaled)[:,1] >= 0.70).astype(bool)

cm = confusion_matrix(y_test, y_pred)
if lrn.classes_[0] == 1:
    cm = np.array([[cm[1,1], cm[1,0]], [cm[0,1], cm[0,0]]])

plot_confusion_matrix(cm, ['0', '1'], )
acc, pr, tpr, fpr = show_data(cm, print_res = 1)

#somehow slightly worse than out of the box GB--<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1426321488335466, 'fraction_numerical': 0.02470715621616301, 'mean_word_length': 4.3010696582311505, 'pattern_counts': {'":': 8, '<': 92, '<?xml version=': 0, '>': 99, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, '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': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '616112', 'license': 'MIT', 'n_tokens_mistral': 7892, 'n_tokens_neox': 7504, 'n_words': 1626, 'path': 'notebooks/6_modeling.ipynb', 'repo_name': 'ryanlrappa/film-profit-prediction', 'size': '20318'}