Bike-Share Demand: Take It Further¶
This notebook continues the Chapter 105 case study. We recap the count model, then work the five extensions from the chapter's Take It Further section: quantifying Poisson vs negative-binomial standard errors, testing the temperature curve, testing on a later year, a weekend-weather interaction, and a week-ahead forecast. Each is visual and explained, closing with what the model says in plain terms and the Step 12 plain-English write-up.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns # seaborn = high-level statistical plots (heatmaps, regplots, pairplots)
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.nonparametric.smoothers_lowess import lowess
from sklearn.linear_model import LinearRegression, Ridge, Lasso, LogisticRegression
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import r2_score, mean_squared_error
ORG="#059669"; DEEP="#047857"; LIGHT="#6ee7b7"; EM="#059669"; EMDEEP="#047857"; INK="#1a2138"; GRID="#e6e9f2"; GREEN="#059669"; RED="#ef4444"; AMBER="#d97706"; BLUE="#2563eb"; PUR="#9333ea"; GREY="#94a3b8"; SLATE="#475569"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
"axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
"axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
sns.set_style("whitegrid")
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
import warnings; warnings.filterwarnings('ignore')
import statsmodels.api as sm, statsmodels.formula.api as smf
from sklearn.model_selection import cross_val_score; from sklearn.linear_model import PoissonRegressor
pd.set_option('display.max_columns', 40)
def load_clean():
try: df = pd.read_csv('../../data/bike_share.csv')
except FileNotFoundError: df = pd.read_csv(BASE + 'bike_share.csv')
df = df.dropna(subset=['rentals']).drop_duplicates(subset='date'); df['rentals']=df['rentals'].astype(int)
norm = lambda s: s.astype(str).str.strip().str.lower()
df['season'] = norm(df['season']).str.capitalize()
df['weather'] = norm(df['weather']).map({'clear':'Clear','mist':'Mist','light rain':'Light Rain','heavy rain':'Heavy Rain'})
df['workingday'] = df['workingday'].astype(str).str.strip().str.upper().map({'1':1,'YES':1,'0':0,'NO':0})
for c in ['humidity','windspeed']: df[c]=pd.to_numeric(df[c],errors='coerce'); df[c]=df[c].fillna(df[c].median())
dt = pd.to_datetime(df['date'], format='mixed', errors='coerce'); df['is_weekend']=(dt.dt.dayofweek>=5).astype(int)
return df
df = load_clean()
FORM = 'rentals ~ temp_c + I((temp_c-14)**2) + humidity + windspeed + C(weather) + C(season) + workingday + holiday'
pois = smf.glm(FORM, df, family=sm.families.Poisson()).fit()
nb = smf.glm(FORM, df, family=sm.families.NegativeBinomial(alpha=1/6)).fit()
print(f'n={len(df)} NB pseudo-R2={1-nb.deviance/nb.null_deviance:.3f} NB AIC={nb.aic:,.0f} Poisson AIC={pois.aic:,.0f}')
n=726 NB pseudo-R2=0.657 NB AIC=11,083 Poisson AIC=182,804
cmp = pd.DataFrame({'rate_ratio': np.exp(nb.params).round(3),
'SE_poisson': pois.bse.round(4), 'SE_negbin': nb.bse.round(4)})
cmp['SE_ratio'] = (nb.bse / pois.bse).round(1) # from the unrounded SEs
print(cmp[['rate_ratio','SE_poisson','SE_negbin','SE_ratio']].to_string())
fig, ax = plt.subplots(figsize=(7.6,4))
terms=[t for t in nb.params.index if t!='Intercept']
ax.barh(range(len(terms)), (nb.bse[terms]/pois.bse[terms]), color=EM)
ax.set_yticks(range(len(terms))); ax.set_yticklabels([t.replace('C(weather)[T.','').replace('C(season)[T.','').replace(']','') for t in terms], fontsize=8)
ax.axvline(1, color=INK, ls='--'); ax.set(title='Negative-binomial SE / Poisson SE (how much Poisson understated)', xlabel='ratio'); plt.tight_layout(); plt.show()
rate_ratio SE_poisson SE_negbin SE_ratio Intercept 1327.226 0.0064 0.0992 15.5 C(weather)[T.Heavy Rain] 0.250 0.0114 0.0927 8.1 C(weather)[T.Light Rain] 0.545 0.0034 0.0454 13.2 C(weather)[T.Mist] 0.898 0.0022 0.0357 16.2 C(season)[T.Spring] 0.813 0.0034 0.0518 15.3 C(season)[T.Summer] 1.029 0.0039 0.0669 17.0 C(season)[T.Winter] 0.862 0.0039 0.0495 12.6 temp_c 1.059 0.0003 0.0037 14.5 I((temp_c - 14) ** 2) 0.999 0.0000 0.0003 14.6 humidity 0.994 0.0001 0.0011 16.2 windspeed 0.988 0.0002 0.0031 15.5 workingday 1.100 0.0022 0.0336 15.5 holiday 0.565 0.0094 0.1205 12.8
What we learn. The rate ratios are essentially the same under both models, but the negative-binomial standard errors are several times larger than Poisson's across the board. Under overdispersion, Poisson's error bars are fiction, they would make weak effects look significant and tight forecasts look certain. The negative binomial reports the honest uncertainty. This is exactly the count-model parallel of using robust standard errors for heteroscedasticity.
nb_lin = smf.glm('rentals ~ temp_c + humidity + windspeed + C(weather) + C(season) + workingday + holiday', df, family=sm.families.NegativeBinomial(alpha=1/6)).fit()
print(f'AIC with temp^2 = {nb.aic:,.0f} without = {nb_lin.aic:,.0f} (lower is better)')
grid = pd.DataFrame({'temp_c':np.linspace(df.temp_c.min(), df.temp_c.max(),60), 'humidity':df.humidity.mean(),
'windspeed':df.windspeed.mean(), 'weather':'Clear','season':'Summer','workingday':1,'holiday':0})
fig, ax = plt.subplots(figsize=(7.6,4))
ax.plot(grid.temp_c, nb.predict(grid), color=EM, lw=2.6, label='with temp squared (curved)')
ax.plot(grid.temp_c, nb_lin.predict(grid), color=GREY, lw=2, ls='--', label='linear temp')
ax.set(title='Fitted demand vs temperature', xlabel='temperature (C)', ylabel='predicted rentals'); ax.legend()
plt.tight_layout(); plt.show()
AIC with temp^2 = 11,083 without = 11,098 (lower is better)
What we learn. The squared term earns its place: the AIC is lower with it, and the curve shows demand rising with warmth, more slowly in the heat (its peak would sit above the warmest day in the data), which the straight-line version misses. This is the flexible-model idea from earlier chapters, applied inside a GLM: a quadratic feature lets a linear model bend.
Why I((temp_c-14)**2) subtracts 14. We center the temperature at the middle of the data (the average day is about 14°C) before squaring. Centering removes most of the correlation between the linear temp_c term and its square (in this data it falls from about +0.97 to roughly 0), keeping both coefficients stable and readable. It does not change the fitted curve or the accuracy (the linear term absorbs the shift); any value near the mean works. The Chapter 105 main notebook has a plot verifying that 14 is the center. I(...) just tells the formula parser to square the number literally rather than read ** as a formula operator.
X = pd.get_dummies(df[['temp_c','humidity','windspeed','weather','season','workingday','holiday','is_weekend']], drop_first=True).astype(float)
X['temp_sq'] = (df['temp_c']-14)**2
# forecasting needs a time-ordered split: train on 2011, test on 2012, and score both with the same R2
year = pd.to_datetime(df['date'], format='mixed', errors='coerce').dt.year
tr, te = df[year==2011], df[year==2012]
fit = smf.glm(FORM, tr, family=sm.families.NegativeBinomial(alpha=1/6)).fit()
r2 = lambda d: 1 - ((d.rentals - fit.predict(d))**2).sum() / ((d.rentals - d.rentals.mean())**2).sum()
print(f'trained on 2011 ({len(tr)} days): R2 = {r2(tr):.3f} tested on 2012 ({len(te)} days): R2 = {r2(te):.3f}')
trained on 2011 (363 days): R2 = 0.548 tested on 2012 (363 days): R2 = 0.557
What we learn. Trained on 2011 and scored on 2012 with the same R-squared, the model does about as well on the new year (0.557) as on the year it learned from (0.548), so it generalizes forward in time. A forecast must be tested on later data than it trained on, because a random split lets the model peek at the future, and a big gap would have warned of overfitting. Always check out-of-sample before trusting a forecast operationally.
mi = smf.glm('rentals ~ temp_c + I((temp_c-14)**2) + humidity + windspeed + C(weather)*is_weekend + workingday + holiday', df, family=sm.families.NegativeBinomial(alpha=1/6)).fit()
ip = mi.pvalues.filter(like=':is_weekend')
print('weather x weekend interaction p-values:'); print(ip.round(3).to_string())
print(f'\nsignificant at 0.05: {(ip<0.05).sum()} of {len(ip)}')
weather x weekend interaction p-values: C(weather)[T.Heavy Rain]:is_weekend 0.272 C(weather)[T.Light Rain]:is_weekend 0.060 C(weather)[T.Mist]:is_weekend 0.984 significant at 0.05: 0 of 3
What we learn. The interaction terms are not clearly significant (the closest, light rain on weekends, has p = 0.06), so the data give no clear evidence that a rainy Saturday is hit harder than a rainy Tuesday in proportional terms. We tested the idea rather than assuming it, and the simpler main-effects model stands.
forecast = pd.DataFrame({
'day':['Mon','Tue','Wed','Thu','Fri','Sat','Sun'],
'temp_c':[24,25,19,17,22,26,23], 'humidity':[55,52,70,80,58,50,60], 'windspeed':[8,10,15,20,9,7,11],
'weather':['Clear','Clear','Light Rain','Heavy Rain','Mist','Clear','Clear'],
'season':['Summer']*7, 'workingday':[1,1,1,1,1,0,0], 'holiday':[0]*7})
pr = nb.get_prediction(forecast).summary_frame(alpha=0.05)
forecast['predicted'] = pr['mean'].round(0).astype(int)
# a 95% range for one actual day: the negative-binomial spread around each mean (weather inputs taken as known)
from scipy.stats import nbinom
size = 6.0 # 1/alpha
forecast['low'] = nbinom.ppf(0.025, size, size/(size+pr['mean'])).astype(int)
forecast['high'] = nbinom.ppf(0.975, size, size/(size+pr['mean'])).astype(int)
print(forecast[['day','weather','temp_c','predicted','low','high']].to_string(index=False))
fig, ax = plt.subplots(figsize=(8,4))
ax.bar(forecast.day, forecast.predicted, color=EM, alpha=0.85)
ax.errorbar(forecast.day, pr['mean'], yerr=[pr['mean']-forecast['low'], forecast['high']-pr['mean']], fmt='none', ecolor=INK, capsize=4)
ax.set(title='Week-ahead forecast with 95% ranges for each day', ylabel='rentals'); plt.tight_layout(); plt.show()
day weather temp_c predicted low high Mon Clear 24 3324 1218 6467 Tue Clear 25 3408 1249 6630 Wed Light Rain 19 1252 458 2437 Thu Heavy Rain 17 462 168 901 Fri Mist 22 2706 991 5265 Sat Clear 26 3343 1225 6503 Sun Clear 23 2732 1001 5315
What we learn. The model becomes a planning tool: a warm, clear day forecasts far more rentals than the cold heavy-rain day, and the mid-week storm is clearly the low point to schedule maintenance. Each bar comes with a 95% range for the actual count on that day, so operations can plan for the likely spread, not a false-precision single number. The range assumes the weather forecast itself is right.
What the model says, in plain terms¶
- Weather is the biggest lever on demand, a rainy day can cut rentals by up to three-quarters, and warmer days are busier (about +33% per 5°C in mild weather, less when it is already warm).
- The calendar matters: working days are busier than weekends, and holidays are quieter than either.
- Daily counts are volatile, far more than a naive model assumes, so we used a negative-binomial model that reports honest confidence, and a test on 2012, after training on 2011, confirmed it forecasts a new year about as well.
- The extensions reinforced the story: the temperature curve is real, the negative binomial's error bars are the trustworthy ones, and the data give no clear evidence that weather hits weekends and weekdays differently.
Step 12 - the plain-English write-up¶
What we did. We took two years of daily rental records, cleaned them, and built a formula that forecasts how many bikes will be rented on a given day from the weather and the calendar, then stress-tested it five ways.
How good is it? It explains a little over half of the day-to-day variance in rentals (R-squared about 0.55) and holds up on a year it has not seen, so it can plan a future week, not just describe the past.
What drives demand: weather first (rain and cold cut it sharply), then temperature, then the calendar.
How sure are we? We used a model built for volatile counts (so the confidence is honest), confirmed the temperature effect curves, checked it generalizes, and produced a week-ahead forecast with ranges, ready to drive tomorrow's bike-rebalancing and staffing plan.
Bottom line: let the weather forecast set the bike plan, scale down for rain and cold, up for warm clear working days, and use the low-demand rainy days for maintenance.