# Install required packages
!pip install --upgrade pip
!pip install -q -r requirements.txt
# sklearn modules
import sklearn
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error
from sklearn.pipeline import make_pipeline
from sklearn.compose import ColumnTransformer, make_column_transformer
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV, KFold
from sklearn.linear_model import Ridge
from sklearn.linear_model import RidgeCV
from sklearn.linear_model import Lasso
# pandas, numpy, and plotting
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as st
from scipy.stats import norm
# schrutepy
from schrutepy import schrutepy
np.random.seed(44)
# Display plots inline
%matplotlib inline
# Plotting defaults
plt.rcParams['figure.figsize'] = (8,5)
plt.rcParams['figure.dpi'] = 80
df = schrutepy.load_schrute()
#Combining number of reviews and average review rating, deleting previous columns
d = pd.read_csv("the_office.csv")
tv = d["total_votes"].mean()
d["Episode Scor"] = 0.7*d["imdb_rating"] + 10*(0.3)*(1-np.exp((d["total_votes"])/np.log(1/2)*tv))
d['Episode Score'] = d['Episode Scor'].apply(lambda x: round(x, 2))
del d['total_votes']
del d['imdb_rating']
del d['Episode Scor']
plt.hist(d['Episode Score'],bins=10,alpha=0.7,color = 'green', density = True)
plt.xlabel('Episode Score')
plt.ylabel('Frequency Density')
plt.title('Histogram of Episode Scores')
score_fit = norm.fit(d['Episode Score'])
x_vals = np.linspace(7.4,10.4,1000)
pdf = norm.pdf(x_vals,score_fit[0],score_fit[1])
plt.plot(x_vals,pdf)
nam = list(d["episode_name"])
cam_1 = ["Hot Girl","Office Olympics","The Fire","The Client","E-Mail Surveilance","Valentine's Day","The Convention","The Return","Launch Party (Parts 1&2)","Business Trip","New Boss","Two Weeks","Dream Team","The Michael Scott Paper Company","Stairmageddon","Promos","Livin' the Dream","Heavy Competition","Broke","Sabre","Mafia","The Delivery (Parts 1&2)","Nepotism","The Sting","The Seminar","Training Day","Michael's Last Dundies","Goodbye, Michael","The Inner Circle","Mrs. California","Turf War","Here Comes Treble","The Boat","Lice","Paper Airplane"]
cam_2 = ["Booze Cruise","Garden Party"]
cam_3 = ["Stress Relief (Parts 1&2)"]
cam_4 = ["Finale"]
cam_5 = ["Search Committee (Parts 1&2)"]
cameos = [0] * 186
for i in cam_1:
cameos[nam.index(i)] = 1
for i in cam_2:
cameos[nam.index(i)] = 2
for i in cam_3:
cameos[nam.index(i)] = 3
for i in cam_4:
cameos[nam.index(i)] = 4
for i in cam_5:
cameos[nam.index(i)] = 5
d["Number of Cameos"] = cameos
d["Episode Length (Minutes)"] = (23,23,22,23,23,23,21,22,22,22,30,22,22,22,23,22,22,22,22,22,22,22,22,21,22,22,22,22,
23,22,22,22,22,22,31,30,30,42,21,22,28,21,21,22,21,30,21,21,29,28,42,42,42,42,42,22,30,30,30,30,30,22,30,42,42,43,22,
22,30,22,22,30,30,22,30,30,22,60,60,60,22,30,30,30,30,21,22,30,22,30,30,30,30,30,30,49,30,30,30,23,30,30,30,30,30,30,
60,30,23,30,30,30,30,30,30,30,30,30,30,30,30,23,30,30,30,60,30,30,30,22,30,30,30,30,20,50,23,23,42,23,23,23,23,23,23,
23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,22,22,22,22,22,22,22,22,22,22,22,22,42,22,22,22,22,42,
43,51,)
d["Episode Location (Binary)"] = (0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,1,0,0,0,1,0,0,0,0,
1,1,1,0,0,0,1,1,0,1,0,0,0,0,1,1,1,1,0,1,0,1,0,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1,
0,0,1,0,0,1,0,0,1,0,0,0,1,0,0,1,0,1,0,1,1,0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,1,1,1,1,1,1,1,1,0,0,1,
0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1)
#finding best director
unique_directors = d['director'].unique().tolist()
directors = list(d["director"])
director_count = []
for i in unique_directors:
dir_count = 0
for j in d.director:
if j == i:
dir_count += 1
director_count.append(dir_count)
count = 0
names = []
for i in range(0,len(unique_directors)):
if director_count[i] >5:
names.append(unique_directors[i])
count+=1
for i in directors:
if i not in names:
directors[directors.index(i)] = "others"
directors1 = pd.get_dummies(directors)
directors1.columns = ['D: Charles McDougall', 'D: David Rogers','D: Greg Daniels','D: Jeffrey Blitz','D: Ken Kwapis',
'D: Ken Whittingham' , 'D: Matt Sohn','D: Paul Feig','D: Paul Lieberstein','D: Randall Einhorn','D: Others']
writers = list(d.writer)
for j in writers:
if j == 'Lee Eisenberg;Gene Stupnitsky':
writers[writers.index(j)] = "Gene Stupnitsky;Lee Eisenberg"
#finding best writers
unique_writers = d['writer'].unique().tolist()
writer_count = []
for i in unique_writers:
wri_count = 0
for j in d.writer:
if j == i:
wri_count += 1
writer_count.append(wri_count)
count = 0
names = []
for i in range(0,len(unique_writers)):
if writer_count[i] >6:
names.append(unique_writers[i])
count+=1
for i in writers:
if i not in names:
writers[writers.index(i)] = "others"
writers1 = pd.get_dummies(writers)
writers1.columns = ['W: B.J. Novak', 'W:Brent Forrester','W:Charlie Grandy','W:Gene Stupnitsky;Lee Eisenberg',
'W:Greg Daniels', 'W:Jennifer Celotta' , 'W:Justin Spitzer','W:Michael Schur','W:Mindy Kaling','W:Paul Lieberstein',
'W: Others']
d = pd.concat([d, writers1], axis=1)
d = pd.concat([d, directors1], axis=1)
# Lineshare (number of lines spoken by characters)
# Remove records that don't contain dialogue
df = df.dropna()
df
episode_names = d.episode_name
# Main 39 characters
main_chars = ["Michael","Dwight","Jim","Pam","Andy","Angela","Kevin","Erin","Oscar","Ryan"]
n_main_chars = len(main_chars) # Number of main characters
n_episodes = len(episode_names) # Number of episodes
array = np.zeros((n_episodes,n_main_chars)) # Create array of zeros with episode number as rows and character names as columns
# Create array of number of lines spoken by each characters
for i in episode_names:
for j in main_chars:
episode_names_list = episode_names.values.tolist()
row = episode_names_list.index(i)
col = main_chars.index(j)
episode = df.loc[df['episode_name'] == i]
character_lines = episode.loc[episode['character'] == j]
number_lines = len(character_lines)
array[row,col] = number_lines
# Create array with proportion of lines spoken by each character
array_proportion = np.zeros((186,10))
for i in range(0,186):
row_tot = sum(array[i,:])
for j in range(0,10):
array_proportion[i,j] = array[i,j]/row_tot
# To produce models using these variables without skewing the dependency on variables, we reduce the number of characters we consider
ten_main_prop = array_proportion
# To show that all columns are indepent, we show ten_main_prop has full rank
print("The rank of ten_main_prop is:",np.linalg.matrix_rank(ten_main_prop))
# Add characters to dataframe
for i in range(0,10):
d[main_chars[i]] = ten_main_prop[:,i]
pd.DataFrame(new_d)
pd.DataFrame(new_d)
new_d = d.drop(columns=['season','episode','director','writer','air_date','main_chars','episode_name'])
new_d_main_vars = new_d[['n_lines','n_directions','n_words','n_speak_char','Episode Score','Number of Cameos','Episode Length (Minutes)']]
plt.figure(figsize=(12, 8))
mask = np.triu(np.ones_like(new_d_main_vars.corr(), dtype=np.bool))
corr = new_d_main_vars.corr()
kot = corr[abs(corr)>=0.5]
ax = plt.axes()
sns.heatmap(kot, mask=mask, vmin=-1, vmax=1, annot=True,ax=ax)
ax.set_title("Correlation Heatmap between Main Variables", fontsize = 16)
plt.show()
new_d_characters = new_d[["Michael","Dwight","Jim","Pam","Andy","Angela","Kevin","Erin","Oscar","Ryan"]]
plt.figure(figsize=(12,8))
mask = np.triu(np.ones_like(new_d_characters.corr(),dtype=np.bool))
corr1 = new_d_characters.corr()
kot1 = corr1[abs(corr1)>=0.35]
ax = plt.axes()
sns.heatmap(kot1, mask=mask, vmin=-1, vmax=1, annot=True, ax = ax)
ax.set_title("Heatmap of correlation between main characters", fontsize = 16)
plt.show()
new_d = new_d.drop(columns=["n_lines","n_words","Dwight","Oscar","Erin","Andy","Kevin"])
# List of non binary and none dropped after heat map correlation test features.
non_binary_column_names = ['n_directions','n_speak_char','Number of Cameos','Episode Length (Minutes)',
'Michael','Jim','Pam','Angela','Ryan']
fig, axes = plt.subplots(3,3, figsize=(12, 12), sharey=True)
degree_range = 5
colours = ['red','green','black','magenta','cyan']
for i in range(3):
for j in range(3):
x = new_d[non_binary_column_names[3*i+j]]
y = new_d['Episode Score']
sns.scatterplot(ax=axes[i,j],x = new_d[non_binary_column_names[3*i+j]],y = new_d['Episode Score'])
for r in range(1,degree_range):
degree = r
X = pd.DataFrame(x)
X_seq = np.linspace(X.min(),X.max(),186).reshape(-1,1)
polyreg=make_pipeline(PolynomialFeatures(degree),LinearRegression())
polyreg.fit(X,y)
axes[i,j].plot(X_seq,polyreg.predict(X_seq),color=colours[r-1], label = 'Degree =' + str(r))
if j == 0:
if i == 0:
axes[i,j].legend()
new_d = new_d.drop(columns=["Jim","Ryan","Pam","Angela"])
#helper functions
def get_coefs(m):
"""Returns the model coefficients from a Scikit-learn model object as an array,
includes the intercept if available.
"""
# If pipeline, use the last step as the model
if (isinstance(m, sklearn.pipeline.Pipeline)):
m = m.steps[-1][1]
if m.intercept_ is None:
return m.coef_
return np.concatenate([[m.intercept_], m.coef_])
def model_fit(m, X, y, title, plot = False):
"""Returns the root mean squared error of a fitted model based on provided X and y values.
Args:
m: sklearn model object
X: model matrix to use for prediction
y: outcome vector to use to calculating rmse and residuals
plot: boolean value, should fit plots be shown
"""
y_hat = m.predict(X)
rmse = mean_squared_error(y, y_hat, squared=False)
res = pd.DataFrame(
data = {'y': y, 'y_hat': y_hat, 'resid': y - y_hat}
)
if plot:
plt.figure(figsize=(12, 6))
plt.subplot(121)
sns.lineplot(x='y', y='y_hat', color="grey", data = pd.DataFrame(data={'y': [min(y),max(y)], 'y_hat': [min(y),max(y)]}))
sns.scatterplot(x='y', y='y_hat', data=res).set_title("Fit plot")
plt.subplot(122)
sns.scatterplot(x='y', y='resid', data=res).set_title("Residual plot")
plt.subplots_adjust(left=0.0)
plt.suptitle(title + str(round(rmse, 4)), fontsize=16)
plt.show()
return rmse
kf = KFold(n_splits = 5, shuffle = True,random_state=1)
result = next(kf.split(new_d), None)
train = new_d.iloc[result[0]]
test = new_d.iloc[result[1]]
x_train = train.drop(columns=['Episode Score'])
y_train = train['Episode Score']
# Non scaled
lm = LinearRegression(fit_intercept = False).fit(x_train, y_train)
#From now on we will only fit to x_train with scaled categorical features
scaled_features = x_train
col_names = ['n_directions', 'n_speak_char', 'Michael', 'Number of Cameos','Episode Length (Minutes)']
features = scaled_features[col_names]
scaler = StandardScaler().fit(features.values)
features = scaler.transform(features.values)
scaled_features[col_names] = features
x_train_scaled = scaled_features
lm_scaled = make_pipeline(
LinearRegression(fit_intercept=False)
).fit(x_train_scaled, y_train)
# Tuning linear ridge model for best alpha
alphas = np.linspace(0.0000001, 0.01, num=200)
lm_scaled_ridge = GridSearchCV(
make_pipeline(
Ridge(fit_intercept = False)
),
param_grid={'ridge__alpha': alphas},
cv=KFold(5, shuffle=True),
scoring="neg_root_mean_squared_error"
).fit(x_train_scaled, y_train)
get_coefs(lm_scaled_ridge.best_estimator_)[1:30]
alphas = np.linspace(0.0000001, 1, num=100)
lm_scaled_lasso = GridSearchCV(
make_pipeline(
Lasso(fit_intercept = False)
),
param_grid={'lasso__alpha': alphas},
cv=KFold(5, shuffle=True),
scoring="neg_root_mean_squared_error"
).fit(x_train_scaled, y_train)
get_coefs(lm_scaled_lasso.best_estimator_)[1:30]
poly_model = make_pipeline(
make_column_transformer(
(PolynomialFeatures(include_bias=False), ['n_directions']),
(PolynomialFeatures(include_bias=False), ['n_speak_char']),
(PolynomialFeatures(include_bias=False), ['Number of Cameos']),
remainder='passthrough'),
LinearRegression(fit_intercept=False))
parameters_poly_model = {
'columntransformer__polynomialfeatures-1__degree': np.arange(1,4,1),
'columntransformer__polynomialfeatures-2__degree': np.arange(1,3,1),
'columntransformer__polynomialfeatures-3__degree': np.arange(1,4,1),
}
kf = KFold(n_splits=5, shuffle=True, random_state=0)
grid_search_poly_model = GridSearchCV(poly_model,
parameters_poly_model,
cv=kf,
scoring="neg_root_mean_squared_error").fit(x_train_scaled, y_train)
get_coefs(grid_search_poly_model.best_estimator_)[1:30]
alphas_poly = np.linspace(0.014, 0.015, num=10)
poly_model_lasso = make_pipeline(
make_column_transformer(
(PolynomialFeatures(include_bias=False), ['n_directions']),
(PolynomialFeatures(include_bias=False), ['n_speak_char']),
(PolynomialFeatures(include_bias=False), ['Number of Cameos']),
remainder = 'passthrough'),
Lasso(fit_intercept = False))
parameters_poly_model_lasso = {
'columntransformer__polynomialfeatures-1__degree': np.arange(1,4,1),
'columntransformer__polynomialfeatures-2__degree': np.arange(1,3,1),
'columntransformer__polynomialfeatures-3__degree': np.arange(1,4,1),
'lasso__alpha': alphas_poly,
}
grid_search_poly_model_lasso = GridSearchCV(poly_model_lasso,
parameters_poly_model_lasso,
cv=kf,
scoring="neg_root_mean_squared_error",
).fit(x_train_scaled, y_train)
get_coefs(grid_search_poly_model_lasso.best_estimator_)[1:30]
alphas = np.linspace(2, 5, num=10)
poly_model_ridge = make_pipeline(
make_column_transformer(
(PolynomialFeatures(include_bias=False), ['n_directions']),
(PolynomialFeatures(include_bias=False), ['n_speak_char']),
(PolynomialFeatures(include_bias=False), ['Number of Cameos']),
remainder = 'passthrough'),
Ridge(fit_intercept = False))
parameters_poly_model_ridge = {
'columntransformer__polynomialfeatures-1__degree': np.arange(1,4,1),
'columntransformer__polynomialfeatures-2__degree': np.arange(1,3,1),
'columntransformer__polynomialfeatures-3__degree': np.arange(1,4,1),
'ridge__alpha': alphas,
}
grid_search_poly_model_ridge = GridSearchCV(poly_model_ridge,
parameters_poly_model_ridge,
cv=kf,
scoring="neg_root_mean_squared_error",
).fit(x_train_scaled, y_train)
get_coefs(grid_search_poly_model_ridge.best_estimator_)[1:30]
x_test = test.drop(columns=['Episode Score'])
y_test =test['Episode Score']
rmse_nslm = model_fit(lm, x_test, y_test, 'Non Scaled Linear Model RMSE: ', plot=False)
rmse_slm = model_fit(lm_scaled, x_test, y_test, 'Scaled Linear Model RMSE: ', plot=False)
rmse_rr_slm = model_fit(lm_scaled_ridge.best_estimator_, x_test, y_test, 'Ridge Regression on Scaled Linear Model RMSE: ', plot=False)
rmse_l_slm = model_fit(lm_scaled_lasso.best_estimator_, x_test, y_test, 'Lasso on Scaled Linear Model RMSE: ', plot=False)
rmse_std_pm = model_fit(grid_search_poly_model.best_estimator_, x_test, y_test, 'Standard Polynomial Model RMSE: ',plot=False)
rmse_l_pm = model_fit(grid_search_poly_model_lasso.best_estimator_, x_test, y_test, 'Lasso Polynomial Model RMSE: ',plot=False)
rmse_rr_pm = model_fit(grid_search_poly_model_ridge.best_estimator_, x_test, y_test, 'Ridge Polynomial Model RMSE: ',plot=False)
rmse_results ={'Method:':['Non-Scaled Linear','Scaled Linear','Ridge Scaled Linear',
'Lasso Scaled Linear','Polynomial','Ridge Polynomial','Lasso Polynomial'],
'RMSE':np.around([rmse_nslm, rmse_slm, rmse_rr_slm, rmse_l_slm, rmse_std_pm, rmse_rr_pm, rmse_l_pm], decimals = 3)
}
rmse_results = pd.DataFrame(rmse_results)
rmse_results.sort_values(by = ['RMSE'])
new_d_variables = new_d.drop(columns=['Episode Score'])
print(poly_model_lasso.best_estimator_.coefs_)
variables = new_d_variables.columns.values
variables_list = list(variables)
variables_list.insert(0, 'intercept')
table1 = {'variable:':variables_list,
'coefficients:': coeffs
}
chosen_model = pd.DataFrame(table1).transpose()
chosen_model
coeffs
array_prop_df = pd.DataFrame(array_proportion)
# Find episode in which Michael has the proportion of lines
index = np.argmax(array_prop_df[[0]])
best_line = {'character:': main_chars,
'proportion' : np.around(array_prop_df.iloc[index,:],2)
}
best_line_df = pd.DataFrame(best_line).transpose()
best_line_df
# Writer
writer_coeffs = lm_scaled_ridge.best_estimator_.named_steps['ridge'].coef_[5:16]
best_writer = variables[np.argmax(writer_coeffs)+5]
best_writer_coeffs = max(writer_coeffs)
# Director
director_coeffs = lm_scaled_ridge.best_estimator_.named_steps['ridge'].coef_[16:27]
best_director = variables[np.argmax(director_coeffs)+16]
best_director_coeffs = max(director_coeffs)
print('The best writer is', best_writer, 'with coefficient', best_writer_coeffs)
print('The best director is', best_director, 'with coefficient', best_director_coeffs)
fig = plt.figure()
ax = plt.axes()
x = np.linspace(0, 60, 1000)
ax.plot(x,chosen_model.iloc[1,0]+chosen_model.iloc[1,1]*x ,
color = 'red',label='n_directions');
ax.plot(x,chosen_model.iloc[1,0]+chosen_model.iloc[1,2]*x ,
color = 'blue',label='n_speak_char');
ax.plot(x,chosen_model.iloc[1,0]+chosen_model.iloc[1,4]*x ,
color = 'green',label='Episode Length');
ax.plot(x,chosen_model.iloc[1,0]+chosen_model.iloc[1,3]*x ,
color = 'cyan',label='Number of Cameos');
plt.ylim(8.6, 10)
plt.xlim(0, max(x))
plt.legend()
plt.title("Continuous Variable in their own dimensional plane")
plt.xlabel("Value")
plt.ylabel("Model Predicted Episode Score")
plt.show()
recc = chosen_model.loc[:, 1:4]
new_int = chosen_model.iloc[1,0]+chosen_model.iloc[1,10]+chosen_model.iloc[1,19]+chosen_model.iloc[1,28]*0.7+chosen_model.iloc[1,29]*0.06+chosen_model.iloc[1,30]*0.04
recc.insert(0, '0', ["Intercept" ,new_int])
recc
# Run the following to render to PDF
!crew install pandoc
!jupyter nbconvert --to pdf project1.ipynb
