{'fraction_non_alphanumeric': 0.029364846827489664, 'fraction_numerical': 0.006673828824429469, 'mean_word_length': 4.724562057398434, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interpretable Machine Learning Models for Modal Split Prediction in Transportation Systems

Modal split prediction in transportation networks has the potential to support network operators in managing traffic congestion and improving transit service reliability. We focus on the problem of hourly prediction of the fraction of travelers choosing one mode of transportation over another using high-dimensional travel time data. We use logistic regression as base model and employ various regularization techniques for variable selection to prevent overfitting and resolve multicollinearity issues. Importantly, we interpret the prediction accuracy results with respect to the inherent variability of modal splits and travelers' aggregate responsiveness to changes in travel time. By visualizing model parameters, we conclude that the subset of segments found important for predictive accuracy changes from hour-to-hour and include segments that are topologically central and/or highly congested. We apply our approach to the San Francisco Bay Area freeway and rapid transit network and demonstrate superior prediction accuracy and interpretability of our method compared to pre-specified variable selection methods.


I. INTRODUCTION
Multi-modal transportation system operations depend crucially on the distribution of aggregate travel demand across travel modes and the way in which this modal split varies over time. Accurate modal split prediction at a fine temporal scale can support transportation authorities in better managing congestion and improving service reliability. This information can be especially valuable in highly uncertain situations when traveler preferences shift significantly from one mode to another. For example, accurate modal split prediction in times when travel demand changes from offpeak to peak periods can help transit operators in scheduling frequent service in high-demand corridors, which in turn can induce further demand shift to transit and reduce freeway congestion during peak periods (see Fig. 1). Additionally, improved understanding of how the modal split is shaped by spatio-temporal variations in travel times can provide guidance on interventions such as control [10] or tolling [15] of specific freeway segments.
In this paper, we propose an interpretable machine learning (ML) approach to (1) predict the modal split in transportation systems on an hourly basis using transportation network travel times; and (2) identify the subset of network segments that are most critical for predictive accuracy. We adopt logistic regression as a base model and investigate the use of different variable selection (regularization) methods to tackle the high dimensionality of travel time data and identify segments that have the most impact in improving the predictive accuracy. We instantiate our approach on data from the San Francisco (SF) Bay Area, namely regions in the East Bay and SF proper, which are served by the Bay Area Rapid Transit (BART) and a congestion-prone freeway system. Our approach directly predicts the hourly aggregate modal split based on raw, segment-level travel time data from a large number of regions where both transit and driving choices are available. Essentially, the logistic hypothesis captures mode choice behavior at an aggregate level while the associated variable selection and regularization methods control the variance of the predicted modal split. This makes our approach well-suited for generating modal split predictions on an hourly basis, in comparison to the classical approach in the travel demand literature that typically relies on survey data to build a discrete-choice model of traveler preferences [2] and estimates demand patterns through extensive microsimulations [1], [3].
Importantly, the resulting models are interpretable in that the selected features correspond to segments impacting modal split predictions due to their location in the network as well as their influence on the overall congestion dynamics. Indeed, previous work has used specific criteria to pre-select segments in freeway networks based on their location or impact on congestion [14], [9]. In contrast, we leverage recent developments in regularized regression [6] and obtain superior predictive performance in comparison to pre-specified variable selection methods that only select segments based on their location in the network (betweenness centrality) or local congestion measures (segment-level speeds). Our computational results demonstrate the effectiveness of variable selection methods in tackling high-dimensional travel time data and identifying a relatively small subset of segments that are most critical for out-of-sample predictive accuracy.
Our paper contributes to the emerging literature that explore the use of ML models such as classification and regression trees (CART) and neural networks for discretechoice modeling [11], [8], [12], [4], [16]; for a more thorough discussion of applications of ML towards choice modeling, we refer readers to [13]. These approaches treat choice modeling as a classification problem and typically focus on out-of-sample predictive accuracy as the main performance measure. Thus, they naturally leverage the ability of ML models to arrive at the desired bias-variance tradeoff by minimizing overfitting. We go a step further by focusing on how variable selection methods can be used to not only build sparse models (similar to [4], [12]), but also to generate new insights by properly interpreting them.
The rest of the paper is organized as follows: In Sec. II, we present our modal split prediction model and discuss the method of regularized maximum likelihood estimation as well as the model selection and evaluation process. In Sec. III, we introduce pre-specified variable selection methods that serve as benchmarks for evaluating our approach and detail how the accuracy of our prediction model can be interpreted in terms of variability of modal split and responsiveness to changes in travel time. Finally in Sec. IV, we apply our modeling and interpretation framework to a case study in the SF Bay Area.

A. Prediction Model
In this section, we present our model for predicting in time interval t the fraction of driving flow among the total flows of driving and transit based on travel time data for segments in a multi-modal transportation network. In our case study, each t corresponds to a one-hour interval, although t can vary in practice depending on the time granularity of available data. We estimate separate models for different times of day since travelers and their demand patterns can vary considerably throughout the day. For the observation corresponding to day i, we denote by q 1 i and q 0 i the driving flow and the transit ridership from region A to region B, respectively. Here, regions A and B correspond respectively to the sets of origin and destination locations that are served by both freeway and transit modes (Fig. 2). We denote by y i ∈ [0, 1] the driving fraction on day i: which is the quantity we aim to predict. Our multi-modal network is comprised of the freeway network and the transit network. We define segments in the freeway network (n driving segments in total) as sections between adjacent entrance or exit ramps and segments in the transit network (n transit segments in total) as sections between two adjacent stations. The total number of segments is n = n driving + n transit .
For a given time interval t and day i, we predict the driving fraction y i using the travel time of all segments in the network during time intervals t, (t − δ), . . . , (t − pδ), where δ > 0 denotes a unit time lag and p the maximum number of time lags incorporated in our model. Our prediction model uses travel time vectors in lagged time intervals because the flows of driving and transit counted in any time interval are induced by travelers' mode choice decisions that may have been made before their trip was counted by the traffic sensor or transit tap-in sensor. As such, some mode choices made in one time interval will influence the flow count for a subsequent time interval. The unit time lag variable δ depends on the granularity of the available travel time data, and p is chosen such that pδ is larger than the maximum travel time from any origin in region A to the sensors that record the flow between A and B in the network. In our case study, we choose p = 6 and δ = 10 minutes -i.e. one hour of travel time data in 10-minute intervals.
Thus, for each time interval t, our goal is to predictŷ i using the travel time vector ..,n is the vector of travel times for all n segments in time interval (t − kδ) on day i. We predict whereβ 0 andβ := (β 0 ,β δ , . . . ,β pδ ) are estimated regression coefficients and x i := (τ 0 i , τ δ i , . . . , τ pδ i ) denotes the combined vector of travel times. Here, σ(z) = 1 1+e −z is the logistic function that constraints the linear outputβ 0 +β x i to (0, 1). 1 Hence,ŷ i can alternatively be interpreted as the probability that a traveler in time interval t on day i chooses to drive rather than take public transit.
We can then estimateβ 0 ,β by maximizing the loglikelihood, where m is the total number of days (i.e. observations) in our dataset. This formulation is equivalent to estimating a logistic regression with m i=1 q 1 i + q 0 i total observations of mode choices such that day i contributes q 1 i driving and q 0 i transit observations, all corresponding to the same travel time vector x i .

B. Estimation Methods
We estimate the constantβ 0 ∈ R and coefficientŝ β ∈ R n(p+1) using the dataset X ∈ R m×n(p+1) , which contains m days of travel time vectors in time periods t, (t − δ), . . . , (t − pδ). For each day i ∈ {1, . . . , m}, the input data for interval t is the vector x i . In many practical settings, constructing a dataset with observations spanning longer durations (many months or years) is not only difficult but also unadvisable given the changing behavior of travelers and network operations (e.g. changes to tolling, transit schedules, road closures, etc.). For these reasons, it is typically the case that n(p + 1) m; in our case study, we have n(p + 1) = 1106 and m = 130. As such, estimating the regression function in (1) suffers from multicollinearity and overfitting. To resolve challenges associated with high dimensionality, we use regularized logistic regression methods for shrinkage and variable selection [6].
For all regularization methods, we standardize the dataset X to have a mean of 0 and standard deviation of 1 before estimation so that regularization is applied evenly across all variables [7].
1) Forward-stepwise selection is a greedy algorithm for selecting a subset of variables with which we estimate a logistic regression [7]. We initialize the algorithm with a constant model and iteratively incorporate the variable that, in combination with all previously selected variables, yields the lowest cross-validation error. This process continues until (i) there is no remaining variable such that adding that variable achieves a lower cross-validation error than that of the previous iteration, or (ii) adding a variable results in a rank-deficient data matrix with more variables than observations.
Forward-stepwise selection is computationally costly as it can require up to 2) Lasso regression controls variance by imposing a penalty on the 1 norm, or sum of absolute values, of the regression coefficients to yield a sparse coefficient estimate in which some coefficients are exactly zero [7]. The lasso estimator maximizes the objective function, where the hyperparameter λ > 0 governs the penalty imposed on the 1 norm of the coefficient estimates. For higher values of λ, lasso automatically selects a smaller subset of variables and sets the coefficients of the remaining variables to zero.
3) Ridge regression controls variance by imposing a penalty on the 2 norm, or sum of squared values, of the regression coefficients to yield coefficients estimates that are shrunk towards zero [7]. The ridge estimator maximizes the objective function, where the hyperparameter λ > 0 governs the penalty imposed on the 2 norm of the coefficient estimates. Higher values of λ yield more shrunken coefficient estimates. 2

C. Model Selection and Evaluation
We select the hyperparameter of each method (i.e. λ for ridge and lasso regression) and the set of predictors in forward-stepwise selection by minimizing the average outof-sample root mean squared error (RMSE) in a K-fold cross-validation. We randomly split the dataset into K folds I 1 , . . . , I K . For each fold I f where f = 1, . . . , K, we estimateβ 0 ,β using data in the remaining folds {1, . . . , m}\I f . We then calculate the out-of-sample root mean squared error (RMSE) for fold I f as The average cross-validation RMSE is an absolute measure of error that describes the standard deviation of out-of-sample prediction errors. Note that this metric has the same units as the driving fraction, i.e. the fraction of travelers choosing the freeway network. Practically, however, the driving fraction shows greater variation, and consequently decreased predictability, in some hours compared to others (see. Fig.  3b). To compare prediction error across hours, we consider the out-of-sample R 2 , which we compute for each fold I f as follows: Here, y is the sample mean of y i for i ∈ I f . The outof-sample R 2 normalizes the mean squared error by the sample variance of the driving fraction for each hour, which allows comparison of performance across hours. Unlike the in-sample R 2 , the out-of-sample R 2 is not lower bounded by 0 and instead ranges from −∞ to 1 with values closer to 1 indicating higher out-of-sample accuracy. For any given hyperparameter(s), we iterate this process over all folds and compute the average cross-validation RMSE and outof-sample R 2 as

III. INTERPRETATION FRAMEWORK A. Pre-specified variable selection
We introduce two approaches to variable selection as benchmarks to evaluate the regularized models described in Sec. II-B. In particular, we compute a metric, either betweenness centrality or average speed, for each segment and regress on a subset of segments that fall within a given range for each metric. Selecting these subsets and comparing their estimated coefficients to those from regularized approaches allows us to identify how topological (i.e. betweenness) and temporally varying demand-based (i.e. average speed) segment characteristics are selected in different time intervals.
1) Selection based on betweenness centrality: We define a betweenness centrality metric for each segment as the number of routes -sequences of segments -that connect an origin in region A to a destination in region B and pass through the given segment. From a topological viewpoint, a segment with higher betweenness centrality affects the travel time of more routes that pass through it. The key hypothesis here is that changes in travel time on segments with a high betweenness centrality will likely have a larger impact on the aggregate driving fraction relative to those with a lower centrality. Consequently, in building a logistic regression, one can select the travel time variables of segments with betweenness centralities above a certain threshold b.
2) Selection based on average speed: We compute the average speed of each segment as the ratio between the segment length and travel time. Segments with low average speeds are considered to be more congested, and consequently changes in travel time on these segments are likely to have a greater impact on the driving fraction y. Thus, an alternative to variable selection based on betweenness centrality is to select the travel time variables of segments with average speeds below a certain threshold v.
For both the betweenness and speed selection methods, we also select the travel time variables of segments observed in time (t − k b δ) and (t − k v δ), respectively. Here, the time lags k b and k v are also hyperparameters.
We estimate the logistic regression coefficients of the travel times on the segments selected based on betweenness centrality or average speed by maximizing the log-likelihood, (β 0 , β | q 1 , q 0 , X), where X contains travel times for the selected segments.

B. Degrees of freedom
For any estimated model, we report the effective degrees of freedom, df < n(p + 1), which can be viewed as a measure of model complexity. For methods that yield sparse coefficient estimates, df is the number of non-zero coefficients. For ridge regression, df can be computed as Trace(X(X X + λI) −1 X ), where I is the identity matrix. Note that df is closely related to the regularization hyperparameter of a given model as greater regularization yields lower values of df.
This quantity provides a behavioral interpretation -higher values of df suggest more travel time information is incorporated in the model prediction, i.e. the aggregate driving fraction can be better predicted with travel time information from a larger set of network segments. However, df is not a measure of responsiveness (which we discuss in the following subsection) as df is the result of estimating a predictive model with a limited number of observations and many correlated predictor variables.

C. Behavioral implications
From a behavioral viewpoint, one can expect that the changes in travel times of certain "critical" segments would impact travelers' preferences for the routes on which these segments lie. Consequently, for any given hour, such segments can have a significant impact on the modal split. We argue that our modeling approach is useful for identifying critical segments in each hour and helps us study how these segments impact the travelers' aggregate responsiveness to travel time changes. In Sec. IV, we show that this behavioral interpretation provides new insights about high variability and low responsiveness of the modal split (and hence difficulty in its prediction) for certain hours of the day. We now formally introduce these quantities below.
1) Variability: For any time interval t, we quantify the inherent variability of the mode share σ 2 in time t by the variance of the driving fraction, i.e.
Time intervals t with higher values of σ 2 correspond to more difficult prediction tasks.
2) Responsiveness: For any time interval t, any freeway segment j, and any time lag k, we denote the correlation ratio between the travel time of segment j in time (t − kδ) with the driving fraction y as Corr(y, τ kδ j ) ∈ (0, 1), which is given by We quantify the responsiveness of mode share to changes in travel time by Corr(y, τ kδ j ), or the correlation between the driving fraction y and the travel times τ kδ j averaged over all segments j and time lags k. The value of the average correlation ratio ρ is in (−1, 1). Negative values suggest that high travel times of freeway segments reduce the driving fraction, while positive values suggest that high driving times increase the driving fraction. Moreover, the average correlation value being closer to −1 indicates that the driving fraction is more responsive to the driving times of the segments.

IV. EMPIRICAL STUDY IN SAN FRANCISCO BAY AREA
The goal of our case study is to predict the fraction of driving flows during each one-hour time interval from 5am to 10pm in the SF Bay Area. Recall that our multimodal network is comprised of the freeway network and the BART network. The origin region A includes the East Bay, and the destination region B includes San Francisco proper and Daly City. The two regions are connected by the San Francisco -Oakland Bay Bridge freeway and the transbay BART tube. The shortest driving route from any origin in A to any destination in B passes through the Bay Bridge, and similarly any BART route that connects stations in the two areas passes through the transbay tube. 3 The network includes 158 freeway segments (determined by neighboring on-ramps and off-ramps) and 38 BART segments (determined by adjacent stations). Regions A and B as well as relevant freeways, loop detectors, and BART stations are depicted in Fig. 2.

A. Data
We use traffic flow data reported by the Caltrans Performance Measurement Systems (PeMS) [5] and transit ridership data reported by the BART authority on workdays between February and September of 2019. We identify and remove any days in which severe incidents have been reported in the California Highway Patrol (CHP) incidents report dataset to have caused lane closures in the two areas as well as any days in which BART operations are severely disrupted. 4 In total, we use data from m = 130 workdays.
Next, we describe how we measure the hourly driving flow, BART ridership, and average travel time of segments.
1) Driving flow. The PeMS data includes the hourly driving flow (i.e. the number of passing vehicles in each hour) recorded by loop detector sensors that are embedded under the freeways. Since any shortest route that connects regions A and B passes over the bridge, the total driving flow from A to B is the flow recorded by the loop detectors on the bridge. 2) BART ridership. The tap-in/tap-out systems at BART stations record the number of completed trips between each pair of stations in each hour. We compute the total ridership from A to B by summing over all station pairs with origins in A and destinations in B. 3) Average travel time of segments. The PeMS dataset records the average speed of vehicles that pass over each loop detector in every 5-minute interval. We compute the average driving time of each segment as the ratio between the segment length and the average speed recorded by the loop detectors located in that segment. 5 The operating time for each BART segment is not recorded by the BART authority, and thus cannot be used in our prediction. 6

Fig. 2: PeMS loop detectors (red) and BART stations (blue)
in the San Francisco Bay Area. Fig. 3a shows the variability of driving travel times from Richmond to the Bay Bridge along the north branch of our network (i.e. I-80W). The highest driving time occurs during morning rush hour and is three times the lowest driving time, which occurs at night. We also observe that the confidence interval of the driving time in each hour is quite wide, which suggests that travelers experience high variability in day-today driving times. Furthermore, Fig. 3b demonstrates that the driving fraction is lowest during morning rush hour and highest at night. We also observe within-day variation in the driving fraction. This suggests that the driving fraction induced by travelers' choices between driving and taking public transit might be sensitive to the changes of driving time in the network over the course of the day.

B. Results and discussion
The results from our case study are summarized in the following figures: (i) Fig. 4 compares the predictive accuracy of various models by reporting the final 5-fold cross validation RMSE for each variable selection method (see Sec. II-B and (a) Driving travel times on I-80W peak in morning rush.
(b) Driving fraction is lowest in morning rush hour. III-A) in each hour. For the sake of comparison, we also include results from applying a random forest approach -a popular nonparametric ML method [7]. (ii) Relative accuracy is reported in Fig. 5a as the average cross-validation out-ofsample R 2 results (see II-C) for each method in each hour. (iii) Model complexity is reported in Fig. 5b as the effective degrees of freedom (see Sec. III-B) selected by our 5-fold cross-validation for each method and each hour.
Based on these results, we can make the following observations: 1) Overall performance: All methods provide accurate prediction of driving fraction with a median RMSE below 0.03 across all hours except for the hours of 10am, 11am, and 10pm.
2) Comparison across hours: Our models achieve their highest predictive performance during morning rush hour (6am-9am) and early afternoon (12pm-4pm). The median RMSE and variance of RMSE across folds in these hours are both lower than in other hours.
One reason for superior predictive performance in morning rush hour is that travelers are more likely to respond to changes in driving time within these hours. Indeed, from Fig.  6b, we can see that 6am-9am are hours in which the driving fraction exhibits higher responsiveness (i.e. more negative correlation) towards the changes in driving time. Additionally, we can see from Fig. 5a that these hours achieve a higher out-of-sample R 2 and also select for high model complexity (i.e. the driving time of more segments are selected in the regression) in all methods except forward-stepwise selection (Fig. 5b). On the other hand, high prediction accuracy in the early afternoon (12pm-4pm) and evening (8pm-9pm) is a result of especially low variability rather than high responsiveness. Finally, we observe that, despite exhibiting moderate responsiveness, prediction accuracy in 10am-11am suffers as a result of high variability. In summary, observations suggest that prediction accuracy in a given hour is restricted by moderate variability unless high responsiveness is also exhibited.
3) Comparison across methods: In Fig. 4, we observe that the forward-stepwise selection approach achieves a lower RMSE than all other methods. This is likely due to the fact that forward-stepwise selection chooses from a larger set of combinations of predictors.
We also observe that the regularized methods generally achieve a lower RMSE than the two pre-specified variable selection methods (see Sec. III-A). This is due to the fact that these regularized methods address model overfitting and multicollinearity in a more flexible manner than variable selection methods that are purely based on betweenness centrality or average speed. However, the performance of either pre-specified variable selection method relative to the other varies across hours. In 8am-9am, selection based on betweenness centrality outperforms selection based on average speed, which coincides with the hours of highest travel times. In almost all other hours, however, selection based on average speed outperforms selection based on betweenness centrality. This might be due to the fact that the highway network is already largely in a congested regime in 8am-9am, and during this time incorporating segments that are more congested is likely less beneficial for prediction than incorporating more central segments. In fact, the heatmaps of model coefficients presented in Sec. IV-C demonstrate that regularized methods are able to capture both congestionbased and spatial characteristics when assigning coefficients to segments.

C. Spatial distribution of critical segments
We compute the weight of each segment j ∈ {1, . . . , n} as the sum of coefficients across all time lags, i.e. w j = p k=0β kδ j . We demonstrate the geospatial distribution of segment weights for ridge regression using heatmaps in Fig. 7.
Lower (more negative) weights indicate that an increase of driving time on the considered segment has a greater impact on reducing the overall driving fraction. Thus, segments marked in red are ones that have a high impact on the driving fraction and thus deemed "critical" for the prediction task. Practically, prior knowledge of these high impact segments can be useful to the transportation authority for effective management of driving and transit demand.
Coefficient heatmaps for all methods demonstrate a similar pattern of how the spatial distribution of high-impact segments changes across hours. In particular, we observe that segments on the northeast branch have a high impact on the modal split in most hours, while segments on the south and southeast branches demonstrate high criticality in 8am-11am.
Finally, we note that the set of high-impact segments identified by lasso, ridge, and forward-stepwise selection  includes segments with high betweenness centrality (i.e. the Bay Bridge and surrounding segments) and segments with low speed (i.e. on the north and northeast branches). This suggests that the superior performance of these regularization methods in comparison to pre-specified selection methods can be explained by their flexibility in selecting for combinations of segments that are topologically central and/or highly  congested. However, segments that are neither topologically central nor highly congested may also be selected; this may be explained by the sensitivity of demand along these segments to occasional congestion.

V. CONCLUSION
In this work, we propose an interpretable machine learning approach for predicting aggregate modal split for travelers in a multi-modal transportation network. We use a regularized logistic regression model, namely with 1 norm and 2 norm regularization as well as forward-stepwise variable selection. We discuss behavioral implications of our results by quantifying the inherent variability of the modal split as well as responsiveness of the modal split to changes in travel time. Finally, we apply our approach to predicting driving fraction in the SF Bay Area and discover that our models achieve their highest accuracy on predictions during: (1) morning rush hour (6am-9am) when travelers are most responsive to travel time, and (2) the early afternoon (12pm-4pm) when variability is lowest.
Visualizing our estimated model parameters as heatmaps allows us to identify critical regions in the freeway network at different hours. As such, our approach provides a class of accurate predictive models along with useful interpretations that may support transportation authorities in anticipating demand shifts in critical regions and making operational decisions accordingly.
Future work may focus on predicting aggregate modal split in a severe incident setting. Preliminary investigation suggests this task is more difficult than predicting aggregate  Fig. 7: Relative values of weights are colored on a spectrum from red to green, where red represents lower weight (higher criticality) values and green represents higher weight (lower criticality) values. Weights are normalized in each hour to convey the geospatial distribution of weights so the weight of any one segment cannot be compared across hours. modal split during nominal conditions. Accurate prediction in incident settings requires the construction of a larger incident event dataset as well as a different approach to analyzing the behavioral adjustments of travelers as they obtain incident information and reroute. Finally, designing resilient network operational decisions based on aggregate modal split predictions is another problem of significant interest.