\section{Introduction}

bike-sharing systems (BSSs) have become increasingly popular in urban areas and have successfully  complemented public transportation systems in dense metropolitan cities. 
In addition to its utility to bike users, the installation of BSSs has been found
to reduce the usage of automobiles \citep{fishman2014bike} and thus traffic congestion and possibly green house emissions \citep{hamilton2018bicycle}.
To fully realize these benefits, efficient allocation of docking stations and bike docks is essential,
which in turn requires understanding the user behaviour and characteristics of the system
based on the abundant data collected on the BSS, in addition to
other urban and environmental factors that are known to influence bike usage.

In line with the increasing popularity of BSSs,
there exists a vast literature on the analysis of bike usage patterns (for and overview, see \cite{shaheen2010bikesharing};
and \cite{fishman2016bikeshare})
Below we provide a brief summary of the literature on quantitative or statistical analysis of BSS usage data which is categorized into two, following \cite{etienne2014model}.
The first branch addresses the problem of clustering stations in a system based on usage patterns according to some measure of similarity \citep{froehlich2009sensing, vogel2011understanding, etienne2014model}.
Regarding the BSS as a network, community detection algorithms have also been adopted for this purpose \citep{austwick2013structure, borgnat2013dynamical, zhou2015understanding}.
\cite{gervini2019exploring} model the demand for bikes as a multivariate temporal point process
and cluster stations based on functional canonical correlations
of log-intensity functions.

The second line of research concerns the problem of 
predicting the station occupancy or the state of the system at a given time.
\cite{faghih2016incorporating} model incoming and outgoing traffic at multiple stations as a panel with variables accounting for spatial and temporal autoregressive structures.
\cite{liu2016rebalancing} model inter-station bike transitioning
for improving the effectiveness of rebalancing operations
by predicting the station drop-off demand.
\cite{torti2021modelling} adopt functional linear regression
to model the directed flow between pairs of administrative divisions that aggregate multiple stations.

In all above, it is well-documented that BSS data show temporal and spatial patterns. 
To address such patterns, some previous works pre-process the datasets by, for example, aggregating stations into administrative regions \citep{torti2021modelling}, partitioning the data using subject-specific knowledge \citep{faghih2016incorporating}, or separately analyzing the data collected on weekdays and at weekends \citep{liu2016rebalancing}.
In complex urban environments, however, it may be difficult to find a single clustering of the data that comprehensively accounts for the patterns in the usage of the BSS, since there exist multiple approaches to produce geographical or temporal divisions according to socioeconomic characteristics, land zones, traffic infrastructure or population composition.
Then, two data points that can belong to the same division for one aspect (for example, according to administrative division) may belong to different divisions in another aspect (for example, business district versus tourist area). 
Besides, collecting in-depth information about the multifaceted nature of a large metropolitan city is typically costly or even impossible.

In this paper, we analyze the hourly bike rental data collected from a BSS in Seoul, South Korea, by adopting a penalized regression modelling approach. The proposed method utilizes the whole dataset without partitioning that still addresses the complex features inherent in the temporal and spatial properties of the data. 
We model the BSS at the granularity of individual stations by including station-specific parameters as well as trends and variables related to air quality and precipitation.
Such a model enjoys considerable flexibility and captures time-dependent usage patterns at individual stations, but it potentially suffers from the risk of overfitting as the number of parameters ($> 4 \times 10^4$) increases linearly with the number of stations. 

To address these issues, we propose a multilayer network fused Lasso penalty which extends the fused Lasso penalty \citep{tibshirani2005}. 
In the absence of a natural ordering among the stations, the proposed penalty imposes the penalization using a multilayer network and promotes fusion of the parameters linked by edges in the network.
In doing so, we view the BSS as a multilayer network with the stations serving as its nodes.
In the network, the layers correspond to different hours of a day and within each layer, a pair of stations are connected based on their geographical distance as a proxy for shared nodal features.
The cross-layer edges come from that the usage patterns tend to vary over the course of a day. 

The model fitted from the penalized regression method adaptively captures spatial and temporal homogeneity in bike usage, without (arbitrarily) partitioning the data which potentially leads to information loss. 
The degree of homogeneity is determined by the data-driven choice of the penalty parameter which does not involve the researcher's subjective decision. 
Our data analysis shows the superior predictive performance of the proposed multilayer fused Lasso penalty over alternative penalization approaches. 
Also, we propose a new network-based model complexity measure which reveals that while similarities exist, the stations exhibit fair amount of heterogeneity.
This conclusion supports that partitioning the stations into a handful of clusters may be inappropriate for such large-scale urban transportation systems.
The proposed idea of treating the underlying structure as a multilayer network, can be useful in problems beyond that of BSS modelling whenever the data exhibits network-like features.

The remainder of this paper is organised as follows.
In Section~\ref{sec:exploratory}, we describe the notable properties of BSSs such as station-specific temporal patterns, spatial homogeneity, and covariate effects,
through an exploratory analysis of 
bike rental data collected in Seoul, South Korea.
Section~\ref{sec:method} introduces our proposed Poisson regression model and
the accompanying penalization strategy that uses of fusion penalties, 
and provides a network-based interpretation of the latter.
Section~\ref{sec:analysis} demonstrates the effectiveness of our proposed approach 
on the examined BSS dataset.
We conclude the paper in Section~\ref{sec:conc}, and describe the algorithms for handling large-scale datasets in Appendix.

\section{Exploratory analysis of BSS data}
\label{sec:exploratory}

In this section, we describe data collected from ``Ddareungi'', a public BSS in Seoul, South Korea.
We use the hourly rental records collected from individuals between April 1, 2019 and May 30, 2019. \footnote{The dataset is available at \url{https://data.seoul.go.kr/}.}
By selecting the temperate months of April and May, we avoid dealing with seasonality or possibly abnormal observations due to extreme weather conditions.
We also exclude three public holidays falling in this period from our analysis.
Thus, our final dataset comprises observations from $T = 57$ days and $S = 1,505$ stations.

Since the launch of the BSS in 2010, the numbers of subscribed users, stations, and available bikes have steadily increased, as has the number of trips.
In addition to this overall trend, the bike usage data exhibits substantial temporal variations across localities, as well as dependence on precipitation and air quality, which is in line with the observations of previous BSSs literature.

\begin{figure}[htbp]
\centering
\subfloat{\includegraphics[width = .425\textwidth]{ST-767.png}}
\subfloat{\includegraphics[width = .425\textwidth]{ST-816.png}}\\
\subfloat{\includegraphics[width = .425\textwidth]{ST-107.png}}
\subfloat{\includegraphics[width = .425\textwidth]{ST-108.png}}
\caption{Hourly bike rentals from four selected stations over eight weeks, adjusted by the capacity (number of docking stations) of individual stations. Here, stations ST-767 and ST-816 are adjacent to one another and so are ST-107 and ST-108. We note that Week~5 contains three public holidays and show some anomalies.} 
\label{fig:temporal}
\end{figure}

Bike rental patterns exhibit substantial heterogeneity across stations. 
In Figure~\ref{fig:temporal}, we plot the hourly bike rental frequencies from selected pairs of stations adjacent to one another over the examined two-month period.
There are clear station-specific patterns over the course of a day and a week, and nearby stations tend to display similar usage patterns. 
Specifically, as stations ST-767 and ST-816 are located in a commercial district with government agencies and large firms, a spike in bike rental frequency is observed during weekday commuting-time.
On the other hand, stations ST-107 and ST-108 are located along a major riverside park in Seoul, and hence, are typically used for leisure activities, as evidenced by the large number of rentals concentrated on weekends. 
Additionally, although daily peaks and troughs in usage can be observed, these temporal patterns do not undergo abrupt changes in the sense that the number of bikes rented out between 9am and 10am is reasonably close to that between 10am and 11am.

\begin{figure}[htbp]
\centering
\subfloat{\resizebox*{4.5cm}{!}{\includegraphics{weekday_mean.png}}}
\subfloat{\resizebox*{4.5cm}{!}{\includegraphics{weekend_mean.png}}
\subfloat{\resizebox*{4.5cm}{!}{\includegraphics{weekend_weekday_diff_mean.png}}}
\caption{Average daily bike rentals at individual stations on weekdays (left), and weekends (middle) and their station-wise differences (right). The differences are taken in log-scale for better visualisation.} 
\label{fig:weekend}
\end{figure}

\begin{figure}[htbp]
\centering
\subfloat{%
\resizebox*{6.5cm}{!}{\includegraphics{rain_x_o_diff_mean.png}}}
\subfloat{%
\resizebox*{6.5cm}{!}{\includegraphics{dust_min_max_diff_mean.png}}}
\caption{Difference in average daily bike rentals (in log-scale) at individual stations due to precipitation (left) and air quality (right).}
\label{fig:raindust}
\end{figure}

Figure~\ref{fig:weekend} plots the daily rental frequencies averaged over weekdays and weekends as well as their differences (in log-scale) at individual stations.
The figure shows that bike usage behavior differs between weekdays and weekends, and that the degree of variation differs across stations.
Additionally, the last panel of Figure~\ref{fig:weekend} exhibits the presence of local clusters that share similar weekday/weekend variations.

Figure~\ref{fig:raindust} shows that bike usage depends on the weather condition and air quality.
We record each day as ``rainy'' if positive precipitation is recorded in any part of the city at any time of the day.
Comparing the average rental frequencies on days with and without precipitation (in log-scale, see the left panel of Figure~\ref{fig:raindust}), precipitation reduces bike usage, as expected.
Similar observations can be made with respect to the air quality: We plot the differences in average rental frequencies (in log-scale) on days when the air quality was ``good'' and ``very bad'', in the right panel of Figure~\ref{fig:raindust}. 
The qualitative categorization into ``good'', ``average'', ``bad,'' and ``very bad'' follows  classification system of the  Korean Ministry of Environment based on PM10 and PM2.5 dust concentrations.
The relative rental frequencies across the BSS display a similar pattern on ``good'' and ``very bad'' days, but the volume of bike rentals is lower on ``very bad'' days.

In summary, the BSS dataset exhibits the following characteristics:
\begin{enumerate}[itemindent=!, label = (C\arabic*)]
\item \label{c:trend} There exists an overall increasing trend.
\item \label{c:temp} Bike rentals show station-specific temporal patterns and the transition is reasonably smooth.
\item \label{c:net} These patterns are shared across stations that are geographically adjacent to one another. 
\item \label{c:cov} Bike rentals are influenced by the weather condition and the air quality.
\end{enumerate}


\section{Model and estimation}
\label{sec:method}

Motivated by the observations made in~\ref{c:trend}--\ref{c:cov} in Section~\ref{sec:exploratory},
we propose a Poisson regression model 
for hourly rental frequencies collected from the entire BSS (Section~\ref{sec:model}).
Section~\ref{sec:penalty} presents the accompanying estimation strategy
and introduces the multilayer fused Lasso penalty,
which is designed to capture the characteristics~\ref{c:temp}--\ref{c:net}.
In Section~\ref{sec:network}, we provide a network interpretation of the proposed penalization technique which aids in understanding and visualizing the penalty.

\subsection{Poisson regression model}
\label{sec:model}

Let $Y_i$ denote the $i$th observation representing the hourly rental frequency
at station $S(i) \in \mathcal S = \{1, \ldots, S\}$
and hour $H(i) \in \mathcal H = \{0, \ldots, 23\}$,
on day $D(i) \in \mathcal D = \{\text{Mo, Tu, \ldots, Su}\}$
for $i = 1, \ldots, n$, with $n = 2,058,840$ denoting the number of observations in the bike usage dataset.
We denote
the time index of the $i$th observation by $t(i) \in \mathcal T = \{0, \ldots, T - 1\}$.
Since we consider the period of two months, as described in Section~\ref{sec:exploratory} 
(excluding three public holidays), we have $T = 57$.
Additionally, $z^{\text{rain}}_i \in \{0, 1\}$ and $\mathbf z^{\text{air}}_i \in \{0, 1\}^4$ 
denote the variables representing
the precipitation and air quality statuses associated with the $i$th observation, respectively.
Consistent with the literature on BSSs and with the observations made in Section~\ref{sec:exploratory},
bike usage depends on the aforementioned variables
and we collect the variables relevant for the $i$th observation
in $\mathbf x_i = (S(i), t(i), D(i), H(i), z^{\text{rain}}_{i}, (\mathbf z^{\text{air}}_{i})^\top)^\top$ (for an overview, see Table~\ref{table:variables}).

\begin{table}[ht!!]
\caption{Variables in the BSS dataset.}
\label{table:variables}
\centering
{
\begin{tabular}{llll}
\toprule
Variable & Range & Description \\
\cmidrule(lr){1-1} \cmidrule(lr){2-2} \cmidrule(lr){3-3}
$S(i)$ & $\mathcal S = \{1, \ldots, S\}$ & Station node for the $i$th observation (categorical) \\
$H(i)$ & $\mathcal H = \{0, \ldots, 23\}$ & Hour of a day for the $i$th observation (categorical)\\
$D(i)$ & $\mathcal D = \{\text{Mo, \ldots, Su}\}$ & Day of a week for the $i$th observation (categorical)\\
$z^{\text{rain}}_i$ & $\{0, 1\}$ & Whether it rained or not for the $i$th observation (categorical)\\
$\mathbf z^{\text{air}}_i$ & $\{0, 1\}^4$ & Whether air quality was very bad, bad, average \\
& & or good for the $i$th observation (categorical) \\
$t(i)$ & $\{0, \ldots, T - 1\}$ & Time point for the $i$th observation (numeric)\\
\bottomrule
\end{tabular}}
\end{table}

Each station in the BSS of Seoul, South Korea, has a fixed number of docks, but unlike other BSSs, this does not determine the capacity of a station, as it is possible to leave bikes even if the docks are fully occupied through chaining them to existing bikes.
Therefore, we adopt a Poisson distribution for modelling $Y_i$, the hourly count of the bikes rented out at station $S(i)$, and model the relationship between $Y_i$ and $\mathbf x_i$ 
via Poisson regression \citep[Chapter~6]{mccullagh1989} as
\begin{align*}
Y_i \vert \mathbf x_i & \sim_{\text{iid}} \text{Poisson}(\mu_i) \quad \text{with} \quad
\mu_i \equiv \mu(\mathbf x_i) = \E(Y_i \vert \mathbf x_i).
\end{align*}
One way to model the link between $\mu_i$ and $\mathbf x_i$ is to use the following linear model:
\begin{align}
\log\l(\frac{\mu_i}{C_{S(i)}}\r) = \theta_{S(i)} + \alpha\, t(i) +
\beta^{\text{rain}}\, z^{\text{rain}}_{i} + \langle \bm\beta^{\text{air}}, \mathbf z^{\text{air}}_{i} \rangle +
\theta^{\text{hod}}_{H(i)} + \theta^{\text{dow}}_{D(i)}
\label{eq:model:no}
\end{align}
where $\bm\beta^{\text{air}} = (\beta^{\text{air}}_j, \, 0 \le j \le 3)^\top$.
Here, $\theta^{\text{hod}}_h$ and $\theta^{\text{dow}}_d$, respectively, contain ``hour of a day'' and ``day of a week'' effects common to all stations.
An offset term $C_s$ relates to the capacity of station~$s$, that is, the number of docks,
so that~\eqref{eq:model:no} can be interpreted as modelling
the expected rental frequency per hour per station capacity.
This model, referred to as the {\it no-interaction} model,
does not permit the temporal effects to be station-specific
and thus is too simple to address~\ref{c:temp}.

Allowing for interactions between the station and temporal effects,
we consider the following {\it full-interaction} model:
\begin{align}
\log\l(\frac{\mu_i}{C_{S(i)}}\r) =& \ \theta_{S(i)} + \alpha\, t(i) 
+ \beta^{\text{rain}}\, z^{\text{rain}}_{i}
+ \langle \bm\beta^{\text{air}}, \mathbf z^{\text{air}}_{i} \rangle 
\nonumber \\
& + \theta^{\text{hod}}_{H(i)} + \theta^{\text{dow}}_{D(i)}
+ \theta^{\text{hod}}_{S(i), H(i)} + \theta^{\text{dow}}_{S(i), D(i)}.
\label{eq:model}
\end{align}
The station-hour and station-day interaction terms $\theta^{\text{hod}}_{s, h}$ and $\theta^{\text{dow}}_{s, d}$ 
allow each station to exhibit individual temporal patterns.
For model identifiability, we set the baseline parameters to zero;
in other words, $\beta^{\text{air}}_0 = \theta^{\text{hod}}_0 = \theta^{\text{dow}}_{\text{Mo}} = \theta^{\text{hod}}_{s, 0} = \theta^{\text{dow}}_{s, \text{Mo}} = \theta^{\text{hod}}_{1, h} = \theta^{\text{dow}}_{1, d} = 0$.

Mmodel~(\ref{eq:model}) accounts for~\ref{c:trend}, \ref{c:cov}, and to a certain extent, \ref{c:temp}
by including the parameters $\alpha$ capturing the overall trend,
$\beta^{\text{rain}}$ and $\bm\beta^{\text{air}}$ capturing the effects of precipitation and air quality,
and $\theta^{\text{hod}}_{s, h}$ and $\theta^{\text{dow}}_{s, d}$ addressing station-specific temporal patterns.
In doing so, we take a different approach from those taken in previous studies in which,
after (arbitrarily) partitioning the dataset according to temporal or spatial variables, or both,
individual partitions are separately modelled (e.g.,\ \cite{austwick2013structure}; \cite{liu2016rebalancing}). 
Instead, by including the interaction terms, we use the full dataset and avoid any information loss from data partitioning.

While the full-interaction model enjoys considerably more flexibility than the no-interaction model in~\eqref{eq:model:no},
it suffers the risk of overfitting the data with the number of parameters to be estimated amounting to
$p = 34 + S + (S-1)\cdot(23 + 6) = 45,155$.
In addition, the characteristic identified in~\ref{c:net},
that is,  stations geographically close to one another tending to exhibit similar usage patterns,
is not adequately accounted for by the model fitted without any constraint.
Such an approach does not benefit from the temporal ordering inherent in the parameters $\theta^{\text{hod}}_{s, h}$, and thus does not fully account for~\ref{c:temp}.
In the next section, we propose a penalized maximum likelihood estimation (MLE) methodology for model in~(\ref{eq:model}) with a multilayer network fused Lasso which explicitly sets out to address these issues.

\subsection{Penalized MLE via multilayer network fused Lasso}
\label{sec:penalty}

We bridge the two models~(\ref{eq:model:no}) and~(\ref{eq:model}) at extreme ends, by adopting a penalized MLE strategy with a Lasso penalty mposed on the interaction parameters, and the fusion penalties designed to capture spatial and temporal homogeneity in bike usage patterns observed in the data, that is, \ref{c:temp} and~\ref{c:net}.

First proposed by \cite{tibshirani2005}, the fused Lasso augments 
the $\ell_1$-penalized least squares estimation method--or Lasso \citep{tibshirani1996regression}--
with a penalty that takes advantage of a meaningful ordering of the variables when such is available.
This fusion penalty is imposed on the $\ell_1$-norm of successive differences in the parameters corresponding to the ordered variables and encourages local constancy therein.
Its use has been extended beyond regression problems
such as trend estimation \citep{tibshirani2014adaptive},
change point detection \citep{harchaoui2010multiple} and graphical modelling \citep{danaher2014joint} among others.
 
Under full-interaction model in~(\ref{eq:model}), 
we partition the parameters into
$\Theta = \{\theta_s, \, s \in \mathcal S\}$, 
$\Theta_H = \{\theta^{\text{hod}}_{h}, \, \theta^{\text{hod}}_{s, h}, \, h \in \mathcal H \setminus \{0\}, \, s \in \mathcal S \setminus \{1\} \}$,
$\Theta_D = \{\theta^{\text{dow}}_{d}, \, \theta^{\text{dow}}_{s, d}, \, d \in \mathcal D \setminus \{\text{Mo}\}, \,  s \in \mathcal S \setminus \{1\} \}$ and
$\Delta = \{\alpha, \beta^{\text{rain}}, \beta^{\text{air}}_1, \beta^{\text{air}}_2, \beta^{\text{air}}_3\}$.
We adopt the fusion penalty to pool the information 
\begin{enumerate}[itemindent=!, label = (\roman*)]
\item \label{goal:one} from adjacent stations for the estimation of $\Theta$, $\Theta_H$ and $\Theta_D$, and
\item \label{goal:two} over the course of a day for the estimation of $\Theta_H$.
\end{enumerate}
We also impose a standard lasso penalty to encourage sparsity. Subsequently, our objective is to minimize the following penalized negative log-likelihood
\begin{align}
- \ell(\Theta, \Theta_H, \Theta_D, \Delta)
&+ \lambda 
\sum_{s \in \mathcal S} \l( \sum_{h \in \mathcal H} \vert \theta^{\text{hod}}_{s, h} \vert
+ \sum_{d \in \mathcal D} \vert \theta^{\text{dow}}_{s, d} \vert \r)
+ \lambda_N \cdot p_N(\Theta, \Theta_H, \Theta_D) 
+ \lambda_H \cdot p_H(\Theta_H)
\label{eq:org_obj_admm}
\end{align}
with respect to $\Theta$, $\Theta_H$, $\Theta_D$ and $\Delta$.
Here, $\ell(\Theta, \Theta_H, \Theta_D, \Delta)$ denotes 
the log-likelihood 
\begin{align*}
\ell(\Theta, \Theta_H, \Theta_D, \Delta)
= -\sum_{i = 1}^n \mu_i + \sum_{i = 1}^n y_i \log(\mu_i) + \text{constant},
\end{align*}
and $\mu_i = \mu_i(\Theta, \Theta_H, \Theta_D, \Delta)$  
is a function of the parameter vectors under~\eqref{eq:model}.
The parameter $\lambda$ controls the degree of the Lasso penalization 
and $\lambda_N$ and $\lambda_H$ control that of the fusion penalization;
below we discuss the construction of
the functions $p_N$ and $p_H$
to achieve the above goals~\ref{goal:one} and~\ref{goal:two}.

To capture the similarities between geographically adjacent stations,
we define a set of neighboring stations for each station $s$
as $\mathcal N_r(s) = \{s^\prime \in \mathcal S \setminus \{s\}: \, d(s, s^\prime) < r\}$ and its cardinality by $\vert \mathcal N_r(s) \vert$,
where $d(s, s^\prime)$ denotes the distance between the two stations $s$ and $s^\prime$
measured as the crow flies, and $r$ denotes a pre-specified distance.
We define 
\begin{align}
& p_N(\Theta, \Theta_H, \Theta_D) = \sum_{s \in \mathcal S} \sqrt{ \vert \mathcal N_r(s) \vert \, \sum_{s^\prime \in \mathcal N_r(s)} \l[ \sum_{h \in \mathcal H} (\phi^{\text{hod}}_{s, h} - \phi^{\text{hod}}_{s', h})^2 +
\sum_{d \in \mathcal D} (\phi^{\text{dow}}_{s, d} - \phi^{\text{dow}}_{s^\prime, d})^2 \r] }
\label{def:admm_net_penalty} 
\end{align}
with $\phi^{\text{hod}}_{s, h} = \theta_s + \theta^{\text{hod}}_{h} + \theta^{\text{hod}}_{s, h}$ and $\phi^{\text{dow}}_{s, d} = \theta_s + \theta^{\text{dow}}_{d} + \theta^{\text{dow}}_{s, d}$
which fully encode the station-specific hourly and daily patterns under~(\ref{eq:model}).
The penalty $p_N$ forces the pairs of parameters $(\phi^{\text{hod}}_{s, h}, \phi^{\text{hod}}_{s', h})$ and $(\phi^{\text{dow}}_{s, d}, \phi^{\text{dow}}_{s^\prime, d})$ to fuse neighboring stations~$s$ and $s'$, 
which encourages them to exhibit similar hourly and daily patterns and thus addresses the behavior noted in~\ref{c:net}.

Through $p_N$, the estimates of $\phi^{\text{hod}}_{s, h}$ and $\phi^{\text{hod}}_{s', h}$ may fuse even when $s' \notin \mathcal N_s(r)$
as long as they are connected via intermediate station nodes,
and the same applies to the estimates of $\phi^{\text{dow}}_{s,d}$ and $\phi^{\text{dow}}_{s',d}$.
This may be interpreted that by adopting $p_N$, we pool information across the BSS to estimate the interaction parameters (for a further discussion, see Section~\ref{sec:network}).
The penalty in~\eqref{def:admm_net_penalty} is related to the spatial fusion penalties \citep{sun2016fused, li2019spatial, sass2021flexible} adopted in the spatial regression literature, but is clearly distinguished as we adopt the group Lasso \citep{yuan2006model} in its formulation.
Every parameter inside the square root is considered to belong to the same group,
and the weighting applied with the size of $\mathcal N_r(s)$ 
follows the convention in the group Lasso literature.
Then, the parameters associated with the stations with a large number of neighbors receive more penalization.

\begin{rem}[Re-parameterization using $\phi^{\text{hod}}_{s, h}$ and $\phi^{\text{dow}}_{s, d}$.]
\label{rem:repar}
We choose to decompose the station-specific hourly and daily effects as $\phi^{\text{hod}}_{s, h} = \theta_s + \theta^{\text{hod}}_h + \theta^{\text{hod}}_{s, h}$ and $\phi^{\text{dow}}_{s, d} = \theta_s + \theta^{\text{dow}}_d + \theta^{\text{dow}}_{s, d}$. 
Then, the proposed method applies the Lasso penalty to the $\ell_1$-norm of $\theta^{\text{hod}}_{s, h}$ and $\theta^{\text{dow}}_{s, d}$ only, which gives it the interpretation of bridging between no-interaction and full-interaction models.
We choose not to impose the Lasso penalization on $\theta_s$ to avoid cancelling out station-specific usage patterns,
or $\theta^{\text{hod}}_h$ and $\theta^{\text{dow}}_d$, which represent overall hourly and daily patterns shared across all stations.

Recall that, for model identifiability, we set $\theta^{\text{hod}}_{s, 0} = 0$.
Thus, we can write
\begin{align*}
\sum_{h \in \mathcal H} (\phi^{\text{hod}}_{s, h} - \phi^{\text{hod}}_{s', h})^2
= (\theta_s - \theta_{s'})^2 + \sum_{h = 1}^{23} (\theta^{\text{hod}}_{s, h} - \theta^{\text{hod}}_{s', h})^2
\end{align*}
and similarly, since $\theta^{\text{dow}}_{s, \text{Mo}} = 0$, we have
\begin{align*}
\sum_{d \in \mathcal D} (\phi^{\text{dow}}_{s, d} - \phi^{\text{dow}}_{s', d})^2
= (\theta_s - \theta_{s'})^2 + \sum_{d \in \mathcal D \setminus \{\text{Mo}\}} (\theta^{\text{dow}}_{s, d} - \theta^{\text{dow}}_{s', d})^2.
\end{align*}
From these observations, we can re-write $p_N$ as
\begin{align*}
& p_N(\Theta, \Theta_H, \Theta_D) =
\sum_{s \in \mathcal S} \sqrt{ \vert \mathcal N_r(s) \vert \, \sum_{s^\prime \in \mathcal N_r(s)} \l[ 2(\theta_s - \theta_{s^\prime})^2 + \sum_{h = 1}^{23} (\theta^{\text{hod}}_{s, h} - \theta^{\text{hod}}_{s', h})^2 +
\sum_{d \in \mathcal D \setminus \{\text{Mo}\}} (\theta^{\text{dow}}_{s, d} - \theta^{\text{dow}}_{s^\prime, d})^2 \r] },
\end{align*}
which shows that $p_N$ implicitly encourages the parameters
$(\theta_s, \theta_{s'})$ to take similar values when stations~$s$ and $s'$ are neighbors.
\end{rem}



There is a natural temporal ordering inherent in $\Theta_H$ that gives rise to
the fusion penalty in its canonical form imposed on $\phi^{\text{hod}}_{s, h}$:
\begin{align}
p_H(\Theta_H) =  \sum_{s \in \mathcal S} 
\sum_{h = 0}^{23} \l\vert \phi^{\text{hod}}_{s, h} - \phi^{\text{hod}}_{s, h + 1} \r\vert 
=  \sum_{s \in \mathcal S} 
\sum_{h = 0}^{23} \l\vert (\theta^{\text{hod}}_{h} +\theta^{\text{hod}}_{s, h})- (\theta^{\text{hod}}_{h+1} + \theta^{\text{hod}}_{s, h + 1}) \r\vert 
\label{def:admm_hour_penalty} 
\end{align}
with $\theta^{\text{hod}}_{s, 24} = \theta^{\text{hod}}_{s, 0}$.
Imposing a penalty on $p_H(\Theta_H)$ encourages the consecutive station-specific hourly effects 
$\phi^{\text{hod}}_{s, h}$ and $\phi^{\text{hod}}_{s, h+1}$ at a given station $s$,
to become close to one another and suppress abrupt changes in usage.
This reflects the fact that rental counts at each station rarely undergo radical shifts
during the course of the day.

Jointly, $p_N$ and $p_H$ comprise the proposed multilayer network fusion penalty.
We devote Section~\ref{sec:network} to its interpretation with the description of the underlying multilayer network.
The impact of the fusion penalty is determined by the sizes of $\lambda_N$ and $\lambda_H$,
which we select via cross-validation, as described in Section~\ref{sec:tuning}.
The Alternating Direction Method of Multipliers (ADMM) algorithm \citep{boyd2011distributed} is employed to solve the convex optimization problem in~\eqref{eq:org_obj_admm}.
Efficient implementation of the algorithm requires careful re-parametrization of  model~(\ref{eq:model}), which makes use of the specific structure of the data, and we discuss this in detail in Appendix~\ref{sec:admm}. 

\subsection{multilayer network interpretation of the fusion penalty}
\label{sec:network}

We provide a network interpretation of the proposed penalized regression method. 
First, we introduce the following networks that underpin the penalty functions $p_N$ and $p_H$:
\begin{align}
\mathfrak{N}_{\text{single}}(r) =& (\mathcal S, \mathcal E_{\text{single}}(r))
\text{ \ with \ } 
\mathcal E_{\text{single}}(r) = \cup_{s \in \mathcal S} \{ (s, s'), \, s' \in \mathcal N_r(s)\},
\label{eq:net:single} \\
\mathfrak{N}_{\text{multi}}(r) =&
(\mathcal S \times \mathcal H, \, \mathcal E_{\text{multi}}(r)) \text{ \ with \ } \nonumber
\\
\mathcal E_{\text{multi}}(r) =&
\l[ \cup_{h \in \mathcal H} \cup_{s \in \mathcal S} \l\{ \l( (s, h), (s', h) \r), \, s' \in \mathcal N_r(s) \r\} \r] \bigcup \l[\cup_{s \in \mathcal S} \l\{ \l( (s, h), (s, h + 1) \r), \, h \in \mathcal H \r\}\r].
\label{eq:net:multi} 
\end{align}
(For an illustrative example of $\mathfrak{N}_{\text{single}}$ and $\mathfrak{N}_{\text{multi}}$, see Figure~\ref{fig:single_multi_network}).
While both $\mathfrak{N}_{\text{single}}$ and $\mathfrak{N}_{\text{multi}}$ as well as their edge sets
depend on the choice of $r$, we suppress this dependency for simplicity when it does not cause any confusion.

\begin{figure}[h]
\centering
\resizebox*{10cm}{!}{\includegraphics{singlenet.png}}
\caption{Illustration of a single-layer network $\mathfrak{N}_{\text{single}}$ (left) and a multilayer network $\mathfrak{N}_{\text{multi}}$ (right). Each dot represents a node (station) and a solid line represents an edge connecting the nodes within each layer and across adjacent layers.}
\label{fig:single_multi_network}
\end{figure}

Network $\mathfrak{N}_{\text{single}}$ is a single-layer, undirected network 
that is solely determined by the sets of neighbors $\mathcal N_r(s)$ of the stations.
In this network, a pair of ``day of a week'' parameters $(\phi^{\text{dow}}_{s, d}, \phi^{\text{dow}}_{s', d})$ for each given day $d$, are encouraged to take values close to one another by the penalty function $p_N$, provided that the stations belong to the same connected component of $\mathfrak{N}_{\text{single}}$.
On the other hand, $\mathfrak{N}_{\text{multi}}$ is a multilayer, undirected network;
we follow the notational convention of \cite{kivela2014multilayer}, in which the set $\mathcal H$ serves as a set of elementary layers for the hourly aspect,
and each edge connects a pair of node-layer tuples $(s, h)$ and $(s', h')$ for some $s, s' \in \mathcal S$ and $h, h' \in \mathcal H$.
Each pair of the ``hour of a day'' parameters 
$\phi^{\text{hod}}_{s, h}$ and $\phi^{\text{hod}}_{s', h'}$
is encouraged to fuse with one another by the penalty functions $p_N$ and $p_H$,
if the corresponding pair of nodes are connected in $\mathfrak{N}_{\text{multi}}$.

Next, we define networks whose edges 
are determined by the coefficient estimates from the penalized MLE.
Let
$\widehat\phi^{\text{hod}}_{s, h}$ and $\widehat\phi^{\text{dow}}_{s, d}$ denote the estimates of the parameters $\phi^{\text{hod}}_{s, h}$ and $\phi^{\text{dow}}_{s, d}$, respectively;
their values depend on the tuning parameters 
$(r, \lambda, \lambda_N, \lambda_H)$; however,
we omit this dependence for simplicity.
Then, the networks associated with these estimates are defined as
\begin{align}
\widehat{\mathfrak{N}}_{D, d} = \l(\mathcal S, \, \mathcal E_{D, d}\r) \text{ \ with \ } 
\mathcal E_{D, d} =& \l\{(s, s'), \, s \ne s': \, \widehat\phi^{\text{dow}}_{s, d} = \widehat\phi^{\text{dow}}_{s', d} \r\} \text{ \ for each \ } d \in \mathcal D, 
\nonumber \\
\widehat{\mathfrak{N}}_H = \l(\mathcal S \times H, \mathcal E_H\r)  \text{ \ with \ }
\mathcal E_H =&  \cup_{h \in \mathcal H} \cup_{s \in \mathcal S} 
\l\{((s, h), (s', h')), \, (s,h) \ne (s',h'): \, \widehat\phi^{\text{hod}}_{s, h} = \widehat\phi^{\text{hod}}_{s', h'} \r\} 
\label{eq:net:hod} 
\end{align}

As with $\mathfrak{N}_{\text{single}}$, the networks $\widehat{\mathfrak{N}}_{D, d}$ are single-layer networks
and an edge joins two station nodes $s$ and $s'$ when their node features (i.e. \ parameter estimates of  $\phi^{\text{dow}}_{s,d}$ and $\phi^{\text{dow}}_{s',d}$ for a given $d \in \mathcal D$) are identical, possibly due to the fusion penalty but not necessarily so.
The network $\widehat{\mathfrak{N}}_H$, as with $\mathfrak{N}_{\text{multi}}$, 
is a multilayer network with the hourly layer given by $\mathcal H$, 
and an edge is formed between a pair of nodes $(s,h)$ and $(s',h')$ when the estimates of $\phi^{\text{hod}}_{s, h}$ and $\phi^{\text{hod}}_{s', h'}$ agree at $(s,h) \ne (s',h')$. 
(For an illustrative example of $\widehat{\mathfrak{N}}_H$, see Figure~\ref{fig:network}).

\begin{figure}[h]
\centering
\resizebox*{13cm}{!}{\includegraphics{multinet.png}}
\caption{Illustration of multilayer networks. 
Left: Each layer of $\mathfrak{N}_{\text{multi}}$ (see~\eqref{eq:net:multi}) embeds the linkages between the stations determined by their geographical distances at a given hour $h$, and the station at layer $h$ is linked to itself at layers $h - 1$ and $h + 1$, 
which underpins how $\phi^{\text{hod}}_{s, h}$ are encouraged to be fused by $p_N$ and $p_H$.
Middle: $\widehat{\mathfrak{N}}_H$ (see~\eqref{eq:net:hod}) is determined by the values of the estimates of $\phi^{\text{hod}}_{s, h}$ with an edge indicating that the connected estimates share the identical values.
Right: A multilayer network formed with its edges obtained as an intersection of the edge sets of $\mathfrak{N}_{\text{multi}}$ and $\widehat{\mathfrak{N}}_H$,
which contains four connected components.} 
\label{fig:network}
\end{figure}

For two networks (either single- or multilayer) 
$\mathfrak{N}_i = (\mathcal V, \mathcal E_i), \, i = 1, 2$, sharing the same node set $\mathcal V$, 
denote by $\mathfrak{N}_1 \cap \mathfrak{N}_2 = (\mathcal V, \mathcal E_1 \cap \mathcal E_2)$ 
the network formed by taking the intersection of their edge sets.
Our proposed penalized regression method takes as an input the observable networks $\mathfrak{N}_{\text{single}}$ and $\mathfrak{N}_{\text{multi}}$, which inform spatial and temporal proximity.
Then, it outputs the networks capturing the homogeneity between the stations nodes, namely, $\mathfrak{N}_{\text{single}} \cap \widehat{\mathfrak{N}}_{D, d}$ (on a given day of a week $d$) and $\mathfrak{N}_{\text{multi}} \cap \widehat{\mathfrak{N}}_H$ (along the hourly layer),
which can provide insights into BSS management and urban planning.
Depending on the choice of penalty parameters, the output networks are not necessarily sparse; in fact, this is the case in our data analysis reported in Section~\ref{sec:analysis}.
This distinguishes our approach from the existing literature on clustering or partitioning the dataset using spatial or temporal variables prior to analysis.


\section{Data analysis}
\label{sec:analysis}

\subsection{Tuning parameter selection}
\label{sec:tuning}


For the selection of $r$ that determines the set of neighbors $\mathcal{N}_r(s)$ for each station $s$, we examine the distance between each station and its $10$ nearest stations,
that is, $d_{s, i}, \, i = 1, \ldots, 10$, for all $s \in \mathcal S$.
The median of the first quartile
of $d_{s, i}$ over $s \in \mathcal S$ is $762$ meters, 
from which we choose to consider $r \in \{750, 1500\}$.
The network $\mathfrak{N}_{\text{single}}(r)$ (see~\eqref{eq:net:single}) with 
$r = 750$ has $34$ connected components, whereas the network with $r = 1500$ has two connected components.
Later, we show that our penalized regression approach
is not sensitive to $r$ on the BSS dataset
through the adaptive selection of the penalty parameters discussed below.

The parameters $\lambda$, $\lambda_N$, and $\lambda_H$ 
control the overall complexity of the fitted model, and 
the latter two, in particular control the degree of spatial and temporal homogeneity
induced by the parameter estimates of $\Theta$, $\Theta_H$ and $\Theta_D$.
We propose the selection of the tuning parameters via cross validation (CV).
In the penalized regression literature, CV is typically performed by randomly partitioning the data into five or ten folds, 
which, in the case of the BSS dataset, ignores the temporal structure therein.
Instead, we make use of the fact that our dataset covers seven weekends, 
and adopt seven-fold CV, in which each fold is ensured to 
include a balanced number of all seven days of the week.

As a CV measure, we adopt the mean squared Pearson residuals (MSPR): 
\begin{align}
\label{eq:cv}
\text{CV}(r, \lambda, \lambda_N, \lambda_H) = 
\frac{1}{7} \sum_{j = 1}^7 \frac{1}{n_j} \sum_{i = 1}^{n_j} 
\frac{(Y^{(j)}_i - \widehat\mu^{(j)}_i(r, \lambda, \lambda_N, \lambda_H))^2}{\widehat\mu^{(j)}_i(r, \lambda, \lambda_N, \lambda_H)}
\end{align}
where for the $j$th fold,
$n_j$ denotes the total number of observations,
$Y_i^{(j)}$ the $i$th observation and
$\widehat\mu_i^{(j)}(r, \lambda, \lambda_N, \lambda_H)$ is the corresponding estimate of the mean from the model fitted to the remaining data with the given tuning parameters. 
We evaluate $\text{CV}(r, \lambda, \lambda_N, \lambda_H)$ on grids of equispaced values on a log-scale for $(\lambda, \lambda_N, \lambda_H)$
and select the combination that returns the smallest CV.


\subsection{Analysis of the BSS data}

\subsubsection{Effects of penalization}

To assess the effect of penalization, we compare our multilayer network fused Lasso-based penalized regression approach (referred to as ``fused Lasso'') with the methods that adopt either the fusion penalty or the Lasso penalty only, which are referred to as ``fusion-only'' and ``Lasso-only''.
We also consider models (\ref{eq:model:no}) and~(\ref{eq:model}) fitted without any penalization,
referred to as ``no-interaction'' and ``full-interaction'', respectively.

\paragraph{Predictive performance}

Due to small sample size and inherent non-exchangeability of the data, we adopt CV measures as an indicator of the predictive performance.
In Table~\ref{tab:mspr} and Figure~\ref{fig:mspr},
we report the fold-wise MSPRs involved in the seven-fold CV in~\eqref{eq:cv} as well as the overall CV error, all evaluated at the penalty parameters selected to minimize the overall CV error for the respective methods.

When comparing no-interaction and full-interaction methods,
the flexibility afforded by allowing for station-specific temporal effects proves useful in enhancing the predictive performance as the latter model consistently attains a considerably smaller MSPR.
We observe further improvement when appropriate penalization is applied to the interaction parameters.
In particular, adopting the proposed fused Lasso penalty returns the minimum CV error across the seven folds compared to the other methods regardless of the choice of $r$.
Between the fusion-only and the Lasso-only methods, the former outperforms the latter, indicating that capturing across-station homogeneity pays off by accounting for the stylized features of the BSS data that neighboring stations exhibit similar usage patterns, (see~\ref{c:net}). 
The MSPR attained with $r = 1500$ is consistently smaller, albeit by a small margin, than that attained with $r = 750$.
Recalling that the choice $r = 1500$ reduces the number of connected components in $\mathfrak{N}_{\text{single}}(r)$ compared to when $r = 750$ (from $54$ to $2$), this result is indicative of the benefit of pooling information across the wider spatial networks of BSSs.

\begin{table}[h!]
\caption{Mean squared Pearson residuals (MSPRs) from each fold used in the seven-fold CV and the overall CV error as measured in~\eqref{eq:cv}, with the penalty parameters selected to minimize the latter for the respective methods where relevant.}
\label{tab:mspr}
\centering
{\small \begin{tabular}{l ccccccc cc}
\toprule
& \multicolumn{7}{c}{Fold-wise CV} & \multicolumn{2}{c}{} \\
Method & 1 & 2 & 3 & 4 & 5 & 6 & 7 & Average  \\
\cmidrule(lr){1-1} \cmidrule(lr){2-8} \cmidrule(lr){9-10}
Fused Lasso (750m) & 1.558 & 1.423 & 1.424 & 1.497 & 1.449 & 1.831 & 1.683 & 1.552 \\
Fused Lasso (1500m) & 1.553 & 1.420 & 1.420 & 1.494 & 1.443 & 1.826 & 1.675 & 1.547 \\
\cmidrule(lr){1-1} \cmidrule(lr){2-8} \cmidrule(lr){9-10}
Fusion-only (750m) & 1.572 & 1.443 & 1.445 & 1.512 & 1.473 & 1.861 & 1.701 & 1.569 \\
Fusion-only (1500m) & 1.568 & 1.440 & 1.438 & 1.511 & 1.463 & 1.847 & 1.700 & 1.567 \\
Lasso-only & 1.609 & 1.479 & 1.473 & 1.537 & 1.493 & 1.902 & 1.728 & 1.603 \\
Full-interaction & 1.729 & 1.601 & 1.556 & 1.629 & 1.596 & 2.016 & 1.841 &  1.710 \\
\cmidrule(lr){1-1} \cmidrule(lr){2-8} \cmidrule(lr){9-10}
No-interaction & 1.991 & 1.857 & 1.791 & 1.918 & 1.838 & 2.264 & 2.118 & 1.969 \\
\bottomrule
\end{tabular}}
\end{table}

\begin{figure}[ht]
\centering
\includegraphics[width = 0.9\textwidth]{mspr.png}
\caption{Mean squared Pearson residuals (MSPRs) from each fold used in the seven-fold CV and the overall CV error as measured in~\eqref{eq:cv},
with the penalty parameters selected to minimize the CV error for the respective methods.}
\label{fig:mspr}
\end{figure}

\paragraph{Sparsity and model complexity}

We examine the reduction in complexity brought about by the Lasso penalty and the fusion penalties $p_N$ and $p_H$ in~\eqref{eq:org_obj_admm}.
The sparsity induced by the Lasso penalty is easily measured by the proportion of non-zero coefficient estimates (see the last column of Table~\ref{tab:mc}).
The model fitted with the fused Lasso penalty is approximately $13$\% (resp.\ $20$\%) less complex than the full-interaction model when $r = 750$ (resp.\ $r = 1500$), 
yet the former has much better predictive power than the latter as  Table~\ref{tab:mspr} shows.
At the same time, further sparsity induced by adopting the Lasso-only penalty does not improve the prediction performance. 
The fusion-only method does not set any parameter estimate to be exactly zero.

\begin{table}[h!]
\caption{Model complexity (MC) of the models fitted the penalty parameters selected to minimize the CV error for the respective methods (where relevant).
We report the overall MC out of all parameters ($\Delta \cup \Theta \cup \Theta_H  \cup \Theta_D$) as well as that out of the interaction  parameters representing hour-of-a-day ($\phi^{\text{hod}}_{s,h}$) and day-of-a-week ($\phi^{\text{dow}}_{s,d}$) effects. Additionally, the proportion of non-zero coefficient estimates is presented.}
\label{tab:mc}
\centering
{\small 
\begin{tabular}{l cccc}
\toprule
& \multicolumn{3}{c}{Parameter set} & Proportion of \\
Method & All & $\Theta_H \setminus \{\theta^{\text{hod}}_h, \, h \in \mathcal H\}$ & $\Theta_D \setminus \{\theta^{\text{dow}}_d, \, d \in \mathcal D\}$ & non-zeros \\
\cmidrule(lr){1-1} \cmidrule(lr){2-2} \cmidrule(lr){3-3} \cmidrule(lr){4-4} \cmidrule(lr){5-5}
Fused Lasso (750m) & 0.828 & 0.790 & 0.977 & 0.867\\
Fused Lasso (1500m) & 0.848 & 0.827 & 0.927 & 0.795 \\
\cmidrule(lr){1-1} \cmidrule(lr){2-2} \cmidrule(lr){3-3} \cmidrule(lr){4-4} \cmidrule(lr){5-5}
Fusion-only (750m) & 0.903 &  0.885 & 0.974 & 1.000\\
Fusion-only (1500m) & 0.853 & 0.836 & 0.916 & 1.000 \\
\cmidrule(lr){1-1} \cmidrule(lr){2-2} \cmidrule(lr){3-3} \cmidrule(lr){4-4} \cmidrule(lr){5-5}
Lasso-only &
-- & -- & -- & 0.632 \\
\bottomrule
\end{tabular}}
\end{table}

We use the networks introduced in Section~\ref{sec:network} to define a model complexity (MC) measure that evaluates the effect of the fusion penalties.
Simply put, the model represents the proportion of coefficient estimates that are not fused by the penalties $p_N$ and $p_H$,
out of the total number of parameters in model in~(\ref{eq:model}), as
\begin{align}
& \text{MC}(r, \lambda, \lambda_N, \lambda_H) 
=  \frac{1}{p}
\l( 34 + \mathcal C(\mathfrak{N}_{\text{multi}}(r) \cap \widehat{\mathfrak{N}}_H) + \sum_{d \in \mathcal D  \setminus \{\text{Mo} \}} \mathcal C(\mathfrak{N}_{\text{single}}(r) \cap \widehat{\mathfrak{N}}_{D, d}) \r)
\label{eq:mc}
\end{align}
where we suppress the dependence on $(r, \lambda, \lambda_N, \lambda_H)$ for notational simplicity.
Recall that $34$ represents the number of parameters that are not penalized, namely $\alpha$, $\beta^{\text{rain}}$, $\{\beta^{\text{air}}_j, \, 1 \le j \le 3\}$, $\{\theta^{\text{hod}}_h, \, h \in \mathcal H\}$ and $\{\theta^{\text{dow}}_d, \, d \in \mathcal D\}$.
Denoting the number of connected components in a network $\mathfrak{N}$ by $\mathcal C(\mathfrak{N})$,
the number of unique parameter estimates of $\phi^{\text{hod}}_{s,h}$, which are not fused by penalization,
is given by $\mathcal C(\mathfrak{N}_{\text{multi}}(r) \cap \widehat{\mathfrak{N}}_H)$.
Figure~\ref{fig:network} shots that the network $\mathfrak{N}_{\text{multi}} \cap \widehat{\mathfrak{N}}_H$ contains four connected components; consequently,
$\mathcal C(\mathfrak{N}_{\text{multi}} \cap \widehat{\mathfrak{N}}_H) = 4$. Similarly,
We can find the number of unique parameter estimates for $\phi^{\text{dow}}_{s,d}$. 
In~\eqref{eq:mc}, we exclude the intersection network $\mathfrak{N}_{\text{single}}(r) \cap \widehat{\mathfrak{N}}_{D, \text{Mo}}$ from the numerator since, due to model identifiability constraints, we have $\phi^{\text{dow}}_{s,\text{Mo}} = \phi^{\text{hod}}_{s, 0} = \theta_{s}$.
That is, the fusion among the station-specific intercept parameters $\theta_s$ has already been accounted for by $\mathfrak{N}_{\text{multi}}(r) \cap \widehat{\mathfrak{N}}_H$ at layer $h = 0$.
We refer to Appendix~\ref{sec:mc} for the efficient calculation of MC which requires some care owing to the multilayered nature of the network in~\eqref{eq:net:multi} underlying the penalties.
We analogously evaluate the complexity of the model fitted using the fusion-only method (for the full results, see Table~\ref{tab:mc}).

The MC ranges between $0$ and $1$, and when its value is closer to $0$, it implies that most stations exhibit homogeneous behavior with their neighboring stations.
The resultant intersection networks
$\mathfrak{N}_{\text{multi}}(r) \cap \widehat{\mathfrak{N}}_H$ and
$\mathfrak{N}_{\text{single}}(r) \cap \widehat{\mathfrak{N}}_{D, d}$
are not highly homogeneous, as evidenced by the MC being close to one, that is,  many stations exhibit individual behavior. 
In other words, station-specific parameters account for a large portion of the variation in bike usage, which supports modelling the data at the individual station level.
The high degree of heterogeneity across the BSS can be attributed to the fact that Seoul is a large city with high population density, so that each station is associated with multiple aspects of usages. 
This indicates that partitioning stations into a handful of clusters may ignore the complex nodal features that drive the usage of bikes at each station.
Instead, our proposed multilayer fused Lasso method imposes fusion penalties through the networks $\mathcal N_{\text{single}}(r)$ and $\mathcal N_{\text{multi}}(r)$, and learns the degree of homogeneity in bike usage patterns across the BSS in a data-driven manner.

Generally, more homogeneity is observed in the hour-of-a-day effects than in the day-of-a-week effects.
Also, since the choice of $r = 1500$ leads to a highly connected $\mathfrak{N}_{\text{single}}(r)$, it leads to a fitted model with a larger MC than the choice of $r = 750$.
Overall, the fused Lasso method tends to further reduce the MC compared to the fusion-only method.
We attribute this to the fact that the Lasso penalty enhances the effect of the fusion penalty by coercing estimated values to be zeros.

\subsubsection{Parameter estimates}

There exist methods for performing inference in high-dimensional generalized linear models \citep{belloni2016post, cai2021statistical}, but they do not easily apply to our setting because of the presence of a fusion penalty.
Instead, we examine the coefficient estimates obtained from different folds used in the cross validation step in~\eqref{eq:cv} (for the estimates of the parameters associated with the overall trend and the effects of the precipitation and the air quality, obtained from each fold as well as from the full data, see Table~\ref{tab:coef}).
We only report the results when $r = 1500$ in the main text because we obtain nearly identical estimates when $r = 750$ as shown in Table~\ref{tab:coef750} in Appendix \ref{sec:aar}).
Table~\ref{tab:coef} shows that, while the values of the estimates vary slightly from one fold to another, their signs and overall magnitude do not change.
This confirms the observations~\ref{c:trend} and~\ref{c:cov} made in Section~\ref{sec:exploratory}, that is, the variables have meaningful effects on overall bike usage across the system.

\begin{table}[htb]
\centering
\caption{Estimated coefficients for the covariate effects by the proposed fused Lasso regression method from each fold used in the seven-fold CV and from the full data when $r = 1500$.
For comparison, we also report the estimates obtained with fusion Lasso-only, full-interaction and no-interaction methods.}
\label{tab:coef}
\resizebox{\columnwidth}{!}
{\small 
\begin{tabular}{c ccccccc c c c c c}
\toprule
& \multicolumn{7}{c}{Fold} & \multicolumn{5}{c}{Full data}  \\
& 1 & 2 & 3 & 4 & 5 & 6 & 7 & Fused & Fusion-only & Lasso-only & Full & No \\
\cmidrule(lr){1-1} \cmidrule(lr){2-8} \cmidrule(lr){9-9} \cmidrule(lr){10-10} \cmidrule(lr){11-11} \cmidrule(lr){12-12} \cmidrule(lr){13-13}
$\alpha$ & 0.057 & 0.052 & 0.051 & 0.050 & 0.050 & 0.054 & 0.050  & 0.052 & 0.052 & 0.052 & 0.052 & 0.052 \\
$\beta^{\text{rain}}$ & -2.301 & -2.356 & -2.380 & -2.376 & -2.238 & -2.182 & -2.382  & -2.324 & -2.325 & -2.321 & -2.323 & -2.310\\ 
$\beta^{\text{air}}_1$ & 0.098 & 0.141 & 0.100 & 0.177 & 0.124 & 0.169 & 0.168 & 0.140 & 0.141 & 0.138 & 0.138 & 0.137 \\
$\beta^{\text{air}}_2$ & 0.117 & 0.174 & 0.099 & 0.177 & 0.133 & 0.159 & 0.188  & 0.149 & 0.150 & 0.145 & 0.144 & 0.146\\ 
$\beta^{\text{air}}_3$ & 0.218 & 0.285 & 0.251 & 0.169 & 0.259 & 0.321 & 0.310 & 0.275 & 0.278 & 0.271 & 0.276 & 0.243\\
\bottomrule
\end{tabular}}
\end{table}

\begin{figure}[htbp]
\centering
\subfloat{\includegraphics[width = .425\textwidth]{hourly_1500.png}}
\subfloat{\includegraphics[width = .425\textwidth]{daily_1500.png}}
\caption{Parameter estimates for $\phi^{\text{hod}}_h, \, h \in \mathcal H$ (left) and $\phi^{\text{dow}}_d, \, d \in \mathcal D$ (right) from each fold used in the seven-fold CV and from the full data when $r = 1500$.}
\label{fig:thetas}
\end{figure}

\begin{figure}[htbp]
\centering
\subfloat{\resizebox*{4.5cm}{!}{\includegraphics{Tue.png}}}
\subfloat{\resizebox*{4.5cm}{!}{\includegraphics{Sun.png}}
\subfloat{\resizebox*{4.5cm}{!}{\includegraphics{Diff.png}}}
\caption{Estimated station-specific bike demands in log-scale (given by $\widehat\theta_s + \widehat\theta^{\text{hod}}_h + \widehat\theta^{\text{dow}}_d + \widehat\theta^{\text{hod}}_{s, h} + \widehat\theta^{\text{dow}}_{s, d}$) from the model fitted with $r = 1500$ at 8am on Tuesdays (left), at 8pm on Sundays (middle) and their differences (right).} 
\label{fig:fitted}
\end{figure}

In Figure~\ref{fig:thetas}, we plot the estimates of $\theta^{\text{hod}}_h$ and $\theta^{\text{dow}}_d$ which are shared by all the stations belonging to the BSS.
We observe that the smooth transition over the course of a day and a week is well captured across the seven folds, and there are peaks corresponding to the high demand by commuters.
In addition, we plot the combined effects of temporal variables on the mean bike demand (in log-sale), namely $\widehat\theta_s + \widehat\theta^{\text{hod}}_h + \widehat\theta^{\text{dow}}_d + \widehat\theta^{\text{hod}}_{s, h} + \widehat\theta^{\text{dow}}_{s, d}$ under~\eqref{eq:model}, for all stations $s \in \mathcal S$ when $(h, d) = (8, \text{Tu})$ (8am on Tuesday) and $(20, \text{Su})$ (8pm on Sunday) (see Figure~\ref{fig:fitted}).
As expected the spatial distribution of bike usage concentration is markedly different when $(h, d) = (8, \text{Tu})$ and $(h, d) = (20, \text{Su})$, as  bikes are primarily used for commuting at 8am on Tuesdays and for leisure activities at 8pm on Sundays.
In Figures~\ref{fig:thetas} and~\ref{fig:fitted}, we report the results obtained with $r = 1500$ only. (For those obtained with $r = 750$, indicating that the parameter estimates are not sensitive to the choice of $r$, 
we refer to Figures~\ref{fig:thetas:750} and~\ref{fig:fitted:750} in Appendix~\ref{sec:aar}.)

\section{Conclusions}
\label{sec:conc}

In this study, we examine the problem of modelling usage data collected 
from a BSS spanning a large metropolitan city.
We model the data at the granularity of individual stations 
by incorporating covariate effects as well as spatial and temporal characteristics commonly observed in bike usage data.
The proposed multilayer fused Lasso penalty is imposed on the networks encoding  the geographical proximity of the stations with an additional hourly layer, and successfully captures the spatial and temporal homogeneity.
Combined with the data-driven choice of penalty parameters, our penalized regression approach strikes a good balance between a simplistic model that does not allow for station-specific behavior,
and a complex model suffering from high dimensionality,
by returning a fitted model with good predictive performance.

We envision that the proposed method is applicable to different datasets with network-like features,
such as those collected from large transportation, communication, or logistic systems.
In particular, when information about the factors driving link homophily 
(such as the nodal features related to land use, slope of terrain, nearby landmarks, and other modes of transportation in the case of BSSs) is not readily available,
our penalized regression approach enables learning of the linkages in an unobservable network
from the fusion of parameters induced by the penalties defined on an observable network
(such as those determined by the distance between pairs of stations).

\subsubsection*{Acknowledgement}

Yunjin Choi was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT, No. NRF-2020R1G1A1A01014088). Haeran Cho was supported by the Leverhulme Trust Research Project Grant (RPG-2019-390).

\bibliographystyle{apalike}
