An Adapted Geographically Weighted Lasso(Ada-GWL) model for estimating metro ridership

TL;DR

This paper introduces Ada-GWL, a novel local regression model that combines Lasso regularization with geographically weighted regression, tailored for accurate metro ridership estimation considering network-specific distances.

Contribution

The study develops and validates Ada-GWL, a new model integrating local model selection and network-based distance metrics for improved metro ridership estimation.

Findings

Ada-GWL outperforms traditional models in estimation accuracy.

The model reveals spatial variation in ridership elasticities.

Clustering analysis aligns with Shenzhen's land use policies.

Abstract

Ridership estimation at station level plays a critical role in metro transportation planning. Among various existing ridership estimation methods, direct demand model has been recognized as an effective approach. However, existing direct demand models including Geographically Weighted Regression (GWR) have rarely included local model selection in ridership estimation. In practice, acquiring insights into metro ridership under multiple influencing factors from a local perspective is important for passenger volume management and transportation planning operations adapting to local conditions. In this study, we propose an Adapted Geographically Weighted Lasso (Ada-GWL) framework for modelling metro ridership, which performs regression-coefficient shrinkage and local model selection. It takes metro network connection intermedia into account and adopts network-based distance metric instead of Euclidean-based distance metric, making it so-called adapted to the context of metro networks. The real-world case of Shenzhen Metro is used to validate the superiority of our proposed model. The results show that the Ada-GWL model performs the best compared with the global model (Ordinary Least Square (OLS), GWR, GWR calibrated with network-based distance metric and GWL in terms of estimation error of the dependent variable and goodness-of-fit. Through understanding the variation of each coefficient across space (elasticities) and variables selection of each station, it provides more realistic conclusions based on local analysis. Besides, through clustering analysis of the stations according to the regression coefficients, clusters' functional characteristics are found to be in compliance with the facts of the functional land use policy of Shenzhen. These results of the proposed Ada-GWL model demonstrate a great spatial explanatory power in transportation planning.

Figures

Tables

Table 1. Table 1: The table of independent variables summary

Categories	Independent variables	Acronym of variables	Source
Land Use	No. of residential units	Residence	Baidu map
	No. of restaurants	Restaurant	Baidu map
	No. retailers/shopping	Shopping	Baidu map
	No. of schools	School	Baidu map
	No. of offices	Offices	Baidu map
	No. of banks	Bank	Baidu map
	No. of hospitals	Hospital	Baidu map
	No. of hotels	Hotel	Baidu map
Network Structure	Distance to the city center	Dis_to_center	Calculated
	Degree centrality	Degree	Calculated
	Betweenness centrality	Betweenness	Calculated
Social Economics	Population	Pop	Worldpop
	Days since opened	Days_open	Wikipedia
Intermodal Traffic Access	No. of bus stations	Bus	Baidu map

Table 2. Table 2: Moran’s I tests for spatial autocorrelation on the dependent variables and some of the candidate predictors.

	Variable	Moran’s I	Expected I	p-value
Dependent variables	Weekly_ridership	0.24116622	-0.00854701	1.21E-06
	Weekday_ridership	0.26389365	-0.00854701	1.44E-07
	Weekend_ridership	0.16612276	-0.00854701	0.0004006
	Evenrush_ridership	0.27729573	-0.00854701	3.90E-08
	Nonrush_ridership	0.24779725	-0.00854701	4.64E-07
Candidate predictors	Residence	0.39365231	-0.00854701	1.35E-14
	Restaurant	0.34171521	-0.00854701	1.69E-11
	Shopping	0.23155195	-0.00854701	3.92E-06
	School	0.28917679	-0.00854701	4.16E-11
	Offices	0.17494683	-0.00854701	1.37E-05
	Bank	0.38831017	-0.00854701	1.98E-14
	Bus	0.31229433	-0.00854701	7.90E-10
	Hospital	0.21127159	-0.00854701	1.26E-05
	Hotel	0.46648095	-0.00854701	2.2e-16
	Dis_to_center	0.95915286	-0.00854701	2.2e-16
	Degree	0.12583388	-0.00854701	0.005471
	Betweenness	0.24105756	-0.00854701	1.67E-06
	Pop	0.71033294	-0.00854701	2.2e-16
	Days_open	0.38318677	-0.00854701	2.11E-13

Table 3. Table 3: Clusters and their functional characteristics

Cluster	Functional zone	Cluster size (Number of stations)	Functional characteristics	Representative stations
1	Employment and educational based	34	The coefficients of offices and schools land-use are relatively larger than those of other variables, so offices and schools can attract more commuting ridership in this zone. The change of offices and schools land-use will have great impact on the ridership of stations in this cluster.	Shenzhen University; Xili; High-tech Park; Daxin; Qianhaiwan;
2	Commercial land	33	The coefficients of shopping land use are relatively larger than other variables, indicating shopping land use has great positive impact on these stations. Surrounding the stations in this cluster, commercial, recreational and other land uses are signally larger than other categories land uses, which attracts a large amount of non- commuting ridership.	Guomao; Huaqiang Rd; Huaqiang North; Shopping Park; Houhai; Convention&Exhibition Center;
3	Traffic hubs land	29	The coefficients of degree and bus stations are relatively larger than those of other variables, meaning that the transfer stations and bus stations have great impact on ridership. The stations of this kind are mainly located in urban external
			transportation hub (including railway station, passenger transport station, cross-border checkpoint, and Special Economic Zone Port), the station ridership is mainly derived from the passenger flow of traffic hub.	Luohu; Shenzhen North Station; Shekou Port; Chiwan; Buji; Xixiang
4	Residential based	22	The coefficients of residential land use are relatively larger than those other variables, indicating the station ridership is mainly affected by the commuting passenger flow from residential areas. stations of this kind are also mainly located in residential based areas (even though the number of residential units surrounding these stations might be not the largest among other kinds of land use).	Baishizhou; Wuhe; Bu Xin; Jixiang;

Equations

$$y_i=\beta +\underset{k=1}{\sum ^p}{\beta }_k{x}_{ik}+{\epsilon }_i$$

$$y_i={\beta }_0\left(u_i,v_i\right)+\underset{k=1}{\sum ^p}{\beta }_k\left(u_i,v_i\right)x_{ik}+{\epsilon }_i$$

$$y\left(i\right)=X\left(i\right)\beta \left(i\right)+\epsilon \left(i\right)$$

$$\hat{\beta }\left(i\right)={\left[X^{T}W\left(i\right)X\right]}^{-1}{X}^{T}W\left(i\right)y$$

$$w_j\left(i\right)=\exp \left[-{\left(\frac{d_{ij}}{b}\right)}^2\right]$$

$${{w}_{j}}\left(i\right)=\left\{\begin{matrix}{{[1-{{({{d}_{ij}}/b)}^{2}}]}^{2}},&if&{{d}_{ij}}<b\\ 0,&if&{{d}_{ij}}\geq b\\ \end{matrix}\right.\$$

$$\hat{y}\left(i\right)=X\left(i\right)\hat{\beta }\left(i\right)$$

$$E{d}_{rs}=R\ast arccos\left(\cos \left(La{t}_s\right)\ast \cos \left(La{t}_r\right)\ast \cos \left(Lo{n}_s-Lo{n}_r\right)+\sin \left(La{t}_r\right)\ast sin(La{t}_s\right))\ast \frac{\pi }{180}$$

$$Dis{t}_{i,j}=R\ast arccos\left(\cos \left(La{t}_j\right)\ast \cos \left(La{t}_i\right)\ast \cos \left(Lo{n}_j-Lo{n}_i\right)+\sin \left(La{t}_i\right)\ast sin(La{t}_j\right))\ast \frac{\pi }{180}$$

$$N{d}_{rs}=\underset{k\in K_{rs}}{\min }\sum _{i=1}^{n_k-1}Dis{t}_{i,i+1}$$

$${\hat{\beta }}^R=arg\underset{\beta }{\min }\underset{i=1}{\sum ^n}{(y_i-{\beta }_0-\underset{k=1}{\sum ^p}{x}_{ik}{\beta }_k)}^2$$

$$\underset{k=1}{\sum ^p}\left|{\beta }_k\right|\le s,$$

$${\hat{\beta }}^R=arg\underset{\beta }{\min }\{\underset{i=1}{\sum ^n}{(y_i-{\beta }_0-\underset{k=1}{\sum ^p}{x}_{ik}{\beta }_k)}^2+\lambda \underset{k=1}{\sum ^p}\left|{\beta }_k\right|\}$$

$$Dis{t}_i=R\ast arccos\left(\cos \left(La{t}_0\right)\ast \cos \left(La{t}_i\right)\ast \cos \left(Lo{n}_0-Lo{n}_i\right)+\sin \left(La{t}_i\right)\ast sin(La{t}_0\right))\ast \frac{\pi }{180}$$

Full text

Keywords— Geographically Weighted Regression, Lasso, Ada-GWL, network-based distance metric, metro ridership, influencing factors

1 Introduction

With the rapid urbanization, metro has become the mainstream of public transportation for its advantages such as fast speed, large capacity, and convenience. At the stage of urban construction and transportation planning, it is crucial to explore multiple urban indicators from a systematic perspective to capture the urban traffic-related phenomenon like metro ridership at station level, which is a critical element for determining the scale of stations and access facilities. Metro ridership at station level is known to be influenced by the interactions of the different components of an urban system (e.g., land-use, socio-economics, etc.). Understanding the influence of these components is vital to accurately estimate travel demand and to effectively make design schemes of urban systems including the identification of which public infrastructures, services, and resources need to be built and deployed. Modelling metro ridership at station level can help to not only estimate and forecast ridership but also analyse the influencing factors on it.

Various methods have been proposed for transit ridership estimation. As one of the best-known models, the four-step (generation, distribution, mode choice, and assignment) model has dominated the history of transport modeling since the 1950s [1]. However, the four-step model has many drawbacks in practice [2], such as limitation in model accuracy, low data precision, insensitivity to land use, institutional barriers, and high expense [3]. The four-step model is generally effective for estimating transit ridership on a regional scale rather than more detailed scales (such as station level) [4]. As an alternative to the four-step model, direct demand models have drawn growing attention for ridership estimation in recent decades. Direct demand models estimate ridership as a function of influencing factors within the Pedestrian Catchment Areas (PCA) via regression analysis, which enable identifying factors that contribute to higher transit ridership[3, 5, 6, 7, 8]. In the models, a PCA is a geographic area for which a station attracts passengers. The size and shape of a catchment area depend on how accessible a station is and how far it is from alternative stations. One can use buffers to create circular catchment areas by a specific distance or use Thiessen polygons to illustrate the area most accessible to each station. The major advantages of direct demand models in travel analysis are simplicity of use, easy interpretation of results, immediate response, and low cost. The direct demand models generally use OLS multiple regression , which assumes parametric stability. With the development of spatial modeling, direct demand models could increase their spatial explanatory power by using GWR, which is designed to model spatial parametric nonstationarity, and variance heterogeneity. However, the issues of GWR models, such as collinearity in the estimated coefficients, and network connection intermedia may influence metro ridership modelling and estimation. Thus, there are still a number of limitations to be promoted in the existing research on transit ridership modeling.

In light of deficiencies of current direct demand models, considering the pros and cons of Geographically Weighted Regression (GWR) for modelling potentially spatially varying relationships, this paper proposes an Adapted Geographically Weighted Lasso (Ada-GWL) framework for modelling the ridership of metro systems, which considers using network-based distance metric rather than Euclidean-based distance metric, and simultaneously shrink regression coefficients and select feature locally in the process of GWR with Least Absolute Shrinkage and Selection Operator (LASSO).The first objective of this study is to identify the association between multiple factors and metro station ridership from a local perspective. The second objective is to estimate ridership at station level accurately. An Ada-GWL model is used to analyze the significance of various local factors that impact metro station ridership. The ridership data of Shenzhen Metro at station-level in the year of 2013 is used to validate the superiority of our proposed model over global model (Ordinary Least Square), GWR, GWR calibrated with network-based distance metric, and GWL model.

Our contributions are three-fold:

(1) This study introduces network-based distance metric during estimating ridership. By implementing the Ada-GWL model on the real world case study of Shenzhen metro systems, we calibrate distance with network-based distance metric instead of Euclidean-based distance metric. The reason is that locations may be related through a certain intermedia system but not in Euclidean space. This paper investigates ridership at metro stations, which are related through railways. In practice, taking network connection intermedia into consideration by using network-based distance metric is important as routine transportation activities normally occur on a determined transportation network but not a continuous space [9]. Nonetheless, few existing direct demand models have included spatial aspects considering network connection intermedia.

(2) We adopt the basic framework of GWR considering spatial autocorrelation of variables to model metro ridership, which outperforms the traditional global regression model OLS in terms of model fitting and spatial explanatory power.

(3) This study enables simultaneous coefficient penalization and model selection during GWR. Even though GWR models have been utilized to model and analyze the transit ridership and its influencing factors, the collinearity of GWR may make interpreting individual coefficients problematic. Also, different influencing factors may exist at different stations, so the feature selection of the model can be conducted from the local instead of global perspective, that is, local penalization and feature selection, which have rarely been done in GWR research.

In general, to the extent of our knowledge, this is the first work that investigates the utility of GWL model calibrated with network-based distance metric for estimating metro ridership and analyzing the influencing factors from a local perspective. It is worth noting that the fundamental framework proposed is not limited to the specific dataset used in this study, but can be extended and generalized to other applications in transportation systems, such as light rail and rapid bus transit systems.

The remainder of this paper is organized as follows. In Section 2, we review the current literature on transit ridership estimation modeling and factors investigated in those models; In Section 3, we present the data description and various influencing factors considered; In Section 4, we describe the methodology for estimating metro transit ridership and identifying significant influencing factors. Section 5 provides the details for model implementation and results analysis and discussion. Finally, we conclude this paper in Section 6.

2 Literature Review

2.1 Transit Ridership Estimation Modelling

This paper summarizes the related studies on direct demand models for ridership estimation (see Figure A1 in A.1). As summarized in Figure A1, the most widely used method is Ordinary Least Squares (OLS) or linear least squares in the study of the influencing factors on transit ridership by using direct demand models, Kuby et al.(2004), Sohn and Shim(2010), Loo et al.(2010), Sung and Oh(2011), Gutierrez et al.( 2011), Thompson et al.(2012), Guerra et al.(2012), Zhao et al.(2013), Chan and Miranda-Moreno (2013), Singhai et al.(2014), Liu et al.(2016) and Pan et al.(2017) all applied OLS regression to model transit ridership and its influencing factors[7, 10, 11, 12, 3, 13, 14, 15, 16, 17, 18, 19]. However, the limitation of OLS models is that they assume the transit ridership is affected by various factors but has nothing to do with the spatial location, without considering the spatial autocorrelation, in other words, OLS models estimate ridership from a global perspective believing that calculated coefficients do not have significant differences in space. Fotheringham (1996) proposed Geographically Weighted Regression (GWR) model which can reveal the spatial relations on the condition of spatial heterogeneity. GWR model has a strong capability of spatial data analysis [20], and it has been widely used to analyse spatial data in economics, geography, ecology, and many other fields of research. At present, Cardozo et al.(2012) compared the performance of OLS and GWR in modelling transit ridership and its influencing factors, and GWR showed the better goodness-of-fit than OLS for forecasting station-level ridership[4]. In addition, Jun et al.(2015) used mixed geographically weighted regression (MGWR) models consisting of local and global independent variables to examine the relationship between the land use characteristics and subway ridership[21]. However, there was still room for improvement in their GWR models, such as considering the impact of network connection intermedia on distance metric and correlation in the estimated coefficients during GWR [22]. Actually, Wheeler and Tiefelsdorf (2005) showed that GWR coefficients can be correlated even though there is no independent variable collinearity, the coefficient correlation increases systematically with increasing collinearity[23]. Therefore, the application of GWR in the estimation of transit ridership needs to be promoted.

The research of GWR with penalized forms (e.g., $\ell 1$ norm/ $\ell 2$ norm) can be found in Wheeler’s studies[24, 22]. Both ridge regression ($\ell 2$ norm) and the Lasso ($\ell 1$ norm) are penalization methods which impose constraints on the regression coefficients. Ridge regression aims to alleviate collinearity effects by making use of the $\ell 2$ penalty to regularize the size of regression coefficients and reducing variance. The Lasso also performs regularization on variables in consideration. It considers the absolute values of the sum of the regression coefficients. It also sets the coefficients to zero thus simultaneously performing regression coefficient penalization and model selection. Wheeler (2007) implemented a ridge-regression version of GWR, called GWRR, and found it could constrain the regression coefficients to reduce local correlation effects in an existing dataset and showed a lower estimation error for the dependent variable than of GWR [24]. The lasso has also been introduced into the GWR framework in Wheeler (2009)’s study. The lasso is considered in the GWR framework since it has capability to simultaneously coefficient regularization and local model selection, also for its potential to reduce prediction and estimation errors for estimating the dependent variable in GWR. Wheeler (2009) implemented a Lasso version of GWR, called Geographically Weighted Lasso (GWL), and demonstrated the benefit of GWL through a comparative analysis with GWR and GWRR.It can stabilize regression coefficients in the presence of collinearity and produce lower estimation error of the dependent variable than does GWR, and also offered a key advantage to GWRR for the identification of insignificant local effects by doing model selection[22]. However, the existing version GWL model didn’t consider the effect of network connection intermedia, which opened up new perspectives that fully justified the pursuit of this research. Thus, we try to design a version of GWL model adapted to metro network that considering network connection intermedia during calculating distance.

A comprehensive review of direct demand models can also be found in the work by Walters and Cervero (2003) and Cardozo et al. (2012) [25, 4].

2.2 Dependent Variables and Independent Variables

With regard to the dependent variables of direct demand models considered, daily ridership was selected in most of the relevant studies (also see Figure A1 in A.1), such as the research of Kuby et al.(2004), Chu (2004), Sohn and Shim (2010) , Loo et al.(2010), Zhao et al.(2013) and Zhang and Wang(2014), which all took the average weekday ridership as dependent variable[7, 8, 10, 11, 15, 9]. In addition, monthly station ridership was considered in the research of Gutierrez et al.(2011)[3]. For shorter time span, Zhao et al.(2013), Singhal et al.(2014), Hu et al.(2016), Li et.al (2016) and Chan and Miranda-Moreno (2013)chose the hourly (eg. Peak hour) ridership as the dependent variable of their models[15, 17, 30, 29, 16].

Concerning independent variables of models, they can be roughly divided into the following categories, Land use variables, social economic variables, accessibility and the network structure. Figure A1 summarized the independent variables considered in existing studies. For the first kind of variables, such as the research of Estupinan and Rodriguez (2008), Loo et al. (2010), Sohn and Shim (2010), Gutierrez et al. (2011), Sung and Oh (2011), Choi et al. (2012), Cardozo et al. (2012) and Zhao et al (2013)[27, 11, 10, 3, 12, 5, 4, 15], which considered commercial, residence, education, entertainment and other mixed land use as independent variables. For the second kind of variables, factors considered mainly are population, employment, and automobile ownership ratio. For example, Chu (2004), Kuby et al. (2004), Loo et al. (2010), Sohn and Shim (2010), Gutierrez et al. (2011), Choi et al. (2012), Cardozo et al. (2012), Guerra et al. (2012), Thompson et al. (2012) and Zhao et al. (2013) analyzed the influence of population and employment on transit ridership [8, 7, 11, 10, 3, 5, 4, 14, 13, 15]. Chu (2004), Loo et al. (2010), Thompson et al (2012), Cardozo et al. (2012) and Zhao et al. (2013) considered the relationship between automobile ownership and the ridership[8, 11, 13, 4, 15]. For the third category, accessibility, Chu (2004), Kuby et al. (2004), Estupinan and Rodriguez (2008), Loo et al. (2010), Sung and Oh (2011), Gutierrez et al. (2011), Choi et al. (2012), Cardozo et al. (2012), Guerra et al. (2012) and Zhao et al. (2013) analyzed the influence of bus feeder system on transit ridership[8, 7, 27, 11, 10, 13, 5, 14, 15]. Moreover, Kuby et al. (2004), Estupinan and Rodriguez (2008), Loo et al. (2010), Sohn and Shim (2010), Thompson et al. (2012), Choi et al. (2012), Guerra et al. (2012) and Zhao et al. (2013) studied on the station accessibility[7, 27, 11, 10, 13, 5, 14, 15]. Finally, with regard to network structure, Kuby et al. (2004) took the influence of transfer station and terminal station on ridership into account, but both of them were regarded as dummy variables in their models[7]. So far, few relevant studies have been carried out from the quantitative perspective by using measurements from the field of complex network, such as quantifying the practical significance (e.g., transfer, terminal station and the importance of stations’ positions that lie on a large number of paths taken by flows) by the calculation of nodes’ degree centrality and betweenness centrality in networks.

3 Methodology

3.1 Geographically Weighted Regression

GWR model is an extension of OLS, which can be formulated as:

[TABLE]

When geographical location factors are introduced into regression parameters to allow local parameter estimation, OLS model evolves into GWR model, which is as following:

[TABLE]

Where ${{y}_{i}}$ and ${{x}_{i1}},{{x}_{i2}},...,{{x}_{ip}}$ are observed values of dependent variable and independent variables ${{x}_{1}},{{x}_{2}},...,{{x}_{p}}$ at the location of $\left({{u}_{i}},{{v}_{i}}\right)$, which is geospatial coordinates of observation point $i~{}\left(i=1,2,\ldots,n\right)$, called a model calibration location. In the context of metro networks, ${{y}_{i}}$ and ${{x}_{i1}},{{x}_{i2}},...,{{x}_{ip}}$ are observed metro ridership and potential influencing factors ${{x}_{1}},{{x}_{2}},...,{{x}_{p}}$ at station $i~{}\left(i=1,2,\ldots,n\right)$ with the geospatial coordinates $\left({{u}_{i}},{{v}_{i}}\right)$. ${{\varepsilon}_{i}}$ is the normally distributed error term (with the expected value 0 and a constant variance). {{\beta}_{k}}\left({{u}_{i}},{{v}_{i}}\right)$$\left(k=1,2,\ldots,p\right) refers to $p$ unknown coefficients associated with spatial position pending for estimation. $p$ is the number of independent variables. The GWR model at location $i$ in matrix term is:

[TABLE]

The vector of estimation regression coefficients at location $i$ is:

[TABLE]

Where $X={{[{{X}^{T}}(1);{{X}^{T}}(2);...;{{X}^{T}}(n)]}^{T}}$ is the matrix of independent variables, which typically includes a column of 1s for the intercept. $W\left(i\right)=diag\left[{{w}_{1}}\left(i\right),\ldots,{{w}_{n}}\left(i\right)\right]$ is the diagonal weights matrix, which contains the geographical weights in its leading diagonal and 0 in its off-diagonal elements. The weights matrix is calculated for each location $i$ and applies weights to observations $j=1,\ldots,n$. Here, more weight generally will be placed more emphasis on closer observations to the model calibration location $i$. $y$ is the $n\times 1$ vector of dependent variable values; and $\hat{\beta}\left(i\right)={{\left({{{\hat{\beta}}}_{i0}},{{{\hat{\beta}}}_{i1}},\ldots,{{{\hat{\beta}}}_{ip}}\right)}^{T}}$ is the vector of $p+1$ local regression coefficients at location $i$ for $p$ independent variables and an intercept term.

The weights matrix $W\left(i\right)$ is calculated from a kernel function $K(j,i)$. There are many kinds of kernel functions. According to the nature of the bandwidth , kernel functions can be divided into two groups: the fixed kernel function and the adaptive kernel function. Gaussian kernel function is a typical example of the fixed kernel, which was expressed as Equation (5).

[TABLE]

where ${d}_{ij}$ is the distance between the calibration location $i$ and location $j$, and $b$ is the kernel bandwidth parameter.The fixed kernel function assumes that the bandwidth $b$ at each location $i$ is a constant across the study area. If the locations $i$ and $j$ satisfy ${d}_{ij}=0$, ${w}_{j}\left(i\right)=1$, whereas ${w}_{j}\left(i\right)$ decreases according to a Gaussian curve as ${d}_{ij}$ increases. However, the weights are nonzero for all data points, no matter how far they are from the location $i$ [20]. If the sample points are distributed as clusters in the study area but not regularly spaced, it is generally appropriate to allow the kernel to adapt this irregularity by increasing its size when the sample points are distributed sparser and decreasing its size when they are denser. A convenient way of implementing this adaptive bandwidth specification is to choose a kernel which allows the same number of sample points for each estimation. A commonly used adaptive kernel function is the bisquare kernel function, as follows:

[TABLE]

where ${w}_{j}\left(i\right)=1$ at the location $i$ (${d}_{ij}=0$) and ${w}_{j}\left(i\right)=0$ at the location $i$ when ${d}_{ij}\geq b$ [31].

Both the Gaussian kernel function and the bi-square kernel function are highly dependent on bandwidth $b$. Thus, the kernel bandwidth $b$ should be estimated firstly to fit the GWR model. The most common method in practice is leave-one-out cross-validation (CV) across all the calibration locations proposed by Cleveland(1979) [32]. CV is an iterative process that finds the kernel bandwidth with the lowest associated prediction error of all the dependent variables $y\left(i\right)$ . For each calibration location $i$, it removes the data for observation $i$ in the model calibration at location $i$ and predicts $y\left(i\right)$ from the other data points leaving $i$ out and the kernel weights associated with the current bandwidth. The Leave-One-Out CV (LOO-CV) function is ${{\min}_{b}}\sum\limits_{i=1}^{n}{[{{y}_{i}}}-{{\hat{y}}_{\neq i}}(b){{]}^{2}}$, where ${{\hat{y}}_{\neq i}}\left(b\right)$ is the fitted value ${y}_{i}$ with location $i$ omitted during the fitting process. The optimal bandwidth $b$ minimizes the CV function.

An alternative to CV in kernel-bandwidth estimation is the Akaike information criterion (AIC), proposed by Fotheringham et al (2002)[33]. CV and the AIC are also methods applied in model selection, and more details about AIC and model selection can be referred to Burnham and Anderson, 2002 [34]. There is little published research on discussing if CV and AIC will calculate an uniform solution or whether one method is preferred under a certain condition.In this paper, we adopt CV method. Actually, the selection of kernel function and bandwidth depends on the specific distribution of sample points in the study area, which has a strong impact on the descriptive and predictive power of GWR models [31]. Finally, with the kernel weights calculated at each calibration location $i$ and the determined bandwidth $b$ from the kernel function $K$, the regression coefficients are estimated at each model calibration location, consequently, the values of dependent variable can be estimated as following:

[TABLE]

3.2 Distance Matrix considering Network Structure

In practice, since the GWR model is applied to metro station ridership estimation, it is necessary to take the network connection intermedia into consideration when calculating the distance, however, existing studies usually calculate the distance matrix only considering coordinates, i.e.,just using Euclidean-based distance metric to calibrate the distance no matter how the locations connected. In practice, in metro networks, one station is connected to another station with railway, so the Euclidean-based distance cannot reflect the actual distance between two stations. The study uses network-based distance to measure the spatial connectivity instead of Euclidean-based distance.The two distance measurement methods can be illustarted as Figure 1.

The Euclidean-based distance $Ed_{rs}$ between stations $r$ and $s$ can be calculated by:

[TABLE]

Here, we calculate the distance $Dis{{t}_{i,j}}$ between two close stations $i$ and $j$ considering the radius of the earth $R$ according to Equation (9):

[TABLE]

Thus, the network-based distance between stations $r$ and $s$ $Nd_{rs}$ can be calculated by:

[TABLE]

Where $k$ denotes the $k$-th path between stations $r$ and $s$,and $k\in{{K}_{rs}}$. ${{K}_{rs}}$ is the set of paths between stations $r$ and $s$. ${n}_{k}$ is the number of stations the $k$-th path passes through. So the network-based distance $Nd_{rs}$ is exactly the distance of the shortest path along the railway lines, which considers the network connection intermedia. In this way, the $n\times n$ interpoint network-based distance matrix $\mathbf{D}$ can be obtained. Therefore, the shortest path length of every two stations considering railway connection will be adopted to calculate the network-based distance matrix. Here, we made use of distances fuction in “igraph” package (see: https://cran.r-project.org/web/packages/igraph/index.html)to calculate the length of all the shortest paths between vertices in the network. As we adjusted the distance measurement of metro networks, the weights calculated from kernel functions will be changed accordingly.

3.3 Least Absolute Shrinkage and Selection Operator (LASSO)

To address the issue of collinearity in the GWR framework, the Least Absolute Shrinkage and Selection Operator (LASSO) is deemed a proper method because it more directly reduce the variance in the regression coefficient while retaining interpretability of covariate effects[23].In this study, we will also consider the lasso version of GWR framework. As a kind of shrinkage method, the lasso not only tends to reduce the variability of the estimates thus can improve the model’s stability, but also can set some of the coefficients to zero, thus allows for variable selection. The lasso makes use of the $\ell 1$ norm. $\ell 1$ penalties are convex and the assumed sparsity can lead to significant computational advantages. The lasso is defined as following:

[TABLE]

subject to

[TABLE]

where $s$ is a parameter that controls the degree of coefficient shrinkage. Tibshirani (1996) proved that the lasso constraint $\underset{k}{\mathop{\sum}}\,\left|{{\beta}_{k}}\right|\leq s$ is equivalent to adding the penalty term $\lambda~{}\underset{k}{\mathop{\sum}}\,\left|{{\beta}_{k}}\right|$ to the Residual Sum of Squares (RSS). Thus a direct relationship between $s$ and $\lambda\geq 0$ which is a complexity parameter that controls the degree of shrinkage of coefficients. Hence, the lasso coefficients can be also expressed as:

[TABLE]

The generally used methods for solving the lasso are standard convex optimizer [35] and least angle regression (LARS) (Efron et al 2004)[36]. In this study, we adopted LARS to solve the lasso. The procedure of the LARS algorithm can be referred to Wheeler, 2009[22].

3.4 Adapted Geographically Weighted Lasso (Ada-GWL) model

To implement the lasso in GWR, Wheeler (2009) developed a framework called Geographically Weighted Lasso (GWL) and wrote an R package “gwrr” including gwl.est function to run GWL models [43]. Based on to the framework of Wheeler (2009), we developed our Ada-GWL model in the context of metro ridership estimation by taking network connection intermedia into consideration[22]. Here, we called lars function in “lars” package [45] to implement lasso methods and GWmodel package to implement GWR framework [42].

For solving the Ada-GWL model we proposed, the algoritms including the sub-routines and the main body are designed. First of all, the sub-routine to find the best local LARS solution at each location by minimizing LOO-CV Root Mean Square Error (RMSE) using LARS is shown in Algorithm 3.4.

Algorithm 1 LocalLARS() algorithm for obtaining the best LARS solution at each location by minizing CV RMSE.

0:   $b$, the bandwidth; $\mathbf{X}$, a set of observed values of independent variables (factors) ${{x}_{1}},{{x}_{2}},...,{{x}_{p}}$ at the location of $\left({{u}_{i}},{{v}_{i}}\right)$; $\mathbf{y}$, a set of observed values of dependent variable(ridership); $\mathbf{D}$, the distance matrix; $n$, the number of observation points(stations); $K$, the kernel function. 

0:   ${{s}_{i}}$, the shrinkage parameter; $\mathbf{z}$, the indicator vector of which variable coefficients are shrunken to zero; $CVerr$, the minimum CV RMSE. 

1:  for each station $i,i=1,\ldots,n$ do

2:     Calculate $\mathbf{W}(i)$ using $\mathbf{D}$ and $b$ according to Equation (5);

3:     ${{\mathbf{W}}^{1/2}}\left(i\right)\leftarrow$ $\sqrt{\mathbf{W}(i)}$, and${{\mathbf{W}}^{1/2}}{{\left(i\right)}_{ii}}\leftarrow 0$;

4:     ${{X}_{w}}={{\mathbf{W}}^{1/2}}\left(i\right)\mathbf{X}$ , and ${{y}_{w}}={{\mathbf{W}}^{1/2}}\left(i\right)\mathbf{y}$ ${{\mathbf{W}}^{1/2}}\left(i\right)$ at station $i$;

5:     Call function lars$({{X}_{w}},{{y}_{w}})$;

6:     return  $Parameters({{s}_{i}},\mathbf{z},CVerr)=\underset{parameters}{\mathop{\arg\min}}\,{{\left({{y}_{i}}-{{{\hat{y}}}_{i}}\right)}^{2}}$

7:  end for

Besides, the sub-routine of bandwidth selection algorithm via binary search is shown in Algorithm 3.4.

Algorithm 2 BandwidthSelector() algorithm for selecting the bandwidth $b$ by the binary search for the minimum error.

0:   $lb$, the lower bound of the initial bandwidth; $ub$, the upper bound of the initial bandwidth; $eps$, loop control value; $\mathbf{X}$, a set of observed values of independent variables (factors) ${{x}_{1}},{{x}_{2}},...,{{x}_{p}}$ at the location of $\left({{u}_{i}},{{v}_{i}}\right)$; $\mathbf{y}$, a set of observed values of dependent variable(ridership); $\mathbf{D}$, the distance matrix; $K$, the kernel function. 

0:   $b$, the bandwidth with the minimum LOO-CV error; $RMSE$, the minimum root mean square error (LOO-CV error). 

1:  $diff\leftarrow ub-lb$, $c\leftarrow(lb+ub)/2$,

2:  while $diff>eps$ do

3:     $l.c\leftarrow(lb+c)/2$, and $c.u\leftarrow(c+ub)/2$;

4:     $l.c.err\leftarrow$ the $CVerr$ returned from LocalLARS$(l.c,\mathbf{X},\mathbf{y},\mathbf{D},n,K)$;

5:     $c.u.err\leftarrow$ the $CVerr$ returned from LocalLARS$(c.u,\mathbf{X},\mathbf{y},\mathbf{D},n,K)$;

6:     if $l.c.err<c.u.err$ then

7:        $ub\leftarrow c.u$;

8:     end if

9:     if $l.c.err>c.u.err$ then

10:        $lb\leftarrow l.c$;

11:     end if

12:     $c\leftarrow(lb+ub)/2$, $diff\leftarrow|ub-lb|$;

13:  end while

14:  Call the function LocalLARS $(lb,\mathbf{X},\mathbf{y},\mathbf{D},n,K)$, LocalLARS $(ub,\mathbf{X},\mathbf{y},\mathbf{D},n,K)$, and LocalLARS $(c,\mathbf{X},\mathbf{y},\mathbf{D},n,K)$, and compare their $RMSE$ values;

15:  return  $b,RMSE$

With the introduction of sub-routines earlier, the main body for implementing Ada-GWL is shown as Algorithm 3.4:

Algorithm 3 Pseudo code to estimate the Ada-GWL solutions.

0:  $G\left(v,e\right)$, the directed graph generated from a metro network; $\left({{u}_{i}},{{v}_{i}}\right)$, the geospatial coordinates of observation station $i~{}\left(i=1,2,\ldots,n\right)$; $\mathbf{X}$, a set of observed values of independent variables (factors) ${{x}_{1}},{{x}_{2}},...,{{x}_{p}}$ at the location of $\left({{u}_{i}},{{v}_{i}}\right)$; $\mathbf{y}$, a set of observed values of dependent variable(ridership); $K$, the kernel function. 

0:   $b$, the kernel bandwidth parameter; $\hat{\beta}\left(i\right)={{\left({{{\hat{\beta}}}_{i0}},{{{\hat{\beta}}}_{i1}},\ldots,{{{\hat{\beta}}}_{ip}}\right)}^{T}}$, the estimated vector of $p+1$ local regression coefficients at station $i$ for $p$ independent variables and an intercept term; ${{\hat{y}}_{i}}$, the fitted value of ridership at each station $i$; $RMSE$, root mean square error; $R-Squared$. 

//Phase 1: Estimate the local scaling GWL parameters when the LOO-CV error is minimum. //

1:  Calculate the network-based distance $\mathbf{D}$ according to Equation (9) and Equation (10) ;

2:  Set the upper and lower bound of the initial bandwidth $ub\leftarrow ceiling(max(\mathbf{D}))$, and $lb\leftarrow min(\mathbf{D})+0.01\times ub$;

3:  $eps\leftarrow RMSE(OLS)\times 0.05$ (or set some other feasible small values) ;

4:  ${b,RMSE}\leftarrow\textsf{BandwidthSelector}(lb,ub,eps,\mathbf{X},\mathbf{y},\mathbf{D},K)$;

5:  ${{s}_{i}},\mathbf{z}\leftarrow\textsf{LocalLARS}(b,\mathbf{X},\mathbf{y},\mathbf{D},n,K)$;

//Phase 2: Estimate the final local scaling GWL solutions using the shrinkage parameter${{s}_{i}}$ and bandwidth $b$ estimated in Phase 1.//

6:  Using the bandwidth $b$ estimated in Phase 1, call the function LocalLARS $(b,\mathbf{X},\mathbf{y},\mathbf{D},n,K)$;

7:  Estimate $\hat{\beta}\left(i\right)$ according to Equation (4), and estimate ${{\hat{y}}_{i}}$ according to Equation (7);

8:  return  $b,\hat{\beta}\left(i\right),{{\hat{y}}_{i}},RMSE,R-Squared$

4 Empirical Study Area and Data

In this paper, we concern investigating factors influencing the transit ridership at metro station level in Shenzhen, which is supported by a relatively large Metro transportation network, consisting of 5 lines (Line 1(Luobao Line), Line 2(Shekou Line), Line 3(Longgang Line), Line 4(Longhua Line) and Line 5 (Huanzhong Line)) and 118 stations in the year of 2013. Figure 2(a) is the schematic map of Shenzhen metro in 2013111Source: http://toursmaps.com/wp-content/uploads/2017/02/shenzhen_metro_map-1.gif. Figure 2(b) shows the spatial distribution of those stations. The population of Shenzhen at the end of 2013 is about 10,628,900, and the population density is 5,323 persons/km2 222Source: http://www.sztj.gov.cn/nj2014/indexeh.htm.

4.1 Data Description

The Shenzhen metro ridership data at station level were aggregated by using the data collected through Automatic Fare Collection (AFC) system of Shenzhen Metro Corporation in China. The dataset includes the total information about entry-exit smart card records. The data used in the research cover a time span of 7 days from October 14th (Mon) to 20th (Sun) in the year of 2013, which includes summed boarding and alighting ridership amounts at two time resolutions with five kinds of time periods: time of the day ((1) rushhour and (2) non-rush hour) and day of the week ((3) average weekday, (4) average weekend and (5) average daily ridership of the whole week). The ridership data at two time resolutions are adopted as dependent variables in the fitting of models. The independent variables represent factors hypothesized to influence station ridership (A detailed description is provided in Table1). These variables can be classified into four categories: (1) Social economic variables; (2) Land use variables; (3) Intermodal traffic access variables; and (4) Network structure variables. Before collecting the data such as land use-related variables and population, a critical first step is to evaluate the walking distance to metro stations, i.e., the size of a PCA, aiming to determine the range of data collection. As the average friendly walking distance is generally assumed to be 500m in large and middle-sized cities according to Kim et al. (2017) [37], we also define the distance of PCA of each Shenzhen metro station as 500m. In our work, we use a buffer to create circular PCA by 500m. Based on the buffer with 500m radius determined, population, all of the land use-related data and the number of bus stations are collected subsequently.

4.2 Dependent variable

This paper aims to identify and analyze different factors influencing the ridership at station level. In all available days(Oct. 14th- Oct. 20th), the count of average daily records of weekdays is about 2,513,330, and that of weekends is about 2,480,449, which is less than that of weekdays, indicating the daily trip frequency of weekdays is higher than that of weekends.

We conduct preliminary statistical analysis using metro AFC data on October 14. Figure3(a) shows the spatial distribution of AFC data records in one day. It presents that the records are most densely distributed at Grand Theatre Station and Laojie Station, closely followed by Huaqiang Road station and Luohu station, and the records of other stations have relatively sparse distribution.

Figure3(b) is about temporal distribution of AFC data records, which shows that the spatial distribution of records has a peak value at both 8:00 and 18:00 on no matter weekdays or weekends, and the records in rush hours on weekends are significantly less than those on weekdays. In addition, it is found that the peak value in the morning rush hour is generally less than that in the evening rush, which may result from people’s travel time (such as commuting time) inconsistency in the morning while consistency in the evening. Besides, the records of other periods except for rush hours of weekends are more than those of weekdays. A possible reason could be that people take more non-commuting trips on weekends than on weekdays.

Additionally, the characteristics of temporal distribution of records on weekdays and weekends are quite similar, which suggests that there are similar metro travel patterns with morning and evening peaks on weekdays and weekends in Shenzhen. Moreover, from the perspective of distribution of records in a single day, it is also found that the maximum evening peak value falls on Friday. It indicates that more people tend to take trips in addition to daily commutes on Friday evening because of the last workday.

According to the analysis above, the travel demands and travel patterns not only differentiate in spatial but also differentiate at different time resolutions, such as the level of time of the day and day of the week.

4.3 Independent variables

The independent variables represent factors hypothesized to influence station ridership (Table1). The variables can be classified into four categories: (1) Land use; (2) Social economic variables; (3) Intermodal traffic access variables; and (4) Network structure variables.

4.3.1 Land use variables

All of the land use-related data within a PCA were collected from Baidu Map with the assistance of API, and land use variables consist of the stations’ nearby residence, entertainment, services, business, education, offices. Specifically, the information covers the numbers of residence, restaurants, schools, offices, hospitals, banks, shopping places, and hotels within 500m PCA.

4.3.2 Social economics

With regard to social economic variables, they consist of the population distribution of Shenzhen in the year of 2013 and operation days since metro stations opened. The information of days since metro lines and stations opened was collected from “Wikipedia” 333Source: https://en.wikipedia.org/wiki/Shenzhen_Metro . The population data were collected from the website of Worldpop 444Source: http://www.worldpop.org.uk/data/get_data/. The format of the population file is Geotiff. The file provides estimated numbers of people per grid square at $8.33\times{{10}^{-4}}$ degrees spatial resolution (approximately 100m at the equator), which can be projected to “GCS-WGS-1984” geographic coordinate system. During the data pre-processing, the population within each buffer can be obtained by summing up the value of the grid falls into the metro station buffer. Figure 4 shows the population distribution of the whole city of Shenzhen in the year of 2013. Meanwhile, the buffers of metro stations with a radius of 500 m were created by using ArcGIS 10.2[39], which are also illustrated in Figure 4. Through the preliminary visualization in Figure 4, it was noted that population is densely distributed near the metro region. The influence of population density within each station buffer on ridership is pending for analysis in the model.

4.3.3 Intermodal traffic access variables

As for intermodal traffic access, here we only considered the feeder bus system, and the number of bus stations within 500m PCA of a metro station was hypothesized to be positively related to station ridership, which was also collected from Baidu Map 555Source: https://map.baidu.com/.

4.3.4 Network structure variables

In this paper, network structure variables comprise the degree centrality and betweenness centrality of the metro network nodes, and spatial characteristics of each station: distance to city center. In the field of complex networks, as degree is a simple centrality measure that counts how many neighbors a node has, and the betweenness centrality for each node refers to the number of shortest paths that pass through the node (Erciyes, 2015) [40], thus they are correlated to the information for transfer stations or terminal stations, and the importance of stations in the aspect of their controlling over flows passing between others of metro networks. As for the distance $Dis{{t}_{i}}$ of each station i to the city center, i.e., Shenzhen Municipal People’s Government, located in Futian District, we calculate it by the following (14) considering the effect of the radius of the earth:

[TABLE]

Where, $R$ is the radius of the earth, $(La{{t}_{0}},~{}Lo{{n}_{0}})$ and $(La{{t}_{i}},~{}Lo{{n}_{i}})$ are the latitude and longitude of the city center and station $i$, respectively. The related geographical data were collected from Google map 666Source: https://maps.google.com.

5 Model Implementation and Results Discussion

In this section, we present the implementation details for metro station ridership modeling using the proposed Ada-GWL. We then use this implementation in our case study for results comparisons (with OLS, GWR, and GWL models), aiming to demonstarte the superiority of Ada-GWL model, and further analyze the results of Ada-GWL combined with the practical meanings.

5.1 Implementation

5.1.1 Data processing

First of all, for the sake of convenience in representing and understanding, we use alphanumeric code instead of Chinese to denote each station name. Here, we define identifiers for station names according to the following rules: (1) non-transfer stations consist of 3 digits, where the first digit denotes the line number, and the rest 2 digits denote the sequential number of station; (2) transfer stations start with character t followed by 3 digits, where the first 2 digits denote the intersection of two lines, and the last digit means the sequential number of intersections between those two lines. For example, “402” represents the 2nd station of Line 4, and “t131” represents the transfer station that is the first intersection of lines 1 and 3. In this way, all of 118 stations can be encoded by such identifiers containing the line and station information literally.

Shenzhen metro implements AFC system in the whole rail network, which brings chance to record accurately the information of passengers’ boarding and alighting points. The information makes it possible to possess the temporal and spatial distribution law of ridership, getting the real-time ridership data of each station. This paper is based on the boarding and alighting ridership (smart card) data of 118 stations of Shenzhen metro covering October 14th to October 20th (a whole week), 2013. In the study, we sum up the boarding (entry) and alighting (exit) records at different time resolutions, such as the level of time of the day and day of the week. Therefore, different regression models with different dependent variables, those are rush-hour(evening:17:00-19:00) ridership, non-rush hour (9:00-17:00, 19:00-23:00) ridership, average weekday ridership, average weekend ridership and average daily ridership of the whole week (the operation times of Shenzhen metro is 6:30-23:00), will be built intending to validate the model’s effectiveness and applicability to different time resolutions, and further to find if the factors influencing the station-level ridership differentiate at different time resolutions 777Average daily ridership of a whole week is the average of total ridership of seven days of week in operation times (6:30-23:00). Rush hour ridership is calculated by the total ridership of evening rush-hours from 17:00 to 19:00 of a whole week divided by 14 hours (multiply 2 hours by 7 days). Non-rush hour ridership is calculated by the total ridership of remaining time (9:00-17:00, 19:00-23:00) except for morning (7:00-9:00) and evening rush-hours of a whole week divided by 84 hours (multiply 12 hours by 7 days)..

For the acquisition of data as independent variables mostly relies on the collection with the assistance of Baidu map API, all of land use variables are the aggregated data collected within 500m PCA, which can be used in the model directly.For the raw population data with Tiff format, they are processed and aggregated within 500m buffers with ArcGIS 10.2 . Figure 5 shows the overall procedures of population data processing with ArcGIS 10.2. The processing procedure contains 3 main steps described as follows.

1. Step 1:

Transform the coordinate information of metro station into points with shapefile format. Create a new point feature layer based on x and y coordinates defined in a source table and geographic coordinate system by using “Make XY Event Layer” tool. 2. Step 2:

Generate 500m buffers of metro stations. Copy features from the input layer to a feature class (shapefile format) by using “Copy features” tool as the layer above is temporary. Then Create buffer polygons around input features to 500 meters by using “buffer” tool. 3. Step 3:

Zonal statistical analysis. Aggregate the population within each buffer. Since the value of each grid represents the number of people per grid square at $8.33\times{{10}^{4}}$ degrees spatial resolution (approximately 100m at the equator), so the population within each buffer can be obtained by summing up the value of the grid falls into the buffer by using “Zonal Statistics as Table” tool.

5.1.2 Spatial autocorrelation test to variables

Before building geographically weighted models, analysis was performed to test whether the candidate variables were spatially autocorrelated. First of all, the spatial distribution of variable candidates is analyzed. The results (Figure 6) indicate that these independent variables do have certain spatial correlation and aggregation to some extent. For example, population is densely distributed in the southeast regions, covering Luohu district and Futian district. Areas with large betweenness centrality generally locate on the closed railway loop of the metro network. Besides, the openning time of stations of Shenzhen metro is associated with the whole line openning time, and the distance to city center is exactly calculated according to the geographical location information, so these two variables have significant spatial correlation.

With the overview of spatial distribution of independent variables, the test of spatial autocorrelation can further detect how strong spatial correlation of variables are, this will provide a theoretical basis for the feasibility of applying a GWR model. Moran’s I is a measure of spatial autocorrelation developed by Patrick Alfred Pierce Moran (1950) [47] . All of the dependent variables and candidate predictors had estimated Moran’s I values higher than the expected E(I) (Table 2), thus, a positive spatial autocorrelation existed . Moran scatter plot can reflect the spatial autocorrelation intuitively. The scatter plot has four quadrants. If the observed value falls to the first and third quadrants, it indicates that there is a strong positive spatial correlation. If it falls to the second and fourth quadrants, it indicates there is a strong negative spatial correlation. Figure 7 shows several variables’ Moran scatter plot.

According to Moran scatter plots (see Figure 7), it can be seen that Moran’s I values of all variables above are not close to 0, which indicates these variables are not randomly distributed in space, and mostly falls to the first and the third quadrants. It is shown that each variable is positively spatial correlated more or less, especially eight of independent variables: population, distance to city center, and days since opened, and the number of residence, restaurants, banks, bus stations, and hotels have strong spatial correlation as Moran’s I is greater than 0.3 (Cressie, 1993) [48]. The result also lays the foundation for the feasibility of follow-up study.

5.1.3 Local collinearity diagnostics

Since strong spatial correlation was found for the variables in the research, it is reasonable to build GWR models to analyze influencing factors on station ridership of Shenzhen Metro. Here, condition number that evaluated collinearity were 32.50 ($<100$), indicating the low degree of multicollinearity among the candidate independent variables. However, through a series of local collinearity diagnostics for the local regression coefficients of GWR model for ridership estimation at different time resolutions, we found there was severe collinearity in local regression coefficients according to local condition numbers of all stations greater than 1000 (see Figure 8) . Multicollinearity is of great concern when condition numbers are greater than 1000, and results become unstable in the presence of strong local collinearity. It just showed that the local regression coefficients were potentially collinear even if the underlying exogenous variables in the data generating process were uncorrelated [23]. Due to the collinearity issue in the GWR model for our case, it is proper to implement Ada-GWL model to the data.

5.1.4 Model implementation

A series of models icluding OLS, GWR, GWL, and the proposed Ada-GWL were built in R software[41] by using the “GWmodel”, “gwrr”, “fields”, “lars” and “igraph packages [42, 43, 44, 45, 46].The algorithm for solving the Ada-GWL model is coded in R programming language and tested on Intel(R) Core(TM) i5-6500 CPU, 3.20 GHz processor. Windows 7 operating system is used for testing. For GWR, GWL, and Ada-GWL models, the selection of kernel function and bandwidth is necessary, which depends on the specific distribution of sample points in the study area, which has a strong impact on the descriptive and predictive power of GWR models [31], it is necessary to consider the distribution of sample points of Shenzhen metro network while selecting the kernel function. Since the stations of Shenzhen metro are relatively regularly spaced, so we adopted Gaussian kernel function to calculate the weight matrix and determined the bandwidth $b$ based on CV.

5.2 Results analysis

5.2.1 Comparative analysis

In the following part, OLS model, GWR model (an spatial extension of OLS), Adapted GWR (Ada-GWR) model for metro network (i.e., GWR model calibrated with network-based distance metric), GWL model (GWR with LASSO enabling simultaneous coefficient penalization and model selection), Ada-GWL model for metro network (GWL model calibrated with network-based distance metric) are implemented to the Shenzhen metro dataset described earlier, the results of models for estimating ridership at different time periods are shown in Figure 9:

The accuracy of the estimated responses is measured by calculating RMSE. The RMSE is the square root of the mean of the squared deviations of the estimates from the true values, and should be small for accurate estimators. R-squared is a statistical measure that represents the proportion of the variance for a dependent variable that’s explained by independent variables. The results of fitting five models to the data provide the error values shown above. In general, the performance ranks of five methods are the same for the five models (different time periods), which are “Ada-GWL$>$GWL$>$Ada-GWR$>$GWR$>$OLS”. First of all, we can see the superiority of four local models (GWR, Ada-GWR, GWL, and Ada-GWL) over global model (OLS) in terms of estimation error of dependent variable and goodness of fit. Secondly, it is noted that the adapted models including both GWR and GWL performed better than original versions of GWR and GWL models, which proved the fitness of the adapted models to the metro network. Thirdly, Ada-GWL for metro network do substantially better than other three models at estimating the dependent variable. Therefore, we can conclude that Ada-GWL model for the metro network we proposed can estimate Shenzhen metro ridership at any time period more accurately.

5.2.2 Analysis of the local regression of Ada-GWL model

1. To see which station a certain variable has the most impact on

Focusing on the result of Ada-GWL model, the local regression coefficients’ distribution for each variable of all stations can be plotted as the bubble plots. Take the modelfor estimating average daily ridership of a whole weel as an example. The spatial distribution of local coefficients is shown in Figure 10.

Through understanding the spatial distribution of local coefficients (elasticities), it is possible to know how relations between the variables vary across space (estimated coefficients) and variables selection of all stations. In above figures, the bubble with the bold outline demonstrates the coefficient for the variable in the station equals to 0, and other bubbles’ colours identify the range of coefficients, the bubble with lighter colours mean the coefficient is larger, and vice versa. Take the population factor as an example, stations with large positive coefficients are mainly distributed in the centre, meaning that more trips per capita were expected in the central south area of the metro network, where commerce, administration and education are concentrated. Besides, we can note that for the factors of degree and hospital, there are a large number of stations with zero coefficient, so these factors are not so important factors for influencing most of Shenzhen metro stations’ ridership.

Besides, in order to verify the applicability of Ada-GWL model to different time resolutions. We apply the models to estimate the ridership of other time resolutions. With regards to the models for non-rush hour ridership, the spatial distribution of local coefficients is illustrated in Figure 11.

According to Figure 11, we can see that more stations with zero coefficient than those of average daily ridership, and the coefficient matrix is very sparse, indicating that for a single station, main factors influencing the hourly ridership at non-rush hour are not manifold.

2. To see which variable has the biggest impact on a certain station

The heat map is applied to demonstrate the coefficient distribution of all variables and all stations (scaled by column (station) at first). For example, for station “313” (i.e., Tianbei Station), population is the factor that has the biggest positive impact on its ridership compared with other factors. For station “t251” (i.e., Huangbeiling Station), the number of shopping malls in its PCA has the biggest positive impact on its ridership compared with other factors.

Furthermore, in order to present the impacts of different factors on different stations more intuitively, based on the coefficients of (scaled by column (station)) of the model for average daily ridership of a whole week, we conduct clustering (Figure13) to further illustrate the association among stations, so that we can see which groups of stations that all variables have similar effects on, that is, for a group of stations, some variables have the same great influence on ridership (coefficients are large) and some variables have the same slight influence on ridership (coefficients are small). K-means clustering is one of the simplest and popular unsupervised machine learning algorithms. The goal of K-means clustering is to find groups in the data, with the number of groups represented by the variable $K$. K-means algorithm works iteratively to assign each data point to one of $K$ groups based on the features that are provided. Data points are clustered based on feature similarity. Here, we choose the k-means clustering method in this study. We determine the number of groups $K$ as 4 through minimizing the difference of Within-Cluster Sum of Squares (WCSS).

According to Figure 13, the four clusters are coloured by four colours. Referring to the distribution of functional zones in Shenzhen 888Source: http://www.szgeoinfo.com:5001/msmap/flex/landuse.html , we can see that the stations of cluster 1(purple) are mainly located in employment and educational based land. Stations in cluster 2 (green) are mainly located in commercial land. Stations in cluster 3 (pink) mainly distribute in traffic hubs land, and cluster 4(blue) are mainly located in residential based areas. To analyze the rationality of the categories as shown in the cluster analysis and the functional characteristics of the catergories, the results are summarized in Table 3.

By this token, Ada-GWL model can offer clear categorizing of functional zones of Shenzhen metro stations belong to, which makes findings more explainable and instructive to metro transportation and urban planning.

Here, we also use heatmap to demonstrate the coefficient distribution of all variables and all stations (scaled by column (station) at first) for non-rush hour ridership (see Figure 14), also, based on the coefficients of stations of the model for non-rush hour ridership, we did clustering to illustrate the association among stations (see Figure 15).

According to Figure 15, stations belong to cluster 1 (green) are mainly located in residential based areas. Stations in cluster 2 (green) are mainly located in commercial land. Stations in cluster 3 (pink) mainly distribute in traffic hubs land, and cluster 4(blue) are mainly located in employment and educational based land. The result is slightly different from that for average daily ridership, indicating different time resolutions will make stations present different characteristics.

6 Conclusion and Impications for Metro Planning and Periphery Development

In summary, this paper proposed an Ada-GWL framework for modelling the ridership of metro systems, which considered the network connection intermedia when calculating the distance matrix instead of Euclidean-based distance, and also simultaneously performed regression-coefficient shrinkage and local model selection. We demonstrated the superiority of our framework through a real world case study of Shenzhen Metro Systems, and meanwhile, we analyzed the influencing factors of Shenzhen metro station-level ridership from a local perspective. Different from previous works, the proposed direct demand model based on Ada-GWL framework engaged local model selection and network-based distance to estimate the station-level ridership accurately.

The results of the case study showed that Ada-GWL model performed the best in terms of goodness-of-fit compared with OLS, the ordinary GWR model, Ada-GWR calibrated with network-based distance, and ordinary GWL model. The implementation of the proposed Ada-GWL model can provide a theoretical basis for interpreting the influencing factors of station-level ridership from a local perspective. First, through understanding the spatial distribution of local coefficients (elasticities), it was possible to know how relations between the variables vary across space (estimated coefficients) and variables selection of all stations, which was helpful to see which station an influencing factor has the most impact on. Second, the coefficient distribution of all variables and all stations made it possible to see which variable has the biggest impact on a given station. Third, through clustering analysis of the stations according to the regression coefficients, clusters’ functional characteristics were found to be in compliance with the facts of the functional land use policy of Shenzhen.

In general, the adapted version of GWL introduced in this paper not only extended the traditional GWR to deal with collinearity issues among regression coefficients and estimation of station-level ridership but also inspired metro planning and periphery development from a local perspective. First, Transit Oriented Development (TOD) planning is suggested to adjust measures to local conditions. Here, local conditions mean more than the geographical factors, and they also incorporate different local effects of influencing factors on metro ridership. Second, different functional zones need to adopt different strategies of diverting passenger flows. For example, in commercial developed areas, enhancing security when diverting passenger flows is vital, while improving efficiency is more important when diverting passenger flows in employment based areas. Third, for stations with unsaturated passenger flows, the attraction to passenger flows can be boosted in combination with the functionality of the region. These findings can also inspire the metro planning and periphery development of other cities.

Appendix A Appendix

A.1 Summary of literature review on direct demand models for ridership estimation

The related studies on direct demand models for ridership estimation are summarized as following.