{'fraction_non_alphanumeric': 0.022914057010173842, 'fraction_numerical': 0.006782560875011457, 'mean_word_length': 5.1168940384974775, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, '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}}
Decoding pedestrian and automated vehicle interactions using immersive virtual reality and interpretable deep learning

To ensure pedestrian friendly streets in the era of automated vehicles, reassessment of current policies, practices, design, rules and regulations of urban areas is of importance. This study investigates pedestrian crossing behaviour, as an important element of urban dynamics that is expected to be affected by the presence of automated vehicles. For this purpose, an interpretable machine learning framework is proposed to explore factors affecting pedestrians' wait time before crossing mid-block crosswalks in the presence of automated vehicles. To collect rich behavioural data, we developed a dynamic and immersive virtual reality experiment, with 180 participants from a heterogeneous population in 4 different locations in the Greater Toronto Area (GTA). Pedestrian wait time behaviour is then analyzed using a data-driven Cox Proportional Hazards (CPH) model, in which the linear combination of the covariates is replaced by a flexible non-linear deep neural network. The proposed model achieved a 5% improvement in goodness of fit, but more importantly, enabled us to incorporate a richer set of covariates. A game theoretic based interpretability method is used to understand the contribution of different covariates to the time pedestrians wait before crossing. Results show that the presence of automated vehicles on roads, wider lane widths, high density on roads, limited sight distance, and lack of walking habits are the main contributing factors to longer wait times. Our study suggested that, to move towards pedestrian-friendly urban areas, national level educational programs for children, enhanced safety measures for seniors, promotion of active modes of transportation, and revised traffic rules and regulations should be considered.


Introduction
With a technological revolution in mobility on the way, unprecedented changes in urban areas are expected to happen. Before switching to Automated Vehicles (AVs), a careful and detailed investigation is required to analyze their impacts on future mobility patterns, travel behaviour and street design. As a part of the human behaviour that is expected to change, interactions between different road users are of importance as they affect safety, efficiency, congestion, quality of service and public trust of transportation systems. In particular, the interaction that we seek to investigate in this study is that of pedestrians, as the most vulnerable road users, and vehicles. Recent instances of AV-pedestrian collisions, e.g. Uber's test AV fatal incident in Tempe, Arizona, and Navya SAS automated bus accident in Vienna, Austria, reveal the vital importance of such investigations. In an urban space dominated by ruleobeying automated vehicles that always stop for pedestrians even at mid-block unsignalized crosswalks, an emphasis on investigating this type of crossing is much needed (Millard-Ball, 2018). Currently, the interaction between drivers and pedestrians is mainly in the form of a silent agreement. Eye contact between the two agents, or head orientations and body movements of the pedestrians observed by the drivers, are some of the behaviours that establish this silent agreement. Drivers then continuously try to predict the next movements of pedestrians based on what they see, while pedestrians decide on their behavioural choices based on what they perceive from the vehicles and drivers.
The question arises, however, that how can we replicate the same interactions in an automated environment? In the absence of a driver, how can we make the automated vehicles socially-aware to be able to anticipate movements and decisions of pedestrians based on the contextual environment? Moreover, what factors, in terms of traffic parameters, rules and regulations, policies, design, and practices need to be taken into account in the transition towards automated urban environment?
A schematic representation of interactions between vehicles and crossing pedestrians is outlined in Figure 1. A pedestrian's behaviour in this context can be roughly abstracted into two parts: a. waiting on a sidewalk, during which a pedestrian intends to cross or decides to wait further until they feel more comfortable to cross, and b. crossing the street, which results in a pedestrian trajectory. The driver/vehicle reacts according to its observation of the environment and prediction of pedestrian's behaviour (c). In this study, we unfold the first part by analyzing pedestrian waiting time, with an emphasis on 1) naturalistic behavioural data 2) data-driven learning of pedestrian behaviour 3) interpreting the model results, and 4) discussing possible practices towards pedestrian friendly future urban areas. Obtaining data for studies involving pedestrian crossing behaviour is often difficult, as it could be dangerous to participate in various crossing scenarios under different conditions. Related studies mainly rely on data extracted from video footage, which often fails to capture details required for policy suggestions and best practices, as it is often collected passively from environments with no or limited form of control on variables of crossing scenarios. To tackle the difficulties of data collection for futuristic studies, and to have more realistic responses than stated preferences surveys, a virtual reality experiment is designed and implemented in this study, involving around 180 participants over a period of 5 months.
To analyze the data, a neural network based survival model is developed, and the results are compared to traditional survival models. To enhance the performance of the model and add to its explainability, a feature selection method and a post-hoc game theoretic based interpretability method are added to the framework, respectively. Various researchers have studied the technology, adaption rate, demands, etc. of automated vehicles. However, according to a survey on recognizing the barriers to cities' AV efforts, the lack of clarity on issues that require city actions is pointed out as one of the biggest obstacles for preparing cities for automated vehicles (Bloomberg and Aspen, 2017). We try to address such issues by providing detailed insights from our experiments, models, and results. This paper builds upon our previous work (Kalatian and Farooq, 2019) where an initial deep learning framework was proposed for modelling waiting times. For the purpose of this study, we enhanced the modelling, added a comprehensive analysis for interpretability of deep learning models, and based on the results, proposed suggestions and practical implications which can be useful for urban planners, AV manufacturers, transportation researchers and policy makers. The rest of this paper is organized as follows: after a brief review over the literature on waiting time studies, neural network-based survival analysis and machine learning interpretability in section 2, the model is described in detail in section 3, along with the methodology implemented for feature selection and model interpretability. The data collection procedure, from experiment design to data summary is then presented in section 4. Model results and analysis are then provided in section 5, followed by a discussion of practical implications and policy insights in section 6. Finally, a conclusion of the research and the future directions of the study will conclude this paper in section 7.

Background
First, a general overview of the research on pedestrian crossing behaviour is presented. As per the scope of the study, we mainly focus on studies on crossing behaviours at unsignalized and mid-block crossings. Survival models have primarily been used in such studies. Thus, a brief review on the background of survival analysis methods and their transformations using deep learning is provided. Finally, we discuss the literature on interpretability of machine learning models, in particular in transportation research. As machine learning models are considered mainly as black boxes, developing a systematic interpretability method is essential in our study to recognize contributing factors to pedestrian wait time and infer insights for practical implications of the results.

Pedestrian unsignalized crossing behaviour
Studies on pedestrian crossing behaviour at unsignalized intersections and mid-block crossings have predominantly focused on the concept of accepted gap time. Das et al. (2005) found that people waiting on medians accept shorter gaps more easily, compared to those waiting on the curb sides. Distribution of gaps was estimated in this study using parametric and non-parametric approaches. Oxley et al. (2005) focused on the effect of age on the ability to choose safe gap time. First, ANOVA was used to compare wait times of participants. According to their results, the primary contributing factor in deciding to cross, which can be interpreted as waiting time, appeared to be the distance to upcoming vehicle and time of arrival. Moreover, elderly participants appeared to select more unsafe gap times, given that their walking speeds were slower. In a relatively more recent work, Kadali et al. (2014) used artificial neural networks to estimate pedestrians' gap acceptance at mid-block unmarked crossings and compared its performance to multiple regression models. Their results showed that ANN has a better prediction performance, being able to consider the effects of more number of variables. However, the authors mentioned the strength of regression models in such cases due to their ability to reflect the significance of contributing variable.
Systematic modelling of pedestrian wait time using cox proportional hazards (CPH) models dates back to early 2000's, where Hamed developed a CPH model to analyze the time pedestrians spend on the sidewalk before initiating a cross on mid-block crosswalks. Data used in this study was collected by observing and manual recordings, and interviews at pedestrian crossings. Among the factors analyzed, previous accident experiences, group size of people crossing, car ownership, trip purpose, gap time between cars and finally, age and gender were appeared to be the contributing factors to wait time of pedestrians in their sample size (Hamed, 2001). Interestingly, parameters related to traffic were not found to be statistically significant in Hamed's study. In another application of survival models in pedestrian wait time analysis, Wang et al. (2011) investigated the effects of personal characteristics on crossing behaviour using a parametric survival analysis. The effect of age, trip purpose, safety awareness and conformity behaviour appeared to be the contributing factors to pedestrian wait time in this study.

Survival models and data-driven machine learning
In survival analysis, the aim is to analyze the time until an event occurs. A collection of statistical modelling methods can be designed for this purpose. The most common methods for survival analysis are the Kaplan-Meier model (Kaplan and Meier, 1958) and Cox Proportional Hazards (CPH) model (Cox, 1972). The Kaplan-Meier model is a non-parametric model for estimating the survival function in homogeneous groups. This model is very easy to implement, but is unable to account for individuals. To take into account covariate vectors and compute survival functions for individuals, CPH model is a conventional solution. CPH model is a semi-parametric model that assumes the time component and the covariate component of hazard function to be proportional. Despite being the common method of survival analysis for years, the assumptions made in Cox Proportional Hazards are not always true and have limitations. It assumes: • A linear combination of covariates within the hazard function.
• A constant hazard function over time.
To address these issues, Yu et al. (2011) introduced the Multi-Task Logistic Regression model. Their model can be interpreted as having different time intervals with different logistic regression models that estimate the probability of occurrence of an event in a interval. Although Multi-Task Logistic Regression models address some of the problems of the previous methods, they still remain linear, thus they cannot capture the nonlinear complexities within the data. Several researchers, particularly in the field of medical sciences, have tried to address this issue by implementing a deep learning approach (Faraggi and Simon, 1995;Katzman et al., 2016;Luck et al., 2017). Faraggi and Simon (1995) first incorporated a feed-forward neural network as a shift to nonlinear proportional hazards model. In their model, they replaced the linear combination of covariates with the output of a neural network with one output node. However, research followed by suggested that their network did not perform better than linear Cox Proportional Hazards models (Mariani et al., 1997). As the mentioned work were done prior to the outburst of modern deep learning algorithms, Katzman et al. (2016) decided to test the performance of more recent deep networks on survival models. They added fully connected and dropout layers and outperformed the performance of previous linear cox proportional hazards models. As a suggested future direction, the authors suggested implementing architectures like Convolutional Neural Networks to enable estimating waiting time directly from medical images.
Strong performances of survival analysis models in different areas, along with their relatively widespread usage in pedestrian wait time modelling, made them a suitable base model for our study. In an earlier study in our lab, we started analyzing crossing behaviour by studying distracted pedestrian's waiting time using Cox Proportional Hazards models, and we addressed various contributing distraction factors (Kalatian et al., 2018). Later, in (Kalatian and Farooq, 2019), we introduced Deepwait, which is a neural network based extension of CPH, empowered with a feature selection algorithm. However, despite the improvements in accuracy we obtained by using Deepwait, lack of a interpretability mechanism makes the application of such models limited, especially as learning the contributing factors are of vital importance. To better understand the problem of interpretability in machine learning models, and to be able to propose policy recommendations based on the results of our model, a review on the existing literature of machine learning interpretability is presented in the following subsection.

Interpretation in machine learning
Widespread application of machine learning methods for decision making and policy planning requires a high level of interpretability for often complex, dense, and deep networks. Despite high prediction accuracies obtained by machine learning models in recent years, interpretability and explainability of the algorithms are yet considered as a cumbersome and in some cases, even unnecessary task. Major arguments against interpretability of machine learning models are the notion that performance of a model is more important than model interpretability and, the burden of understanding the model can prevent adoption of machine learning models (Lipton, 2016). Nevertheless, the trade off between model complexity and explainability may weigh in favor of using traditional explainable models in some areas where gaining information on why? and how? the model works is as important as the predictivity of the model.
To fill the existing void in machine learning models, various approaches to make machine learning models explainable have emerged in recent years. The black-box nature of machine learning algorithms can be addressed by model-agnostic methods, which are post-hoc interpretability methods that rely solely on input and output of the models, disregarding their structure (Molnar, 2019). Mainly practiced model-agnostic methods include sensitivity analysis plots (Friedman, 2001;Goldstein et al., 2015), feature permutation methods (Fisher et al., 2018), surrogate local interpretable estimators (Ribeiro et al., 2016), and game theoretic based approaches (Štrumbelj and Kononenko, 2014;Lundberg and Lee, 2017). In transportation, not so many studies have addressed the interpretability of the models despite the relatively large number of works on machine learning and deep learning applications in the field. In the interpretability of choice analysis for instance, Hagenauer and Helbich (2017) used a permutation based variable importance method upon their DNN model to analyze the performance of the model when the values of a variable changes randomly. In (Wang et al., 2018), by considering the inputs to the Softmax layer as the utility functions, authors compute economic information using DNNs. They conclude that economical information aggregated over the population can be reliably derived from DNN. However, disaggregate information for individuals still lacks the reliability acquired by discrete choice models. Wong and Farooq (2020) used Jacobian determinant of generative models to calculate the elasticity of mode choice with respect to different explanatory variables. Jacobian matrix in this study was generated for each instance of all the conditional outputs, and density of elasticises was estimated across the data points. Use of Jacobian matrix, however, is computationally expensive and cannot guarantee to observe the interactive effects of variables. For prediction and estimation tasks, Das and Tsapakis (2019) used partial dependency plots to intuitively show the effects of different variables and detect contributing factors on predicting traffic volumes. Most of the existing work in transportation rely on intuitive and illustration-based interpretability frameworks without considering the interactive effects of variables in a systematic manner. To address this issue, more recent and advance methods of game-theoretic based interpretability methods will be used in this study to find the contributing factors to pedestrians' waiting time before crossing.

Methodology
Pedestrian wait time before they start crossing can be modelled as the time it take before an event occurs. Thus survival analysis is a popular and powerful tool used for this purpose. In the following, first traditional survival models, and in particular, Cox Proportional Hazards (CPH) models are explained. We used CPH as the base model for comparison in this study. Then, the data-driven deep neural network version of CPH is proposed as the framework used in this study, which enables incorporating more covariates while increasing goodness of fit. Finally, interpretability mechanism used to explain model results and obtain practical implications is outlined.

Survival analysis
Survival analysis techniques are used in cases the time to occurrence of an event is of interest. Although originally developed for applications in medical science where the event is treatment of a disease or the death of the patients, survival methods have found their way in other fields, including transportation. By defining survival function S (t) as the probability that the event occurs after time t, and hazard function h(t) as the probability of the event occurring at time t, the following relationship can be established between the two: Different methods have been proposed for estimating hazard function, parametric, semi-parametric and non-parametric. Among them, Cox Proportional Hazards (CPH) model is probably the most widely used method, which enables considering the effects of covariates on the hazard function (Cox, 1972). Equation 2 shows the hazard function formula in CPH.
In the above equation, h 0 (t) is the baseline hazard, which is the hazard function regardless of the values of the covariates and is a function of time, and β i and Z i are the coefficient and value of covariate i, respectively (Therneau and Grambsch, 2013). In the terminology of survival analysis, e i β i Z i is called partial hazard function or risk function, which is representative of the effects of values of covariates on the hazard function and is independent of time. As it can be seen in equation 2, the log-partial hazard function in CPH is in a linear combination of covariates. To find the coefficients in CPH, the partial likelihood shown in equation 3 is maximized.
Where, index j in the denominator selects the instances that are still at risk when event of instance k occurs. In this study, we use CPH with linear log-partial hazard function as the base model, in which an event is defined as initiation of a cross by a participant. The simplistic assumption of linearity of log-partial hazard function is what drives us to improve and update the model to account for today's complex data and advanced modelling tools available.

Data-driven survival analysis
Rapid and ubiquitous emergence of data-driven approaches in recent years provides an opportunity to further improve the models for survival analysis. Due to the availability of complex and high-dimensional data today, traditional methods may not be capable of capturing nonlinearities within data. Thus, we build upon the classic CPH (Equation 2) and replace the linear log-partial hazard function ( i β i Z i ) with a deep neural network, g(w) (Kalatian and Farooq, 2019). The proposed framework used in this study is outlined in Figure 2, and the model components are explained in the next two subsections.

Feature importance ranking
High-dimensionality of modern data sources has led researchers and machine learning engineers to incorporate feature selection methods to prevent overfitting and performance degradation (Alelyani et al., 2018). Although deep neural networks are known for their capability of detecting correlation among features, a Variance Inflation Factor (VIF) analysis and feature selection algorithm is implemented in the preprocessing step to 1. develop a base case cox proportional hazards model, and 2. compare the effectiveness of deep survival models with and without data preprocessing.
Initially proposed by Kononenko et al. (1997) for classification, Relief family of algorithms consider the ability of variables to distinguish among instances that are close to each other,and rank their quality based upon this ability. RRELIEFF (Robnik-Šikonja and , which is an adaptation of Relief for using in regression, is used in this study to prioritize covariates, and reduce the dimensionality of the dataset collected. Based on Relief in its simplest form, the importance of a covariate in a classification task increases, if its value is different for instances from different classes, and decreases if its value is different for instances of the same class. RreliefF, which is the regression version of the family, replaces the classes with a probability defined based on the relative distance between the target values of the instances. In mathematical terms, importance weight of covariate f , W f , is calculated based on the following equation (Robnik-Šikonja and Kononenko, 2003): Where: • P f = P(different value of f in nearest instances) • P v = P(different target value in nearest instances)

Hazard function
The formulation of modified hazard function in our proposed model is shown in equation 5: where g w is the dense neural network with weights w. Top n covariates (Z n ), based on the feature importance ranking, will be imported into the network's input layer. After a number of fully connected hidden layers, each followed by dropout layers and batch normalization to prevent overfitting to the training data, the output of the network g(Z n ) will be calculated as the value of the last one-node layer, and will be used as our estimation of log-partial hazard. In the deep network described, network architecture parameters i.e. number of covariates to be selected, number of hidden layers, number of nodes in each layer and dropout rate as well as optimization parameters, including learning rate decay and momentum, are used as the hyperparameters to be found. The loss function to be minimized in network training is derived from the partial likelihood function of CPH. As shown in equation 6, network's loss function would be the average negative logarithm of CPH's partial likelihood: where N is the total number of events or instances, and Z nk is the is the value of top n covariates for instance k.

Model interpretability
Transition from regression models to neural based models comes with the cost of losing useful information on the effects of covariates on the model output. Developing traditional CPH models, the effect of each covariate on the log-partial hazard function can be estimated by looking into covariate coefficients. A greater coefficient for a covariate means an increase in log-partial hazard function, and subsequently in hazard function. In other words, a positive coefficient for a covariate means that greater value of that covariate results in higher probability of the event to occur at any time, thus a shorter wait time. By using z-scores and p-values, the significance of variables can be measured in traditional based CPH models. The case is different for neural networks, however, as it is impossible to track the nonlinearities involving multiple variables at the same time and in different layers.
To address this issue, SHAP (SHapley Additive exPlanations) is used in this study. Introduced in (Lundberg and Lee, 2017), SHAP is a game theoretic model explanation method based on Shapley Values. In game theory, Shapley value is a method to fairly distribute the payoff of a job done by a coalition of players with different skills, to the players. By replacing the payoff with model prediction and players with features, SHAP calculates the average marginal contribution of each feature to model prediction for each instance. In mathematical terms, Shapley value of each feature i is calculated by equation 7 as follows (Lipovetsky and Conklin, 2001): In which S is a subset of all features F, g S ∪{i} is the model trained using a subset with feature {i} present, and g S is the model trained with the feature i missing. One of the advantages of using Shapley values over other interpretation methods is that it considers the fact that the contribution of a feature depends on the values of other features, thus the contributions are computed for all possible subsets. In their paper, Lundberg and Lee (2017) show that SHAP solution satisfies required conditions to fairly distribute the payoff to players. By using Shapley values, importance of the features can be estimated as the average absolute Shapley values of features over the whole instances.

Controlled experiment and data
4.1. VIRE Virtual Immersive Reality Experiment (VIRE) is used in this study as the tool to develop the controlled experiments and collect the required data. We selected virtual reality, as our experiments involve situations that are either dangerous in reality and may result in disastrous outcomes, or are futuristic scenarios that are not yet feasible to implement. On the other hand, stated preferences surveys, which are the conventional tool of data collection in similar studies, cannot provide naturalistic datasets as participants have no prior exposure to automated vehicles.
To study perception of people, pictures, maps and videos have been used in the past. Recent developments in Virtual Reality (VR) environments have made investigation of human perception and behaviours in controlled conditions possible. Combination of the physical environment and virtual elements, i.e. information or images, broadens the opportunities for content delivery Jennett et al., 2008;Animesh et al., 2011;Nah et al., 2011;Faiola et al., 2013). Experiments in VR environments have successfully been conducted in different fields in cognitive studies . In the field of transportation, researchers have used VR in pedestrian route choice studies and their reaction to information in evacuation scenarios. However, these studies mainly lack the interactive and dynamic potential of VR. More specifically, user's actions in the VR environment need to get responses from the elements of the environment, and vice versa.
VIRE is an in-house developed virtual reality framework designed for experiments in controlled environments . The unique feature of VIRE is that the virtual objects are designed such that they react to the participant's actions. For instance, if a participant is crossing the road, based on their location, the approaching vehicle may slow down or completely stop, thus giving a dynamic, immersive, and realistic experience. A scenario is generated based on a set of variables and then it is projected onto the eyes of participants through an immersive head mounted display. Participants, finding themselves on a simulated 3D sidewalk, are asked to start crossing when they feel it is comfortable to do so. During the experiments, data on the movements of pedestrians is recorded using motion sensors and the reactions of virtual objects are computed in real-time.
Each scenario is defined by 9 controlled variables, which are selected based on related literature on future urban streets (Bloomberg and Aspen, 2017), while considering the available facilities. These variables are categorized as follows: • Rules and regulations: speed limit, minimum allowed gap time between vehicles Multiple levels are defined for each of the controlled variables, as presented in Table 1. The standards of speed limit, suggested minimum gap between cars and lane width at the time of experiment in a typical road in downtown Toronto were 50 km/hr, 2 sec, and 3 m, respectively. To take into account future possible modifications to the current standards, one of the levels for each of these three variables are set in our experiment to have the same value of a typical Toronto street standard, with the other two levels considering possible changes. Number of braking levels, which is only defined for automated vehicles, represents the level of smoothness of braking system of the AVs. Having a greater number of levels shows a more smooth deceleration when AVs face a barrier. Simulated traffic on the road can either be fully human driven, fully automated or a mixture of automated vehicles and human driven cars. In the experiments, users distinguish human driven and automated vehicles by observing the absence of a driver in the driving seat in AVs, and a box-shaped LiDAR sensor on top of them. Moreover, AVs and human driven cars are different in braking systems, which might not be realized by the participants, unless they do the experiment as such that a vehicle has to brake for them. As for the demand on the road, values are extracted from flow rates on a typical road in downtown Toronto during off hour and peak hour hourly average (530 and 1100 veh/hr/lane respectively) and the mean value of these flow rates (750 veh/hr/lane). Two flow variables are also derived based on the controlled variables and used in the analysis: distance between cars and traffic density. Finally, the environmental conditions in the simulation are designed as such that participants experience day and night time of the day, as well as clear and snowy weather, in which the sight distance of pedestrians is affected. While doing the experiments, coordinates and head orientations of participants are recorded in intervals of 100 milliseconds, as well as the coordinates of simulated vehicles. Scenarios are then generated based on the above-mentioned controlled variables. Each possible combination of different levels of controlled variables is a potential scenario. To select the scenarios for the experiments, a design of experiment is conducted, which is described in detail in the next section.

Design of experiment
To study the effects of various different attributes on pedestrians' crossing behaviour, a design of experiment (DoE) needs to be conducted. The first and simplest method for DoE is full factorial design, in which after assigning different levels to each attribute, all possible combinations of different attributes' levels are considered for the experiment. In our case, however, this is an impossible approach to take as there would be a total of 8,747 possible combinations based on the 9 attributes we defined and their levels. To overcome this problem, several Fractional Factorial Designs have been developed based on the objectives of the design, models to be developed and the amount of prior information available (Cavalcante and Roorda, 2011). In fractional factorial designs, a subset of choice tasks are selected. The most widely used strategy for selecting a subset of tasks is orthogonal design, in which tasks are selected in a way to produce zero correlation between attributes (Ortuzar and Willumsen, 1994). However, orthogonality may limit the experiment design by putting constraints in terms of the number of runs required and possible settings for factor levels. Recent studies on the topic have led to Optimal Designs, which focus on the efficiency of the experiment, rather than orthogonality (Rose and Bliemer, 2009). In optimal and efficient designs researchers try to find the efficient design in terms of a selected measure of quality for parameters' estimates. In the most widely used optimal design, D-Optimal design, for instance, variances and covariances of parameter estimates determine attribute level combinations (Atkin-son et al., 2007). For the purpose of data collection for this study, we use an optimal design of experiment, as it allows experiments to be conducted in a more flexible manner in terms of number of runs required and tasks to be observed.
In D-Optimal designs, a number of assumptions and parameter estimates for the model are required in order to calculate covariance matrix, as this matrix varies based on the model used. In order to deal with parameter priors required in these methods, two approaches can be taken. First approach is to consider parameter priors to be zero, which is called the null hypothesis (Street et al., 2005). This assumption indirectly results in the orthogonality of null hypothesis D-Optimal designs. In addition, it is assumed in this type of methods that the model used is Multinomial Logit (MNL). Not requiring a priori knowledge on parameter estimates makes this method useful when no prior information is available. The second approach, on the other hand, assumes that some estimates are available on the values of parameters. Rose and Bliemer suggested that a priori parameter estimates can be obtained from: the literature, pilot studies, focus groups and expert judgments (Rose et al., 2008). The main advantage of this approach is that it enable researchers to use any model type, in contrary to the former approach that was applicable only to MNL.
In this study, we obtained estimates on a priori parameters based on a trial experiment on a limited number of participants. To do so, a demo of the data collection was conducted with 5 participants, each going through 30 randomly selected scenarios to get an initial idea on the model parameters. As the objective of the data collection was to obtain required data to model pedestrian wait time, we assumed a cox proportional hazards model (CPH) for waiting time 2 . Thus, a design of experiment was developed based on the CPH model on the trial data, and top scenarios based on importance weights were selected for our experiment. In total, 90 different scenarios were selected. The design of experiment formulation for wait time modelling, which is the topic of interest in this study is presented in the next section.

D-optimal design for Cox Proportional Hazards (CPH) model
Schmidt and Schwabe (2015) proposed the D-Optimal design for a CPH with three covariates. In this section, we extend and generalize their formulation for cases where the number of covariates in the Hazards model are more than 3. Cox proportional hazards model with fixed termination time c, is specified as: Where Y i is the survival time for respondent (participant) i, and β is the vector of parameters associated with k independent variables. Now if the data is type I censored, meaning that the experiment stops at a predetermined time which is the case of our experiment, the Fisher information matrix is given by: in which x is the set of independent variables, and for a linear specification, we have: Suppose that we need to generate m scenarios for the experiment, then In which x m represents the m th combination of k independent variables and (0 ≤ ω m ≤ 1) is the associated weight, where m ω m = 1. For such model the information matrix M(ξ, β) will be: Given β, x, and m, we need to find a design ξ * β that maximizes det(M(ξ * β , β)). For a high number of independent variables, we can use Markov Chain Monte Carlo simulation for finding ξ * β . In particular, Simulated Annealing (SA) is a great candidate for such problems (Kirkpatrick et al., 1983).

Data collection campaign
Data collection campaign for this study started in April 2018, as a trial experiment with limited number of people, and then continued with main experiments until September 2018. In total, 180 people participated in the experiments, 160 of which succeeded to complete the experiment. The completed participants consisted of 113 adults, and 47 kids and teenagers. To make the data as inclusive and heterogeneous as possible, we conducted our experiment in four different locations. Started at Ryerson University with post-secondary students, we then moved to Toronto City Hall and North York Civic Center to include professionals in the field who were more familiar with the nature of the experiments and the questions we were trying to answer. We repeated the experiments in Markham City Library, in which we included general public. Finally, we had our experiments with Maxim City summer school, which included kids and teenagers between 9 to 15 years old. In Figure 3, a view of a sample scenario and the experiment setup with participants is shown. A video of a participant, while doing the experiment is also attached in the supplementary materials.
The process of an experiment for a participant is as follows: The participant is first asked to fill a questionnaire on his/her sociodemographic information, travel patterns and previous experiences with virtual reality. The participant is then familiarized with the VR environment for 5 minutes. After familiarization, 15 scenarios are randomly drawn from the 90 selected scenarios by the experiment design. For an adult participant, each scenario is conducted two times, accumulating to 30 total experiments with a total time of about 30 minutes. Our trial study showed that continuing experiments for over 30 minutes causes fatigue among participants and thus the results may be affected. For adolescent participants, each scenario is conducted once, with a total time of 15 minutes for the whole experiment. It should be noted that data from child participants are not used in this particular study as their behaviour when facing VR environment and AVs mostly involved unpredicted reactions and crossing behaviour, and we believe a different data cleaning is required on their data. A video of a young participant, while doing the experiment is provided in the supplementary materials.

Data summary
Before applying feature selection algorithms, the first step in data preprocessing would be to remove highly correlated variables from the dataset. To do so, Variable Inflation Factor of covariates are calculated, which for each covariate represent the multicollinearity of that covariate. In general, a higher VIF for a covariate indicates the higher ability to predict that covariate using linear regression based on other covariates on the dataset. After removing experiments where pedestrians could not complete the crossing, or could not follow the experiments' procedures, a total of 2,291 responses were remained to be studied. Collected data is then used to estimate models for predicting pedestrian waiting time. In Figure 4, the frequencies of crosses with different wait times are provided. As it can be inferred from this figure, a majority of the crosses (54%) occurred in the first two seconds, decreasing afterwards with a relatively constant rate. According to the figure, Only 6% of the crossings took more than 20 seconds of wait time. All quartiles of waiting time for participants, as well as the average waiting time, are depicted thorough box plots in Figure 5 for different levels of all covariates. As shown in the figure, more participants in our experiment tended to wait longer when exposed to automated vehicles, and in mixed traffic conditions participants even took longer to cross, compared to fully automated conditions. The reason might be the effect of having two different types of vehicles in these scenarios, which increases uncertainty among participants. In the second plot, waiting time over different levels of speed limit, is not as expected. Participants in our sample data waited longer in slower traffics. The reason for this behaviour can be that participants waited for a number of vehicles in scene to pass first before crossing, which takes longer for slow vehicles. To address this issue, other traffic parameters like density and flow are analyzed to replace speed limit. The next box plot reveals that in wider lane widths, participants in our experiment waited longer before crossing. No observable difference is identified in the median of minimum gap time allowed, with 1 second gap times having a slightly higher wait times compared to 2 seconds gap times. According to the arrival rate box plot, participants in more congested areas with higher traffic flows waited longer before crossing. Higher values for derived variable of vehicles' density also had a positive correlation with waiting time. Values of density are continuous in the experiment, and are categorized to three equal intervals in the figure. Three defined braking types for automated vehicles are shown next. As participants are not informed about the type of braking level before the experiments, it seems that a particular trend does not exist with them and waiting time. Two-way roads with median seemed to result in shorter waiting times compared to two-way roads with no median. Regarding environmental variables, no notable pattern based on weather conditions and time of the day is seen in the plots. In the age categories analyzed, two extremes of the analyzed age groups have greater waiting time, and in gender, females have slightly longer wait times. Five next plots are relevant to travel habits of participants. Almost all of these plots show longer wait times for participants who tend to have less tendency to have active modes of transportation. Finally, participants who have prior experience in using virtual reality wait shorter on the sidewalk, as they probably feel more comfortable using the devices. It should be noted that a detailed analysis of covariates require developing survival models, which are presented in detail in section 5.

Modelling results and analysis
Three models are developed and compared in this section: a traditional CPH model with linear log-partial hazard function, a deep neural network-based CPH function (DCPH1), and a deep neural network-based CPH function empowered by VIF and RreliefF (DCPH2).

Cox Proportional Hazards model
Applying VIF before developing the model is essential to remove highly correlated variables, as the models consisting of them fails to converge because of non-invertibility of singular matrices. The model is then developed in Python's lifeline package (Davidson-Pilon et al., 2019) starting with all the remained variables and iteratively removing insignificant variables. Models are developed over a training dataset, consisting of 80% of the data, and evaluated on a test dataset, both of which selected randomly. To measure the goodness of fit of the model, Concordance Index (C-index) is used as the global index for comparing survival times. C-index is defined as the fraction of pairs of observations, in which the observation with higher predicted value for log-partial hazard function has shorter waiting time (higher risk of crossing). Developing a CPH model over the train dataset using 10-fold cross validation, a C-Index of 0.57 is achieved on the test set. Details of the performance of CPH model are presented in Table 3.

Deep CPH model (DCPH)
After applying VIF to remove highly correlated variables, variables are ranked based on their importance according to their RReliefF importance weights, before training the network. Hyperparameter n is introduced as the number of covariates to be used in training the network. Other hyperparameters include: number of hidden layers, number of nodes in hidden layers, learning rate and learning rate decay for optimization, whether to use batch normalization, and dropout layer and its rate. To find the best model, random search for optimizing hyperparameters is utilized (Bergstra and Bengio, 2012). In Table 2, the optimum values for hyperparameters based on a 10-fold cross-validation over 100 epochs on training dataset are presented. Same test dataset as the one used for CPH is used again to evaluate the performance of the network over 100 epochs. Using a framework with n =19 top covariates as inputs, three hidden layers of 90 units each, batch normalization and a dropout rate of 0.1, a C-index of 0.64 in cross-validation and 0.62 for test set is achieved, showing an improvement of 5 percent compared to CPH.
To further assess the performance of the proposed framework, another neural-network based survival model is developed without incorporation of RreliefF and VIF (DCPH1). C-indices of all three models, along with the number of covariates incorporated in the models are presented in Table 3. It can be observed that DCPH2 outperforms two other models, with 5% and 2% improvement over Linear CPH and DCPH1, respectively. Comparing number of  covariates used, using DCPH2, the process of removing insignificant covariates and retraining the models to find the best combination of covariates manually is no longer required. Compared to DCPH1, DCPH2 results in a better C-index with less number of covariates, meaning a less computationally expensive prediction in DCPH2.

Cox Proportional Hazards model
Covariates used for developing CPH are listed in Table 4, along with their coefficients, Hazard Ratios and P-Values as a measure of significance. P-Value of 0.1 (90% confidence interval) is set as the threshold of significance for selecting the covariates. Third column of the table provides Hazard Ratio, which is defined as the ratio of the hazard rates of two conditions of a covariate. A hazard ratio of 0.84 for binary covariate Age Over 50, for instance, means that the hazard of a pedestrian older than 50, over the hazard of other pedestrians equals 0.84, implying lower probability of crossing (longer waiting times) for pedestrians aged over 50. Hazard ratio is calculated as the exponential of coefficients in the second column. In other words, a higher coefficient, implies higher hazard ratio, and a hazard ratio greater than 1 (positive coefficient) means that the probability of a cross is higher for instances having greater values of that covariate, leading to shorter waiting times. Based on the results in Table 4, Participants who walk for shopping, have previous VR experience, aged between 30 and 39, or have no cars in the household wait shorter in our sample of collected data compared to their counterparts with different values for each of the mentioned covariates. Moreover, experiments on two-way roads with median take shorter wait time compared to two-way roads with no median as the baseline variable of this covariate. On the other hand, participants who indicated their main mode of travel as cars, are aged over 50, or are females, have longer wait times. In addition, experiments in higher traffic densities and higher lane widths makes pedestrians wait longer in our sample data.
Going deeper into the effects of covariates, it can be inferred that participants who tend to have the walking habits in their life-style, i.e. walk for doing their shops and have no cars in their households, feel more comfortable crossing the streets, and have shorter waiting times. On the other hand, participants who use cars as their main mode of transportation, are expected to be less familiar with the crossing conditions and have to wait longer before initiating a cross. Considering gender and age covariates, it can be observed that females in our sample size waited longer before crossing, which is in line with some of the previous studies in the literature. The eldest age group of participants, who are over the age of 50, also tend to spend more time on the sidewalk before crossing the street. One the other hand, participants in the second youngest age group wait on the sidewalk for shorter times compared to the baseline category for age: 40 to 49. The developed model does not consider age category of between 18 and 29 as a significant variable, which can be resulted by the high heterogeneity of the participants in this age category. Regarding variables related to geometry of the road, participants tend to wait longer before crossings on roads with larger lane width. Significance of this variable can lead to the promotion of larger sizes for sidewalks and narrower lane widths. As stated earlier, out of the three road types, two-way roads with no median are settled as the baseline category of road type. Compared to the baseline category, medians help pedestrians cross with shorter waiting times. Other significant covariate set by the scenario is traffic density on road, for which higher values lead to longer waiting times by the participants, due to the time required to find the appropriate safe gaps for crossing. Finally, having previous VR experience helps participants of our experiment cross in shorter times. Despite having a 5-minute familiarization session to the equipment and environment before the experiments, it seems that further tutorials can be helpful to participants with no prior VR experience in future studies. An interesting observation on the results of CPH is the insignificance of the level of vehicles' automation, as well as other scenario-related variables, on participants' wait time. As different pedestrian's behaviours was expected for different scenarios, it can be inferred that linear CPH is incapable of capturing complex relationships among the covariates and wait time.

Deep CPH (DCPH) model
SHAP in its core, estimates the contribution of covariates to the waiting time for each instance separately. In Figure 6, SHAP Values for the covariates for all the instances are depicted. Each SHAP value, corresponds to the contribution of the covariate to log-partial hazard values, compared to when the covariate has a baseline value, called background value in SHAP terminology (Molnar, 2019). Default baseline values for a covariate in the original paper are set to be equal to the average value of that covariate. In this study, we used the default values for continues covariates. However, the question arises for binary variables, in which having a real number value does not have a meaningful interpretation outside of the model. To address this issue, we set the baseline value for binary variables to zero, thus comparing the contribution of each binary covariate of an instance to cases where the same covariate has a value of zero. By doing so, SHAP value of a covariate in instances that the value of that covariate is zero is not calculated (blue dots on the central axis in Figure 6) and contribution of a covariate is assessed based on SHAP values of the covariate in other instances. To have a numerical understanding of SHAP values, average and standard deviation of SHAP values for all 17 covariates used in training the network, and over all instances, are presented in Table 5, sorted by their importance according to SHAP values. As a rule of thumb, we consider the overall contribution of a covariate interpretable, if it has a relatively uniform effect on a major part of instances. In other words, if for a covariate, the absolute mean of SHAP values is greater than the standard deviation of the SHAP values, over all instances, we call the effect of that covariate uniform. It should be noted that non-uniformity of SHAP values does not mean insignificance of that covariate. Non-uniformity reveals that the interpretability method used could not capture the complex non-linear multi-level relationship among covariates. Intuitively, a covariate has non-uniform SHAP values when it has a wide range of SHAP values, extending over negative and positive values, as in these cases, the effects can vary based on the instance. Interactions with other covariates may be investigated to further analyze the effects of such covariates in future studies. In this study, however, we avoid multi-level analysis for the sake of avoiding complexities involved and numerous possible combinations available for investigation.
Looking into SHAP values in Figure 6, it can be observed that unlike CPH that was not able to capture the effects, automation level covariates have generally a negative contribution to log-partial hazard, leading to longer waiting times. Relatively high mean and low standard deviation of fully automated and mixed traffic conditions in Table 5 confirms the significant negative effect of these two covariates on log-partial hazard. As pedestrians feel more unfamiliar and concerned with cars with no drivers, the effects were expected, but could not be captured by linear CPH. Positive contribution captured for previous VR experience is in line with the results obtained in CPH, which confirms the necessity of stronger VR tutorials before experiments. Similar to CPH, wider lane width and higher traffic density both have negative effect on log-partial hazard, implying longer waiting times. Their mean and standard deviation confirms the uniformity of SHAP values for these two covariates. Six of the used covariates are relevant to walking habits: Walk to work and shopping, main transportation mode of cars and active modes, having driving licence and having more than 1 car in the house. Among them, only walking to work appears to be interpretable based on mean and standard deviation of SHAP values, which shows shorter waiting time for participants who indicated that they walk to commute to their works. Environmental covariates, on the other hand, are added to the model and appear to have uniform effect. Results showed that participants crossing in simulated snowy scenarios, have longer wait times, which may be traced to poor sight distance. On the other hand, in the night scenarios, participants tended to have shorter wait times, which is not consistent with the supposed poor sight distance in the simulation. This may be due to the fact that in night scenarios simulated, unlike snowy simulated scenarios, the sight observed by participants was not limited and changes were only made to the color of the sky to project the mental effect of night. This can be addressed in future data collection campaigns. Finally, the effect of road type covariate did not appear to be interpretable over all instances, and vehicle arrival rate, although not one of CPH covariates, appeared to be one if the selected covariates by DCPH2, without captured uniform contribution for all instances.

Discussion and policy implications
There are several key practical implications and useful policy recommendations that can be derived from the virtual reality experiments, models, and results of this study, which can be of use to urban planners, policy makers, manufacturers of automated vehicles, and researchers seeking to implement virtual reality and data-driven modelling in their work.
Regarding participants' experience in VR, extra attention should be paid to children and senior participants in terms of training and educational programs, and their comfort during the experiments. Moreover, providing facilities for people with special needs should be taken into account in future research to make experiments more inclusive. Out of the 180 participants recruited for the experiments, 20 who were children and senior participants, failed to complete the tasks. Although we had one person using wheelchair as a participant, more accessibility services should be provided in future to increase the participation of people with different disabilities. As per our observations, some children tended to see the virtual reality tools as gaming platforms. Irrational crossing behaviours to explore the 3D environment, followed by their loss of interest in continuing the experiment after observation of different elements of the environment, led to their failure in completing the experiment tasks. The virtual AVs were designed such that they could foresee if the pedestrian is in danger and would apply braking when needed. Some teenage participants realized that and started to play with the virtual AV. They would come in front of an AV so it would stop. Participant would then move back, and virtual AV would sense this action and start moving, at which point, the participant would move in front of the AV again causing it to stop, repeating such actions several times. This behaviour is expected to happen when the fully automated vehicles are operating on urban roads, especially with low speeds, resulting in disruption of the upstream traffic. To avoid such behaviour, VR based education and training is needed. Municipalities can also look into proactive regulations to minimize such disruptions. As for senior participants, their inconvenience while using VR headsets, inability to ambulate independently without mobility devices such as wheelchair, and feeling nauseated and fatigued while interacting with 3D immersive environment were identified in our observations as the main cause of them not completing the tasks. Interpretation of our model results can be useful for practical implications for urban decision makers and car manufacturers. Our results show that participants tend to be more conservative and cautious in the presence of automated vehicles. Either in fully automated or mixed automated and human driven conditions of traffic, participants waited longer compared to solely human driven conditions. Considering the unfamiliarity of pedestrians with automated vehicles in the first years of introduction of AVs, this observation should be taken into account by city planners and decision makers to make required modifications to enhance pedestrians' crossing experience. Nationwide educational training programs should be practiced before the transition to automated environments, to familiarize pedestrians with new dynamics of the city. Immersive and dynamic virtual reality technology can play an important role in such programs, to ensure providing safe and naturalistic experience to users. Manufacturers should consider alternative ways to replicate the current silent social interactions between the driver and pedestrians in automated environments. Some manufacturers, like drive.ai, have introduced screens on their vehicles that can give visual cues to pedestrians. Such new treatments can be systematically optimized and their effectiveness can be studied in VIRE environment. Based on results of our study, narrower lane widths, lower traffic densities, and better sight distances are also revealed to be affecting parameters in pedestrians' crossing behaviors, leading to shorter wait times. Wider and more comfortable sidewalks, narrower lane widths, enhanced lighting equipment, and incorporation of pedestrian-to-vehicle communication technologies are some of the solutions that can be implemented before diving into future automated urban areas. Our study reveals that having frequent walking habits positively affects the crossing experience. Promoting active modes and developing more pedestrian friendly infrastructure could lead to better crossing experience on the streets of the future.

Conclusion and future direction
In this paper, we investigated factors affecting pedestrians' wait time before crossing unsignalized crosswalks. With the upcoming revolutionary technologies on the roads, it is of vital importance to re-think and re-asses pedestrian behaviour in the presence of these unprecedented technologies. Studying new technologies require new tools, thus, a virtual reality based immersive and dynamic experiment was introduced in this study to obtain high-dimensional data from customizable scenarios in a safe way. Over a period of five months, a total of 180 people, selected in a heterogeneous and inclusive way, participated in our experiments. Participants were asked to cross an unsignalized crosswalk in various scenarios, incorporating speed limit, lane width, vehicles' arrival rate, road type, level of automation, braking system of cars, minimum gap between vehicles, weather conditions, and time of the day. Moreover, demographic information and travel habits of participants were collected using questionnaires before each experiment.
To analyze the effects of different covariates on pedestrians' wait time, survival analysis models were developed in this study. A traditional Cox Proportional Hazards (CPH) model as the baseline model, and a neural network based CPH model empowered by feature selection and interpretability methods were developed and trained in this study. Using an interpretable framework helps us capture nonlinearities among high-dimensional data, resulting in better goodness of fit achieved compared to the baseline model. Moreover, our framework outperforms simple neural-network based CPH models by following an embedded algorithm for systematic feature selection. Using a game theoretic-based interpretability method, we seek to replace traditional CPH methods by fulfilling their power in capturing the effects of covariates on baseline hazards function.
Based on the interpretability results of our study, we suggested some key practical implications and policy recommendations. Our results showed that pedestrians wait longer before crossing in the presence of automated vehicles. Wider lane widths, higher traffic densities, poor sight distances, and lack of walking habits in pedestrians found to be the other contributing factors to pedestrian wait time. Widespread educational campaigns before introduction of AVs on streets, enhanced safety measures on AVs, promoting active transportation modes, incorporating pedestrianfriendly infrastructures on streets and using pedestrian-to-vehicle communications are some of the possible solutions for future urban areas. In conclusion, this study tries to contribute to the current literature in three general aspects: • Utilizing immersive virtual reality tools for relatively large-scale data collection, backed by a systematic design of experiment to optimize the information that can be inferred from the data. A D-Optimal design for CPH models with more than three covariates was introduced in this paper for the first time. To the best of our knowledge, our VR data collection campaign is one of the biggest such campaigns for pedestrians studies.
• Developing a neural-network-based survival model, to analyze the effects of different parameters on pedestrians' wait time. Incorporating neural network within a CPH model, we achieved an improvement in accuracy by 5% compared to linear CPH models.
• Using SHAP, as a modern game theoretic-based approach for neural network interpretation. Interpreting neural networks is gaining more popularity in recent years, in attempts to transfer black-box models to explainable models that can be used for policy and decision making. Limited number of studies have touched model interpretability in transportation.
Our study was not without limitations. In terms of the data collection, as stated in the paper, participants training procedure can be enhanced to decrease the effect of previous VR experience on the performance. Providing accessibility services for people with disabilities to have more inclusive data collection is another important direction that needs to be followed. Other methods of design of experiment may be also tried to compare the performance of the proposed D-Optimal design introduced in this study. Regarding the model, other deep networks can be utilized within or independent of the CPH model. Predicting pedestrian wait time using developed hazard functions can be addressed in future studies, which requires developing a time-dependant baseline hazard function. Current findings can be useful in the development of new audio, visual or direct vehicle-to-pedestrian communication methods that can make the pedestrian crossing more convenient. The effectiveness of such methods can also be tested in the virtual reality environment. Finally, our paper is a part of an on-going research on interactions of humans and automated vehicle, which aims to develop prediction and training tools to develop socially-aware automated vehicles.