Bike rental: data and features

Examples

We want to predict how many bicycles are rented in a given hour in Washington DC, from the weather, the time of the day, the day of the week and the season.

Loading

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml

bikes = fetch_openml(data_id=42712, as_frame=True, parser="auto").frame
print(bikes.shape)
bikes.head()

(17379, 13)

	season	month	hour	holiday	weekday	workingday	weather	temp	feel_temp	humidity	count
0	spring	1	0	False	6	False	clear	9.84	14.395	0.81	16
1	spring	1	1	False	6	False	clear	9.02	13.635	0.80	40
2	spring	1	2	False	6	False	clear	9.02	13.635	0.80	32
3	spring	1	3	False	6	False	clear	9.84	14.395	0.75	13
4	spring	1	4	False	6	False	clear	9.84	14.395	0.75	1

Every row is one hour of the years 2011 (year 0) and 2012 (year 1).

Looking at the raw data

Code

print(bikes["weather"].value_counts())

weather
clear         11413
misty          4544
rain           1419
heavy_rain        3
Name: count, dtype: int64

Heavy rain occurs in three hours only. No method can learn anything about it from three examples, so we drop these rows.

bikes = bikes[bikes["weather"] != "heavy_rain"]
bikes["weather"] = bikes["weather"].cat.remove_unused_categories()

drawing the mean rentals per hour and per month

fig, axes = plt.subplots(1, 2, figsize=(9, 3.2))
for wd, label in [("True", "working day"), ("False", "weekend or holiday")]:
    m = bikes[bikes["workingday"].astype(str) == wd].groupby("hour")["count"].mean()
    axes[0].plot(m.index, m.values, "o-", ms=3, label=label)
axes[0].set(xlabel="hour of day", ylabel="mean rentals")
axes[0].legend()
for year in [0, 1]:
    m = bikes[bikes["year"] == year].groupby("month")["count"].mean()
    axes[1].plot(m.index, m.values, "o-", ms=3, label=str(2011 + year))
axes[1].set(xlabel="month", ylabel="mean rentals")
axes[1].legend()
plt.tight_layout()
plt.show()

Figure 45.1: Left, the mean rentals per hour of the day, on working days and on other days. Right, the mean rentals per month, for the two years.

On working days there are two commuter peaks, at 8 and at 17 to 18 hours. On the other days there is one broad peak in the afternoon. The relation with the hour is neither linear nor the same on all days. There are more rentals in summer than in winter, and more in 2012 than in 2011.

Code

print(bikes[["temp", "feel_temp", "humidity", "windspeed"]].corr().round(2))

           temp  feel_temp  humidity  windspeed
temp       1.00       0.99     -0.07      -0.02
feel_temp  0.99       1.00     -0.05      -0.06
humidity  -0.07      -0.05      1.00      -0.29
windspeed -0.02      -0.06     -0.29       1.00

The temperature and the felt temperature carry nearly the same information.

Encoding the hour

Hour 23 is next to hour 0, but far from it on the number line. We can give the model the hour as a number, as \(\sin(2\pi h/24)\), \(\cos(2\pi h/24)\) and the same at twice the frequency, as 23 one-hot columns, or as 23 one-hot columns separately for working days and for the other days.

A Poisson regression on the temperature, the year, the working day indicator and one of four encodings of the hour, fitted to all rows. The points are the observed mean rentals per hour, the line the mean of the fitted rate over the same hours. Only with one-hot columns per day type can the fit have two peaks on working days and one on the other days.

A ladder of models

We start with models that are hard to get wrong and add flexibility step by step. All are scored by the same five-fold cross-validation, with the Poisson deviance, the loss that a Poisson regression minimizes.

from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import PoissonRegressor
from sklearn.model_selection import KFold, cross_val_score

y = bikes["count"].to_numpy()
X = bikes.drop(columns=["count", "feel_temp"])
X["hour_daytype"] = X["hour"].astype(str) + "_" + X["workingday"].astype(str)
numeric = ["temp", "humidity", "windspeed"]


def glm(onehot, raw=[]):
    pre = ColumnTransformer([
        ("onehot", OneHotEncoder(drop="first"), onehot),
        ("scale", StandardScaler(), numeric + raw)])
    return Pipeline([("pre", pre), ("glm", PoissonRegressor(alpha=1e-6, max_iter=5000))])


cv = KFold(5)
ladder = {
    "hour as a number": glm(["holiday", "workingday", "weather", "season"],
                            raw=["hour", "weekday", "month", "year"]),
    "one-hot hour": glm(["hour", "weekday", "month", "year", "holiday",
                         "workingday", "weather", "season"]),
    "one-hot hour per day type": glm(["hour_daytype", "weekday", "month", "year",
                                      "holiday", "weather", "season"]),
}
for name, model in ladder.items():
    dev = -cross_val_score(model, X, y, cv=cv, scoring="neg_mean_poisson_deviance")
    print(f"{name:27s} {dev.mean():6.1f}")

1: A new categorical column with 48 values, one per hour and day type. Its one-hot code lets a linear model give each hour its own effect on working days and on the other days.
2: The scaler is part of the pipeline, so in every fold of the cross-validation it is fitted on the training part only.
3: Without shuffling, every fold is a contiguous block of hours. Neighbouring hours are similar, so shuffled folds would put the hour before and the hour after every test hour into the training set, and the estimate would be too optimistic.

hour as a number              98.1
one-hot hour                  38.1
one-hot hour per day type     18.7

Code

from sklearn.metrics import mean_poisson_deviance
print(f"constant rate: {mean_poisson_deviance(y, np.full(len(y), y.mean())):.1f}")

constant rate: 166.4

The representation matters most. With the hour as one-hot columns the deviance drops from about 98 to about 38, and with an interaction between the hour and the day type to about 19, with the same method each time. The interaction we had to find ourselves, by looking at the data.