Bayesian Nonparametric Approach to Average Annual Daily Traffic Estimation for Bridges

Abstract

Every year, the U.S. government provides several billions of dollars in the form of federal funding for transportation services in the U.S.A. Decision making with regard to the use of these funds largely depends on performance indicators like average annual daily traffic (AADT). In this paper, Bayesian nonparametric models are developed through machine learning for the estimation of AADT on bridges. The effect of hyperparameter choice on the accuracy of estimations produced by Bayesian nonparametric models is also assessed. The predictions produced using the Bayesian nonparametric approach are then compared with predictions from a popular Frequentist approach for the selected bridges. Evaluation metrics like the mean absolute percentage error are subsequently employed in model evaluation. Based on the results, the best methods for AADT forecasting for the selected bridges are recommended.

Every year, the U.S. government provides several billions of dollars in the form of federal funding for transportation services in the U.S.A. ( 1 ). For the fiscal year 2019, the Federal Highway Administration (FHWA) requested $46.0 billion for activities that included increasing roadway safety and mobility, while improving the condition of bridges and highways ( 2 ). Decision making with regard to the use of these funds largely depends on several indicators, and for professionals in transportation, traffic flow plays an important role in this decision-making process. Traffic counts are performed to assess traffic flow. The data obtained from performing traffic counts can be used in long-term planning, network renewal and development, customer service, and operational activities ( 3 ).

Under the long-term planning category, traffic counts can be used in traffic modeling and simulation, safety studies and crash analysis, setting and measuring levels of service, and policy development ( 3 ). Activities like deriving traffic loads for bridge and pavement design, project planning and network comparison can be undertaken when using traffic counts for the purpose of network renewal and development ( 3 ). Operational uses include traffic management, network monitoring, and the derivation of traffic loads for planning and design of maintenance interventions ( 3 ). Road users can also receive responses to traffic related enquiries with the aid of these counts ( 3 ). Traffic counts are rarely used in their raw form. Performance indicators like annual average daily traffic (AADT), average daily traffic, monthly average daily traffic (MADT) and peak-hour factor are computed from the raw traffic information.

AADT is key in network planning, renewal, and development. This performance measure indicates the activity level of a road or a bridge. It estimates the average traffic volume for every day of a year at a specified roadway location. There are traditional methods of computing AADT including the simple average method, the American Association of State Highway Transportation Officials (AASHTO) method (average of averages method), and the FHWA AADT method. Several such models have been created using various predictors for the estimation of future AADTs for different roads. This research, however, focuses on AADT estimation for bridges.

Because funding is generally limited for the maintenance and retrofitting of bridges in some areas ( 4 ), it is imperative to develop AADT models which avoid dependency on predictors which may not be readily available to aid in quick and precise decision making with regard to funding allocation. Consequently, this research aims to develop Bayesian nonparametric models through machine learning for the estimation of bridge AADT. This research also seeks to explore the effect of hyperparameter value choice on the accuracy of estimations produced by Bayesian nonparametric models. Another objective of this research is to evaluate the predictions produced by both Bayesian and Frequentist approaches.

Background

Traffic data is collected by offices in state departments of transportation which may be in charge of research, or by subdivisions which may not be directly involved in research or traffic engineering ( 5 ). Traffic data is collected for maintenance, operations, forecasting, and many other purposes. Traffic counts represent a section of traffic data collected. Some types of traffic counts are pedestrian counts and total volume counts ( 6 ). In an area designated for a study, the pedestrian traffic throughout the study period is termed a pedestrian count. Total volume counts for a selected roadway are mainly used in the calculation of average daily traffic. Various methods of data collection are employed in obtaining the different traffic counts mentioned. Selected “in-situ” technologies can be used to collect traffic data with the help of detecting equipment placed alongside the roadway. The methods of data collection can be categorized as intrusive and non-intrusive ( 7 ).

A method is classified as intrusive if the main tools used are a recorder and sensor placed on or in the road. The major intrusive methods that have been employed over the years are pneumatic road tubes, piezoelectric sensors, and magnetic loops. Non-intrusive methods are done by observing traffic remotely. These methods range from very traditional methods to new technologies. Examples of non-intrusive methods include, but are not limited to: manual counts, passive and active infra-red, passive magnetic, microwave radar, ultrasonic and passive acoustic and video image detection ( 8 ). Traffic data collected using the various intrusive and non-intrusive methods are used to compute essential performance measures. AADT is one such performance measure used by transportation professionals. It is the yearly average of the number of vehicles passing a point per day. It is also simply referred to as the vehicle flow over a road section on an average day of the year ( 7 ).

Three methods that can be used to compute AADT are: simple average, AASHTO, and AASHTO with day of week, month of year adjustment factors. The simple average method can be used to compute AADT if traffic data are available for all days of the year for a specified roadway. Using the simple average method ( 9 ), the formula for calculating AADT values for a vehicle class, c, is given by Equation 1:

AAD T_{C} = \frac{1}{n} \sum_{i = 1}^{n} VO L_{i, c}

(1)

where

VOLi = total traffic on ith day of year,

n = number of days in a particular year, and

c = FHWA vehicle class.

Traffic data for each day of the year may not always be available to use in the simple average method. For this reason, AASHTO developed a second method to estimate AADT while reducing potential errors in its calculation induced by missing days. This AASHTO method makes use of any known periodicity of traffic volume by month in the year, and day of week. The AASHTO method can be used to calculate the AADT on condition that at least one day of each day of week has data for each month ( 10 ). This second method is given by Equation 2:

AAD T_{AASHTO} = \frac{1}{12} \sum_{m = 1}^{12} MAD T_{AASHT O_{m}}

(2)

MAD T_{AASHT O_{m}} = \frac{1}{7} \sum_{j = 1}^{7} [\frac{1}{n_{jm}} \sum_{i = 1}^{n_{jm}} VO L_{ijm}]

where

VOL_ijm = total traffic volume for ith occurrence of the jth day of week within the mth month

i = occurrence of a particular day of the week in a particular month (i = 1, …n_jm) for which traffic volume is available

j = day of the week (j = 1,2, …7)

m = month (m = 1,…12)

n_jm = the count of the jth day of week during the mth month for which traffic volume is available (n_jm ranges from 1 to 5 depending on day of week, month, and data availability)

From the method formulation, the AASHTO method is seen to have a mathematical bias since it gives equal weights to each day of the week in each month and then gives equal weights to all months. As it gives equal weighting to all MADT values, the AASHTO method overweights the months with fewer than 31 days and underweights the 31-day months. For this reason, even if traffic volume were available every day in each year, the AASHTO method calculation would be less likely to match the simple average ( 10 ). This bias is fixed by using the AASHTO method with day of week, month of year adjustment factors which is given by Equation 3:

AAD T_{AASHTO (MOD)} = \frac{\sum_{m = 1}^{12} d_{m} * MAD T_{AASHTO {(MOD)}_{m}}}{\sum_{m = 1}^{12} d_{m}}

(3)

MAD T_{AASHTO {(MOD)}_{m}} = \frac{\sum_{j = 1}^{7} w_{jm} [\frac{1}{n_{jm}} \sum_{i = 1}^{n_{jm}} VO L_{ijm}]}{\sum_{j = 1}^{7} w_{jm}}

where

VOL_ijm = total traffic volume for ith occurrence of the jth day of week within the mth month

i = occurrence of a particular day of the week in a particular month (i = 1, …n_wm) for which traffic volume is available

j = day of the week (j = 1, 2, …7)

m = month (m = 1, …12)

n_jm = the count of the jth day of week during the mth month for which traffic volume is available (n_jm ranges from 1 to 5 depending on day of week, month, and data availability)

w_jm = the weighting for the number of times the jth day of week occurs during the mth month (either 4 of 5); the sum of the weights in the denominator is the number of calendar days in the month (i.e., 28, 29, 30, or 31)

d_m = the weighting for the number of days (i.e., 28, 29, 30, or 31) for the mth month in the particular year.

The areas of research in relation to AADT estimation include expanding short-duration counts to yearly values, predicting AADTs, and spatially extrapolating counts from one location to another ( 11 ). This research focuses on the second area mentioned, which is AADT prediction/forecasting. Regression analysis is the most common method explored by researchers in estimating AADT (12 –15). This method has been used to predict AADT for minor roads at intersections ( 15 ), rural and urban roads ( 5 ), and some classes of highways ( 16 ). Recently, more supervised learning techniques have been explored and employed in the prediction of AADT. The use of random forest (RF) ( 15 ), and neural networks (15, 17, 18) has become quite popular in this area. Researchers have also explored AADT prediction using support vector regression with data-dependent parameters ( 5 ). Geographically weighted regression (GWR) has been used to predict traffic counts in unmeasured locations ( 19 ) and the results of published research indicate that, as compared with ordinary linear regression models, the GWR models were more accurate ( 14 ).

An efficient way of performing variable selection for regression has been employed in AADT prediction called the smoothly clipped absolute deviation penalty (SCAD). It is able to select significant predictors and estimate unknown regression coefficients simultaneously. Yang et al. ( 20 ) employed SCAD in their research and concluded from the results that it further improved the local AADT estimation when satellite information was incorporated. In addition to these methods, spatial techniques have been found to be quite useful in obtaining predictions ( 21 ). AADT predictions were made for selected roads in Texas using universal kriging ( 19 ). Universal kriging was found to reduce errors in statistically significant ways although its performance was low at sites with low traffic counts. The use of nonparametric regression models has proved to be beneficial in AADT forecasting. These regression methods do not require prior training and make predictions with reference to a group of similar cases located around the current input state at the time of prediction ( 18 ). Researchers who used nonparametric regression methods revealed that it performed well with AADTs calculated from yearly traffic counts with no missing data ( 18 ).

Gaussian Process Regression

The Bayesian approach typically starts with a parametric model that describes a phenomenon. Then a prior distribution is obtained for the unknown parameters of the model which represent our past knowledge or belief about the phenomenon before observing any data. After observing new data that is generated by the model, assumptions or beliefs can be updated ( 22 ). Bayesian methods are rooted in Bayes’ theorem to obtain a posterior probability density for the unknown parameters (Equation 4):

P (θ | y) = \frac{P (y | θ) P (θ)}{\int P (y | θ) P (θ) d θ}

(4)

In this equation, $θ$ is the vector of unknown parameters that form the model, $P (θ)$ is the prior sampling density function of $θ$ , and y is a sample drawn from the real distribution with sampling density function $P (y | θ)$ (likelihood or evidence). The normalizing constant is $\int P (y | θ) P (θ) d θ$ . Consequently, the posterior distribution for parameter θ considers the prior distribution for $θ$ and the observed data $y$ . The objective in Bayesian inference is to obtain the posterior distribution, which is the conditional distribution of θ given y. Bayesian approaches differ from Frequentist methods in this way: Bayesian inference is built on the belief that there is uncertainty in the estimation of parameter values given a fixed number of observations and this is expressed in the posterior distributions of parameters. Sampling processes like Gibbs sampling and Markov Chain Monte Carlo are extensively employed to obtain posterior distributions. Models created could either be parametric or nonparametric.

Parametric Bayesian models assume that the set of parameters is finite. Given the parameters, model predictions are independent of the observed data. The parameters are able to provide all the information about the observed data ( 22 ). The complexity of the model is therefore bounded although the amount of data is unbounded. For this reason, parametric models are known for their lack of flexibility. Examples of parametric models are: polynomial regression, mixture models, k-means, hidden Markov models, factor analysis, and logistic regression. A nonparametric Bayesian model is a Bayesian model whose parameter space is said to have an infinite dimension. In other words, it is a model that assumes a potentially infinite number of parameters. This provides high flexibility as compared with the better modeling of parametric methods. Nonparametric Bayesian methods are advantageous because they provide a simple framework for modeling complex data which has features that vary as the dataset grows ( 23 ).

To define nonparametric Bayesian models, the probability distribution of the prior in an infinite-dimensional space needs to be specified ( 24 ). A distribution on an infinite-dimensional space is a stochastic process with paths in that infinite-dimensional space. Such distributions are typically harder to define than distributions on real coordinate spaces with known dimensions. Fortunately, there are many tools obtained from stochastic process theory and applied probability that can be used to obtain the needed distributions ( 24 ). Nonparametric models include infinite latent factor models, infinite hidden Markov models, Dirichlet process mixtures, Gaussian process classifiers, and Gaussian processes. These models are applied in feature discovery, time series, clustering, classification and function approximations, respectively ( 23 ).

The nonparametric Bayesian approach selected to be utilized in this research is Gaussian process regression. Traditional regression methods are used when a particular underlying function is expected. For instance, least square methods could be used in linear regression when the underlying function is suspected to be linear. Model selection can also be used to select the best model among different polynomial regression methods based on assumptions of polynomial functions ( 25 ). Gaussian process regression provides an alternate regression approach. Instead of assuming that the underlying function f(x) relates to a specific model (e.g., f(x) = mx + c), a Gaussian process rather represents f(x) obliquely, albeit rigorously, and allows the data to “speak” clearly for itself. This means that it combines the basic structure of Bayesian inference and interpretable parameters with an ability to approximate an infinite number of functions ( 26 ). Gaussian processes have been employed to solve a wide array of engineering problems ( 27 ).

Gaussian process regression (GPR) can be classified as a supervised learning method ( 25 ). Gaussian processes expand multivariate Gaussian distributions to infinite dimensions. Generally, Gaussian processes generate data found in a domain so that any finite subset of the range should follow a multivariate Gaussian distribution ( 25 ). The mean of a Gaussian process is usually assumed to be zero everywhere. If this assumption is made, it can also be stated that a covariance function, k(x, x′) relates the observations in the dataset to each other ( 25 ). A well-known covariance function is the “squared exponential” in Equation 5,

k (x, x') = σ_{f}^{2} \exp [\frac{- {(x - x')}^{2}}{2 l^{2}}]

(5)

where

$σ_{f}^{2} = the maximum allowable covariance$

l = length parameter

k(x, x′) = covariance function.

In Equation 5, if x≈x′, then k(x, x′) approaches this maximum, implying that f(x) is almost perfectly correlated with f(x′). This is required for smooth functions where neighbors are alike. On the other hand, if x is farther from x′, then k(x, x′) ≈ 0, therefore the points will likely not “see” one other ( 25 ). This can be seen during interpolation where distant observations have little effect on new x values. The degree of separation is determined by the length parameter, introducing flexibility into the covariance function ( 25 ). During data collection and processing, errors could be introduced from the method of measurement and processing techniques to make the data noisy. The relationship between an observation y in the data and the underlying function f(x) can be represented in a Gaussian noise model ( 25 ) (Equation 6):

y = f (x) + N (0, σ_{n}^{2})

(6)

This noise can be integrated into k(x, x′), in the form (Equation 7):

k (x, x') = σ_{f}^{2} \exp [\frac{- {(x - x')}^{2}}{2 l^{2}}] + σ_{n}^{2} δ (x, x')

(7)

where δ(x, x′) = the Kronecker delta function.

Before initiating GPR, the covariance function is computed among all possible combinations of data points. The findings are then written in three matrices (Equations 8–10) ( 25 ).

K = [\begin{matrix} k (x_{1}, x_{1}) & k (x_{1}, x_{2}) & \dots & k (x_{1}, x_{n}) \\ k (x_{2}, x_{1}) & k (x_{2}, x_{2}) & \dots & k (x_{2}, x_{n}) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ k (x_{n}, x_{1}) & k (x_{n}, x_{2}) & \dots & k (x_{n}, n) \end{matrix}]

(8)

K_{*} = [\begin{matrix} k (x_{*}, x_{1}) & k (x_{*}, x_{2}) & \dots & k (x_{n}, x_{n}) \end{matrix}]

(9)

K_{* *} = k (x_{*}, x_{*})

(10)

An important assumption in Gaussian process modeling is that the data can be representative of a sample obtained from a multivariate Gaussian distribution ( 25 ). This is shown by Equation 11:

[\begin{matrix} y \\ y_{*} \end{matrix}] ~ N (0, [\begin{matrix} K & K_{*}^{T} \\ K_{*} & K_{* *} \end{matrix}])

(11)

where T = matrix transposition.

Given that there are n observations in an arbitrary data set, y = {y₁,…, y_n}, the conditional probability of interest is p(y_∗| y ). Simply explained, we must compute how likely a certain prediction for y_∗ is, given the data. This probability follows a Gaussian distribution (Equation 12) ( 25 ):

y_{*} | y ~ N (K_{*} K^{- 1} y, K_{* *} - K_{*} K^{- 1} K_{*}^{T})

(12)

The best estimate for y_∗ is given by the mean of the following distribution (Equation 13): ( 25 ):

\bar{y_{*}} = K_{*} K^{- 1} y

(13)

and the uncertainty of the computed estimate is captured by the variance in the distribution (Equation 14) ( 25 ):

var (y_{*}) = K_{* *} - K_{*} K^{- 1} K_{*}^{T}

(14)

Methodology

In this section, the four-step research method used in this paper is outlined, as shown in Figure 1. These steps are further expanded in the following subsections.

Figure 1.

Research method steps.

Data Collection/Description

As part of a new initiative by New York State Governor Andrew M. Cuomo to provide access to government data and information, the Open NY initiative was launched in March 2013 by signed Executive Order 95 ( 28 ). State agencies were directed to identify, catalog, and publish their data on the state’s open data website under this initiative, which has been administered by the Office of Information Technology Services since March 2013. The data used in this research was traffic count data for selected bridges which was made available under the Open NY initiative ( 29 ). The selected bridges in this dataset were the Rip Van Winkle, Kingston-Rhinecliff, Mid-Hudson, Newburgh-Beacon, and Bear Mountain bridges. The roads hosted by the bridges and their functional classes are listed below:

Rip Van Winkle Bridge—NY23—Urban Principal Arterial (Non-Interstate)

Kingston-Rhinecliff Bridge—NY199—Rural Principal Arterial (Non-Interstate)

Mid-Hudson Bridge—US44—Urban Principal Arterial Expressway

Newburgh-Beacon—I84—Urban Principal Arterial Interstate

Bear Mountain Bridge—US202—Urban Principal Arterial (Non-Interstate)

Annual traffic counts were provided by bridge from 1933 to 2018 for each of the selected bridges under the New York State Bridge Authority system.

Data Processing and Exploration

The traffic count data was transformed into AADT data by using the simple average method given by Equation 1. Descriptive statistics for the traffic counts and AADT of each bridge were then computed and the results are presented in Table 1. As seen in the table, the Newburgh-Beacon Bridge recorded the highest and lowest annual volumes of traffic over the study period. Thus, it also registered the highest standard deviation among all the selected bridges. The data used in the actual analysis was from 1963 to 2018 because traffic volume data was not available for all the bridges before this period.

Table 1.

Descriptive Statistics for Data

Bridge	Rip Van Winkle		Kingston-Rhinecliff		Mid-Hudson		Newburgh-Beacon		Bear Mountain
	Count	AADT	Count	AADT	Count	AADT	Count	AADT	Count	AADT
Maximum	5,968,688	16,353	8,221,926	22,526	14,586,706	39,964	26,922,078	73,759	7,905,754	21,660
Minimum	1,411,114	3,866	1,306,393	3,579	4,598,638	12,599	284,359	779	1,809,561	4,958
Mean	3,763,820	10,312	5,181,090	14,195	10,338,234	28,324	16,692,768	45,734	4,796,192	13,140
SD	1,463,544	4,010	2,237,342	6,130	3,497,573	9,582	8,068,693	22,106	1,861,892	5,101

Note: AADT = average annual daily traffic; SD = standard deviation.

To visualize the changes in traffic volume over the years, a time series heat map of all traffic volumes was plotted, as shown in Figure 2. In this plot, and with reference to the color bar at the side of the plot, it can be seen that the Rip Van Winkle, Kingston-Rhinecliff, and Bear Mountain bridges have much lower annual volumes than the Mid-Hudson and Newburgh-Beacon bridges. The actual trends over the years are not very clear for any of the bridges except the Newburgh-Beacon Bridge. To visualize trends more clearly, the data for each bridge were standardized. This process made sure the data from each bridge were comparable with the data from the other bridges, making the data internally consistent. The process of standardization creates compatibility and similarity and keeps measurement errors to a minimum ( 30 ). Standardization was also done to improve the numerical stability of the model to be created and to decrease training time. Gaussian processes assume that the underlying data is normally distributed ( 31 ). The mean of the prior distribution is often set to zero in GPR ( 32 ). Since standardization transforms data to have a mean of zero, it was found to be more appropriate than normalization for this research. The time series heat map of the standardized data is presented in Figure 3. In this figure, we can clearly see increasing traffic volumes over the years for all the bridges. The data was then split for model training (80%) and testing (20%).

Figure 2.

Time series heat map of traffic volumes from 1963 to 2018.

Figure 3.

Plot of standardized traffic volumes from 1963 to 2018.

Model Fitting

Gaussian Process Models

GPR modeling was implemented in Scikit-learn library for machine learning in Python ( 33 ). Gaussian processes are implemented in the “Gaussian Process Regressor” and “Gaussian Process Classifier” in Scikit-learn. “Gaussian Process Regressor” was employed for regression in this research ( 34 ). The prior was specified in the initial stages of implementation in Equation 15:

f (x) ~ GP (m (x), k (x, x'))

(15)

where

m(x) = mean function

k(x, x') = covariance function.

The prior mean was assumed to be the mean of the training data. The covariance function or kernel, $k (x, x')$ , was then specified by passing kernel objects. The kernels used in this research were ( 34 ):

Radial-basis function (RBF) kernel in Equation 16:

k (x_{i}, x_{j}) = \exp (- \frac{d {(x_{i}, x_{j})}^{2}}{2 l^{2}})

(16)

where

$d (\cdot, \cdot) = Euclidean distance$

$l = length - scale parameter (l > 0)$ .

Matérn kernel in Equation 17:

k (x_{i}, x_{j}) = \frac{1}{Γ (ν) 2^{ν - 1}} {(\frac{\sqrt{2 ν}}{l} d (x_{i}, x_{j}))}^{ν} K_{ν} (\frac{\sqrt{2 ν}}{l} d (x_{i}, x_{j}))

(17)

where

$K_{ν} (\cdot) = modified Bessel function$

$Γ (\cdot) = Gamma function$ .

ν = \frac{1}{2}, ν = \frac{3}{2} or ν = \frac{5}{2}

Rational quadratic kernel in Equation 18:

k (x_{i}, x_{j}) = {(1 + \frac{d {(x_{i}, x_{j})}^{2}}{2 α l^{2}})}^{- α}

(18)

where

$l = length - scale parameter (l > 0)$

$α = scale mixture parameter$ .

Exponential-sine-squared (Exp-Sine-Squared) kernel in Equation 19:

k (x_{i}, x_{j}) = \exp (- \frac{2 si n^{2} (π d (x_{i}, x_{j}) / p)}{l^{2}})

(19)

$l = length - scale parameter (l > 0)$

$p$ = periodicity parameter $(p > 0)$ .

Dot-product kernel in Equation 20:

k (x_{i}, x_{j}) = σ_{0}^{2} + x_{i} \cdot x_{j}

(20)

where $σ$ = parameter which controls the inhomogenity of the kernel.

The hyperparameters of each kernel were adjusted and then optimized during fitting of Gaussian Process Regressor. For instance $l$ and $α$ were the hyperparameters adjusted for the rational quadratic kernel (Equation 18). The hyperparameters adjusted for each kernel can be seen in Equations 16–20. Optimization was performed by maximizing the log-marginal-likelihood based on the passed optimizer. The optimizer used was “fmin_l_bfgs_b” which is also known as the limited-memory BFGS optimization algorithm. This algorithm gives an approximation of the Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS). This method of implementing Gaussian processes was to help avoid the computationally-intensive procedure of cross-validation that is used in choosing optimal hyperparameters for the model ( 35 ). Models were fit using each kernel mentioned for each of the five bridges selected for this study. During model fitting using each kernel, sample sizes were varied and the effect on predictions was recorded. It has been stated in the literature that moderate noise level could be useful in dealing with numeric issues during fitting ( 34 ), therefore noise levels for each kernel were varied and evaluated. The models for each bridge were then plotted with a 95% or two sigma confidence band.

Autoregressive Integrated Moving Average Models

The Frequentist models to be compared with the Gaussian process models were Box–Jenkins autoregressive integrated moving average (ARIMA) models. ARIMA is an interactive process which involves four stages: (a) identification, (b) estimation, (c) diagnostic checking, and (d) forecasting of time series ( 36 ). The general mathematical equation for ARIMA(p, d, q) models is shown in Equation 21:

ARIMA (p, d, q) : \hat{y_{t}} = α + β_{1} y_{t - 1} + \dots + β_{p} y_{t - p} + θ_{1} e_{t - 1} + \dots + θ_{q} e_{t - q}

(21)

where

p = number of autoregressive terms,

d = number of nonseasonal differences needed for stationarity, and

q = number of lagged forecast errors in the prediction equation.

θ = moving average parameter

y ₌the dth difference of Y;

therefore if d = 1, y_t = Y_t–Y_t–1.

Models were developed for each bridge with different combinations of p, d, and q values from 1 to 4. This resulted in the creation of over 40 models. The best models for each bridge were selected based on an evaluation metric described in the next section and compared with the GPR models.

The Akaike’s information criterion (AIC) and the Bayesian information criterion (BIC) were then computed for each model using Equations 22 and 23 ( 37 ). The selected model for each bridge was the one with the least AIC and BIC values (37, 38).

AIC = - 2 \ln (maximum likelihood) + 2 m

(22)

BIC = - 2 \ln (maximum likelihood) + m \ln (n)

(23)

where

m = number of the estimated parameters

n = number of the observations.

Model Evaluation

The entire dataset was initially split into two sub-datasets for each bridge—training (80%) and testing (20%)—by sampling randomly. The testing dataset was set aside to evaluate how close the predictions were to the actual values of AADT. The same training and testing datasets were used to develop models using ARIMA and GPR. The best ARIMA models were selected based on their AIC and BIC values. The selected model for each bridge was the one with the least AIC and BIC values.

The evaluation metrics used in this research to evaluate and compare the GPR and ARIMA models were the mean square error (MSE), root mean square error (RMSE), and mean absolute percentage error (MAPE). MSE measures the average of the squared difference between the predicted values and the actual value. It cannot be used to compare accuracy across time series with different scales. RMSE is the standard deviation of the prediction errors. MAPE measures the mean or average of the absolute percentage errors of predictions. It is commonly used because it is relatively easier to interpret and explain. Formulas for computing MSE, RMSE, and MAPE are shown in Equations 24–26, respectively:

MSE = \frac{1}{n} \sum {(y - \hat{y})}^{2}

(24)

RMSE = \sqrt{\frac{1}{n} \sum {(y - \hat{y})}^{2}}

(25)

M A P E = \frac{100 %}{n} \sum | \frac{y - \hat{y}}{y} |

(26)

where

y = actual value

$\hat{y}$ = predicted value

n = number of observations.

The best method for modeling this data was then selected based on Lewis’ scale of interpretation of estimation accuracy ( 34 ). On this scale, MAPE of less than 10% is a highly accurate forecast, 11%–20% is a good forecast, 21%–50% is a reasonable forecast, and 51% and above is an inaccurate forecast. This scale was used to inform model selection.

Results and Discussion

The Rip Van Winkle, Kingston-Rhinecliff, Mid-Hudson, Newburgh-Beacon and Bear Mountain bridges are referred to as Winkle, Kingston, Hudson, Beacon and Bear, respectively, in parts of the discussion and the tables. The priors for some kernels that were specified at the initial stages of GPR have been visualized in Figures 4 –6 showing the mean which was assumed to be the mean of the training data. The plots shown are for the Kingston-Rhinecliff Bridge.

Figure 4.

Gaussian process regression model for the Kingston-Rhinecliff Bridge using the dot product kernel.

Figure 5.

Gaussian process regression model for the Kingston-Rhinecliff Bridge using the rational quadratic kernel.

Figure 6.

Gaussian process regression model for the Kingston-Rhinecliff Bridge using the radial-basis function (RBF) kernel.

Five kernels were used to fit models using the testing dataset for each bridge, resulting in 25 models (Table 2). For each of these models, the testing data was reduced by 20% to observe the effect on the MAPE. It was observed that the overall MAPE% was decreased by 36% when this was done, highlighting the importance of collecting as much data as possible for model development. Increasing noise levels (alpha) from 0.0015 to 0.15 resulted in an overall MAPE increase of 26%. Model evaluation indicated that the MSE and MAPE values were closest for the Bear Mountain Bridge models. The dot-product kernel produced the lowest MSE value (0.025) among the Bear Mountain Bridge models. On the other hand, the MAPE value was lowest for the model created with the rational quadratic kernel (9%). The worst model created with data from the Mid-Hudson Bridge was fit using the dot-product kernel (MSE = 0.032, MAPE = 16%).

Table 2.

Evaluation of Models Created with Different Kernels

Kernel	Bear		Hudson		Winkle		Kingston		Beacon
	MSE	MAPE (%)	MSE	MAPE (%)	MSE	MAPE (%)	MSE	MAPE (%)	MSE	MAPE (%)
Radial-basis function	0.031	11	0.009	1	0.014	3	0.033	38	0.007	2
Rational quadratic	0.031	9	0.004	1	0.006	2	0.003	11	0.003	2
Exp-sine-squared	0.032	12	0.010	0	0.014	3	0.657	37	0.284	16
Dot-product	0.025	10	0.032	16	0.022	2	0.026	17	0.022	11
Matérn	0.031	10	0.002	1	0.004	1	0.001	5	0.002	2

Note: MSE = mean square error; MAPE = mean absolute percentage error.

Out of the 25 models, the one fitted with the exp-sine-squared kernel for the Mid-Hudson Bridge gave the lowest MAPE (approximately zero) indicating high accuracy. The RBF, rational quadratic and Matérn kernels all gave a MAPE of 1% for the Hudson Bridge models. The Rip Van Winkle Bridge models had higher MSE values that the Mid-Hudson Bridge models and MAPE values which were between 1% and 3%. The best Rip Van Winkle Bridge model was fitted with the Matérn kernel (MSE = 0.004, MAPE = 1%). Although the MSE values for the Kingston-Rhinecliff Bridge models were comparable with the MSE values of the other models, the MAPE metric revealed that they were the worst performing models with MAPE values between 5% and 38%. The Newburgh-Beacon Bridge models created with the RBF, rational quadratic and Matérn kernels had the lowest MSE values and had the same MAPE value of 2%. Overall the Matérn kernel produced an average MSE (0.008) which was the least value among all models and the rational quadratic kernel produced the least average MAPE (5%) for the bridge models. The exp-sine-squared model produced the highest average MAPE and MSE values across the bridge models.

AIC values of between 108 and 236 were obtained for the 14 ARIMA models generated. Their BIC values were between 115 and 296. The best ARIMA models selected for each bridge were those with the least AIC and BIC values. The best models for the Kingston-Rhinecliff, Rip Van Winkle, Mid-Hudson, Newburgh-Beacon and Bear Mountain bridges had the following configurations: MA(4), ARIMA(2, 1, 2), ARMA(2, 2), ARMA(1, 1) and ARMA(3, 1), respectively. Their respective AIC values were 115, 115, 108, 114, and 114. Their corresponding BIC values were 125, 121, 120, 115, and 135. The Mid-Hudson Bridge model was the best performing model (MAPE = 23%), and the Bear Mountain Bridge model was the worst performing model (MAPE = 76%) among the selected models. Further details about the models can be found in Table 3.

Table 3.

Best ARIMA Models for Each Bridge

Kingston - MA(4)
	Estimate	Standard error	Lower CI	Upper CI
Intercept	−0.07	0.07	−0.21	0.07
ma.L1	−0.16	0.16	−0.49	0.16
ma.L2	−0.19	0.17	−0.53	0.14
ma.L3	−0.35	0.19	−0.72	0.01
ma.L4	0.17	0.2	−0.22	0.55
Bear—ARIMA (3, 1)
	Estimate	Standard error	Lower CI	Upper CI
Intercept	0.02	0.01	0.01	0.04
ar.L1	1.25	0.5	0.27	2.23
ar.L2	0.21	0.94	−1.63	2.04
ar.L3	−0.58	0.47	−1.5	0.35
ma.L1	−1.58	3.04	−7.54	4.37
Hudson—ARIMA (2, 2)
	Estimate	Standard error	Lower CI	Upper CI
Intercept	0.01	0.02	−0.03	0.05
ar.L1	−0.11	0.2	−0.49	0.27
ar.L2	0.08	0.16	−0.23	0.39
ma.L1	−0.6	42.35	−83.6	82.4
ma.L2	−0.78	19.8	−39.6	38.03
Beacon ARIMA (1, 1)
	Estimate	Standard error	Lower CI	Upper CI
Intercept	0.04	0.02	0.02	0.07
ar.L1	0.78	0.08	0.62	0.94
ma.L1	−1.13	44100	−86400	86400
Winkle—ARIMA (2, 1, 2)
	Estimate	Standard error	Lower CI	Upper CI
Intercept	0.03	0.06	−0.09	0.15
ar.L1	−1.74	0.15	−2.03	−1.45
ar.L2	−0.87	0.12	−1.11	−0.63
ma.L1	0.87	6511.32	−12800	12800
ma.L2	−0.71	1538.01	−3015.14	3013.73

Note: ARIMA = autoregressive integrated moving average; CI = Confidence Interval; ar = autogregressive; ma = moving average; L = lag.

The RMSE, MAPE, and MSE values for the best performing GPR and ARIMA models can be found in Table 4. The results show that the RMSE and MSE values for the GPR models were far lower than the values for the ARIMA models. This difference could also be seen in the MAPE values for both models. K-fold cross-validation procedure was undertaken to divide the dataset into five non-overlapping folds. Each of the five folds was used as a held back test set, while all other folds were collectively used as a training dataset. A total of five models were fit and evaluated on the hold-out test sets and the mean performance was reported. The results did not differ by much from the original results obtained before cross-validation (Table 4).

Table 4.

Comparison of GPR and ARIMA Models

Bridge	GPR			ARIMA
	MSE	RMSE	MAPE (%)	MSE	RMSE	MAPE (%)
Bear	0.0311	0.1764	9	1.0112	1.0056	76
Hudson	0.0023	0.0480	1	0.7134	0.8446	23
Winkle	0.0042	0.0648	1	0.9031	0.9503	33
Kingston	0.0011	0.0332	5	0.8221	0.9067	26
Beacon	0.0024	0.0490	2	0.9203	0.9593	62

Note: ARIMA = autoregressive integrated moving average; GPR = Gaussian process regression; MSE = mean square error; RMSE = root mean square error; MAPE = mean absolute percentage error.

Overall, the Mid-Hudson Bridge models from both methods performed best (GPR: 1%, ARIMA: 23%). The worst performing models were the Bear Mountain Bridge models (GPR: 9%, ARIMA: 76%). The best performing ARIMA model had a 22% higher MAPE than the best GPR model. All GPR models had MAPE of below 10%. On Lewis’ scale of interpretation of estimation accuracy, this indicates highly accurate forecasting power. The MAPEs of the Kingston-Rhinecliff, Rip Van Winkle, and Mid-Hudson ARIMA models fell between the Lewis’ predetermined range of 20%–50%, indicating reasonable forecasting power. The MAPEs of the Newburgh-Beacon and Bear Mountain bridge ARIMA models were greater than 50%, classifying them as models with weak forecasting power. The highest MAPE recorded for both models (76%) was far less than what was recorded in a similar study conducted by Wu and Xu ( 15 ), whose highest MAPE was 291.1%. The MAPE range for the GPR models (1%–9%) was also much lower than was recorded by Xia et al. ( 12 ) which was 1.31%–57%. When all models are considered, the MAPE values obtained in this research had a wider range than the aforementioned study (1%–76%). The overall MAPE for the GPR models was 3.6% and that of the ARIMA models was 44%. The data used for the models was collected from 1963 to 2018. Additional predictions were, however, made for the years 2019 and 2020 using the best models for each bridge. The predictions were made with respect to the last AADT recorded (2018). These predictions were made ceteris paribus. The changes in AADT for 2019 and 2020 respectively were as follows: Kingston-Rhinecliff Bridge (1.9% and 3.4%), Rip Van Winkle Bridge (3.7% and 7%), Mid-Hudson Bridge (2% and 2.7%), Bear Mountain Bridge (–26% and –61%), Newburgh-Beacon Bridge (3% and 5.5%).

Concluding Remarks

This study developed models to estimate AADT of selected bridges in New York using GPR and ARIMA methods. The GPR and ARIMA models were assessed and it was determined that the GPR models provided better forecasts than the ARIMA models. The researchers therefore recommend the use of GPR when estimating AADT for these bridges, especially if a high forecasting accuracy needs to be achieved. The advantage in using the proposed method is that extra predictors like functional class, lane configuration, and so forth that are needed for other methods whose predictions rely on the use of multiple variables, are not needed to estimate AADT. The major contribution of this paper is the development of GPR AADT forecasting models exclusively for these selected bridges which, to the best of the authors’ knowledge, has not been undertaken by any other research study. Based on the performance of the GPR model, the authors cautiously conclude that this method is adequate for AADT estimation for bridges that host principal arterial roads. This is because of the small sample size of the bridges studies (n = 5). Although good results were obtained using GPR, the authors suggest that future studies should be conducted on data from additional bridge sites with larger sample sizes to ascertain the potency of this method in AADT estimation. Bridges that host roads of other functional classes should be included in these studies. Also, future studies should be conducted using AADTs computed with methods other than the simple average method.

Footnotes

Author Contributions

The authors confirm contribution to the paper as follows: study conception and design: Grace Ashley, Nii Attoh-Okine; data collection: Grace Ashley; analysis and interpretation of results: Grace Ashley, Nii Attoh-Okine; draft manuscript preparation: Grace Ashley, Nii Attoh-Okine. Both authors reviewed the results and approved the final version of the manuscript.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

U.S. Census Bureau. 2020Census Count Guides Funding of New Roads and Bridges. https://www.census.gov/library/stories/2019/12/2020-census-count-guides-funding-of-new-roads-and-bridges.html. Accessed July 7, 2020.

Federal Highway Administration. Federal Highway Administration (FHWA) Fiscal Year 2019 Budget. U.S. Department of Transportation. https://www.fhwa.dot.gov/cfo/fhwa-fy-2019-cj-final.pdf. Accessed July 9, 2020.

Rodwell

Guidelines for Traffic Counting. Roading Infrastructure Management Support Group (RIMS) Body of Knowledge, 2013, pp. 9–11.

MnDOT Research Services. Local Bridge Removal Policies and Programs. Transportation Research Services, Minnesota Department of Transportation, 2017.

Castro-Neto

Jeong

M. K.

Han

L. D.

AADT Prediction Using Support Vector Regression with Data-Dependent Parameters. Expert Systems with Applications, Vol. 36, No. 2, 2009, pp. 2979–2986. https://doi.org/10.1016/j.eswa.2008.01.073.

Mooney

S. J.

DiMaggio

C. J.

Lovasi

G. S.

Neckerman

K. M

Bader

M. D. M.

Teitler

J. O.

Sheehan

D. M.

Jack

D. W.

Rundle

A. G.

Use of Google Street View to Assess Environmental Contributions to Pedestrian Injury. American Journal of Public Health, Vol. 106, No. 3, 2016, pp. 462–469. https://doi.org/10.2105/AJPH.2015.302978.

Leduc

Road Traffic Data: Collection Methods and Applications. Working Papers on Energy, Transport and Climate Change, Vol. 1, No. 55, 2008, pp. 1–55.

Federal Highway Administration. Highway Performance Monitoring System (HPMS). https://www.fhwa.dot.gov/policyinformation/hpms/volumeroutes/ch4.cfm. Accessed July 24, 2020.

Grande

Wood

Ominski

Regehr

J. D.

Evaluating Annual Average Daily Traffic Calculation Methods with Continuous Truck Traffic Data. Transportation Research Record: Journal of the Transportation Research Board, 2017. 2644: 30–38.

10.

Federal Highway Administration. Assessing Roadway Traffic Count Duration and Frequency Impacts on Annual Average Daily Traffic (AADT) Estimation: Assessing Accuracy Issues with Current Known Methods in AADT Estimation from Continuous Traffic Monitoring Data. FHWA, U.S. Department of Transportation, Washington, D.C., 2014.

11.

Lowry

Dixon

GIS Tools to Estimate Average Annual Daily Traffic. National Institute for Advanced Transportation Technology, Moscow, ID, 2012.

12.

Xia

Zhao

Chen

Shen

L. D.

Ospina

Estimation of Annual Average Daily Traffic for Nonstate Roads in a Florida County. Transportation Research Record: Journal of the Transportation Research Board, 1999. 1660: 32–40.

13.

Zhao

Chung

Estimation of Annual Average Daily Traffic in a Florida County Using GIS and Regression. Transportation Research Record: Journal of the Transportation Research Board, 2001. 1769: 113–122.

14.

Zhao

Park

Using Geographically Weighted Regression Models to Estimate Annual Average Daily Traffic. Transportation Research Record: Journal of the Transportation Research Board, 2004. 1879: 99–107.

15.

Annual Average Daily Traffic Prediction Model for Minor Roads at Intersections. Journal of Transportation Engineering, Part A: Systems, Vol. 145, No. 10, 2019, p. 04019041. https://doi.org/10.1061/JTEPBS.0000262.

16.

Wang

Kockelman

K. M.

Forecasting Network Data Spatial Interpolation of Traffic Counts from Texas Data. Transportation Research Record: Journal of the Transportation Research Board, 2009. 2105: 100–108.

17.

Khan

M. Z.

Khan

S. M.

Chowdhury

Dey

Development and Evaluation of Recurrent Neural Network Based Models for Hourly Traffic Volume and AADT Prediction. Transportation Research Record: Journal of the Transportation Research Board, 2018. 2673: 489–503.

18.

Tang

Y. F.

Lam

W. H. K.

P. L. P.

Comparison of Four Modeling Techniques for Short-Term AADT Forecasting in Hong Kong. Journal of Transportation Engineering, Vol. 129, No. 3, 2003, pp. 271–277. https://doi.org/10.1061/(ASCE)0733-947X(2003)129:3(271).

19.

Selby

Kockelman

K.M.

Spatial Prediction of Traffic Levels in Unmeasured Locations: Applications of Universal Kriging and Geographically Weighted Regression. Journal of Transport Geography, Vol. 29, 2013, pp. 24–32.

20.

Yang

Wang

S. G.

Bao

Efficient Local AADT Estimation via SCAD Variable Selection Based on Regression Models. Proc., 2011 Chinese Control and Decision Conference (CCDC), Mianyang, China, IEEE, New York, 2011, pp. 1898–1902.

21.

Klatko

T. J.

Saeed

T. U.

Volovski

Labi

Fricker

J. D.

Sinha

K. C.

Addressing the Local-Road VMT Estimation Problem Using Spatial Interpolation Techniques. Journal of Transportation Engineering, Part A: Systems, Vol. 143, No. 8, 2017, p. 04017038. https://doi.org/10.1061/jtepbs.0000064.

22.

Rocca

. Bayesian Inference Problem, MCMC and Variational Inference. Towards Data Science. https://towardsdatascience.com/bayesian-inference-problem-mcmc-and-variational-inference-25a8aa9bce29. Accessed July 11, 2020.

23.

Ghahramani

Parametric vs Nonparametric Models. Tübingen, Germany, 2015.

24.

Orbanz

Lecture Notes on Bayesian Nonparametrics. Journal of Mathematical Psychology, Vol. 56, 2014, pp. 1–12.

25.

Ebden

Gaussian Processes for Regression: A Quick Introduction. Robotics Research Group, University of Oxford, 2008.

26.

Belyaev

Burnaev

Kapushev

Gaussian Process Regression for Structured Data Sets. In Proc., Lecture Notes in Computer Science: International Symposium on Statistical Learning and Data Sciences ( A.

Gammerman

Vovk

Papadopoulos

, eds.), Vol. 9047, Egham, April 20–23, 2015, Springer, Cham, pp. 106–115.

27.

Williams

Gaussian Processes for Machine Learning. ANC, Harrison, NY, 2007.

28.

Office of Information Technology Services. New York State Open Data Program Overview. NY Open Data, pp. 1–6. https://data.ny.gov/dataset/NYS-Open-Data-Program-Overview/7ruf-yihe. Accessed July 11, 2020.

29.

State of New York. Traffic Counts for Select Bridges: Beginning1933. https://catalog.data.gov/dataset/traffic-counts-for-select-bridges-beginning-1933. Accessed July 9, 2020.

30.

Muralidharan

Sigma

A Note on Transformation, Standardization and Normalization. International Journal of Operations and Quantitative Management, Vol. 9, No. 1&2, 2016, pp. 116–122.

31.

Bilogur

. Gaussian Process Regression and Classification. Kaggle. https://www.kaggle.com/residentmario/gaussian-process-regression-and-classification. Accessed October 19, 2020.

32.

Schulz

Speekenbrink

Krause

A Tutorial on Gaussian Process Regression: Modelling, Exploring, and Exploiting Functions. Journal of Mathematical Psychology, Vol. 85, 2018, pp. 1–16. https://doi.org/10.1016/j.jmp.2018.03.001.

33.

Pedregosa

Varoquaux

Gramfort

Michel

Thirion

Grisel

Blondel

Prettenhofer

Weiss

Dubourg

Vanderplas

Passos

Cournapeau

Brucher

Perrot

Duchesnay

Scikit-Learn: Machine Learning in Python. Journal of Machine Learning Research, Vol. 12, 2011, pp. 2825–2830.

34.

Scikit-learn. Gaussian Processes. scikit-learn 0.23.1. https://scikit-learn.org/stable/modules/gaussian_process.html. Accessed July 13, 2020.

35.

Fonnesbeck

. Fitting Gaussian Process Models in Python. Data Science Blog by Domino. https://blog.dominodatalab.com/fitting-gaussian-process-models-python/. Accessed July 13, 2020.

36.

Omar

Hoang

V. H.

Liu

D.-R.

A Hybrid Neural Network Model for Sales Forecasting Based on ARIMA and Search Popularity of Article Titles. Computational Intelligence and Neuroscience, Vol. 2016, 2016. https://doi.org/10.1155/2016/9656453.

37.

Wang

Liu

Comparison of Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) in Selection of Stock-Recruitment Relationships. Fisheries Research, Vol. 77, 2006, pp. 220–225. https://doi.org/10.1016/j.fishres.2005.08.011.

38.

Bozdogan

Model Selection and Akaike’s Information Criterion (AIC): The General Theory and its Analytical Extensions. Psychometrika, Vol. 52, No. 3, 1987, pp. 345–370. https://doi.org/10.1007/BF02294361.