<jupyter_start><jupyter_text>Modeling<jupyter_code>import pandas as pd
import numpy as np
import pickle
import os
from datetime import datetime
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

from sklearn.linear_model import LinearRegression, Lasso, LassoCV, Ridge, RidgeCV
from sklearn.model_selection import train_test_split
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn import preprocessing

import statsmodels.api as sm
import statsmodels.formula.api as smf
import scipy.stats as stats
import patsy<jupyter_output><empty_output><jupyter_text>Read pickle:<jupyter_code># read pickles
movies_df = pd.read_pickle('data/movies_df_no_outliers.pkl')
movies_genre_df = pd.read_pickle('data/movies_genre_df_no_outliers.pkl')<jupyter_output><empty_output><jupyter_text>BaselineLook at one feature (you might get lucky and this is THE feature needed to explain your observations). You can see the fit immediately from the seaborn plot. Run statsmodels or sklearn to look at R^2, etc Baseline model & OLS regression summary stats Let's pick production budget as the single feature for the baseline model:<jupyter_code>X = movies_df.loc[:,['production_budget']]
y = movies_df['domestic_gross']
X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                                    test_size=0.2, random_state=1)                       
print(X_train.shape, X_test.shape)
model = sm.OLS(y_train, sm.add_constant(X_train))
results = model.fit()

results.summary()<jupyter_output><empty_output><jupyter_text>Residual plot:<jupyter_code>y_pred = results.predict(sm.add_constant(X_test))

plt.scatter(y_pred,y_test-y_pred,alpha=.1)
plt.plot(np.linspace(0,600000,1000),np.linspace(0,0,1000))
plt.title('Residuals')
plt.xlabel('Predicted Value')
plt.ylabel('Residual');
def diagnostic_plot(x, y):
    plt.figure(figsize=(20,5))
    
    rgr = LinearRegression()
    rgr.fit(x,y)
    pred = rgr.predict(x)

    # Plot fit 
    plt.subplot(1, 3, 1)
    plt.scatter(x,y, alpha=0.1)
    plt.plot(x, pred, color='blue',linewidth=1)
    plt.title("Regression fit")
    plt.xlabel("x")
    plt.ylabel("y")
    
    # Plot residuals
    plt.subplot(1, 3, 2)
    res = y - pred
    plt.scatter(pred, res, alpha=0.1)
    plt.title("Residual plot")
    plt.xlabel("prediction")
    plt.ylabel("residuals")
    
    # Plot Q-Q plot
    plt.subplot(1, 3, 3)
    #Generates a probability plot of sample data against the quantiles of a 
    # specified theoretical distribution 
    stats.probplot(res, dist="norm", plot=plt)
    plt.title("Normal Q-Q plot")
diagnostic_plot(X_test, y_test)<jupyter_output><empty_output><jupyter_text>Appears to be some level of heteroskedasticity per residual plot (middle), with residuals getting larger as the prediction values get larger, giving a cone-shape spread 5-fold cross validation to validate R^2 score:<jupyter_code>lm = LinearRegression()

cross_val_score(lm, X_train, y_train, # estimator, features, target
                cv=5, # number of folds 
                scoring='r2') # scoring metric
print(np.mean(cross_val_score(lm, X_train, y_train, cv=5, scoring='r2')))<jupyter_output>0.3917427199688304<jupyter_text>A model with just the budget feature results in R^2 of 0.391<jupyter_code>#Root Mean Squared Error (RMSE)
from sklearn.metrics import mean_squared_error

def RMSE(actual, predicted):
    return np.sqrt(mean_squared_error(actual,predicted))   
#Mean Absolute Error (MAE)
def mae(y_true, y_pred):
    return np.mean(np.abs(y_pred - y_true))<jupyter_output><empty_output><jupyter_text>Expand modelRun it with a few features (3 or 4), so it is still humanly interpretable. Plot distribution of target variable, errors vs. predicted value to see if non-linearity or heteroskedasticity are an issue. Check fit of the model. Possibly run statsmodels to calculate p-values and see if they correspond to your intuition Expand model to include director power, star power, star appearances and star points, for a total of five features:<jupyter_code>expanded_df = movies_df[['production_budget', 'director_power', 
                         'star_power', 'star_appearances', 'star_points', 
                         'domestic_gross']] 
X = expanded_df.loc[:,['production_budget', 'director_power', 
                      'star_power', 'star_appearances', 'star_points']] 

y = expanded_df['domestic_gross']
X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                                    test_size=0.2, random_state=1)

#X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, 
#                                                  test_size=0.25, random_state=1) # 0.25 x 0.8 = 0.2                                                 
print(X_train.shape, X_test.shape)
model = sm.OLS(y_train, sm.add_constant(X_train))
results = model.fit()

results.summary()
y_pred = results.predict(sm.add_constant(X_test))

plt.scatter(y_pred,y_test-y_pred,alpha=.1)
plt.plot(np.linspace(0,600000,1000),np.linspace(0,0,1000))
plt.title('Residuals')
plt.xlabel('Predicted Value')
plt.ylabel('Residual');<jupyter_output><empty_output><jupyter_text>5-fold cross validation to validate R^2 score:<jupyter_code>lm = LinearRegression()

cross_val_score(lm, X_train, y_train, # estimator, features, target
                cv=5, # number of folds 
                scoring='r2') # scoring metric
print(np.mean(cross_val_score(lm, X_train, y_train, cv=5, scoring='r2')))<jupyter_output>0.4041984336841802<jupyter_text>Adding in those four features increased out R^2 from 0.391 to 0.404 **Now let's add the remaining features, for a total of 9 features:**.<jupyter_code>movies_df.shape
complete_df = movies_df[[
    'runtime', 'year', 'genre_count', 'title_length', 'production_budget',
    'director_power', 'star_power', 'star_appearances', 'star_points',
    'domestic_gross'
]]
X = complete_df.loc[:, [
    'runtime', 'year', 'genre_count', 'title_length', 'production_budget',
    'director_power', 'star_power', 'star_appearances', 'star_points'
]]

y = expanded_df['domestic_gross']
X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                                    test_size=0.2, random_state=1)                       
print(X_train.shape, X_test.shape)
complete_df.shape
model = sm.OLS(y_train, sm.add_constant(X_train))
results = model.fit()

results.summary()<jupyter_output><empty_output><jupyter_text>**Runtime, genre count and title length don't appear to be statistically significant**<jupyter_code>y_pred = results.predict(sm.add_constant(X_test))

plt.scatter(y_pred,y_test-y_pred,alpha=.1)
plt.plot(np.linspace(0,600000,1000),np.linspace(0,0,1000))
plt.title('Residuals')
plt.xlabel('Predicted Value')
plt.ylabel('Residual');<jupyter_output><empty_output><jupyter_text>5-fold cross validation to validate R^2 score:<jupyter_code>lm = LinearRegression()

cross_val_score(lm, X_train, y_train, # estimator, features, target
                cv=5, # number of folds 
                scoring='r2') # scoring metric
print(np.mean(cross_val_score(lm, X_train, y_train, cv=5, scoring='r2')))<jupyter_output>0.40754127244234645<jupyter_text>This model evaluates to an R^2 of 0.407, a 0.003 increase from 0.404   Categorical Features Incorporate categorical features in genre, certificate and month via one-hot encoding:<jupyter_code>movies_matrix = pd.get_dummies(movies_genre_df, prefix=['genre'], columns=['genre'])
movies_matrix = pd.get_dummies(movies_matrix, prefix=['certificate', 'month'], columns=['certificate', 'month'], drop_first=True)
movies_matrix.head(3)
movies_genre_df.shape
movies_matrix.shape<jupyter_output><empty_output><jupyter_text>Aggregate the dataframe to one row per movie by summing the encoded genre features:<jupyter_code>groupby_cols = [
    'runtime',
    'year',
    'domestic_gross',
    'production_budget',
    'genre_count',
    'director',
    'title_length',
    'director_power',
    'star_power',
    'star_points',
    'star_appearances',
    'certificate_NC-17',
    'certificate_Not Rated',  #'certificate_G' 
    'certificate_PG',
    'certificate_PG-13',
    'certificate_R',
    'month_Sep',
    'month_Jul',
    'month_Oct',  #'month_Apr'
    'month_Aug',
    'month_Dec',
    'month_Feb',
    'month_Jan',
    'month_Jun',
    'month_Mar',
    'month_May',
    'month_Nov'
]

encoded_cols = [
    'genre_Action', 'genre_Adventure', 'genre_Animation', 'genre_Biography',
    'genre_Comedy', 'genre_Crime', 'genre_Documentary', 'genre_Drama',
    'genre_Family', 'genre_Fantasy', 'genre_History', 'genre_Horror',
    'genre_Music', 'genre_Musical', 'genre_Mystery', 'genre_News',
    'genre_Romance', 'genre_Sci-Fi', 'genre_Sport', 'genre_Thriller',
    'genre_War', 'genre_Western'
]

movies_matrix = movies_matrix.groupby(
    groupby_cols, as_index=False)[encoded_cols].sum().drop_duplicates(subset=[
        'production_budget', 'domestic_gross', 'title_length', 'director',
        'year', 'runtime'
    ])
movies_matrix.shape
movies_matrix.head()<jupyter_output><empty_output><jupyter_text>Save movies matrix to pickle and csv:<jupyter_code>movies_matrix.info()
movies_matrix.columns
all_df = movies_matrix[[
    'runtime',
    'year',
    'genre_count',
    'title_length',
    'production_budget',
    'domestic_gross',  #director
    'director_power',
    'star_power',
    'star_points',
    'star_appearances',
    'certificate_NC-17',
    'certificate_Not Rated',  #'certificate_G' 
    'certificate_PG',
    'certificate_PG-13',
    'certificate_R',
    'month_Sep',
    'month_Jul',
    'month_Oct',  #'month_Apr'
    'month_Aug',
    'month_Dec',
    'month_Feb',
    'month_Jan',
    'month_Jun',
    'month_Mar',
    'month_May',
    'month_Nov',
    'genre_Action',
    'genre_Adventure',
    'genre_Animation',
    'genre_Biography',
    'genre_Comedy',
    'genre_Crime',
    'genre_Documentary',
    'genre_Drama',
    'genre_Family',
    'genre_Fantasy',
    'genre_History',
    'genre_Horror',
    'genre_Music',
    'genre_Musical',
    'genre_Mystery',
    'genre_News',
    'genre_Romance',
    'genre_Sci-Fi',
    'genre_Sport',
    'genre_Thriller',
    'genre_War',
    'genre_Western'
]]
X = all_df.loc[:, [
    'runtime',
    'year',
    'genre_count',
    'title_length',  #director
    'production_budget',
    'director_power',
    'star_power',
    'star_points',
    'star_appearances',
    'certificate_NC-17',
    'certificate_Not Rated',  #'certificate_G' 
    'certificate_PG',
    'certificate_PG-13',
    'certificate_R',
    'month_Sep',
    'month_Jul',
    'month_Oct',  #'month_Apr'
    'month_Aug',
    'month_Dec',
    'month_Feb',
    'month_Jan',
    'month_Jun',
    'month_Mar',
    'month_May',
    'month_Nov',
    'genre_Action',
    'genre_Adventure',
    'genre_Animation',
    'genre_Biography',
    'genre_Comedy',
    'genre_Crime',
    'genre_Documentary',
    'genre_Drama',
    'genre_Family',
    'genre_Fantasy',
    'genre_History',
    'genre_Horror',
    'genre_Music',
    'genre_Musical',
    'genre_Mystery',
    'genre_News',
    'genre_Romance',
    'genre_Sci-Fi',
    'genre_Sport',
    'genre_Thriller',
    'genre_War',
    'genre_Western'
]]

y = all_df['domestic_gross']
X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                                    test_size=0.2, random_state=1)
print(X_train.shape, X_test.shape)
model = sm.OLS(y_train, sm.add_constant(X_train)) #, X.astype(float))
results = model.fit()

results.summary()<jupyter_output><empty_output><jupyter_text>5-fold cross validation to validate R^2 score:<jupyter_code>lm = LinearRegression()

cross_val_score(lm, X_train, y_train, # estimator, features, target
                cv=5, # number of folds 
                scoring='r2') # scoring metric
print(np.mean(cross_val_score(lm, X_train, y_train, cv=5, scoring='r2')))<jupyter_output>0.43953745973231956<jupyter_text>R^2 increased from 0.407 to 0.439  Challenge modelConsidering the bias vs. variance tradeoff, let's add complexity by adding non-linearity through polynomial features:<jupyter_code>X_train.shape
lm = LinearRegression()
lm.fit(X_train, y_train)
lm.score(X_train, y_train)
poly = PolynomialFeatures(degree=2)
x_train_poly = poly.fit_transform(X_train)
lm.fit(x_train_poly, y_train)
x_test_poly = poly.transform(X_test)
lm.score(x_train_poly, y_train) #score result on train dataset
lm.score(x_test_poly, y_test) #score result on test dataset
cross_val_score(lm, x_train_poly, y_train, # estimator, features, target
                cv=5, # number of folds 
                scoring='r2') # scoring metric
print(np.mean(cross_val_score(lm, x_train_poly, y_train, cv=5, scoring='r2'))) #score results across test sets
p.get_feature_names()
m.coef_
x_train_poly.shape<jupyter_output><empty_output><jupyter_text>RegularizationAdd or remove complexity (using a more sophisticated model or making your model more conservative by using regularization). Compare with unrefined model (error vs. y_pred, R^2, MAE, RMSE). Use regularization to fine tune model (Ridge, Lasso; Lasso may eliminate coefficients that you included at first **Fit the Standard Scaler to the training data:**<jupyter_code>x_train_poly = p.fit_transform(X_train)
std = StandardScaler()
x_train_poly_scale = std.fit_transform(x_train_poly)
x_test_poly_scale = std.transform(x_test_poly) # added
lm.fit(x_train_poly_scale, y_train)
type(X_train), type(x_train_poly_scale)
plt.hist(x_train_poly_scale[:,3])
x_train_poly_scale.shape
x_train_poly_scale<jupyter_output><empty_output><jupyter_text>LASSO<jupyter_code>## Fit a LASSO model on the standardized data
lasso_model = Lasso(alpha = 100)
lasso_model.fit(x_train_poly_scale, y_train)
## Note that now we can meaningful compare the importance of
## different features, since they're on the same scale

## But it's now difficult to interpret the coefficients
## We would need to translate back to the original feature scales by dividing
## each coefficient by the original column's standard deviation

lasso_coef_list = list(zip(p.get_feature_names(input_features=X_train.columns), lasso_model.coef_)) 
lasso_coef_list.sort(key=lambda x:x[1])
lasso_coef_list
# Make predictions on the test set using the new model
test_set_pred = lasso_model.predict(x_test_poly_scale)
# Find the MAE and R^2 on the test set using this model
print(mae(y_test, test_set_pred), r2_score(y_test, test_set_pred))<jupyter_output>42870240.82496922 0.21584775202840223<jupyter_text>**Automated Regularization Strength Tuning via Cross-validation:**<jupyter_code>## Scale the data as before
# std = StandardScaler()
# std.fit(X_train.values)

x_train_poly = poly.fit_transform(X_train)
std = StandardScaler()
x_train_poly_scale = std.fit_transform(x_train_poly)
#m.fit(x_train_poly_scale, y_train)
## Scale the Predictors on both the train and test set
#X_tr = std.transform(X_train.values)
#X_te = std.transform(X_test.values)
x_test_poly_scale = std.transform(x_test_poly)
print(len(x_train_poly_scale), len(x_test_poly_scale))
x_test_poly_scale
# Run the cross validation, find the best alpha, refit the model on all the data with that alpha
alphavec = 10**np.linspace(-2,2,200)

lasso_model = LassoCV(alphas = alphavec, cv=5)
lasso_model.fit(x_train_poly_scale, y_train)
# This is the best alpha value it found - not far from the value selected using simple validation
lasso_model.alpha_
# These are the (standardized) coefficients found when it refit using that best alpha
lasso_coef_listf = list(zip(poly.get_feature_names(input_features=X_train.columns), lasso_model.coef_))
lasso_coef_listf.sort(key=lambda x:x[1])
lasso_coef_listf
# Make predictions on the test set using the new model
test_set_pred = lasso_model.predict(x_test_poly_scale)
# Find the MAE and R^2 on the test set using this model
mae(y_test, test_set_pred)
r2_score(y_test, test_set_pred)<jupyter_output><empty_output><jupyter_text>Ridge **Apply RidgeCV**<jupyter_code>alphas=[1e-3, 1e-2, 1e-1, 1, 10]

for alpha in alphas:
    ridge_model = RidgeCV(alphas=[alpha]).fit(x_train_poly_scale, y_train)
    print(ridge_model.score(x_train_poly_scale, y_train), ridge_model.score(x_test_poly_scale, y_test))
ridge_model.alpha_
test_set_pred = ridge_model.predict(x_test_poly_scale)
mae(y_test, test_set_pred)
r2_score(y_test, test_set_pred)<jupyter_output><empty_output><jupyter_text>Random Forest<jupyter_code>X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, 
                                                  test_size=0.25, random_state=1) # 0.25 x 0.8 = 0.2     
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor

rf = RandomForestRegressor(n_estimators=1800, max_features=3)
rf.fit(X_train, y_train)
rf.score(X_val, y_val)
r2_score(y_val, rf.predict(X_val))
RMSE(rf.predict(X_val), y_val)
mae(rf.predict(X_val), y_val)<jupyter_output><empty_output><jupyter_text>Gradient Boosting<jupyter_code>gbm = GradientBoostingRegressor(n_estimators=1600, max_depth=3, learning_rate=.01)
gbm.fit(X_train, y_train)
gbm.score(X_val, y_val)
r2_score(y_val, gbm.predict(X_val))
RMSE(gbm.predict(X_val),y_val)
mae(gbm.predict(X_val), y_val)
feat_imps = list(zip(X_train.columns, rf.feature_importances_))
feat_imps = sorted(feat_imps, key = lambda x: x[1], reverse=False)
feat_imps = pd.DataFrame(feat_imps, columns=['feature','importance'])
feat_imps[26:46].plot(x='feature',y='importance',kind='barh') 

plt.title('Feature Importance')
plt.ylabel('Feature') 
plt.savefig("feature_importance.png") ;<jupyter_output><empty_output><jupyter_text>Choosing Random Forest Tree as the final model due to the higher R^2 and lower RMSE relative to GBM  Final model Combine back the X_train and X_val, and Y_train and Y_val, into the original 80% training set, to re-train the random forest tree model:<jupyter_code>X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, 
                                                  test_size=0.25, random_state=1) # 0.25 x 0.8 = 0.2     
X_train_final = pd.concat([X_train, X_val], ignore_index=True)
y_train_final = pd.concat([y_train, y_val], ignore_index=True)
rf = RandomForestRegressor(n_estimators=1800, max_features=3)
rf.fit(X_train_final, y_train_final)
rf.score(X_test, y_test)
r2_score(y_test, rf.predict(X_test))
RMSE(rf.predict(X_test), y_test)
mae(rf.predict(X_test), y_test)<jupyter_output><empty_output><jupyter_text>Let's see the results look between actual and predicted:<jupyter_code>y_pred = rf.predict(X_test)

plt.scatter(y_pred, y_test, alpha=.2)
plt.plot(np.linspace(0,600000,1000),np.linspace(0,0,1000))
plt.title('Actual vs Predicted Domestic Gross')
plt.ticklabel_format(style='plain', axis='y')
plt.ticklabel_format(style='plain', axis='x')
plt.xlabel('Predicted Value')
plt.ylabel('Actual Value') ;
y_pred = rf.predict(X_test)

plt.scatter(y_pred, y_test-y_pred, alpha=.2)
plt.plot(np.linspace(0,600000,1000),np.linspace(0,0,1000))
plt.title('Residuals')
plt.xlabel('Predicted Value')
plt.ylabel('Residual');
y_test-y_pred<jupyter_output><empty_output><jupyter_text>Error Analysis<jupyter_code>y_test.head()
X_test.head()<jupyter_output><empty_output><jupyter_text>Merge y_test with X_test:<jupyter_code>merge_test = X_test.merge(y_test.rename('y_test'), left_index=True, right_index=True)
merge_test.head()<jupyter_output><empty_output><jupyter_text>Now merge with y_pred:<jupyter_code>merge_test = merge_test.merge(y_pred.rename('y_pred'), left_index=True, right_index=True) 
merge_test.columns
merge_test.y_test = merge_test.y_test.astype('float64')
merge_test.y_test.unique()
merge_test['y_residual'] = merge_test.y_test - merge_test.y_pred
all_df.iloc[1483]
movies_matrix.iloc[1483]
movies_genre_df.iloc[1483]
merge_test
# merge_test = pd.merge([merge_test, all_df[]], axis=1)

merge_test = pd.merge(merge_test,
                      movies_matrix[[
                          'runtime', 'year', 'title_length',
                          'production_budget', 'director_power', 'director'
                      ]],
                      how='left',
                      left_on=[
                          'runtime', 'year', 'title_length',
                          'production_budget', 'director_power'
                      ],
                      right_on=[
                          'runtime', 'year', 'title_length',
                          'production_budget', 'director_power'
                      ])
merge_test = pd.merge(merge_test,
                      movies_df[[
                          'runtime', 'year', 'title_length',
                          'production_budget', 'director_power', 'director',
                          'movie'
                      ]],
                      how='left',
                      left_on=[
                          'runtime', 'year', 'title_length',
                          'production_budget', 'director_power', 'director',
                      ],
                      right_on=[
                          'runtime', 'year', 'title_length',
                          'production_budget', 'director_power', 'director',
                      ])
merge_test[merge_test.director == '<NAME>']
movies_genre_df[movies_genre_df.director == '<NAME>an']
merge_test.columns
merge_test['absolute_residual'] = merge_test.y_residual.abs()
merge_test[merge_test.year.isin([2018, 2019])].sort_values(by='absolute_residual', ascending=False)
merge_test[merge_test.year.isin([2018, 2019])].sort_values(by='absolute_residual', ascending=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11767401365531778, 'fraction_numerical': 0.019566512395179397, 'mean_word_length': 4.343780607247797, 'pattern_counts': {'":': 0, '<': 106, '<?xml version=': 0, '>': 106, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '537477', 'license': 'MIT', 'n_tokens_mistral': 8095, 'n_tokens_neox': 7695, 'n_words': 1509, 'path': 'movies_modeling.ipynb', 'repo_name': 'syntheticjohn/predict_movie_domestic_gross', 'size': '21843'}