APPLYING MULTIVARIATE GEOSTATISTICS FOR TRANSIT RIDERSHIP MODELING AT THE BUS STOP LEVEL

: Travel demand models have been developed and refined over the years to consider a characteristic normally found in travel data: spatial autocorrelation. Another important feature of travel demand data is its multivariate nature. However, regarding the public transportation demand, there is a lack of multivariate spatial models that consider the scarce nature of travel data, which generally are expensive to collect, and also need an appropriate level of detail. Thus, the main aim of this study was to estimate the Boarding variable along a bus line from the city of São Paulo - Brazil, by means of a multivariate geostatistical modeling at the bus stop level. As specific objectives, a comparative analysis conducted by applying Universal Kriging, Ordinary Kriging and Ordinary Least Squares Regression for the same travel demand variable was proposed. From goodness-of-fit measures, the results indicated that Geostatistics is a competitive tool comparing to classical modeling, emphasizing the multivariate interpolator Universal Kriging. Therefore, three main contributions can be highlighted: (1) the methodological advance of using a multivariate geostatistical approach, at the bus stop level, on public transportation demand modeling; (2) the benefits provided by the models regarding the land use and bus network planning; and (3) resource savings of field surveys for collecting travel data.


Introduction and Background
Increasing concern about the environment and a discussion about sustainability have strongly influenced public policies around the world.In Brazil, law 12,587/2012, known as the Urban Mobility Law, points out that non-motorized and public transportation modes should be which Simple Kriging results are compared to those of Ordinary Kriging (Daya and Bejari 2015;Taharin and Roslee 2017;Viswanathan et al. 2015), in which Ordinary Kriging is compared to Universal Kriging (Hiemstra et al. 2010;Kiš 2016;Liu et al. 2015;Mubarak et al. 2015;Nalder and Wein 1998;Wang and Zhu 2016), and in which the three techniques are simultaneously compared (Asa et al. 2012;Seo et al. 2015).In short, since UK includes explanatory variables in its formulation, it normally outperforms the other interpolators, especially when there is some largescale trend present in the interest variable structure.Afterwards, OK, which assumes that the interest variable mean is unknown and varies locally, demonstrates the best results compared to SK, whose mean is global, constant and known.
In spite of several comparative studies already developed, the conclusion reported in these studies is not consensual.In the aforementioned articles, the interpolators´ performance varied substantially according to the type of data under analysis.Regarding the travel demand, not many studies were observed that compare the performance of geostatistical interpolators.In the case study proposed by Shamo, Asa and Membah (2015), the interest variable (Annual Average Daily Traffic) refers only to rural highway segments, which does not offer, a priori, a contribution to the urban public transportation planning.Besides this, the authors themselves reinforced the idea that the best kriging technique and semivariogram can only be obtained from the structure present in the available information about the interest variable.
Regarding urban bus transportation planning, which is highly important to the supply and demand balance of the PT system, passenger flow along the bus lines is a valuable information and, often, hard to acquire.Marques and Pitombo (2021), Marques and Pitombo (2019) and Marques (2019) proved that Geostatistics, more specifically Ordinary Kriging, demonstrates an excellent potential in estimating the three variables, collected from a Boarding and Alighting counts survey, that express the passenger demand along a bus route.They are: Boardings and Alightings (number of passengers entering and leaving the bus line at each bus stop, respectively) and Loading (passenger volume inside the bus at each line segment contained between two consecutive bus stops).Since this survey demands high resources, the results found by those authors suggest that it is possible to perform the Boarding and Alighting counts only in some bus line segments and, by kriging, estimate, with relative accuracy, the demand variable for non-sampled bus stops and segments.This study, however, did not make any comparison between OK and other geostatistical interpolators to verify which one of them could best fit the passenger volume estimate along a public transport line.
It is worthwhile mentioning that the spatial modeling of public transportation passengers at the bus stop level and train, metro or bus station is the most detailed treatment that can be applied to PT network planning.Due to this, this approach is the most recent among the techniques that seek to program supply and understand transportation and land use relationships.In the scientific literature, several studies of this kind can be found, most at the station level (Blainey and Mulley 2013;Blainey and Preston 2010;Cardozo et al. 2012;Chakour and Eluru 2013;2016;Chiou, Jou and Yang 2015;Choi et al. 2012;Chow et al. 2006;Gutiérrez et al. 2011;Sun et al. 2016) and a few at the bus stop level (Chu 2004;Dill et al. 2013;Kerkman, Martens and Meurs 2015;Pulugurtha and Agurla 2012;Ryan and Frank 2009).However, due to the difficulty in acquiring the variables to be modeled (Boardings and Alightings), in the case of the bus stop level, to the best of the authors´ knowledge, these studies have still not provided a spatial approach of ridership until the present moment.Even in the station level cases, the studies retrieved basically focus on applying Geographically Weighted Regression and generalized linear models to ridership data.Only the station level study of Zhang and Wang (2014), which applies Universal Kriging to the Boarding Bulletin of Geodetic Sciences.27(2): e2021019, 2021.
variable, was found so far, meaning that approaches based on multivariate Geostatistics at the bus stop level were not yet observed.
Thus, the aim of this study is to estimate a public transportation demand variable, along a bus line, by means of a multivariate geostatistical modeling at the bus stop level.As specific objectives, a comparative analysis conducted by applying Universal Kriging, Ordinary Kriging and Ordinary Least Squares Regression for the same variable under analysis is proposed.
Finally, the following main research gaps associated to this study can be enumerated: (1) Multivariate modeling of public transportation demand at the bus stop level by means of a geostatistical approach; (2) Lack of spatial approaches of transit ridership at the bus stop level; (3) The need for assessing the improvement, in goodness-of-fit measures, caused by the inclusion of explanatory variables to the geostatistical modeling; and (4) Passenger volume modeling at the bus stop level as they are the most appropriate elements for performing this analysis.This article contains 5 sections, including this introduction.The next section summarizes the few studies that perform ridership modeling at the bus stop level.Section 3 introduces the materials used in the case study and the method applied to them.Then, the results, as well as discussions about them, are presented in Section 4. Lastly, Section 5 draws the conclusions and also proposes suggestions for future research.

Ridership models at the bus stop level
While the traditional transportation planning (Ortúzar and Willumsen 2011) is done by means of Traffic Analysis Zones and continues as the most popular method for mobility diagnosis and solution proposal, Cervero (2006) argues that ridership modeling at the local level can provide demand estimates quickly and economically.Moreover, in spite of a regional approach, which uses averaged values of data for each Traffic Analysis Zone, boarding and alighting modeling per bus stop, train, metro or bus station can capture the effect of transit-oriented development on public transport demand, i.e., the influence of built environment variables on transit usage.
From smart card data, boarding and alighting per train or metro station are readily available.On the other hand, bus ridership at the stop level is not easy to collect.Concerning this, cities often depend on expensive surveys, such as boarding and alighting surveys, or automatic counters, which are not widely popularized yet.It may be possible to obtain boarding and alighting per bus stop from smart card data and GPS information, but some assumptions have to be made that affect the accuracy of the results, especially in the case of Alighting.Therefore, boarding and alighting surveys remain the only way to collect ridership at the bus stop level accurately.Table 1 shows studies that perform ridership modeling at the bus stop level.
From Table 1, it can be seen that the models used are limited to ordinary least squares regressions with logarithmic transformation to correct the asymmetry of the interest variable.Models for count data were also applied, but none of them present a spatial approach of bus ridership.Pulugurtha and Agurla (2012) tried to include spatial dependence of boarding through a weighting function, but only in a deterministic way.
Moreover, explanatory variables used in the boarding and alighting modeling can be divided into two groups: demand and supply variables.Demand independent variables intend to capture the effect of sociodemographic and land use features around bus stops on ridership.On the other hand, infrastructure and public transport service characteristics are addressed by the supply independent variables.In order to minimize the amount of information needed for the spatial modeling, the present study proposed a simple method for selecting the best predictors, as described in Section 3.

Materials and Method
The dataset used in this case study refers to the Boarding per bus stop data (number of passengers entering the bus line at each bus stop) over line 856R-10 from the city of São Paulo -Brazil.The results, from a Boarding and Alighting count survey performed along this line on a typical day Bulletin of Geodetic Sciences.27(2): e2021019, 2021.
In order to compose the group of explanatory variables to be included in Universal Kriging and Ordinary Least Squares Regression, both features from bus stops themselves and from their influence area were collected.From a catchment area of radius 400m centered in the bus stops (Zhao et al. 2003), the following variables were calculated: population (inhabitants) and population density (inhabitants per hectare), based on the 2017 Origin and Destination Survey (Metrô 2019) shapefile, which is given in Traffic Analysis Zones; and averaged values of household income and car ownership, female (%), population with no complete higher education (%), households with no private vehicles (%), percent of people aged up to 14, up to 17, aged between 18 and 22, 18 and 29, 18 and 39 and above 60 years old.These data were obtained from the sampled households of the 2017 O/D Survey that were within the catchment area; area, in hectares, of the 16 predominant land use classes according to the shapefile of predominant land use in 2016 (GeoSampa), which is disaggregated at the block level; and number of roads and intersections, length (meters) and road density (meters per hectare) inside each catchment area, based on the São Paulo road system (Open Street Map) shapefile.The number of points of interest (POI), also given by OSM shapefile, inside each influence area, was also considered.Overlapping catchment areas were prevented by using Thiessen polygons, similar to the method adopted by Zhang and Wang (2014) and Sun et al. (2016), in a GIS environment.
Besides the road system variables collected from Open Street Map, other indicators were adopted as a proxy of accessibility as well.Together with the Boarding/Alighting count survey results, SPTrans also made the General Transit Feed Specification (GTFS) data, from the São Paulo PT network, available.Knowing the code of the 57 bus stops covered by line 856R-10, the following was calculated from GTFS data: the number of bus lines that passed by each of these stops, and the average frequency of those lines; Euclidean and network distance between each bus stop and the nearest bus terminal, nearest metro station and nearest train station.Two intermodal proximity measures considering the shortest Euclidean and network distance between each bus stop and the nearest metro or train station were also included.While Euclidean distance is based on a straight line, network distance is calculated along the road system.These distance measures were obtained from the 57 bus stop shapefiles along with the São Paulo bus terminals, metro stations and train stations shapefiles, and Open Street Map road system.Versions of the populational, road system and accessibility variables, transformed by the natural logarithm, were also considered, and, in the cases where the raw data contained zeros, it was added to 1 before applying the transformation (Bartlett 1947).In order to include only the attributes encompassed by the bus stops´ influence area, the attributes of the original shapefiles went through an aerial interpolation.As stated in Table 1, the data collected for the modeling procedure covers both supply and demand independent variables.
Afterwards, dependent and independent variables were selected using a joint analysis of linear correlation and spatial autocorrelation.In order to choose the variable of interest, the Moran index (Moran 1948) was calculated for the Boarding and Alighting data in the six time bands mentioned above.After that, the degree of association between the cases with the highest and statistically significant values of Moran's index and all explanatory variables was tested by the Pearson linear correlation coefficient (R).In order to eliminate multicollinearity, at this stage, the Bulletin of Geodetic Sciences.27(2): e2021019, 2021.
R value between two potential predictors was limited to 0.60.Therefore, when a pair of independent variables had a high correlation with the variable of interest, but R with each other above 0.60, the variable with the least correlation with the dependent variable was discarded.This threshold was considered acceptable to avoid the omitted variable bias as well, since a pair of highly correlated variables does not always represent a cause-effect relationship.Other criteria for choosing dependent and independent variables were: expected correlation signal and presence of independent variables from both supply and demand groups.Thus, the number of Boardings, transformed by the natural logarithm in the 5th time band, also known as Night Peak (NP, from 20h to 23h59), was chosen as the dependent variable.As potential predictors, the following variables were kept: population, number of POIs, number of road intersections, road length, number of other bus lines, mean household income and average frequency of other bus lines in the same time band as Boardings, all transformed by natural logarithm; also, population with no complete higher education (%), residential, commercial and services area (ha), and network distance, in meters, between each bus stop and the nearest metro station, were considered.
The modeling step started by initially calibrating a linear regression model.To select the best predictors among those considered, a stepwise method was applied, in which only three independent variables remained.Regarding the modeling area, in general, there is a trade-off between the prediction power of the technique and the number of explanatory variables used in the model, whose data source might be hard to access.The desirable scenario is to have a minimum number of explanatory variables (that are preferably easy to acquire) associated to a satisfactory performance of the model.Based on this, the following procedure was adopted: initially, a simple linear regression model was calibrated with each one of the three explanatory variables, separately; then, three linear regressions were estimated using two predictor combinations; afterwards, a third model considering the three variables as predictors was generated.This approach was repeated in the geostatistical modeling by means of UK as this estimator also includes explanatory variables in its formulation.The purpose of this analysis was to verify whether the models with the least explanatory variable are also competitive in terms of minimizing errors between real and estimated values, and how much the spatial approach improves bus ridership estimates compared to traditional linear regression.
All linear regression models were calibrated using the Ordinary Least Squares method (Yan and Su 2009).Considering only the cases in which all predictors were statistically significant in linear regression (p < 0.10), the geostatistical modeling steps were performed.They are: (1) Empirical semivariogram calculation and model fitting; (2) Cross validation; and (3) Estimation by OK and UK.
The semivariogram (ℎ), or variogram 2(ℎ), is the main graphical tool of Geostatistics as it visualizes the spatial structure of the variable under analysis.The calculation of the empirical, or experimental, semivariogram is given by Equation (1) (Cressie 1993;Matheron 1971).
(1) (): value of the Regionalized Variable Z in the sampled geographical position ; : number of pairs situated at distance ℎ.
Equation 1 refers to Ordinary Kriging, in which the semivariogram is calculated straight based on the RV information.Concerning UK, this calculation is applied to the residual term, in which a spatial structure is assumed.Then, a theoretical model is adjusted to the empirical semivariogram values.The process of fitting a well-defined function to the empirical semivariogram points consists of obtaining three main parameters, the nugget effect, partial sill and range, from a preestablished method (Cressie 1993).In the present case study, geostatistical modeling was performed by means of the three main theoretical semivariogram models: Exponential (Exp), Gaussian (Gau) and Spherical (Sph) (Olea 2006), in order to verify if one of them demonstrates a much better adjustment compared to the others.
The process of kriging a Regionalized Variable basically consists of obtaining the optimum weights for the linear combination of weights and neighboring values that results in a continuous surface of estimated points, which also covers the non-sampled locations.The kriging estimator is given by Equation (2) (Cressie 1993;Matheron 1971). (2) * ( 0 ): estimated value of Regionalized Variable at the geographic position  0 ;   : optimum weight assigned by kriging to the neighbor  value.
Although both OK and UK are linear combinations, the first one assumes a constant and local, but unknown mean () of the dependent variable observations (Equation ( 3)), while the latter relaxes this assumption by considering the presence of a large-scale trend over the response variable structure (Equation ( 4)).

𝑍 = 𝜇 + 𝜀
(3) in which  is the error term of the model,   represents the explanatory variables, and  +1 expresses the linear function parameters to be calibrated.Thus, Universal Kriging assumes that the Regionalized Variable values are affected not only by their neighbors (small range variation), but also that there is a systematic component in their structure, caused by the influence of the built environment around the treatment elements, which are, in this case, the bus stops.Besides this, UK allows this large-scale variation to be modeled through the inclusion of explanatory variables to the kriging estimator.Thus, instead of considering the errors completely as white noise, it is assumed that the RV spatial structure is present in the residual term oscillation, where the semivariogram function is calculated (Cressie 1993).
Ordinary Kriging weights   are obtained from a matrix operation, represented in Equation ( 5).
The matrix on the left corresponds to the theoretical semivariance between sample points []; vector [] in the middle contains the kriging weights; and the vector on the right expresses the theoretical semivariance between the sample points and the point to be estimated [𝑀] .Therefore, OK weights are calculated according to Equation ( 6) for each point to be estimated.
On the other hand, Universal Kriging formulation deals with parameters in linear function, which is similar to classical regression, and residual semivariogram.Therefore, its calibration process is complex and must be performed in an iterative way.First, the linear model is calibrated and, after the residual term is calculated, the nugget effect, partial sill and range are obtained.Other values for these parameters, nearby the original ones, are tested until there is some convergence to an optimum error between the observed and estimated value criteria (Cressie 1993; Selby and Kockelman 2013; Zhang and Wang 2014).In short, UK estimates are given by Equation ( 7).
Where  0 is the matrix of explanatory variable observations of point  0 ,  is the vector of linear parameter estimates,   0 represents the vector of estimated covariances between sample points and point  0 , while   expresses the matrix of estimated covariances between sample points.It is worth remembering that covariance () and semivariogram () functions are related according to Equation (8).
Where  0 and  1 stand out, respectively, for the nugget effect and partial sill parameters from the theoretical semivariogram.
Concerning geostatistical estimates, cross validation is performed by the leave-one-out method (Cressie 1993).This technique consists of removing the database points one by one and calculating their value from the remaining points and theoretical semivariogram parameters (and also the linear function, when it refers to UK).Therefore, from the observed value at the points and respective estimated value, several goodness-of-fit measures can be established to assess the performance of the applied spatial statistics tool.Regarding the linear regression, the estimate considered in this study was the number of Boarding predicted by the model equation.Thus, some of the goodness-of-fit measures suggested by Hollander and Liu (2008) were calculated, which are: Mean Absolute Error (MAE), Root Mean Square Error (RMSE), Mean Absolute Percentage Error (MAPE) and Pearson linear correlation coefficient between the observed and predicted values (R).
The cited goodness-of-fit measures were applied to the results of each estimate and, hence, it was possible to assess and compare the accuracy of results found from such techniques, and to select those that demonstrated the best performance.In the UK cases, results from the semivariogram that provided the smallest errors were selected to compare them with the respective linear regression estimates.The computational resources that gave support to the method stages were: ArcGIS 10.1, QGIS 3.0.3and GRASS GIS 7.4.0(Bundala, Bergenheim and Metz 2014) to collect the potential predictors; GeoDa (Anselin 2004

Results and Discussion
Figure 1 shows thematic maps for the dependent variable, Boardings, i.e., the number of passengers entering each bus stop on a typical day (Tuesday) in 2017, in the aggregated set of bus trips made from 20:00 to 23:59; and for the three explanatory variables selected by the stepwise method.They are: (1) natural logarithm of population (lnpop); (2) residential, commercial and services area, in hectares (res_com_serv area); and (3) network distance, in meters, between each bus stop and the nearest metro station (metrodist_net).As lnpop and res_com_serv area belong to the demand variable group, and metrodist_net to the supply one, this result was deemed satisfactory.
As expected, bus stops located at regions with more inhabitants tend to have a higher number of Boardings.This pattern can also be noted in the case of residential, commercial and service area, meaning that the higher the land use mixture, the higher Boardings will be.Pearson's correlation coefficient between ln_boarding and ln_pop and between ln_boarding and res_com_serv area was, respectively, 0.68 and 0.45.On the other hand, despite some bus stops located near metro stations are showing less passenger flow, there are many points nearer metro stations that do present a high number of Boardings.This relationship resulted in a R value of -0.26 between ln_boarding and metrostation_net.Thus, it can be stated that most 856R-10 line users, in the period from 20:00 to 23:59, come from metro lines, probably returning from work to home.
Figure 1 also reveals that the number of Boardings per bus stop in line 856R-10 shows, in general, five volume peaks: the first one is next to the beginning of the route, the second and third are halfway, and the last two are near the end of the line.Such peaks interlay with lower passenger flow points, starting at the first bus stops of the line, which present a reduced number of Boardings.This pattern resulted in a Moran's index of about 0.26, which increased to 0.48 with the logarithmic transformation.In both cases, the index value was statistically significant (pseudo p-value < 0.05), proving the presence of spatial dependence in Boardings per bus stop data.
Descriptive statistics of dependent and independent variables are presented in Table 2. Travel demand variables, in general, are given as count data and show asymmetry very often.Thus, their relationship with explanatory variables may not be linear.In this case, logarithmic transformations contribute to linearizing the model equation, addressing the real nature of the data and, hence, improving results.
As shown in Table 2, mean and median measures for ln_boardings and ln_pop are similar, given their normality.Standard deviation for all variables, as well as minimum and maximum values, reveal the presence of a wide range of values, meaning the inclusion of more diversified data in the modeling, thus making it possible to use the models to estimate ridership for various conditions.Moreover, it is important to mention that Boardings and res_com_serv_area were zero for three and five bus stops, respectively.In the case of Boardings, some points at the end of the route did not have any passengers entering the bus line in the period from 20:00 to 23:59, probably because at this time most users are returning home from work and, hence, at the end of the line, most passengers are leaving the vehicle rather than entering it.Res_Com_Serv Area (ha)    The exponential model had a good fit to the experimental ln_boardings semivariogram.On the other hand, the experimental semivariogram seems to increase without bound as the lag distance increases, which could indicate the presence of a large-scale trend in the interest variable that is not being modeled (Oliver and Webster 2015).This might be the reason why Ordinary Kriging estimates are almost twice the observed values, given a Mean Absolute Percentage Error of 96%.It is worth remembering that, although geostatistical and traditional modeling were performed based on the values of Boardings with logarithmic transformation, goodness-of-fit measures were calculated using the estimates with inverse transformation, so they could be directly compared to the real values.
Table 4 shows the resulting parameters from Universal Kriging and linear regression.As for Ordinary Kriging, the best semivariogram model, i.e., the theoretical semivariogram that yielded the best goodness-of-fit measures, in all predictor combination cases, was the exponential one.Therefore, for the sake of brevity, Universal Kriging results shown in Table 4 correspond only to those from the exponential model.As expected, from the linear correlation analysis, population and residential, commercial and service area have a positive effect on ridership.Although the signal of metrodist_net is negative, it means that the closer a bus stop is from a metro station, the higher the number of Boardings at it will be.Moreover, it should be noted that all parameter estimates show little variation across the models (except for the intercept in the second model), which suggests that some factors, such as multicollinearity, that could cause misunderstanding in the coefficient's values, are not present.
Based on statistical significance, one can assume that the order of importance of predictors used might be: ln_pop, res_com_serv area and metrodist_net, which was also the sequence of predictors entering in the stepwise selection method.The percentage of the nugget effect in relation to the sill increased compared to the univariate case.In spite of that, in two of the five models, this parameter remains below 50%.According to Cambardella et al. (1994), variables with nugget-to-sill ratio of 25% up to 75% can still be considered as spatially dependent, in a moderate way.Conversely, range was significantly reduced, showing values from 1.2km to 2.2km, approximately.
Table 5 presents the goodness-of-fit measures applied to models shown in Table 4. Ordinary Kriging results, based on exponential semivariogram, are also displayed.Based on the goodness-of-fit measures, Universal Kriging models can be ranked, from the best to the worst, as follows: 1.1, 2.2, 3, 1.2 and 2.1.The best models for linear regression, in turn, were 3 and 2.1, followed by 1.1, 2.2 and 1.2.Comparing all eleven models simultaneously, UK estimates outperformed all other models, meaning that even the UK cases with only one or two predictors showed better results than linear regression with three predictors.Ordinary Kriging, which is a univariate technique, presented a MAE and RMSE lower than those of linear regression with res_com_serv area as the predictor.
Although models with more predictors may better explain the variance of interest variable, estimates can show no or little improvement when a new explanatory variable is added to the model, even a statistically significant one.The best results, from both Universal Kriging and linear regression, are highlighted in bold in Table 5.In the case of Universal Kriging, the model with only ln_pop as the predictor yielded the best estimates, while for linear regression, the best results are C.
those from models 2.1, which use ln_pop and res_com_serv area, and 3, which uses all three predictors.
The reason for that could be the fact that when multiple predictors are added to the linear combination part of UK, spatial structure of residuals starts to get blurred.As shown in Table 4, the nugget effect of cases 2.1 and 3 are the highest ones, corresponding to 70% of sill, approximately.However, even in these cases, estimates can still be improved through geostatistical modeling as Universal Kriging do not overlook the remaining spatial dependence on residuals.
Table 5 also proves that kriging estimates can, in fact, be improved by including explanatory variables in geostatistical modeling.Comparing Ordinary Kriging results with those of the UK last ranked case (model 2.1), there is a reduction in MAE, RMSE and MAPE of about 27%, 26% and 50%, respectively, while R increased 131%.Considering the best model of UK (1.1), these numbers increase to 37%, 36%, 56% and 170%, respectively.Moreover, ridership estimates can also be significantly improved by geostatistical modeling compared to linear regression: the most subtle improvements were for model 2.1, which showed reductions of 18%, 12% and 23% in MAE, RMSE and MAPE, respectively, and an increase of 17% in R. On the other hand, MAE and RMSE reduced 34% and 31%, respectively, in model 1.2, and R increase reached 128%.The best MAPE improvement corresponded, in turn, to model 1.1, with a reduction of 42%.These results indicate that not only geostatistical modeling can provide the best ridership estimates, but also that improvements will depend on what predictors are being used.
Finally, linear regression models from Table 1 exhibited the following adjusted coefficients of determination (adjusted R²): 0.328 and 0.330 (Ryan and Frank 2009), 0.69, 0.62 and 0.53 (Dill et al. 2013), and 0.772 and 0.762 (Kerkman, Martens and Meurs 2015).Meanwhile, adjusted R² for linear regression models in Table 5 was: (1.1) 0.453, (1.2) 0.188, (2.1) 0.545, (2.2) 0.518, and (3) 0.572.It should be noted that despite using much less information, some linear regression results obtained in the present study, which were outperformed by UK, are similar or slightly better than the first two, which suggests that the three predictors used were correctly specified, as they can explain a significant part of the ridership variance, show little variation when a new predictor is added to the model, and are statistically significant.
In order to provide a disaggregated analysis of errors and allow a comparison between models, Figure 3 shows maps of error ratios for Ordinary Kriging (a); Linear regression and Universal Kriging, both with all predictors, which was considered the best result of linear regression (b and c, respectively); Linear regression and Universal Kriging, both with the ln_pop as the predictor, which is the best result of UK (d and e, respectively).
Bulletin of Geodetic Sciences.27(2): e2021019, 2021.Three bus stops had an observed Boarding value equal to zero.Therefore, the error ratio could not be calculated for these cases, which is the reason why they do not appear in Figure 3. From the minimum and maximum error ratios, as well as the limits for each error group and the amount of bus stops in each group, the following conclusion can be drawn: the best estimates come from UK with ln_pop, then UK with all the predictors, followed by linear regression with three predictors, linear regression with ln_pop, and lastly Ordinary Kriging.
Despite the fact that Ordinary Kriging showed some very high errors, a detailed analysis of percentiles reveals that OK and linear regression with all predictors had the same amount of bus stops with an error ratio between -30% and 30%, approximately, which corresponds to 37% of the total data.Linear regression with ln_pop, UK with three predictors and UK with ln_pop showed, respectively, 34%, 45% and 50% of bus stops with an error rate between -30% to 30%, which was considered a satisfactory range of error.
As Ordinary Kriging assumes the interest variable mean is a constant, OK modeling of variables that present a wide range of variation usually yields high errors.Conversely, the same amount of bus stops showed error ratios ranging from -30% and 30% in both OK, which is a univariate technique, and linear regression with all predictors.On the other hand, as it does not include any explanatory variable, OK can only be applied to short-term public transportation planning, in which all built environment and transportation system variables are assumed to remain constant.
Following the bus stop sequence from top to bottom in Figure 3, extreme error ratios occurred at bus stops 32, in all cases, and 42, in Ordinary Kriging estimates.The main reason for that might be the size of catchment areas devoted to these points, which are the smallest ones due to high proximity to neighboring bus stops.This problem could be solved by running an alternative modeling in which all catchment areas would have the same size, overlapping each other, and then include some explanatory variable that could control the occurrence of competitive bus stops, as performed by Kerkman, Martens and Meurs (2015).

Conclusions and Final Remarks
Public Transportation plays an important role in the sustainable development of cities and social inclusion.In order to promote the proper functioning of this system, travel demand models have been developed and refined over the years, seeking to consider a characteristic normally found in travel data: spatial autocorrelation.Another important feature of travel demand data is its multivariate nature.However, regarding the bus transit demand, there is a lack of multivariate spatial models that consider the scarce nature of travel data, which are expensive to collect, and also need an appropriate level of detail.Thus, the main aim of this study was to estimate the Boarding variable along a bus line from the city of Sao Paulo -Brazil, by means of a multivariate geostatistical modeling at the bus stop level.As specific objectives, a comparative analysis conducted by applying Universal Kriging, Ordinary Kriging and Ordinary Least Squares Regression to the same travel demand variable was proposed.
In general, results showed that the inclusion of explanatory variables to the kriging estimator contributes, in fact, to increasing the prediction power of the technique.However, the performance of the models with only one predictor did not follow the same pattern in both geostatistical and traditional modeling.This reinforces the opportunity to investigate what would be the best predictors to be used in transportation demand spatial approaches to avoid those that would not bring significant improvements, but whose acquisition would require additional costs.
Results also suggested that Ordinary Kriging, which does not require additional information about explanatory variables, can be competitive to linear regression with only one predictor.This comes, probably, from the fact that OK already considers the spatial autocorrelation present in the Boarding variable.However, this interpolator has the disadvantage of not being able, from only the available data about the interest variable, to predict its values for other scenarios, including future ones.This capacity is observed only in Universal Kriging and Linear Regression.In addition, estimates from all geostatistical cases revealed a better adjustment of exponential semivariograms to Boarding data.
Although the results from Universal Kriging may suggest that the lower the number of predictors, the better the estimates will be, we do not encourage ignoring additional information when it is available and contributes, in fact, to explaining interest variables.However, when detailed data is not provided, which is the case of various cities, in development countries, especially the small and medium-sized ones, spatial models with little information available could also yield good estimates.In general, results showed that traditional modeling can always be improved by geostatistical multivariate interpolators, not only in cases where there is only one predictor, but also when a large amount of information is used.Best results from UK showed 50% of bus stops with error between -30% and 30%.In turn, regarding the best results from linear regression, only 37% of bus stops had errors within this range.
Therefore, three main contributions are highlighted: the methodological advance of using a detailed geostatistical approach, the bus stop level, on bus ridership modeling; the benefits provided by the models regarding the land use and bus network planning; and resource savings of field surveys for collecting travel data.In order to compare the achieved results with another spatial method that, similar to the geostatistical interpolators, also creates a surface of estimated values, Geographically Weighted Regression is recommended for the same dataset used in the present study.Nevertheless, it is opportune to compare the OK and UK results to those of generalized linear models (Poisson and Negative Binomial regressions), which consider the Sciences.27(2): e2021019, 2021.

Figure 1 :
Figure 1: Patterns of (from top to bottom) Boardings; Population; Residential, commercial and services area; and distance to the nearest metro station along the bus line 856R-10

Figure 3 :
Figure 3: Error ratios of (from left to right) Ordinary Kriging (a), Linear Regression and Universal Kriging with all three predictors (b and c), and Linear Regression and Universal Kriging with ln_pop (d and e).

Table 3 .
In spite of the low percentages of nugget effect relative to the sill (nugget effect plus partial sill), goodness-of-fit measures are not quite satisfactory.Comparing the three theoretical models, the exponential one provided the best estimates.Experimental semivariogram for ln_boarding and the fitted exponential model are shown in Figure2.

Table 3 :
Ordinary Kriging results * statistically significant at the 0.05 level (one-tailed).MAE, RMSE, MAPE and R are, respectively, Mean Absolute Error, Root Mean Square Error, Mean Absolute Percentage Error and Pearson Linear Correlation Coefficient between predicted and observed values.