IPython Notebooks for StatLearning Exercises

Earlier this year, I attended the StatLearning: Statistical Learning course, a free online course taught by Stanford University professors Trevor Hastie and Rob Tibshirani. They are also the authors of The Elements of Statistical Learning (ESL) and co-authors of its less math-heavy sibling: An Introduction to Statistical Learning (ISL). The course was based on the ISL book. Each week's videos were accompanied by some hands-on exercises in R.

I personally find it easier to work with Python than R. R seems to have grown organically with very little central oversight, so function and package names are often non-intuitive, and often have duplicate or overlapping functionality. In general, an educated guess about an R function has about the same likelihood of being right as a completely random one - unless you know the function or package, your chances are 50-50. On the other hand, with Python, an educated guess has a 40-90 percent chance of being right, depending on the library and how educated your guess was. So while the good profs were patiently explaining the R code, I was mostly busy fantasizing about writing all of it in Python some day.

At the time, I had worked a bit with scikit-learn and NumPy. I had heard about Pandas and knew it was the Python implementation of DataFrames, but hadn't actually worked with it. Over the past couple of months, I have had the opportunity to work with Pandas and IPython Notebooks for a project I did with my kids, and as a result I now quite enjoy the power and expressivity that these libraries provide.

So I decided to apply my newly acquired skills to do this rewrite. One of my incentives for doing this was the chance to get a fairly comprehensive guided tour of scikit-learn algorithms that I wouldn't normally use. Of course, the tour depends a lot on the guide, and the course is taught from the point of view of a statistician than a machine learning person. Since my toolchain (scikit-learn, NumPy, SciPy, Pandas, MatplotLib and a bit of statsmodels) is more focused towards Machine Learning, there were times when I wasn't able to replicate the functionality completely and accurately.

There are 9 notebooks listed below, corresponding to the exercises for Chapters 2-10 of the course. The notebooks and data can be found on my GitHub in the project statlearning-notebooks. You can also read the notebooks directly on the nbviewer.ipython.org via the links in the README.md file.

• Chapter 2: Basic Operations
• Chapter 3: Linear Regression
• Chapter 4: Classification
• Chapter 5: Cross-Validation and Bootstrap
• Chapter 6: Feature Selection
• Chapter 7: Nonlinear Models
• Chapter 8: Decision Trees
• Chapter 9: Support Vector Machines
• Chapter 10: Unsupervised Methods

This exercise introduced me to a lot of scikit-learn algorithms that I had not used before. Since there are quite a few functionality mismatches between R and scikit-learn, trying to match it often led me to novel ideas described on sites like StackOverflow and Cross-Validated, some of which I implemented (and others I have linked to). I also learned quite a bit about plotting with matplotlib, since the original exercises use R's rich plotting features as a matter of course, some of which require additional work in Python.

Overall, I found that the group of Python libraries were more than adequate for most tasks in the exercises, and (at least in my eyes) resulted in cleaner, more readable code. Take a look at these pages to get an overview of what scikit-learn and Pandas, my two top level libraries, can do. However, R also offers lots of functionality - there is lot of overlap, but in some cases R provides algorithms that scikit-learn doesn't. However, scikit-learn has many more algorithms compared to R. So it makes sense to learn and use both as needed.

If you are considering using my group of Python libraries for data analysis, then the notebooks should be useful as examples. For more advanced programmers, if you think there are better ways to do something than what I have done, I would appreciate hearing from you (or since its on GitHub, a pull request would be good too!).

Sentiment Analysis using Classification

At the Introduction to Data Science course I took last year at Coursera, one of our Programming Assignments was to do sentiment analysis by aggregating the positivity and negativity of words in the text against the AFINN word list, a list of words manually annotated with positive and negative valences representing the sentiment indicated by the word.

At the time I wondered if perhaps the word list approach was not too labor intensive, since one must go through a manual process for each domain to identify and score positive and negative words. I figured it may be better to just treat it as a classification problem - manually identifying documents (instead of words) as positive or negative, then use that to train a classifier that can predict the sentiment of unseen documents. But then I got busy with other things and forgot about this until a few days ago, when I came across this post where it describes using classification for sentiment analysis.

The author, Andy Bromberg, describes using NLTK and Python to classify movie reviews as positive or negative. He also refers to a previous attempt using R and the AFINN polarity wordlist, similar to the Programming Assignment I described earlier. In addition, the post describes how feature selection was used to increase the accuracy of the classifier.

As a learning exercise, I decided to do something similar with Scikit-Learn. I used the review training data from the Yelp Recruiting Competition on Kaggle, which I had entered as part of the Peer Assessments in the Intro to Data Science course. Part of the data consisted of 229,907 restaurant reviews in JSON format which had votes by users to indicate usefulness, funnyness and coolness of the review. I used the text as a bag of words and consider a review to be useful, funny and cool respectively if they have more than 0 votes for that attribute. This is used to train 3 binary classifiers that can predict these attributes in new reviews. Following along with Andy's post, I then used the Chi-squared metric to find the most useful features and measured accuracy, precision and recall of the classifiers for different feature sizes.

The code to build and test each classifier using 10-fold cross validation is shown below. We first read the review files, parsing each line into a JSON object, then extracting the text and the useful, funny and cool votes. We then convert the text into a sparse matrix where each word is a feature. In our first pass, we use every word (116,713 unique words) in our text. For each of the useful, funny or cool attribute, we use the matrix and the binarized vote vector for that attribute to construct a Naive Bayes classifier. We then test the classifier and compute accuracy, precision and recall. Next we calculate the most informative features for the attribute using the Chi-squared test, and build models for 1000, 3000, 10000, 30000, and 100000 top features and calculate their accuracy, precision and recall.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76

# Source: src/yelp_ufc/build_classifier.py from sklearn.cross_validation import KFold from sklearn.feature_extraction.text import CountVectorizer from sklearn.feature_selection import chi2 from sklearn.metrics import accuracy_score, precision_score, recall_score from sklearn.naive_bayes import MultinomialNB import json import numpy as np import operator def read_data(fname): f = open(fname, 'rb') texts = [] ys = [] for line in f: rec = json.loads(line.strip()) texts.append(rec["text"]) ys.append([ 1 if int(rec["votes"]["useful"]) > 0 else 0, 1 if int(rec["votes"]["funny"]) > 0 else 0, 1 if int(rec["votes"]["cool"]) > 0 else 0]) f.close() return texts, np.matrix(ys) def vectorize(texts, vocab=[]): vectorizer = CountVectorizer(min_df=0, stop_words="english") if len(vocab) > 0: vectorizer = CountVectorizer(min_df=0, stop_words="english", vocabulary=vocab) X = vectorizer.fit_transform(texts) return vectorizer.vocabulary_, X def cross_validate(ufc_val, X, y, nfeats): nrows = X.shape[0] kfold = KFold(nrows, 10) scores = [] for train, test in kfold: Xtrain, Xtest, ytrain, ytest = X[train], X[test], y[train], y[test] clf = MultinomialNB() clf.fit(Xtrain, ytrain) ypred = clf.predict(Xtest) accuracy = accuracy_score(ytest, ypred) precision = precision_score(ytest, ypred) recall = recall_score(ytest, ypred) scores.append((accuracy, precision, recall)) print ",".join([ufc_val, str(nfeats), str(np.mean([x[0] for x in scores])), str(np.mean([x[1] for x in scores])), str(np.mean([x[2] for x in scores]))]) def sorted_features(ufc_val, V, X, y, topN): iv = {v:k for k, v in V.items()} chi2_scores = chi2(X, y)[0] top_features = [(x[1], iv[x[0]], x[0]) for x in sorted(enumerate(chi2_scores), key=operator.itemgetter(1), reverse=True)] print "TOP 10 FEATURES FOR:", ufc_val for top_feature in top_features[0:10]: print "%7.3f %s (%d)" % (top_feature[0], top_feature[1], top_feature[2]) return [x[1] for x in top_features] def main(): ufc = {0:"useful", 1:"funny", 2:"cool"} texts, ys = read_data("../../data/yelp_ufc/yelp_training_set_review.json") print ",".join(["attrtype", "nfeats", "accuracy", "precision", "recall"]) for ufc_idx, ufc_val in ufc.items(): y = ys[:, ufc_idx].A1 V, X = vectorize(texts) cross_validate(ufc_val, X, y, -1) sorted_feats = sorted_features(ufc_val, V, X, y, 10) for nfeats in [1000, 3000, 10000, 30000, 100000]: V, X = vectorize(texts, sorted_feats[0:nfeats]) cross_validate(ufc_val, X, y, nfeats) if __name__ == "__main__": main()

The top 10 features for each classifier (ie the words that have highest "polarity" for that particular attribute) are shown below. The first column is the Chi-squared score for the word, the second column is the word itself, and the third column is the index of the word in the sparse matrix.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

TOP 10 FEATURES FOR: useful 5170.064 like (60636) 4835.884 just (56649) 2595.147 don (32684) 2456.476 know (58199) 2346.778 really (84130) 2083.032 time (104423) 2063.618 people (76776) 2039.718 place (78659) 1873.081 think (103835) 1858.230 little (61092) TOP 10 FEATURES FOR: funny 9087.141 like (60636) 6049.875 just (56649) 4848.157 know (58199) 4664.542 don (32684) 3361.983 people (76776) 2649.594 think (103835) 2505.478 oh (72420) 2415.325 ll (61174) 2349.312 really (84130) 2345.851 bar (11472) TOP 10 FEATURES FOR: cool 6675.123 like (60636) 4616.683 just (56649) 3173.775 know (58199) 3010.526 really (84130) 2847.494 bar (11472) 2715.794 little (61092) 2670.838 don (32684) 2300.151 people (76776) 2217.659 place (78659) 2216.888 ve (110157)

We also plot some graphs for each classifier showing how the accuracy, precision and recall vary with the number of features. The horizontal lines represent the accuracy, precision and recall achieved using the full data set. As can be seen, the metrics improve as more features are added but tend to flatten out eventually.

The code to build these graphs out of the metrics printed out by our classifier training code uses Pandas dataframe plotting functionality and is shown below:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

# Source: src/yelp_ufc/plot_results.py import matplotlib.pyplot as plt import pandas as pd import sys def main(): assert(len(sys.argv) == 2) df = pd.read_csv("all.csv") adf = df.ix[df.attrtype == sys.argv[1]] adf_all = adf.ix[adf.nfeats < 0] adf_rest = adf.ix[adf.nfeats > 0] print adf_all print adf_rest adf_rest = adf_rest.drop("attrtype", 1) adf_rest = adf_rest.set_index("nfeats") adf_rest["accuracy_all"] = adf_all[["accuracy"]].values[0][0] adf_rest["precision_all"] = adf_all[["precision"]].values[0][0] adf_rest["recall_all"] = adf_all[["recall"]].values[0][0] adf_rest.plot(title=sys.argv[1]) plt.show() if __name__ == "__main__": main()

Thats all I have for today. Many thanks to Andy Bromberg for posting his analysis, without which my analysis would not have happened. The code for this blog can also be found on my GitHub.