Teaching a Machine to Read Pitchfork

#import the artists and albums from 'meddulla's github page
pitchfork_data = pd.read_csv(
    'https://gist.githubusercontent.com/meddulla/8277139/raw/f1595d50cada910d082bc1dfd8ef47ff49910cb3/pitchfork-reviews.csv',
           sep = ',')
pitchfork_data.head()

# a note on scraping pitchfork: according to their robots.txt file, scraping is allowed,
# as long as users aren't scraping any of their search pages.

#input: urls in a pandas series format
#output: review text for each of those links as an array of strings
from bs4 import BeautifulSoup

def get_multiple_reviews(urls):
    if isinstance(urls, pd.Series):
        review_text = []

        for link in urls:
            #pull url html into a beautifulsoup object        
            try:
                r = requests.get(link)
                soup = BeautifulSoup(r.text, 'lxml')
                #pull all paragraph tags/contents with soup.findAll, and then put the text into one string
                review_text.append('\n\n'.join(map(str, [i.text.lower() for i in soup.findAll('p')])))
            except Exception:
                review_text.append(np.nan)
        return review_text
    
    else:
        print('Check inputs - confirm they\'re in pd.Series format')
        return

##this code executes the review scraping function. since scraping 5000 reviews takes a while, I've saved it as a csv file. 
#sample_df = pitchfork_data.sample(5000)
#sample_df['review'] = sample_reviews_merged
#sample_df.to_csv('C:\\Users\\Matt\\Documents\\Python_Scripts\\CS109_Zips\\2013-homework\\HW Practice\\pitchfork_reviews.csv')

#pull in scraped data saved as a csv for convenience
sample_df = pd.read_csv(
    'C:\\Users\\Matt\\Documents\\Python_Scripts\\CS109_Zips\\2013-homework\\HW Practice\\pitchfork_reviews.csv', 
     sep=',', encoding = 'latin-1', index_col = 'Unnamed: 0')
sample_df.head()

print('min/mean/max: ',pitchfork_data.score.min(), '/', pitchfork_data.score.mean(), '/', pitchfork_data.score.max())
print('# reviews: ',len(sample_df))

print('NA values before cleaning: ',sum(sample_df.review.isnull()))
sample_df = sample_df.reset_index(drop=True)
sample_df.drop(sample_df.index[sample_df[sample_df.review.isnull()].index.values], axis=0, inplace=True)
print('NA values now: ',sum(sample_df.review.isnull()))

print('# of reviews: ', len(sample_df))
print('# of reviewers: ', len(sample_df.reviewer.unique()))

#Histogram to look at the spread of reviewers
plt.title('How Prolific are the Reviewers?')
plt.xlabel('Quantity of Reviews')
plt.ylabel('Number of Reviewers')
plt.hist(sample_df.groupby('reviewer').review.count(), bins=40)
plt.show()

top_5_reviewers = sample_df.groupby('reviewer').review.count().nlargest(5)
print('The 5 most prolific reviewers:\n', top_5_reviewers)

#pull the names of all reviewers with >n reviews
#use that to get average reviews for relatively prolific reviewers
floor_score = 10
df_count = sample_df.groupby('reviewer').count()
df_mean = sample_df.groupby('reviewer').mean()
plt.title('Average score for reviewers with >%i reviews' %floor_score)
plt.ylabel('# Reviewers')
plt.xlabel('Score')
plt.hist(df_mean[sample_df.groupby('reviewer').count().review>floor_score].score)
plt.show()

#import vectorizer and fit it to our data
from sklearn.feature_extraction.text import CountVectorizer
vectorizer = CountVectorizer()
vectorizer.fit(sample_df.review)

#assign vectorizer data to an array
x = vectorizer.transform(sample_df.review)
x = x.toarray()
print('number of words in the dataframe: ', len(vectorizer.get_feature_names()))
words = pd.DataFrame({'word':vectorizer.get_feature_names(),'freq':x.sum(axis=0)})
print('20 highest frequency words: ', [i for i in words.nlargest(20, 'freq').word.values])

print('# words that appear once: ', words[words.freq == 1].freq.count())
print('The most common word appears ', words.freq.max(),' times')
print('Most common word: ', words[words.freq == words.freq.max()].values[0][1])

vectorizer = CountVectorizer(min_df = 1e-03, max_df = 0.1)
vectorizer.fit(sample_df.review)
x = vectorizer.transform(sample_df.review)
x = x.toarray()
len(vectorizer.get_feature_names())
words = pd.DataFrame({'word':vectorizer.get_feature_names(),'freq':x.sum(axis=0)})
print('20 highest frequency words: ', [i for i in words.nlargest(20, 'freq').word.values])

print('# words that appear once: ', words[words.freq == 1].freq.count())
print('The most common word appears ', words.freq.max(),' times')
print('Most common word: ', words[words.freq == words.freq.max()].values[0][1])

plt.hist(words.freq, bins=400)
plt.title('Stop words eliminated')
plt.xlabel('Frequency')
plt.ylabel('Number of words that appear F times')
plt.axis([0,500,0,10000])
plt.show()

print('number of words appearing more than once: ', words[words.freq>1].freq.count())
print('number of words appearing once: ', words[words.freq==1].freq.count())

cutoff_score = 8
print('Percent of reviews that are \'best new music\' (>%i.0): ' %cutoff_score, sum(sample_df.score>cutoff_score)
      /sum(sample_df.score>0)*100)

#We'll use this function in our grid search to reconstruct the vectorizer
#critics is the dataframe of reviews, etc.
from sklearn.feature_extraction.text import CountVectorizer

def make_xy(critics, mind=1, maxd=0.99):
    #create y matrix, where y=1 if it's greater than our cutoff score
    bnm = critics.score>cutoff_score
    y = bnm.astype(int).values
    
    #set up the vectorizer, with terms for min_df and max_df
    vectorizer = CountVectorizer(min_df = mind, max_df = maxd)
    vectorizer.fit(critics.review)
    x = vectorizer.transform(critics.review)
    x = x.toarray()
    
    return x, y, vectorizer

X, Y, v = make_xy(sample_df, 1, 0.99)

print('Number of words analyzed: ', len(v.get_feature_names()))
print('Average words per article per the min/max_df above: ', np.array([i.sum() for i in X]).mean())

from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import MultinomialNB
#Select best parameters (min_df, max_df, alpha)

#cs109 parameter selection. gridsearch above
import sklearn.model_selection
#we'll use the area under the roc curve as a metric
from sklearn.metrics import roc_auc_score

#the grid of parameters to search over
alphas = [.1, 0.5, 1, 3, 5, 10, 50]
min_dfs = [1e-4, 1e-3, 1e-2]
max_dfs = [1e-1, 0.2, 0.5]

#Find the best value for alpha and min_df, and the best classifier
best_alpha = None
best_min_df = None
best_max_df = None
max_auc = -np.inf

for a in alphas:
    for min_df in min_dfs:   
        for max_df in max_dfs:
            #configure vectorizer
            #vectorizer = CountVectorizer(min_df = min_df, max_df = max_df)       
            X, Y, v = make_xy(sample_df, min_df, max_df)

            #split up testing/training
            xtrain, xtest, ytrain, ytest = train_test_split(X, Y, test_size = 0.2)
            
            #fit the naive bayes model
            nb = MultinomialNB()
            nb.set_params(alpha=a)
            nb.fit(xtrain, ytrain)
            
            #function below implements kfold validation to find best auc
            y_score = nb.predict_proba(xtest)[:,1]
            nll = roc_auc_score(ytest, y_score)
            
            #set best params
            if nll > max_auc:
                max_auc = nll
                best_alpha = a
                best_min_df = min_df
                best_max_df = max_df

#show best parameters
print(max_auc, '|', best_alpha, '|', best_min_df, '|', best_max_df)

from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import MultinomialNB

#split our data into training data (80%) and testing (20%)
X, Y, v = make_xy(sample_df, best_min_df, best_max_df)
#provide indices so we can examine later
tts_indices = list(range(len(X)))
xtrain, xtest, ytrain, ytest, indices_train, indices_test = train_test_split(
    X, Y, tts_indices, test_size = 0.2, random_state = 10)

#fit the naive bayes model
nb = MultinomialNB(alpha=best_alpha)
nb.fit(xtrain, ytrain)

#test training prediction
print('training accuracy: ', sum(nb.predict(xtrain)==ytrain)/len(ytrain))
print('testing accuracy: ', sum(nb.predict(xtest)==ytest)/len(ytest))
#well this model is most definitely overfit - training performance is much better than testing
#also with just 15% of music as bnm, I'm concerned with 81% accuracy

#and that scorer is area under the receiver operating characteristic curve! (aucroc)
from sklearn.metrics import roc_auc_score
y_score = nb.predict_proba(xtest)[:,1]
roc_auc_score(ytest, y_score)

#ROC Curve
from sklearn.metrics import roc_curve, auc
fpr, tpr, thresh = roc_curve(ytest, y_score)
plt.plot(fpr, tpr)
plt.title('ROC Curve')
plt.xlabel('False positive rate (100-specificity)')
plt.ylabel('True positive rate (sensitivity)')
plt.plot(np.linspace(0,1,100),np.linspace(0,1,100))

#calculate youden's j and identify the threshold at which it's maximized
youden_j = tpr + (1-fpr) -1
ideal_thresh = thresh[np.where(youden_j == youden_j.max())]
plt.plot(np.linspace(1, len(thresh),len(thresh)),thresh)
plt.plot(np.linspace(1, len(thresh), len(thresh)),youden_j)
plt.plot(np.where(youden_j == youden_j.max()),ideal_thresh,'o', c='r')
print('ideal threshold: ',ideal_thresh)
plt.show()

#let's assess whether our model is well-calibrated (i.e. whether its confidence % is about equal to its accuracy %)
hst = plt.hist(nb.predict_proba(xtest)[:,1], bins=20)

bin_counts = hst[0]
probabilities = nb.predict_proba(xtest)
digits = np.digitize(probabilities[:,1], hst[1])

#next step is to figure out what % of each bin in our test data is fresh
#so we'll need to assign ytest to bins per digits
binframe = pd.DataFrame({'y':ytest, 'bin':digits})
bnm_prob = binframe.groupby('bin').sum()/binframe.groupby('bin').count()

#number of reviews in each bin divided by the total number of (test) reviews
pct_df = binframe.groupby('bin').count()/binframe.bin.count()

#pct_df won't have the right # of items if the bins aren't nicely distributed, so we'll pad it
ind = list(range(1,len(hst[0])+2))
indy_df = pd.DataFrame({'b':np.zeros(len(ind))}, index=ind)
indy_df = indy_df.join(pct_df)
indy_df = indy_df.max(axis=1)

#plot the cumulative distribution function
plt.xlabel('P(fresh)')
plt.ylabel('P(predicted fresh)')
plt.title('The model is overconfident at both extremes, typical of Naive Bayes models')

#plt.plot(hst[1][:-1], np.cumsum(pct_df.values),'.')
plt.plot(hst[1], np.cumsum(indy_df.values),'.')
plt.plot(hst[1][:-1],hst[1][:-1], alpha=.5)
plt.legend(['Distribution of predicted values', 'True distribution'])

#create list with the indices and test results (w), and split into correct and incorrect
all_pred = list(zip(indices_test,nb.predict(xtest)==ytest))
correct_pred = [i[0] for i in all_pred if i[1] == True]
incorrect_pred = [i[0] for i in all_pred if i[1] == False]

#Where false, the prediction didn't agree with reality
print(sample_df.iloc[incorrect_pred[0],:],'\n\n', 
      sample_df.iloc[incorrect_pred[1],:],'\n\n', sample_df.iloc[incorrect_pred[2],:])

sample_df.iloc[incorrect_pred,2].mean()
#the reviews it missed had an average score of 7.64, which is fairly close to our dividing line.

print(sample_df.iloc[incorrect_pred[0],5],'\n\n')
      #sample_df.iloc[incorrect_pred[1],:],'\n\n', sample_df.iloc[incorrect_pred[2],:])

#testing to see which words indicated good results vs. poor
words = np.array(v.get_feature_names())
x = np.eye(xtest.shape[1])
probs = nb.predict_log_proba(x)[:, 0]
ind = np.argsort(probs)

print('good words: ',words[ind[:10]])
print('bad words: ',words[ind[-10:]])

#testing out review snippets

#reviews comprised of words that seem like they would be good/bad, but did not appear in best/worst 100 words
print('good sounding words: ', nb.predict_proba(v.transform(
    ['meticulous and timeless masterpiece, overwhelmingly excellent'])))
print('bad sounding words: ', nb.predict_proba(v.transform(
    ['disgusting overused meaningless drivel'])))

#reviews based on the best/worst 100 words of our vectorizer
print('best words: ', nb.predict_proba(v.transform(
    ['this otherwordly quintet is absolute gas, a terrific achievement for 20th century france'])))
print('worst words: ', nb.predict_proba(v.transform(
    ['an aimless, jerky and redundant shoegaze group from the 2000s trades principle for cliches'])))