# RAMP on predicting cyclist traffic in Paris¶

Authors: Roman Yurchak (Symerio); also partially inspired by the air_passengers starting kit.

## Introduction¶

The dataset was collected with cyclist counters installed by Paris city council in multiple locations. It contains hourly information about cyclist traffic, as well as the following features,

- counter name
- counter site name
- date
- counter installation date
- latitude and longitude

Available features are quite scarce. However, we can also use any external data that can help us to predict the target variable.

In [3]:

    from pathlib import Path

    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import sklearn

# Loading the data with pandas¶

First, download the data files,

- train.parquet
- test.parquet

and put them to into the data folder.

Data is stored in Parquet format, an efficient columnar data format. We can load the train set with pandas,

In [4]:

    data = pd.read_parquet(Path("data") / "train.parquet")

In [5]:

    data.head()

Out[5]:

| counter_id | counter_name | site_id | site_name | bike_count | date | counter_installation_date | counter_technical_id | latitude | longitude | log_bike_count
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
48321 | 100007049-102007049 | 28 boulevard Diderot E-O | 100007049 | 28 boulevard Diderot | 0.0 | 2020-09-01 02:00:00 | 2013-01-18 | Y2H15027244 | 48.846028 | 2.375429 | 0.000000
48324 | 100007049-102007049 | 28 boulevard Diderot E-O | 100007049 | 28 boulevard Diderot | 1.0 | 2020-09-01 03:00:00 | 2013-01-18 | Y2H15027244 | 48.846028 | 2.375429 | 0.693147
48327 | 100007049-102007049 | 28 boulevard Diderot E-O | 100007049 | 28 boulevard Diderot | 0.0 | 2020-09-01 04:00:00 | 2013-01-18 | Y2H15027244 | 48.846028 | 2.375429 | 0.000000
48330 | 100007049-102007049 | 28 boulevard Diderot E-O | 100007049 | 28 boulevard Diderot | 4.0 | 2020-09-01 15:00:00 | 2013-01-18 | Y2H15027244 | 48.846028 | 2.375429 | 1.609438
48333 | 100007049-102007049 | 28 boulevard Diderot E-O | 100007049 | 28 boulevard Diderot | 9.0 | 2020-09-01 18:00:00 | 2013-01-18 | Y2H15027244 | 48.846028 | 2.375429 | 2.302585

We can check general information about different columns,

In [6]:

    data.info()

    <class 'pandas.core.frame.DataFrame'>
    Int64Index: 455163 entries, 48321 to 928462
    Data columns (total 11 columns):
     #   Column                     Non-Null Count   Dtype
    ---  ------                     --------------   -----
     0   counter_id                 455163 non-null  category
     1   counter_name               455163 non-null  category
     2   site_id                    455163 non-null  int64
     3   site_name                  455163 non-null  category
     4   bike_count                 455163 non-null  float64
     5   date                       455163 non-null  datetime64[ns]
     6   counter_installation_date  455163 non-null  datetime64[ns]
     7   counter_technical_id       455163 non-null  category
     8   latitude                   455163 non-null  float64
     9   longitude                  455163 non-null  float64
     10  log_bike_count             455163 non-null  float64
    dtypes: category(4), datetime64[ns](2), float64(4), int64(1)
    memory usage: 29.5 MB

and in particular the number of unique entries in each column,

In [7]:

    data.nunique(axis=0)

Out[7]:

    counter_id                     56
    counter_name                   56
    site_id                        30
    site_name                      30
    bike_count                    977
    date                         8230
    counter_installation_date      22
    counter_technical_id           30
    latitude                       30
    longitude                      30
    log_bike_count                977
    dtype: int64

We have a 30 counting sites where sometimes multiple counters are installed per location. Let's look at the most frequented stations,

In [8]:

    data.groupby(["site_name", "counter_name"])["bike_count"].sum().sort_values(
        ascending=False
    ).head(10).to_frame()

Out[8]:

|  | bike_count
--- | --- | ---
site_name | counter_name |
Totem 73 boulevard de Sébastopol | Totem 73 boulevard de Sébastopol S-N | 1809231.0
Totem 64 Rue de Rivoli | Totem 64 Rue de Rivoli O-E | 1406900.0
Totem 73 boulevard de Sébastopol | Totem 73 boulevard de Sébastopol N-S | 1357868.0
67 boulevard Voltaire SE-NO | 67 boulevard Voltaire SE-NO | 1036575.0
Totem 64 Rue de Rivoli | Totem 64 Rue de Rivoli E-O | 914089.0
27 quai de la Tournelle | 27 quai de la Tournelle SE-NO | 888717.0
Quai d'Orsay | Quai d'Orsay E-O | 849724.0
Totem Cours la Reine | Totem Cours la Reine O-E | 806149.0
Face au 48 quai de la marne | Face au 48 quai de la marne SO-NE | 806071.0
Face au 48 quai de la marne NE-SO | 759194.0

# Visualizing the data¶

Let's visualize the data, starting from the spatial distribution of counters on the map

In [9]:

    import folium

    m = folium.Map(location=data[["latitude", "longitude"]].mean(axis=0), zoom_start=13)

    for _, row in (
        data[["counter_name", "latitude", "longitude"]]
        .drop_duplicates("counter_name")
        .iterrows()
    ):
        folium.Marker(
            row[["latitude", "longitude"]].values.tolist(), popup=row["counter_name"]
        ).add_to(m)

    m

Out[9]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Note that in this RAMP problem we consider only the 30 most frequented counting sites, to limit data size.

Next we will look into the temporal distribution of the most frequented bike counter. If we plot it directly we will not see much because there are half a million data points,

In [10]:

    mask = data["counter_name"] == "Totem 73 boulevard de Sébastopol S-N"

    data[mask].plot(x="date", y="bike_count")

Out[10]:

    <AxesSubplot: xlabel='date'>

Instead we aggregate the data, for instance, by week to have a clearer overall picture,

In [11]:

    mask = data["counter_name"] == "Totem 73 boulevard de Sébastopol S-N"

    data[mask].groupby(pd.Grouper(freq="1w", key="date"))[["bike_count"]].sum().plot()

Out[11]:

    <AxesSubplot: xlabel='date'>

While at the same time, we can zoom on a week in particular for a more short-term visualization,

In [12]:

    fig, ax = plt.subplots(figsize=(10, 4))

    mask = (
        (data["counter_name"] == "Totem 73 boulevard de Sébastopol S-N")
        & (data["date"] > pd.to_datetime("2021/03/01"))
        & (data["date"] < pd.to_datetime("2021/03/08"))
    )

    data[mask].plot(x="date", y="bike_count", ax=ax)

Out[12]:

    <AxesSubplot: xlabel='date'>

The hourly pattern has a clear variation between work days and weekends (7 and 8 March 2021).

If we look at the distribution of the target variable it skewed and non normal,

In [13]:

    import seaborn as sns

    ax = sns.histplot(data, x="bike_count", kde=True, bins=50)

Least square loss would not be appropriate to model it since it is designed for normal error distributions. One way to precede would be to transform the variable with a logarithmic transformation,

    data['log_bike_count'] = np.log(1 + data['bike_count'])

In [14]:

    ax = sns.histplot(data, x="log_bike_count", kde=True, bins=50)

which has a more pronounced central mode, but is still non symmetric. In the following, we use `log_bike_count` as the target variable as otherwise `bike_count` ranges over 3 orders of magnitude and least square loss would be dominated by the few large values.

## Feature extraction¶

To account for the temporal aspects of the data, we cannot input the `date` field directly into the model. Instead we extract the features on different time-scales from the `date` field,

In [15]:

    def _encode_dates(X):
        X = X.copy()  # modify a copy of X
        # Encode the date information from the DateOfDeparture columns
        X.loc[:, "year"] = X["date"].dt.year
        X.loc[:, "month"] = X["date"].dt.month
        X.loc[:, "day"] = X["date"].dt.day
        X.loc[:, "weekday"] = X["date"].dt.weekday
        X.loc[:, "hour"] = X["date"].dt.hour

        # Finally we can drop the original columns from the dataframe
        return X.drop(columns=["date"])

In [16]:

    data["date"].head()

Out[16]:

    48321   2020-09-01 02:00:00
    48324   2020-09-01 03:00:00
    48327   2020-09-01 04:00:00
    48330   2020-09-01 15:00:00
    48333   2020-09-01 18:00:00
    Name: date, dtype: datetime64[ns]

In [17]:

    _encode_dates(data[["date"]].head())

Out[17]:

| year | month | day | weekday | hour
--- | --- | --- | --- | --- | ---
48321 | 2020 | 9 | 1 | 1 | 2
48324 | 2020 | 9 | 1 | 1 | 3
48327 | 2020 | 9 | 1 | 1 | 4
48330 | 2020 | 9 | 1 | 1 | 15
48333 | 2020 | 9 | 1 | 1 | 18

To use this function with scikit-learn estimators we wrap it with FunctionTransformer,

In [18]:

    from sklearn.preprocessing import FunctionTransformer

    date_encoder = FunctionTransformer(_encode_dates, validate=False)
    date_encoder.fit_transform(data[["date"]]).head()

Out[18]:

| year | month | day | weekday | hour
--- | --- | --- | --- | --- | ---
48321 | 2020 | 9 | 1 | 1 | 2
48324 | 2020 | 9 | 1 | 1 | 3
48327 | 2020 | 9 | 1 | 1 | 4
48330 | 2020 | 9 | 1 | 1 | 15
48333 | 2020 | 9 | 1 | 1 | 18

Since it is unlikely that, for instance, that `hour` is linearly correlated with the target variable, we would need to additionally encode categorical features for linear models. This is classically done with OneHotEncoder, though other encoding strategies exist.

In [19]:

    from sklearn.preprocessing import OneHotEncoder

    enc = OneHotEncoder(sparse=False)

    enc.fit_transform(_encode_dates(data[["date"]])[["hour"]].head())

Out[19]:

    array([[1., 0., 0., 0., 0.],
           [0., 1., 0., 0., 0.],
           [0., 0., 1., 0., 0.],
           [0., 0., 0., 1., 0.],
           [0., 0., 0., 0., 1.]])

## Linear model¶

Let's now construct our first linear model with Ridge. We use a few helper functions defined in `problem.py` of the starting kit to load the public train and test data:

In [20]:

    import problem

    X_train, y_train = problem.get_train_data()
    X_test, y_test = problem.get_test_data()

In [21]:

    X_train.head(2)

Out[21]:

| counter_id | counter_name | site_id | site_name | date | counter_installation_date | counter_technical_id | latitude | longitude
--- | --- | --- | --- | --- | --- | --- | --- | --- | ---
400125 | 100049407-353255860 | 152 boulevard du Montparnasse E-O | 100049407 | 152 boulevard du Montparnasse | 2020-09-01 01:00:00 | 2018-12-07 | Y2H19070373 | 48.840801 | 2.333233
408305 | 100049407-353255859 | 152 boulevard du Montparnasse O-E | 100049407 | 152 boulevard du Montparnasse | 2020-09-01 01:00:00 | 2018-12-07 | Y2H19070373 | 48.840801 | 2.333233

and

In [22]:

    y_train

Out[22]:

    array([1.60943791, 1.38629436, 0.        , ..., 2.48490665, 1.60943791,
           1.38629436])

Where `y` contains the `log_bike_count` variable.

The test set is in the future as compared to the train set,

In [23]:

    print(
        f'Train: n_samples={X_train.shape[0]},  {X_train["date"].min()} to {X_train["date"].max()}'
    )
    print(
        f'Test: n_samples={X_test.shape[0]},  {X_test["date"].min()} to {X_test["date"].max()}'
    )

    Train: n_samples=455163,  2020-09-01 01:00:00 to 2021-08-09 23:00:00
    Test: n_samples=41608,  2021-08-10 01:00:00 to 2021-09-09 23:00:00

In [24]:

    _encode_dates(X_train[["date"]]).columns.tolist()

Out[24]:

    ['year', 'month', 'day', 'weekday', 'hour']

In [25]:

    from sklearn.compose import ColumnTransformer
    from sklearn.linear_model import Ridge
    from sklearn.pipeline import make_pipeline

    date_encoder = FunctionTransformer(_encode_dates)
    date_cols = _encode_dates(X_train[["date"]]).columns.tolist()

    categorical_encoder = OneHotEncoder(handle_unknown="ignore")
    categorical_cols = ["counter_name", "site_name"]

    preprocessor = ColumnTransformer(
        [
            ("date", OneHotEncoder(handle_unknown="ignore"), date_cols),
            ("cat", categorical_encoder, categorical_cols),
        ]
    )

    regressor = Ridge()

    pipe = make_pipeline(date_encoder, preprocessor, regressor)
    pipe.fit(X_train, y_train)

Out[25]:

    Pipeline(steps=[('functiontransformer',
                     FunctionTransformer(func=<function _encode_dates at 0x7f5e41c58550>)),
                    ('columntransformer',
                     ColumnTransformer(transformers=[('date',
                                                      OneHotEncoder(handle_unknown='ignore'),
                                                      ['year', 'month', 'day',
                                                       'weekday', 'hour']),
                                                     ('cat',
                                                      OneHotEncoder(handle_unknown='ignore'),
                                                      ['counter_name',
                                                       'site_name'])])),
                    ('ridge', Ridge())])

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
Pipeline

    Pipeline(steps=[('functiontransformer',
                     FunctionTransformer(func=<function _encode_dates at 0x7f5e41c58550>)),
                    ('columntransformer',
                     ColumnTransformer(transformers=[('date',
                                                      OneHotEncoder(handle_unknown='ignore'),
                                                      ['year', 'month', 'day',
                                                       'weekday', 'hour']),
                                                     ('cat',
                                                      OneHotEncoder(handle_unknown='ignore'),
                                                      ['counter_name',
                                                       'site_name'])])),
                    ('ridge', Ridge())])

FunctionTransformer

    FunctionTransformer(func=<function _encode_dates at 0x7f5e41c58550>)

columntransformer: ColumnTransformer

    ColumnTransformer(transformers=[('date', OneHotEncoder(handle_unknown='ignore'),
                                     ['year', 'month', 'day', 'weekday', 'hour']),
                                    ('cat', OneHotEncoder(handle_unknown='ignore'),
                                     ['counter_name', 'site_name'])])

date

    ['year', 'month', 'day', 'weekday', 'hour']

OneHotEncoder

    OneHotEncoder(handle_unknown='ignore')

cat

    ['counter_name', 'site_name']

OneHotEncoder

    OneHotEncoder(handle_unknown='ignore')

Ridge

    Ridge()

We then evaluate this model with the RMSE metric,

In [26]:

    from sklearn.metrics import mean_squared_error

    print(
        f"Train set, RMSE={mean_squared_error(y_train, pipe.predict(X_train), squared=False):.2f}"
    )
    print(
        f"Test set, RMSE={mean_squared_error(y_test, pipe.predict(X_test), squared=False):.2f}"
    )

    Train set, RMSE=0.80
    Test set, RMSE=0.73

The model doesn't have enough capacity to generalize on the train set, since we have lots of data with relatively few parameters. However it happened to work somewhat better on the test set. We can compare these results with the baseline predicting the mean value,

In [27]:

    print("Baseline mean prediction.")
    print(
        f"Train set, RMSE={mean_squared_error(y_train, np.full(y_train.shape, y_train.mean()), squared=False):.2f}"
    )
    print(
        f"Test set, RMSE={mean_squared_error(y_test, np.full(y_test.shape, y_test.mean()), squared=False):.2f}"
    )

    Baseline mean prediction.
    Train set, RMSE=1.68
    Test set, RMSE=1.44

which illustrates that we are performing better than the baseline.

Let's visualize the predictions for one of the stations,

In [28]:

    mask = (
        (X_test["counter_name"] == "Totem 73 boulevard de Sébastopol S-N")
        & (X_test["date"] > pd.to_datetime("2021/09/01"))
        & (X_test["date"] < pd.to_datetime("2021/09/08"))
    )

    df_viz = X_test.loc[mask].copy()
    df_viz["bike_count"] = np.exp(y_test[mask.values]) - 1
    df_viz["bike_count (predicted)"] = np.exp(pipe.predict(X_test[mask])) - 1

In [29]:

    fig, ax = plt.subplots(figsize=(12, 4))

    df_viz.plot(x="date", y="bike_count", ax=ax)
    df_viz.plot(x="date", y="bike_count (predicted)", ax=ax, ls="--")
    ax.set_title("Predictions with Ridge")
    ax.set_ylabel("bike_count")

Out[29]:

    Text(0, 0.5, 'bike_count')

So we start to see the daily trend, and some of the week day differences are accounted for, however we still miss the details and the spikes in the evening are under-estimated.

A useful way to visualize the error is to plot `y_pred` as a function of `y_true`,

In [30]:

    fig, ax = plt.subplots()

    df_viz = pd.DataFrame({"y_true": y_test, "y_pred": pipe.predict(X_test)}).sample(
        10000, random_state=0
    )

    df_viz.plot.scatter(x="y_true", y="y_pred", s=8, alpha=0.1, ax=ax)

Out[30]:

    <AxesSubplot: xlabel='y_true', ylabel='y_pred'>

It is recommended to use cross-validation for hyper-parameter tuning with GridSearchCV or more reliable model evaluation with cross_val_score. In this case, because we want the test data to always be in the future as compared to the train data, we can use TimeSeriesSplit,

The disadvantage, is that we can either have the training set size be different for each fold which is not ideal for hyper-parameter tuning (current figure), or have constant sized small training set which is also not ideal given the data periodicity. This explains that generally we will have worse cross-validation scores than test scores,

In [31]:

    from sklearn.model_selection import TimeSeriesSplit, cross_val_score

    cv = TimeSeriesSplit(n_splits=6)

    # When using a scorer in scikit-learn it always needs to be better when smaller, hence the minus sign.
    scores = cross_val_score(
        pipe, X_train, y_train, cv=cv, scoring="neg_root_mean_squared_error"
    )
    print("RMSE: ", scores)
    print(f"RMSE (all folds): {-scores.mean():.3} ± {(-scores).std():.3}")

    RMSE:  [-0.96320864 -0.87311411 -0.85383867 -0.87000261 -1.0607788  -0.97297134]
    RMSE (all folds): 0.932 ± 0.0738

In [ ]:

In [ ]: