Forecasting Daily Bike-Share Demand: a count-regression case study¶
The third end-to-end case study, and the one where the outcome is a count: how many bikes will be rented today? We follow the same 12-step method, on a messy daily-operations export (three date formats, categories in every spelling, missing weather readings, duplicate days, a missing-target row), and meet the signature hazard of count models: overdispersion. Library-first with pandas, seaborn, and statsmodels.
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 statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.model_selection import cross_val_score; from sklearn.linear_model import PoissonRegressor
pd.set_option('display.max_columns', 40)
Objective. Predict the daily rental count and identify what drives demand, so operations can pre-position bikes, schedule staff, and plan maintenance.
Why a count model? Rentals are non-negative whole numbers with no upper cap, so ordinary linear regression is the wrong tool (it can predict negative demand and assumes constant variance). The right family is Poisson, with a log link, and we will see it needs one upgrade.
try: raw = pd.read_csv('../../data/bike_share.csv')
except FileNotFoundError: raw = pd.read_csv(BASE + 'bike_share.csv')
print(raw.shape); raw.head(4)
(745, 9)
| date | season | holiday | workingday | weather | temp_c | humidity | windspeed | rentals | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 2011-06-04 | spring | 0 | 0 | Mist | 23.5 | 65.0 | 15.1 | 2486.0 |
| 1 | 02/03/2011 | Winter | 0 | 1 | Mist | 4.3 | 46.0 | 4.0 | 1343.0 |
| 2 | 2011-05-20 | spring | 0 | 1 | Mist | 24.9 | 39.0 | 18.8 | 1481.0 |
| 3 | 17-Mar-2012 | winter | 0 | 0 | Mist | 10.6 | 32.0 | 22.7 | 907.0 |
raw.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 745 entries, 0 to 744 Data columns (total 9 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 date 745 non-null object 1 season 745 non-null object 2 holiday 745 non-null int64 3 workingday 745 non-null object 4 weather 745 non-null object 5 temp_c 745 non-null float64 6 humidity 727 non-null float64 7 windspeed 736 non-null float64 8 rentals 740 non-null float64 dtypes: float64(4), int64(1), object(4) memory usage: 52.5+ KB
print('season spellings :', raw.season.unique()[:8])
print('weather spellings :', raw.weather.unique())
print('workingday values :', raw.workingday.unique())
print('date formats :', raw.date.head(4).tolist())
print('duplicate dates :', raw.date.duplicated().sum(), ' | rentals dtype:', raw.rentals.dtype)
print('missing:'); print(raw.isna().sum()[lambda s: s>0].to_string())
season spellings : ['spring' 'Winter' 'winter' 'Summer' 'fall' 'FALL' 'SPRING' 'summer'] weather spellings : ['Mist' 'Clear' 'Heavy Rain' 'Light Rain' 'clear' 'mist' 'light rain' 'heavy rain'] workingday values : ['0' '1' 'Yes' 'No'] date formats : ['2011-06-04', '02/03/2011', '2011-05-20', '17-Mar-2012'] duplicate dates : 14 | rentals dtype: float64 missing: humidity 18 windspeed 9 rentals 5
The usual real-world mess: season and weather in mixed case, workingday as a blend of 0/1 and Yes/No, three different date formats, 14 duplicate days, a handful of missing humidity and windspeed readings, and a few rows with no rentals (our target).
df = raw.dropna(subset=['rentals']).drop_duplicates(subset='date').copy()
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') # parse all three formats at once
df['dt'] = dt; df['month'] = dt.dt.month; df['is_weekend'] = (dt.dt.dayofweek >= 5).astype(int)
print(f'rows {len(raw)} -> {len(df)} (dropped 5 missing-target + 14 duplicates)')
print('seasons:', sorted(df.season.unique()), '| weather:', sorted(df.weather.dropna().unique()))
rows 745 -> 726 (dropped 5 missing-target + 14 duplicates) seasons: ['Fall', 'Spring', 'Summer', 'Winter'] | weather: ['Clear', 'Heavy Rain', 'Light Rain', 'Mist']
fig, ax = plt.subplots(1, 3, figsize=(14,4))
ax[0].hist(df.rentals, bins=35, color=EM, alpha=0.8); ax[0].set(title='Daily rentals (a count)', xlabel='rentals', ylabel='days')
df.groupby('weather').rentals.mean().reindex(['Clear','Mist','Light Rain','Heavy Rain']).plot(kind='bar', color=[EM,AMBER,RED,'#7f1d1d'], ax=ax[1])
ax[1].set(title='Rentals by weather', ylabel='avg rentals'); ax[1].tick_params(axis='x', rotation=20)
ax[2].scatter(df.temp_c, df.rentals, s=10, color=EM, alpha=0.4); ax[2].set(title='Rentals vs temperature', xlabel='temp (C)', ylabel='rentals')
plt.tight_layout(); plt.show()
# rentals over time: the strong seasonal rhythm across two years
ts = df.sort_values('dt')
fig, ax = plt.subplots(figsize=(12,3.6))
ax.plot(ts.dt, ts.rentals, color=EM, lw=0.9, alpha=0.85)
ax.set(title='Daily rentals over two years: a clear seasonal cycle', xlabel='date', ylabel='rentals')
plt.tight_layout(); plt.show()
Demand rises steeply with good weather (Clear days average far more than Heavy Rain days) and with temperature, up to a point. The count distribution is right-skewed and, crucially, its spread grows with its level, our first hint of overdispersion.
We already engineered month and is_weekend from the date. Weather and season are nominal, so they enter as dummy variables via C(...). Because demand rises with temperature and then flattens off, we add a squared temperature term written I((temp_c-14)**2) (a touch of the flexible-model idea from earlier chapters). And the outcome is a count, so the model family is Poisson with a log link, not ordinary least squares.
Why subtract 14 inside the square? We center the temperature at the middle of the data (the average day is about 14°C) before squaring. Centering strips out most of the correlation between the linear temp_c term and its square, which keeps both coefficients stable and easy to read. Note that the 14 does not change the fitted curve or the accuracy at all (the linear term absorbs the shift, algebraically (temp_c-14)**2 just adds a constant and a multiple of temp_c); it is purely a choice for numerical stability and interpretation, and any value near the mean would serve. The I(...) wrapper tells the formula parser to treat **2 as literal arithmetic (square it) rather than as a formula operator. The cell below verifies, with numbers and a picture, that 14 really is the center.
# WHERE DOES 14 COME FROM? verify it centers the temperature
mid = (df.temp_c.min() + df.temp_c.max()) / 2
print(f'temp_c: mean {df.temp_c.mean():.1f} C | median {df.temp_c.median():.1f} C | range midpoint {mid:.1f} C -> all about 14')
r_raw = np.corrcoef(df.temp_c, df.temp_c**2)[0,1]
r_ctr = np.corrcoef(df.temp_c, (df.temp_c-14)**2)[0,1]
print(f'corr(temp, temp^2) = {r_raw:+.2f} (uncentered square climbs with temp -> collinear)')
print(f'corr(temp, (temp-14)^2) = {r_ctr:+.2f} (centered square is symmetric -> not collinear)')
fig, ax = plt.subplots(1, 2, figsize=(12,4.2))
ax[0].hist(df.temp_c, bins=30, color=EM, alpha=0.8)
ax[0].axvline(df.temp_c.mean(), color=INK, ls='--', lw=2, label=f'mean = {df.temp_c.mean():.1f} C')
ax[0].axvline(14, color=RED, lw=2, label='centering value = 14')
ax[0].set(title='14 sits at the center of the temperature data', xlabel='temp (C)', ylabel='days'); ax[0].legend()
grid = np.linspace(df.temp_c.min(), df.temp_c.max(), 120)
ax[1].plot(grid, grid**2/(grid**2).max(), color=GREY, ls='--', lw=2.2, label=f'temp^2 (rescaled): climbs with temp, r={r_raw:+.2f}')
ax[1].plot(grid, (grid-14)**2/((grid-14)**2).max(), color=EM, lw=2.6, label=f'(temp-14)^2 (rescaled): symmetric U at 14, r={r_ctr:+.2f}')
ax[1].axvline(14, color=RED, lw=2)
ax[1].set(title='Centering turns the squared term into a U centered on 14', xlabel='temp (C)', ylabel='squared term (rescaled 0-1)'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
temp_c: mean 13.8 C | median 13.2 C | range midpoint 13.3 C -> all about 14 corr(temp, temp^2) = +0.97 (uncentered square climbs with temp -> collinear) corr(temp, (temp-14)^2) = -0.06 (centered square is symmetric -> not collinear)
Reading the check. Left: the mean, median, and range midpoint of temperature all land at about 14°C, so 14 genuinely is the center of the data (the red line sits right on the average). Right: subtracting 14 before squaring turns the term into a symmetric U whose low point is exactly 14, while the raw temp^2 just climbs with temperature. That difference is the whole point: the correlation between temperature and its square falls from about +0.97 (raw, badly collinear) to roughly 0 (centered), so the straight-line piece and the curved piece no longer fight each other for the same signal. Any center near the mean gives the same benefit, 14 is simply the clean round number at the middle of this data.
print('mean rentals =', round(df.rentals.mean()), ' variance =', round(df.rentals.var()))
print('variance / mean =', round(df.rentals.var()/df.rentals.mean(),1), ' (Poisson expects this to be ~1)')
print('\navg rentals by weather:'); print(df.groupby('weather').rentals.mean().round(0).to_string())
mean rentals = 1509 variance = 1124117 variance / mean = 744.9 (Poisson expects this to be ~1) avg rentals by weather: weather Clear 1711.0 Heavy Rain 369.0 Light Rain 945.0 Mist 1515.0
The variance-to-mean ratio is in the hundreds, not near 1. A pure Poisson model assumes they are equal, so it will badly understate the uncertainty. We fit Poisson first (to see the problem), then fix it.
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()
print(f'Poisson: pseudo-R2 = {1-pois.deviance/pois.null_deviance:.3f}')
print(f'Pearson chi2 / df = {pois.pearson_chi2/pois.df_resid:.1f} (should be ~1 if Poisson fits)')
Poisson: pseudo-R2 = 0.634 Pearson chi2 / df = 255.0 (should be ~1 if Poisson fits)
The fit explains the data well, but that Pearson chi-square / df in the hundreds is a screaming overdispersion warning: the Poisson standard errors are far too small. Step 9 fixes it.
nb = smf.glm(form, df, family=sm.families.NegativeBinomial(alpha=1/6)).fit()
print(f'AIC Poisson = {pois.aic:,.0f} Negative Binomial = {nb.aic:,.0f} (much lower is much better)')
print(f'NB pseudo-R2 = {1-nb.deviance/nb.null_deviance:.3f}')
Xv = sm.add_constant(df[['temp_c','humidity','windspeed']])
print('VIF:', {c: round(variance_inflation_factor(Xv.values,i),1) for i,c in enumerate(Xv.columns) if c!='const'})
AIC Poisson = 182,804 Negative Binomial = 11,083 (much lower is much better)
NB pseudo-R2 = 0.657
VIF: {'temp_c': np.float64(1.0), 'humidity': np.float64(1.0), 'windspeed': np.float64(1.0)}
# overdispersion, made visual, and the model's forecast accuracy
df['mu'] = pois.mu
g = df.assign(b=pd.qcut(df.mu,10,duplicates='drop')).groupby('b', observed=True).rentals.agg(['mean','var'])
fig, ax = plt.subplots(1, 2, figsize=(12,4.4))
ax[0].scatter(g['mean'], g['var'], s=55, color=EM, zorder=3); hi=max(g['mean'].max(), g['var'].max())
ax[0].plot([0,hi],[0,hi], ls='--', color=GREY, label='Poisson: var = mean'); ax[0].legend()
ax[0].set(title='Overdispersion: variance >> mean', xlabel='mean', ylabel='variance')
ax[1].scatter(df.rentals, nb.fittedvalues, s=12, color=EM, alpha=0.4)
lims=[df.rentals.min(), df.rentals.max()]; ax[1].plot(lims, lims, ls='--', color=INK)
ax[1].set(title='Predicted vs actual (negative binomial)', xlabel='actual', ylabel='predicted')
plt.tight_layout(); plt.show()
Validation verdict. The variance sits far above the mean (left), confirming overdispersion, so we switch to the negative binomial, which adds a dispersion parameter and slashes the AIC by orders of magnitude. VIFs for temperature, humidity and wind are about 1 (no multicollinearity among the weather inputs). The predicted-vs-actual plot (right) tracks typical days well but underpredicts the busiest ones (deviance pseudo-R-squared about 0.66). Same rate ratios as Poisson, but now with honest standard errors.
RR = np.exp(nb.params)
print(f"Heavy rain vs clear : x{RR['C(weather)[T.Heavy Rain]']:.2f} (demand collapses)")
print(f"Light rain vs clear : x{RR['C(weather)[T.Light Rain]']:.2f}")
print(f"Each +5 C, 11.5->16.5 C: x{np.exp(nb.params['temp_c']*5):.2f} (the squared term bends it: smaller when already warm)")
print(f"Working day : x{RR['workingday']:.2f}")
print(f"Holiday : x{RR['holiday']:.2f}")
ex = pd.DataFrame({'temp_c':[22],'humidity':[55],'windspeed':[10],'weather':['Clear'],'season':['Summer'],'workingday':[1],'holiday':[0]})
print(f"\nForecast, warm clear working day: {nb.predict(ex).iloc[0]:,.0f} rentals")
Heavy rain vs clear : x0.25 (demand collapses) Light rain vs clear : x0.54 Each +5 C, 11.5->16.5 C: x1.33 (the squared term bends it: smaller when already warm) Working day : x1.10 Holiday : x0.57 Forecast, warm clear working day: 3,033 rentals
# rate-ratio forest: each factor's multiplier on demand
terms = ['C(weather)[T.Heavy Rain]','C(weather)[T.Light Rain]','holiday','C(weather)[T.Mist]','workingday']
labels = ['Heavy rain','Light rain','Holiday','Mist','Working day']
rr = np.exp(nb.params[terms]); lo = np.exp(nb.conf_int().loc[terms,0]); hi = np.exp(nb.conf_int().loc[terms,1])
order = rr.sort_values().index; yy = np.arange(len(order))
fig, ax = plt.subplots(figsize=(7.8,4))
for i,t in enumerate(order):
col = RED if hi[t]<1 else (EM if lo[t]>1 else GREY)
ax.plot([lo[t],hi[t]],[i,i], color=col, lw=2.5); ax.scatter(rr[t], i, color=col, s=55, zorder=3)
ax.axvline(1, color=INK, ls='--'); ax.set_yticks(yy); ax.set_yticklabels([labels[list(rr.index).index(t)] for t in order])
ax.set(title='Rate ratios: demand multiplier vs a clear, non-working day', xlabel='rate ratio'); plt.tight_layout(); plt.show()
What the model says. Weather rules demand: a heavy-rain day sees only about a quarter of clear-day rentals, light rain about half. Warmth helps (+33% for 5°C near the typical 14°C, less in warm weather), working days run about 10% busier than weekends, and holidays about 43% quieter than ordinary working days.
Actionable: pre-position fewer bikes and staff on forecast rain and cold; scale up for warm, clear working days; plan holidays below even a weekend. The forecast feeds directly into the next day's rebalancing plan.
- Persist the fitted model and cleaning as one pipeline (
joblib). - Score daily: each evening, feed tomorrow's weather forecast (temp, conditions, wind) plus the calendar (weekday, holiday) to get a predicted rental count, and a range.
- Feed operations: the forecast drives bike rebalancing, staff scheduling, and maintenance windows (do it on low-demand rainy days).
- Monitor and retrain: track forecast error, and retrain as ridership grows or the system expands; watch for a changing base level (a growing city).
- Guardrail: report the prediction interval, not just a point count, and do not extrapolate to weather far outside the training range.
What we did. We took two years of daily rental records, cleaned them (removed duplicate days and a few unusable rows, fixed inconsistent labels, filled in a handful of missing weather readings, and standardized the dates), and built a formula that predicts how many bikes will be rented on a given day.
How good is it? It explains about two-thirds of the deviance in daily demand (a likelihood-based measure that runs above the ordinary R-squared). It tracks typical days well and underpredicts the busiest ones, so a forecast should come with a range, as the Take It Further notebook builds.
What drives demand: weather first (a rainy day can cut rentals by 75%), then temperature (warmer is busier), then the calendar (working days busier than weekends and holidays).
A note on the method. Daily counts bounce around far more than a simple model expects, so we used a model built for that (a 'negative binomial'), which gives honest confidence in the forecast.
Bottom line: let tomorrow's weather forecast set tomorrow's bike plan, scale down for rain and cold, scale up for warm clear working days.