<jupyter_start><jupyter_text>Lede Algorithms 2018 Week 2 class 2 - Logistic RegressionThis notebook builds up to logistic regression, that is, fitting a model that predicts the output where there are only two choices. But before that, we're going to do many comparisons and visualizations to talk about the underlying ideas that regression is getting at -- especially the idea of "controlling for" another factor. Throughout we'll be using data on citations for speeding from Massechusets from 2001. This data was originally analyzed by the Boston Globe for their series [Speed Trap: Who gets a ticket, who gets a break?](http://archive.boston.com/globe/metro/packages/tickets/). You can read about their analysis in this [methodology document](http://archive.boston.com/globe/metro/packages/tickets/study.pdf), and get the original data [here](http://archive.boston.com/globe/metro/packages/tickets/warnings.zip). I've added a header row the CSV file, but otherwise `tickets-warnings.csv` in this directory is identical.Why such old data? Well, for one thing it's a classic example of multi-variable regression in journalism. But also, for a brief period of April and May 2001, the state government entered both tickets and warnings from the paper citations. In Massechusets at that time, both were written on the same form, with only a checkbox indicating which it was. A ticket means a fine and a raises your insurance premium. A warning means nothing. Having data on both means we can ask who gets a ticket vs. a warning for the same circumstances.We have a ton of data about each stop, including the speed, age, race, gender, and home zip code of the driver. The journalistic question is, are some groups more likely to get tickets, and not warnings, than others?<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
tw_raw = pd.read_csv('tickets-warnings.csv')
tw_raw<jupyter_output><empty_output><jupyter_text>For a complete description of all 57 fields, see the [methodology](http://archive.boston.com/globe/metro/packages/tickets/study.pdf). For our purposes, we are only intersted in the citations for speeding, and most interested in the fields `TYPE` (ticket or warning), `MINORITY`, `SEX`, `AGE`, and `MPHOVER`.<jupyter_code>tw = tw_raw[tw_raw['DESCRIPT']=='SPEEDING'].reset_index(drop=True)
tw = tw[['DATE', 'TYPE', 'AGENCY3', 'SEX','MINORITY','AGE','MPH','MPHOVER','INTOWN']]
tw.head()<jupyter_output><empty_output><jupyter_text>Let's start with the total number of minority vs. white drivers<jupyter_code>tw.groupby('MINORITY').count()<jupyter_output><empty_output><jupyter_text>What fraction of citations were given to minorities?<jupyter_code>12836/(12836+79275)<jupyter_output><empty_output><jupyter_text>Is 13% a lot or a little? We could compare to the minority population of the towns where these stops happened, but what about people driving through? And there's no guarantee drivers in both groups have the same distribution of speeds. This is often called the "denominator problem." To deal with this issue, we're going to compare the probability of getting a warning, as opposed to a ticket, between the groups.First, let's clean up the data a little but dropping the (relatively small) number of records where the minority status is unknown.<jupyter_code># drop the unknown race entries
tw = tw[tw['MINORITY'] != 'U'].reset_index(drop=True)  # renumber the rows, and don't create a new column
sz = tw.groupby(['MINORITY','TYPE']).size()
sz
pd.crosstab(tw.MINORITY,tw.TYPE)
pd.crosstab(tw.MINORITY,tw.TYPE).plot(kind='bar')
# percent of minority citations which were tickets
sz['M']['T'] / sum(sz['M'])
# percent of non-minority citations which were tickets
sz['W']['T'] / sum(sz['W'])<jupyter_output><empty_output><jupyter_text>We can make a nice little chart of this...<jupyter_code>pd.DataFrame([sz['M']['T'] / sum(sz['M']), sz['W']['T'] / sum(sz['W'])], index=['Minority','White']).plot(kind='bar')<jupyter_output><empty_output><jupyter_text>But how do we know that the drivers in these two groups are really the same? We imagine that there are more tickets given to faster drivers (`MPHOVER` is how much over the speed limit the driver was). Legend also has it that traffic cops are harsher to drivers from out of town (`INTOWN` records where the driver was from.)To get a sense of this, let's compare histograms of the two groups.<jupyter_code>tw['MPHOVER'].hist()<jupyter_output><empty_output><jupyter_text>We've used `MPHOVER` because speed over the limit is probably more interesting than how fast in absolute terms. But checking raw `MPH`, we see a problem with the data: some of the entries don't record the speed! Fortunately, not very many.<jupyter_code>tw['MPH'].hist()<jupyter_output><empty_output><jupyter_text>To see if the distributions are different between white and minorities, let's try plotting them on the same chart.<jupyter_code>w_over = tw[tw['MINORITY']=='W']['MPHOVER'].dropna()  # need dropna because some rows don't record speed
m_over = tw[tw['MINORITY']=='M']['MPHOVER'].dropna()

bins = np.linspace(0, 100, 20)
plt.hist(w_over,bins,alpha=0.5)
plt.hist(m_over,bins,alpha=0.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Still very hard to see, because the scaling is so different. Let's normalize the histograms to display percentages instead of counts.<jupyter_code>plt.hist(w_over,bins,alpha=0.5,density=True)
plt.hist(m_over,bins,alpha=0.5,density=True)<jupyter_output><empty_output><jupyter_text>Hmm, they look pretty close. The minorities plot has a slightly fatter tail, so that's some evidence that they're driving faster, but it doesn't seem to be a huge difference. Let's take a look at the averages, just to get another sense.<jupyter_code>w_over.mean()
m_over.mean()<jupyter_output><empty_output><jupyter_text>So, about half a mile an hour on average difference. Still, this could make the difference if, for example, it pushes a lot of drivers greater than some threshold like 20 mph over the limit. We're going to need to untangle the effects here.In fact there could be all sorts of things going on here. Age and sex could be factors. In town or out of town drivers could matter. The `AGENCY` column records whether the stop was made by Massechusets state police, Boston police, or all other police departments, and that could matter both in terms of how the different agencies operate and in terms of who is driving within their jurisdiction. For our purposes, bias means that minorities get a ticket more often than white people *when everything else is equal.* Is everything else equal? Let's plot a pile of histograms to take a look.<jupyter_code>def comparehist(colname):
    w = tw[tw['MINORITY']=='W'][colname].dropna()
    m = tw[tw['MINORITY']=='M'][colname].dropna()
    
    if tw[colname].dtype == np.float64:
        # it's numeric, plot a histogram
        bins = np.linspace(0, 100, 20)
        plt.hist(w,bins,alpha=0.5,density=True)
        plt.hist(m,bins,alpha=0.5,density=True)
    else:
        # if not numerican, assume it's categorical and plot a bar chart of value counts normalized to percentages
        wc = w.value_counts()
        wc /= sum(wc)
        plt.bar(wc.index, wc, alpha=0.5)
        mc = m.value_counts()
        mc /= sum(mc)
        plt.bar(mc.index, mc, alpha=0.5)
        
    plt.title(colname)
    plt.show()

comparehist('AGE')
comparehist('SEX')
comparehist('MPH')
comparehist('MPHOVER')
comparehist('AGENCY3')
comparehist('INTOWN')<jupyter_output><empty_output><jupyter_text>What can we say here? The minority drivers in this sample have more people in their 20s, skew male, are going slightly faster when they're pulled over, are less likely to be stopped by local police and more likely to be stopped by Boston police, and tend to be out of town more often. Any or all of these effects could account for the difference in ticketing rates. There could even be interactions here: suppose Boston police are harder on people who drive more than 10 mph over the limit. These minority drivers are both more commonly in Boston and going faster on average.  Enter logistic regressionWe're going to use regression to try to analyze the relationships between the variables and whether someone gets a ticket. But everything is going to be a little bit different than last class because the variable we are going to try to predict is binary: it's a yes/no, did they get a ticket or not, instead of a numeric value.<jupyter_code>from sklearn.linear_model import LinearRegression,LogisticRegression<jupyter_output><empty_output><jupyter_text>To set up the regression, let's look at a graph of MPH over the limit on the x axis, vs. whether the driver got a ticket on the y axis, represented as 0 or 1.<jupyter_code># add a column that is 1 if they got a ticket (not a warning)
tw['BUSTED'] = tw['TYPE'].replace({'W':0, 'T':1})

# Remove rows where MPH is unknown
tw = tw[~tw['MPHOVER'].isnull()].reset_index()
# sort by MPHOVER (makes the plotting below simpler)
tw = tw.sort_values(by='MPHOVER')
x = tw[['MPHOVER']].values
y = tw[['BUSTED']].values
tw.plot(kind='scatter',x='MPHOVER',y='BUSTED')<jupyter_output><empty_output><jupyter_text>This has a general "upward slope." Let's use linear regression to try to fit a trend line to it.<jupyter_code>lm = LinearRegression()
lm.fit(x,y)
tw.plot(kind='scatter',x='MPHOVER',y='BUSTED')
plt.plot(x,lm.predict(x),'-',color='m')<jupyter_output>/Users/jonathanstray/anaconda/lib/python3.6/site-packages/sklearn/linear_model/base.py:509: RuntimeWarning: internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver.
  linalg.lstsq(X, y)<jupyter_text>This is a little strange, because we're trying to predict binary outcome with a straight line. If you like, you can imagine that whenever this line is 0.5 or greater, we guess "got a ticket." But really we should use a model that always predicts a value between 0 and 1 -- a probability. To do this we fit a "logistic" function, not a line. This is a function which is an S-shape, and turns any number into a number from 0 to 1<jupyter_code>def logistic(x):
    return 1 / (1 + np.exp(-x))

linear_range = np.linspace(-10,10)
plt.plot(linear_range, logistic(linear_range))<jupyter_output><empty_output><jupyter_text>The exact center of this function, at x=0, has a probability of 0.5, meaning "I'm indifferent guessing between these two options."<jupyter_code>logistic(0)<jupyter_output><empty_output><jupyter_text>Let's apply logistic regression and plot the answer. To interpret the output, we take our normal `y=mx+b` line but run it through the logistic function.<jupyter_code># Run the logistic regression on our data
lm = LogisticRegression()
lm.fit(x,y)

# plot original data
tw.plot(kind='scatter',x='MPHOVER',y='BUSTED')

# Now plot the prediction of the fitted model
m = lm.coef_[0]
b = lm.intercept_
line_y = m*x+b
prediction = logistic(line_y) 
plt.plot(x, prediction, '-', color='m')<jupyter_output>/Users/jonathanstray/anaconda/lib/python3.6/site-packages/sklearn/utils/validation.py:578: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
  y = column_or_1d(y, warn=True)<jupyter_text>We can do the same thing with the `predict_proba` function (though it returns two columns: one is the probability of  FALSE and the other is the probability of TRUE, so we plot only the second column)<jupyter_code>tw.plot(kind='scatter',x='MPHOVER',y='BUSTED')
plt.plot(x, lm.predict_proba(x)[:,1], '-', color='m')<jupyter_output><empty_output><jupyter_text>Odds ratiosHere's a useful fact: the "slope" coefficients of a logistic regression means something when the indepenent variable is binary: it's the *odds ratio*.Suppose we have a dinner party where each person eats chicken or fish. But oh no! Some people got sick afterwards. To see if it was related to what they ate, we can calculate the odds ratio.<jupyter_code># the fish column is 0 if they ate chicken, 1 if fish
dinner = pd.DataFrame([[0,1],[0,0],[0,0],[1,0],[1,1],[1,1]], columns=['fish','sick'])
dinner<jupyter_output><empty_output><jupyter_text>The "odds" of getting sick if you ate chicken are just the number of people who ate chicken and got sick, divided by the number of people who ate chicken and did not get sick.<jupyter_code># useful little function
ct = pd.crosstab(dinner.fish, dinner.sick)
ct
# odds of getting sick if you ate chicken
chicken_odds = ct[0][1] / ct[0][0]
chicken_odds
# odds of getting sick if you ate fish
fish_odds = ct[1][1] / ct[1][0]
fish_odds<jupyter_output><empty_output><jupyter_text>Then "odds ratio" is, yup, the ratio of these sets of odds. It measures how your odds of getting sick *changed* if you ate the fish.<jupyter_code>fish_odds / chicken_odds<jupyter_output><empty_output><jupyter_text>Now let's do a logistic regression to predict who gets sick<jupyter_code>x=dinner[['fish']].values
y=dinner[['sick']].values
lm = LogisticRegression(C=1e5) # this C thing... we just need some large value here, it's a long story
lm.fit(x,y)
# To get the odds ratio, we use the following formula:
np.exp(lm.coef_)<jupyter_output><empty_output><jupyter_text>Controlling for other variablesWe can feed as many independent variables as we like into linear regression. With two variables, instead of fitting a plane to the data we are fitting a curved surface. But here's the cool thing: just like linear regression each coefficient tells us how much that independent variable effects the probability of getting a ticket *if all other variables stay the same.*<jupyter_code># Code these variables a binary
tw.SEX = tw.SEX.replace({'M':0,'F':1, 'U':0})
tw.MINORITY = tw.MINORITY.replace({'W':0,'M':1})
tw.INTOWN = tw.INTOWN.replace({'N':0,'Y':1, 'U':0})

# remove unknown ages
tw = tw[~tw.AGE.isnull()].reset_index()
# Do a logistic regression on many fields at the same time
fields = ['MPHOVER','MINORITY','SEX', 'AGE', 'INTOWN']
x = tw[fields].values
y = tw[['BUSTED']].values
lm = LogisticRegression() 
lm.fit(x,y)
       
# print out the odds ratios for each variable
results = pd.DataFrame(np.exp(lm.coef_), columns=fields)
results<jupyter_output>/Users/jonathanstray/anaconda/lib/python3.6/site-packages/sklearn/utils/validation.py:578: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
  y = column_or_1d(y, warn=True)
{'fraction_non_alphanumeric': 0.10194042088002186, 'fraction_numerical': 0.010795299262093468, 'mean_word_length': 5.352864583333333, 'pattern_counts': {'":': 0, '<': 94, '<?xml version=': 0, '>': 94, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '611484', 'license': 'CC-BY-3.0', 'n_tokens_mistral': 4655, 'n_tokens_neox': 4399, 'n_words': 1774, 'path': 'week-3/week-3-2-class.ipynb', 'repo_name': 'xeophin/lede-algorithms', 'size': '14636'}