Analyzing the NYC Subway Dataset¶
Student: Ada Lee
Date: 2015-10-2
Section 0. References¶
[1] http://stackoverflow.com/questions/3541713/how-to-plot-two-histograms-together-in-r
[2] https://en.wikipedia.org/wiki/Dummy_variable_(statistics)
[3] http://statsmodels.sourceforge.net/devel/generated/statsmodels.regression.linear_model.OLS.html
[4] http://pandas.pydata.org/pandas-docs/stable/generated/pandas.get_dummies.html
[5] http://pandas.pydata.org/pandas-docs/stable/merging.html
[6] http://blog.minitab.com/blog/adventures-in-statistics/regression-analysis-how-do-i-interpret-r-squared-and-assess-the-goodness-of-fit
[7] http://docs.ggplot2.org/current/geom_abline.html
Section 1. Statistical Test¶
import pandas as pd
import numpy as np
from ggplot import *
turnstile_weather = pd.read_csv('data/turnstile_weather_v2.csv')
entries_with_rain = turnstile_weather['ENTRIESn_hourly'][turnstile_weather['rain'] == 1]
entries_without_rain = turnstile_weather['ENTRIESn_hourly'][turnstile_weather['rain'] == 0]
avg_entries_with_rain = np.mean(entries_with_rain)
avg_entries_without_rain = np.mean(entries_without_rain)
print "average entries with and without rain:{}, {}\n".format(avg_entries_with_rain, avg_entries_without_rain)
%matplotlib inline
ggplot(aes(x='ENTRIESn_hourly', fill='rain'), data=turnstile_weather) + geom_density(alpha = 0.2) + ylab('entries density') + ggtitle('Subways Entries Density for Rainy Days and Non-rainy Days')
Choose statistical test¶
From the graph we can see that both distribution are right-skewed, so we could not use z-test or t-test. We also don't know distribution of the entries population for rainy days and non-raniy days. The Mann Whitney U test is a non-parametric test that can be used to test, for two populations with unknown distributions. So Mann Whitney U test would be appropriate for our choice. I would prefer to use two-tailed test.
Hypothesis¶
Giving random draws x from entries population in rainy days, y from entries population in non-rainy days, we have below null hypothesis and alternative hypothesis
$$H_0: P(x > y) = 0.5$$
$$H_a: P ( x > y ) \neq 0.5$$
Calculate p value¶
import scipy.stats
U, one_tailed_p = scipy.stats.mannwhitneyu(entries_with_rain, entries_without_rain)
p = one_tailed_p * 2
print "p value for two-tailed test is", p
print "\nmean entries in rainy days is", avg_entries_with_rain
print "mean entries in non-rainy days is", avg_entries_without_rain
Significance and interpretation¶
We use two-tailed test and get p value $p=5.48\times 10^{-6},\; p << 0.01$, so we reject the null hypothesis, thus we accept the alternative hypothesis $P ( x > y ) \neq 0.5$, in other words probability of entries in rainy days larger than enties in non-rainy days is not equal to 0.5. This means that rain do affect subway entries.
Section 2. Linear Regression¶
I use OLS for linear regression, and try to add information to features as more as possible, so I would add all numeric variables about the weather, they are precipi, pressurei, tempi, wspdi, meanprecipi, meanpressurei, meantempi, meanwspdi. For dummy variables, UNIT and station are similar, and I find that only use UNIT as feature would result in $R^2$ 0.375, only use station as feature would result in $R^2$ 0.324, so I prefer to add UNIT as one of the dummy variables. For dummy variables about dates, datetime's information already cover DATEn, TIMEn, hour, day_week, weekday, so I would only add datetime as one of dummy variables. Finally I would add dummy variable conds, because adding it increases $R^2$ for my model. For weather_lat and weather_lon, I don't think it would affect entries, and acturally I also try to add these, but it does not improve $R^2$, so I would not contain it to my final model.
import statsmodels.api as sm
def linear_regression(features, values):
"""
Perform linear regression given a data set with an arbitrary number of features.
"""
features = sm.add_constant(features)
model = sm.OLS(values, features)
results = model.fit()
return results
values = turnstile_weather['ENTRIESn_hourly']
features = pd.get_dummies(turnstile_weather['UNIT'], prefix='unit')
print "R square for model that only contain UNIT as feature:", linear_regression(features, values).rsquared
features = pd.get_dummies(turnstile_weather['station'], prefix='station')
print "R square for model that only contain station as feature:", linear_regression(features, values).rsquared
## regression without cond
numerics = ['precipi', 'pressurei','tempi', 'wspdi',
'meanprecipi', 'meanpressurei', 'meantempi', 'meanwspdi']
dummies = ['fog', 'rain']
features = turnstile_weather[numerics + dummies]
dummy_units = pd.get_dummies(turnstile_weather['UNIT'], prefix='unit')
features = features.join(dummy_units)
dummy_units = pd.get_dummies(turnstile_weather['datetime'])
features = features.join(dummy_units)
print "R square for model that does not contain cond as dummy feature:", linear_regression(features, values).rsquared
## add cond
dummy_units = pd.get_dummies(turnstile_weather['conds'])
features = features.join(dummy_units)
model_results = linear_regression(features, values) # This is my final model, so I would like to keep it
print "R square for model that add cond as dummy feature:", model_results.rsquared
## add weather_lat and weather_lon
features.add(turnstile_weather[['weather_lat', 'weather_lon']])
print "R square for model that add weather_lat and weather_lon:", linear_regression(features, values).rsquared
Visualzie to see how the model works¶
intercept = model_results.params[0]
params = model_results.params[1:]
predictions = intercept + np.dot(features, params)
values_and_predictions = pd.DataFrame({'entries':values, 'predictions':predictions})
ggplot(aes(x='entries', y='predictions'), data=values_and_predictions) + geom_point() + geom_abline(slope=1, intercept=0, color='red') + ggtitle('Predictions vs. Real Entries')
Look like our model is not very good, we will discuss later.
Parameters of the non-dummy features in the linear regression model and parameter for some dummy variables.¶
from pprint import pprint
print "parameters for non-dummy features:\n"
pprint(params[numerics])
print "\n\nparameter for some dummy variables:\n"
pprint(params[dummies])
$R^2$¶
The definition of R-squared is fairly straight-forward; it is the percentage of the response variable variation that is explained by a linear model.
rsquared = model_results.rsquared
print "rsquared:", rsquared
$R^2$ is roughly 0.575, that means that 57.5% variance of entries are expalined by our linear model.
Plot residuals¶
residuals = pd.DataFrame({'residuals': np.array(predictions) - np.array(values)})
ggplot(aes(x='residuals'), data=residuals) + geom_density() + ggtitle('Residuals of Predictions on Entries') + ylab('density')
We can see that some residuals are very large, and residuals are not symmetric, this make us to question linearity of features and values. Linear model look like not look very well for this dataset.
ggplot(aes(x='ENTRIESn_hourly'), data=turnstile_weather[turnstile_weather['rain'] == 0]) + geom_histogram(fill='blue', binwidth=1000) + ylim(0, 20000) + ggtitle('Non-rainy Day Subway Entries') + ylab('frequency')
ggplot(aes(x='ENTRIESn_hourly'), data=turnstile_weather[turnstile_weather['rain'] == 1]) + geom_histogram(fill='green', binwidth=1000) + ylim(0, 20000) + ggtitle('Rainy Day Subway Entries') + ylab('frequency')
import collections as c
c.Counter(turnstile_weather['rain'])
We see that both plots are right-skewed, and there are a lot more data for rainy days. Actually there are 33064 non-rainy datas, and 9585 rainy datas.
Ridership by time-of-day¶
from pandasql import sqldf
q = '''
select hour, avg(ENTRIESn_hourly) as avg_entries
from turnstile_weather
group by hour'''
avg_entries_for_hours = sqldf(q.lower(), locals())
avg_entries_for_hours
ggplot(aes(x='hour', y='avg_entries'), data=avg_entries_for_hours) + geom_line() + geom_point() + xlim(-.5, 20.5) + ggtitle('Average Ridership by Time-of-day') + scale_x_continuous(breaks=range(0, 21, 4))
We can see that in hour 20 people take subway most, and in hour 12 people take subway the second most, in hour 4 people take subway least frequent.
Section 4. Conclusion¶
params['rain']
From our Mann Whitney U-test, we see that rain do affect subways entries. From our linear model, we can coefficient of rain is 9.64, that means that when it is rainy, more people tend to take subway. In conclusion, rain significantly increase NYC subway entries.
Section 5. Reflection¶
All our data is from year 2011, month May, so there should be 31 days of date for each unit, but among 42649 units, there are 8 units which does not have date on all 31 days in May, specially unit 'R459' ouly have 16 days of data. Among hours 0, 4, 8, 12, 16, 20, hour 8 has a lot less data than other hours. These two kinds of data missing a significant proportion is shortcoming of data I found. Linear regression seems not be very good in this question since we only get R2 around 0.5.
unit_vs_date = {key:set() for key in turnstile_weather['UNIT']}
for unit, date in zip(turnstile_weather['UNIT'], turnstile_weather['DATEn']):
unit_vs_date[unit].add(date)
unit_vs_days = [(unit, len(unit_vs_date[unit])) for unit in sorted(unit_vs_date.keys())]
print 'total units:', len(turnstile_weather['UNIT'])
print '\n units that has less than 31 days date:'
pprint(filter(lambda x: x[1] != 31, unit_vs_days))
print '\ndata count for hours:'
pprint(sorted(c.Counter(turnstile_weather['hour']).items()))
unit_vs_date['R459']