Random-Parameter Multivariate Negative Binomial Regression for Modeling Impacts of Contributing Factors on the Crash Frequency by Crash Types

Highways provide the basis for safe and eﬃcient driving. Road geometry plays a critical role in dynamic driving systems. Contributing factors such as plane, longitudinal alignment, and traﬃc volume, as well as drivers’ sight characteristics, determine the safe operating speed of cars and trucks. In turn, the operating speed inﬂuences the frequency and type of crashes on the highways. Methods . Independent negative binomial and Poisson models are considered as the base approaches to modeling in this study. However, random-parameter models reduce unobserved heterogeneity and obtain higher dimensions. Therefore, we propose the random-parameter multivariate negative binomial (RPMNB) model to analyze the inﬂuence of the traﬃc, speed, road geometry, and sight characteristics on the rear-end, bumping-guardrail, other, noncasualty, and casualty crashes. Subsequently, we compute the goodness-of-ﬁt and predictive measures to conﬁrm the superiority of the proposed model. Finally, we also calculate the elasticity eﬀects to augment the comparison. Results . Among the signiﬁcant variables, black spots, average annual daily traﬃc volume (AADT), operating speed of cars, speed diﬀerence of cars, and length of the present plane curve positively inﬂuence the crash risk, whereas the speed diﬀerence of trucks, length of the longitudinal slope corresponding to the minimum grade, and stopping sight distance negatively inﬂuence the crash risk. Based on the results, several practical and eﬃcient measures can be taken to promote safety during the road design and operating processes. Moreover, the goodness-of-ﬁt and predictive measures clearly highlight the greater performance of the RPMNB model compared to standard models. The elasticity eﬀects across all the models show comparable performance with the RPMNB model. Thus, the RPMNB model reduces the unobserved heterogeneity and yields better performance in terms of precision, with more consistent explanatory power compared to the traditional models.


Introduction
e highway is a dynamic system that includes drivers, vehicles, roads, and environment. Drivers play important roles, and inattention, distraction, and misjudgment are the main causes of highway crashes. Nowadays, research studies mainly focus on drivers. Meanwhile, the road characteristics, which serve as an important basis of the dynamic system and influence safety, are often neglected [1][2][3]. As the basis of driving, road geometry has a very important impact on safety. Almost 40% of all crashes attributed to driver mistake or vehicle failure are fundamentally caused by road geometry [4]. Babukov et al. studied the internal relationship between plane, profile parameters, and road safety in the design of highways and made significant suggestions for improving safety [4][5][6][7]. McGee et al. established statistical models between the curve length, curvature, traffic volume, and the crash frequency [8].
After studying 141,812 crashes, Pei and Ma focused on the relationship between road design factors (plane, vertical, cross section, and intersection) and crashes and put forward corresponding countermeasures [9]. Many types of research have evaluated the safety of road geometry by visualization and computer simulation. Hassan et al. simulated the influence of the 3D environment on driver vision and analyzed the influence of a combination of the horizontal and vertical curve on road safety [10][11][12][13][14]. en, the Federal Highway Administration developed the Interactive Highway Safety Design Model (IHSDM). Software used the lane width, lane number, horizontal alignment, vertical alignment, cross section, and superelevation as the evaluation parameters to check the consistency of the design [15][16][17]. Drivers often attempt to improve the efficiency of transportation when road and traffic conditions allow speeding, which causes many crashes on the highway. e Solomon and the Monash University Accident Research Center modeled the relationship between crash rate and vehicle speed and found that the higher the difference between the operating speed and the average speed, the higher the crash rate [18,19]. e average speed is often exceeded because of the favorable geometry and environmental conditions of modern highways. Meanwhile, the speed difference between cars and trucks increases owing to the relatively poor power and braking performance of trucks, which lead to frequent rearend crashes. e initial road alignment design was based on the design speed, which was first proposed by the American Association of State Highway and Transportation Officials (AASHTO) as the key factor in road design. At the same time, designers and researchers also used design speed to evaluate the consistency of horizontal and longitudinal profiles. e expected speed would greatly exceed the design speed during favorable road conditions, and the traditional design method based on the design speed led to poor continuity and imbalance in actual driving. In addition, the drivers' visual perception expressed poor consistency with actual vehicle control [20]. Due to these limitations, AASHTO later adopted the operating speed to dominate the design, which is defined as the driving speed at the 85th percentile selected by the actual measurement of mediumskilled drivers under a free-flow state with good weather and road conditions [21]. Lamm et al. studied the influence of curve on the operating speed based on the design data, vehicle speed, traffic volume, and crash data of 261 intercontinental highway sections and determined the model of the curve radius and operating speed [22,23]. Anderson et al. investigated the influence of curve radius on the operating speed according to the speed attenuation while passing the curve [24,25]. Krammes et al. established the operation speed model by collecting the operating speed and design data of 138 plane curve sections based on 1,126 curve sections as samples for a preliminary evaluation to determine the causes of accidents in curve sections and select improvement methods [26]. Collins et al. evaluated the freeflow condition by measuring the actual driving speed at the midpoint of the curve [27][28][29][30][31][32][33]. Sil and Maji measured the operating speed through the curve on three sections and established a model between the speed at the curve and the combined curves on the four-lane highway [34]. Based on the existing research abroad, the operating speed model of different vehicle types has been adopted in research in China. In 2004, the Guidelines for Safety Evaluation of Highway Projects proposed the highway speed model based on many research datasets. Specifications for the Highway Safety Audit published in 2015 optimized and improved the speed model of low-grade highways [35]. e evaluated sections are divided into tangent sections, longitudinal slope sections, horizontal curve sections, curved slope combination sections, tunnel sections, interchange sections, and other sections for calculation according to the radius of the horizontal curve and the slope of the profile. A sufficient driving sight distance is directly related to the safety and speed on the highway and is an important index of evaluation. e stopping sight distance enables drivers to stop, meet, and overtake smoothly, and the transverse clear distance measures the safe zone allowed for lateral offset during driving. Both distances are measurable representations of sight characteristics [35].
In summary, previous research studies focused on the contributing factors of road geometry or operating speed characteristics separately, and the sight characteristics determined by the road geometry alignments were rarely studied.
is study takes the geometric, speed, and sight characteristics into consideration simultaneously under a single condition and then adopts the random-parameter approach to reduce the interactions between various variables to acquire a higher dimension and discover the contributing factors. e traditional negative binomial and Poisson models used in previous research studies become complicated when analyzing crash types. In contrast, this random-parameter multivariate negative binomial (RPMNB) model allows for a more simplified framework and reduces the unobserved heterogeneity with higher accuracy. In this study, we collect three-year (2009-2011) statistical crash data of the Beijing-Shanghai Highway as a sample to analyze whether the road geometry, operating speed, traffic volume, and sight characteristics affect the crash types (rear-end, bumping-guardrail, roll-over, etc.).

Methodology
Many methodological approaches, such as the multivariate Poisson (lognormal) model, zero-inflated negative binomial model, and Poisson lognormal spatial and/or temporal model typically address the crash rate considering the number of crashes occurring over a roadway segment or at an intersection [36][37][38].
Generally, some of the significant factors affecting the crash rate are missed in the collected data or difficult to analyze. ese factors (called unobserved heterogeneity) make a variation in the influence on the observed factors, which may lead to misspecified parameters and erroneous inferences. To reduce the influence of the unobserved 2 Discrete Dynamics in Nature and Society heterogeneity, random parameters, latent-class (finitemixture) models, and Markov switching models are considered in the multivariate models [39].

Poisson (Negative Binomial) Regression Model.
In statistics researches on crash frequency, Poisson regression is a generalized linear model analysis used to model crash data. e response variable Y is assumed yielding to Poisson distribution, with the expected value modeled by a linear combination of unknown parameters. In Poisson regression, the Poisson incidence rate μ is determined by X k (the regressor variables) [40][41][42]: e fundamental Poisson regression model (PRM) for an observation i is written as where t is the exposure of a period of time and β 1 , β 2 , . . . , β k are coefficients estimated as unknown parameters. Negative binomial (NB) regression is a common generalization of Poisson regression including a gamma noise variable [43]. is model is very popular because it models the Poisson heterogeneity with a gamma distribution, and the variance is not equal to the mean restrictively. e negative binomial regression model (NBRM) meets the equation When α is significantly different from 0, the negative binomial regression is appropriate. Otherwise, the Poisson model is better.

Zero-Inflated Poisson Model.
When the crash data contain excess zero-count values in the model, the wellknown zero-inflated Poisson model is adopted.
where y i is the nonnegative integer value; λ is the i th expected Poisson count; and π is the probability of extra zeros.

Random Parameters for Unobserved
Heterogeneity. e PRM restricts the mean to be equal to the variance (E � VAR), and the PRM model does not fit well in some cases. When the model does not hold the equality, the data may be overdispersed (E < VAR) or underdispersed (E > VAR), and the standard errors of the estimated parameter of the PRM will be incorrect. To account for overdispersion in the crash count data, PRM is promoted and derived.
where e ε ik are error terms following the gamma distribution with mean 1 and variance α.
In response to the nonconstant explanatory variables in the models, we developed the random parameters in each estimated parameter to account for unobserved heterogeneity [39].
where β lk denotes the l th explanatory variable for observation l; b k are the mean parameter estimates; and φ lk is a randomly distributed term capturing unobserved heterogeneity.

Multivariate Negative Binomial Model.
In recent research studies, the general framework of the random-parameter multivariate negative binomial (RPMNB) model is mostly proposed with the expected number of crashes [44,45], λ, in the i th road segment and k th crash types (in this paper, these types are rear-end, bumping-guardrail, and crash): where X ik � (1, X 1k , X 2k , . . . , X Nk ) is the independent variable vector; β lk � (β 0k , β 1k , . . . , β Nk ) is the coefficient vector; and ε ik is the multivariate error term distributed with zero mean and variance θ 2 . In the RPMNB model, the correlation ρ follows the unstructured correlation covariance matrix, which represents the correlation between ε ik of models for crash type a and b: Discrete Dynamics in Nature and Society

Model Comparison and Evaluation
2.4.1. Goodness of Fit. When analyzing the structure of the crash data model across different crash types, it is important and necessary to confirm the structure of the unobserved parameters. e likelihood ratio test is used to assess the models [46].
where β R and β U denote the log-likelihood at the convergence of the restricted and unrestricted model, respectively. Bayesian information criterion (BIC) is also used for model comparison, which is a generalized version of the Akaike information criterion (AIC) considering the Bayesian equivalent [47]: where k and L denote the number of model parameters and the likelihood function, respectively. AIC introduces the penalty term to minimize the parameters of the model, which helps to reduce the possibility of overfitting and promote the degree of model fitting (maximum likelihood). BIC considers the number of samples, leading to a larger penalty term than AIC. When the number of samples is too large, BIC can effectively prevent the excessively high complexity caused by the precision of the model. So, the models with smaller AIC and BIC values perform better.

Prediction Accuracy.
Other than BIC, we used root mean square error (RMSE) to evaluate the accuracy of the models. Similar to BIC, smaller RMSE values mean the model predicts more accurately.
At the same time, the mean absolute error (MAE) and mean absolute percentage error (MAPE) are also used to estimate the accuracy of models: where O j is the j th observation value; P j is the j th predicted value; and n 0 is the number of observations.

Data Description
Yearly crash data from the Beijing-Shanghai Highway from Xinyi to Jiangdu, from 2009 to 2011, were collected from the traffic management department. Originally, the data were used to record the emergency vehicle on the highway. e crash data contained 3,293 crashes in detail, including crash type, location, time, vehicle type, climate, road surface condition, and casualty condition. Rear-end and bumpingguardrail crashes made up 52.9% and 30.3% of the crash data, respectively. e rear-end and bumping-guardrail crash type are the focus of this research, but the other crash types (roll-over, fire, scrub, and others) are also utilized in the models to better explore the correlations to the first two crash types. Table 1 summarizes descriptive statistics of the five crash types in each section on the highway.
Besides, roadway geometric data are collected by road design and construction drawings, which include road geometric features. e undivided four-lane highway from Xinyi to Jiangdu in this study is 259.5 km in total, with the design speed of 120 km/h, and includes 27 interchanges and rest areas. We divide the road into 426 different sections according to different horizontal alignment, vertical alignment, and interchange both-way. We obtain the AADT of 426 sections reported by roadway management agencies.
We adopt the modified crash frequency method to identify the crash black spots and calculate the average crash number of each section: e critical value of the crash number is determined with the confidence level of 95%: where λ is the average crash number; m i is the crash number in the i th section; n is the total number of sections; and R is the critical value of the crash number. e actual crash number of each section is compared with the critical value, and the section is identified as crash black spots if the actual number is greater than R.
Based on the models in the Specifications for Highway Safety Audit published in 2015, we calculated the operating speeds of cars and trucks by segments according to different geometric features [35]. e stopping sight distance means the shortest distance required for one ordinary driver to react and slow down or stop when encountering obstacles while driving at a certain speed. Based on the Guidelines for Design of Highway Grade-separated Intersections [48], the stopping sight distances of the car and truck are equations (16) and (17), respectively.
e truck drivers can see the vertical plane of the obstacles at a considerable distance from the point of view with the low speed, but this advantage is not enough to make up 4 Discrete Dynamics in Nature and Society for the poor braking performance. Despite the high viewpoint, truck drivers also lose their sight in places with limited lateral line of sight, especially where S car and S truck are the stopping sight distance of the car and truck, respectively; v 85 is the the operating speed (km/h); t is the reaction time, takes 2.5 s generally (judging time as 1.5 s and running time as 1.0 s); g is the gravitational acceleration, takes 9.8 m/s 2 ; i is the longitudinal grade; and f is the longitudinal friction coefficient between truck tires and the road surface which takes 0.17 generally. Corrugated beam guardrails are set in the middle and beside the road across all sections, which affect the driverss ight. We consider the biggest transverse clear distance to confirm the sight safety, which means the distance between the curve of sight and the track. When the plane curve is sharp, the transverse clear distance should be judged on the inside lane. We calculate the required stopping sight distance of each section for safety following the equation where H is the biggest transverse clear distance; R s is the radius of the inside lane; and c is the central angle of line of sight.
As shown in Table 2, we summarize the traffic, speed, geometric, and sight characteristics as the independent variables.  Tables 3 and 4, and Appendix provides the results of corresponding traditional models.

Results and Discussion
Subsequently, each base effect is estimated for common exogenous variables across five crash types, and we estimate the deviation of variables versus the base for each crash type. e corresponding t-statistic will present statistically significant if the deviation term offers a difference from the base effect. Based on the t-statistic, the parameter does not reveal differential sensitivity for the base crash type if the variable is statistically insignificant.

Model Estimation Results.
In this section, we conduct a detailed discussion of the significant factors affecting crash count components on different crash types. e model estimation results for the RPMP model and RPMNB model are shown in Tables 3 and 4, respectively.
Similar to the traditional approach, we presented the individual effect of each exogenous variable accommodating to the crash propensity. Taken the constants estimated in various crash type propensity equations as an example, the effect of the casualty crash in Table 3 can be computed as 1.867 (base of the rear-end crash) ±1.000 (casualty crash deviation) � 0.867. At the same time, we identify the number of significant parameters estimated across the five crash types and estimate one individual parameter across 5 crash types. e positive value of the variable in Table 3 indicates there will be more crashes with the increase of the variable, and less crashes, otherwise. As shown in Tables 3 and 4, the significant factors in the two models are not exactly the same. en, we focus on the results of the RPMNB model while analyzing the effects of factors, and the differences in the RPMP model are appended.

Constant.
e constant represents the intercept of the crash type with exogenous variables, without any substantive interpretation.

Traffic Characteristics.
Whether there is an interchange in the section is found to negatively influence other, noncasualty, and casualty crashes indicating a lower likelihood of crash propensity for these three crash types in interchange sections. Compared with the general sections, more vehicles accelerate, decelerate, and change lanes leading to more interweaving areas in the interchange sections, getting more attention from drivers and causing fewer crashes to some extent. e results are consistent with the earlier research on the highway interchange [48]. And the impact is also found not significantly different for rear-   Discrete Dynamics in Nature and Society end and bumping-guardrail crashes, which shows that these two crash types are not associated with the interchange areas. e corresponding parameter of black spots offers a positive impact on crash occurrence for rear-end, bumpingguardrail, other, noncasualty, and casualty crashes revealing a higher likelihood of crashes with the increased proportion of black spots. More crashes across all the crash types occurred in the black spots than other sections. Similarly, the AADT offers a positive influence across the five crash types in the RPMNB model indicating more crashes occur with more traffic. Greater traffic volume results in higher crash risk [49]. en, the AADT indicates no significance on the casualty crashes in the RPMP model.

Speed Characteristics.
In terms of operating speed, the estimated results of both RPMP model and RPMNB model show that the higher operating speed of cars is likely to result in less crash risk across five crash types. e finding is expected because drivers of cars are likely to drive at a higher speed under good traffic conditions in the sections with great geometric alignment. e estimated results show that the speed difference of cars with the adjacent segment shows the same positive effect on all the crash types. e higher speed difference of cars between adjacent segments leads to higher crash risk, and this finding is corresponding to the Solomon curve [18]. Oppositely, the speed difference of trucks with the adjacent segment shows the same negative effect on all the crash types except the casualty crash. e result of trucks is contrary, and the estimated variable of the speed difference of trucks shows a negative influence on the rear-end, bumping-guardrail, other, and noncasualty crashes and an insignificant effect on the casualty crashes. e reasonable explanation is that the poor power and braking performance of trucks limit the operating speed. When the design speed reaches a certain value, the operating speed will not increase with the improvement of the geometric alignment. And the operating speed restricted by the poor geometric alignment of the sections reduces the crash risk. e operating speed of trucks has a positive influence on the risk of rear-end crashes and a negative influence on the noncasualty crashes. e higher operating speed of trucks leads to more rear-end crashes and less casualty in the crashes.
Relative to the cars, the trucks show poor power and braking performance. e greater operating speed results in a higher propensity of rear-end crashes and a smaller difference between cars and trucks, which lead to less collision impact and damage in the crashes.

Geometric Characteristics.
e length of the present plane curve presents a positive effect on the five crash types, which is perhaps indicating that the longer length of the plane curve results in more crash risk propensity. e reasonable explanation for this result is that the drivers will adapt to the curvature of the curve gradually, causing the fatigue or distraction of drivers. Otherwise, this adaption causes the drivers to accelerate to the speed not suitable for the geometric alignment, leading to more crashes. e Discrete Dynamics in Nature and Society positive effect of the radius of the present plane curve on the rear-end and casualty crashes shows the same tendency. Concerning the radius of the plane curve of the front section, the variable is found to have a positive influence on the rear-end and noncasualty crashes with relatively little influence coefficients (0.000000359 and 0.000000365) in the RPMNB model. However, the estimated results are not found in the RPMP model. e vehicles will run at an uncoordinated speed when transiting to a smaller radius and increasing the rear-end crash propensity. Similarly, the coefficient of the length of the plane curve of the front section shows a positive influence on the casualty crashes. e length of the longitudinal slope corresponding to the minimum grade is also found to have a negative influence on the probability of rear-end, noncasualty, and casualty crashes in the RPMNB model. And in the RPMP model, only the rear-end and noncasualty crashes show significance. Among the road sections with the small longitudinal grade, the trucks show greater speed with less speed difference from cars, which promotes safety. In the RPMNB model, the length of the longitudinal slope corresponding to the maximum grade also presents a negative influence on the noncasualty crashes. However, the results are not present in the RPMP model.

Visibility Characteristics.
As expected, the coefficient associated with the stopping sight distance of trucks offers a negative effect on the crash risk of all the crash types. e likelihood of all single-vehicle crashes will reduce with higher stopping sight distance of trucks. e visibility field and distance in front of the vehicle are important for safe and effective driving on the highway. e speed and direction of the vehicles depend on the visibility of the ahead road and surrounding environment. erefore, the higher stopping sight distance of the truck facilitates the drivers to control the direction more accurately, leading to fewer crashes. Meanwhile, the transverse clear distance of cars indicates negative effects on the rear-end, bumpingguardrail, and noncasualty crashes in the RPMP model. And the transverse clear distance of cars only shows a negative influence on bumping-guardrail and noncasualty crashes in the RPMNB model. e more transverse safe zone allows for a higher tolerance for vehicles, which leads to fewer crashes.

Model Comparisons.
Statistical measures such as loglikelihood, AIC, and BIC values are calculated to measure the goodness of fit in the estimated models, and the results are shown in Table 5.
Based on the comparison of AIC and BIC values, we firstly find the models with zero-inflated effects are not suitable for the crash data in this study. Secondly, the models considering effects caused by unobserved heterogeneity perform better than the independent models. Finally, the proposed models estimate efficient and accurate parameters in the parsimonious model systems. To assess the predictive accuracy of the estimated models, we used RMSE (root mean square error), MAE (mean absolute error), and MAPE (mean absolute percentage error) for discussion. Table 6 presents the values of these measures.
Despite the difference in the total number of parameters between the two models (56 vs. 52), the performance of the two models across five crash types is quite similar. From the total values of the three measures, we can observe that the RPMNB model performs better than the RPMP model. From the perspective of different crash types, the RPMP model performs marginally better than the RPMNB model for the measures with respect to bumping-guardrail and noncasualty crashes, and RPMNB performs better in terms of rear-end, other, noncasualty, and casualty crashes. e difference in the number of parameters across the two models makes no effects on the deviation measures.

Elasticity Effects.
To provide more insight and explain the marginal effects of the exogenous variables, the elasticity effect is computed for RPMNB and RPMNB across all the crash types. Elasticity effect denotes an estimate of the effect of a variable on the expected frequency assuming all the other variables take the average values. e elasticity effect is the effect on the expected frequency λ i of a 1% change in the variable, which is as follows: where x ik is the value of the k th independent variable for observation i; β k is the estimated parameter for the k th independent variable; and λ i is the expected frequency for observation i. As shown in Figure 1, there is not any large difference in the elasticity effects of the two models across five crash types. And almost half of the number of the variables make very little effects (41 out of 79).
We observe that there are significant differences in the elasticity effects for different crash types across different models. With respect to the length of the present plane curve for the rear-end crashes, the RPMNB model predicts a 27.4% increase, while a 52.7% increase from the RPMP. Most likely, the nonlinearity of the two models results in such differences.

Conclusions
e empirical analysis was based on the traffic crash data from the Beijing-Shanghai Highway for the year 2009-2011 and the road design and construction drawings including the road geometric features. e road geometric alignment, operating speed, traffic volume, and sight distance are considered as the exogenous variables across the traditional and proposed models. e traditional negative binomial and Poisson models are commonly adopted in the analysis on crash count frequency to analyze the impacts of various factors. To reduce the influence of the unobserved heterogeneity and get higher dimensions, we adopt the random-parameter multivariate NB (RPMNB) model and random-parameter multivariate Poisson (RPMP) model to analyze the impacts on different 8 Discrete Dynamics in Nature and Society crash types including rear-end, bumping-guardrails, noncasualty, and casualty crashes, with the traditional models as the base. Viewing on the traffic characteristic, more rear-end, bumping-guardrail, other, and noncasualty crashes occurred with more black spots and AADT, while other, noncasualty, and casualty crashes happened in interchange sections less likely. For speed characteristics, the crash risk of five types is greater with the higher operating speed of cars, while greater speed difference of cars results in more crashes across five types. Oppositely, the greater speed difference of trucks results in less rear-end, bumping-guardrail, other, and noncasualty crashes. As for the road geometry, more crashes occurred on the longer length of the present plane curve, while the greater length of the longitudinal slope corresponding to the minimum grade reduces the crash propensity. Considering the sight characteristics, truck drivers with greater stopping sight distance perform safer on the highways.
From the perspective of road design, to avoid the long length of a plane curve with big curvature and reduce the length of the longitudinal slope with large grade are  Figure 1: Elasticity effects across two models (RPMNB and RPMP) for 5 crash types.  With respect to the data fit, the comparison exercise including log-likelihood, AIC, and BIC value highlights the superiority of the multivariate models over the traditional models. Based on the exercise, the RPMNB model shows identical performance even with more numbers of significant parameters compared with the RPMP model (56 vs. 52).  To quantify the predictive performance measure, we calculate the elasticity effects of significant variables including interchange, black spots, AADT, V O−car , Δ V O−car , ΔV O−truck , L present , S truck , and H car for both RPMNB and RPMP models. Significant differences in the elasticity effects exist for different crash types across the two models, which are attributed to the nonlinearity of the two models.
Really, more efforts can be done on further research. With the respect to the variables, more detailed information can be added to reduce the parameter explosion, such as the proportion of trucks in each section and distinguishing the rear-end crash between different vehicle types. Otherwise, the spatial and temporal terms of the crashes are beneficial for the models to reduce the potential unobserved effects,   which can be considered in the future. Furthermore, the other methodological approaches such as latent-class, finitemixture, and Markov switching models can be used to confirm the finding of this study.