IMPROVING TRAUMA TRIAGE MODELS FOR MOTOR VEHICLE CRASHES USING EVENT DATA RECORDERS AND FUNCTIONAL DATA ANALYSIS.

: Quick identification of severe injury crashes can help Emergency Medical Services (EMS) better allocate their scarce resources to improve the survival of severely injured crash victims by providing them with a fast and timely response. Data broadcast from a vehicle’s Event Data Recorder (EDR) provide an opportunity to capture crash information and send them to EMS near real-time. A key feature of EDR data is a longitudinal measure of crash deceleration. We used functional data analysis (FDA) to ascertain key features of the deceleration trajectories (absolute integral, absolute in- tegral of its slope, and residual variance) to develop and verify a risk predic- tion model for serious (AIS 3+) injuries. We used data from the 2002-2012 EDR reports and the National Highway and National Automotive Sampling System (NASS) Crashworthiness Data System (CDS) datasets available on the National Transportation Safety Administration (NHTSA) website. We consider a variety of approaches to model deceleration data, including non- penalized and penalized splines and a variable selection method, ultimately obtaining a model with a weighted AUC of 0.93. A novel feature of our approach is the use of residual variance as a measure of predictive risk. Our model can be viewed as an important first step towards developing a real- time prediction model capable of predicting the risk of severe injury in any motor vehicle crash.


Introduction
Vehicle crashes are responsible for over 35,000 deaths per year in the United States and are among the leading cause of death for people ages 1-44 (Murphy   638 IMPROVING TRAUMA TRIAGE MODELS FOR MOTOR VEHICLE CRASHES USING EVENT DATA  RECORDERS AND FUNCTIONAL DATA ANALYSIS   et al., 2013). Among those fatally injured, most do not die immediately, resulting in the need for a fast medical response to prevent a fatal outcome (Clark et al., 2013). The "Golden Hour" is a term used to indicate the need for fast medical response to prevent fatalities in those who have suffered an acute event. The evidence for the golden hour is mixed (Newgard et al., 2010), but an analysis of crash data by Clark et al. (2013) showed a small improvement in survival with Emergency Medical Services (EMS) response within 30 minutes and an even larger improvement with early hospital intervention.
Hence, to improve the survival of victims with severe injury from motor ve-hicle crashes, the ability to identify these crashes and relay this information to EMS quickly is very important. The capability and increasing presence of Event Data Recorders (EDRs) in vehicles provides such an opportunity. The EDR is a "black box" that records deceleration and other signals in the event of a crash. (Note: EDRs technically report acceleration. However, deceleration is more natural when describing this trajectory.) This information can be trans-mitted immediately to a dispatcher and then to EMS after a crash, potentially resulting in faster deployment of EMS vehicles.
Although transmission of the information about the crash itself can reduce delays, having a prediction model that can identify crashes with a high risk of severe injury would help EMS to allocate scarce resources more efficiently. A number of researchers have developed such algorithms to predict which crashes have a higher probability of resulting in serious injury to one or more occupants. For example, Kononen et al. (2011) describe a triage algorithm based on EDR-available time independent data elements that predicts the presence of severe injury. They also included two additional predictors obtained through contact with a vehicle occupant. Augenstein et al. (2003) also developed a triage algo-rithm to predict serious injury, but their algorithm uses a number of predictors that cannot be obtained from EDRs.
Although the algorithms proposed by Kononen et al. (2011) and Augenstein et al. (2003) were designed to be used with EDR data, they were not directly developed from such data because there was not enough publicly available EDR data at the time of development. In order to compensate for this mismatch, both algorithms used an estimate of the change in velocity of the vehicle due to the crash (delta-v) as a primary predictor of injury outcome rather than accelera-tion histories measured by the EDR. This approach is problematic for several reasons. Delta-v estimates in crash databases have been shown to underestimate EDR-based delta-v by as much as 23% (Hampton and Gabler, 2010). Second, pa-rameters derived from deceleration histories (e.g. average deceleration) have been shown to have better linear correlation with the occurrence of injury compared to the database estimated delta-v (Mizuno et al., 2014). In addition, aspects of the deceleration history, for example the rapid changes in the deceleration, may provide additional information about the seriousness of the injury.
In this paper, we present a new triage algorithm, developed using EDR data rather than data from crash databases. We use functional data analysis (FDA) techniques (Ramsay et al., 2009) designed for analyzing time histories, to de-velop a prediction model of serious injury in crashes. We applied FDA to the EDR-reported deceleration and tested four FDA approaches to estimate this de-celeration. Because deceleration trajectories have different duration and time stamp intervals, we summarized each trajectory using a measure for overall de-celeration, a measure for the deceleration fluctuation, a measure of the residual variability in deceleration after removing the mean trend, and a measure of the length of time the force of the crash affected the vehicle. We then combined these statistics with the limited set of time-invariant variables available from the EDR to obtain our prediction model. We also considered the degree to which the use of deceleration statistics improves the prediction capability of our model compared to a model with only time-invariant variables.
The remainder of this article is organized as follows. Section 2 summarizes the crash history data and EDR data available in the National Automotive Sampling System Crashworthiness Data System (NASS-CDS). Section 3 describes the vari-ous forms of functional data analysis (FDA) we used to estimate the deceleration data and the analyses we conducted. In Section 4 we present the results of our analyses while Section 5 discusses the strengths and limitations of our model and future extensions.

NASS-CDS data
NASS-CDS (National Highway and Traffic Safety Administration, 2008) is an annual representative three-stage probability sample of passenger vehicle crashes sponsored by the National Highway and Transportation Safety Authority (NHTSA). To be eligible for inclusion, a crash must: (1) be police reported, (2) involve a harmful event (property damage and/or personal injury) resulting from a crash, and (3) involve at least one towed passenger car or light truck or van in trans-port on a traffic way. When a crash is selected, NASS-CDS investigators obtain police reports and conduct interviews with the occupants to collect information such as driver's age and sex, vehicle curb weight, type of vehicle, severity of in-jury measured using the Abbreviated Injury Scale (AIS) score (Association for the Advancement of Automotive Medicine, 1990), and the principal direction of impact from the crash.
Because sampling in NASS-CDS was a multi-stage design with probability of selection based on model year and the severity of injury, subjects in less serious 640 IMPROVING TRAUMA TRIAGE MODELS FOR MOTOR VEHICLE CRASHES USING EVENT DATA  RECORDERS AND FUNCTIONAL DATA ANALYSIS crashes were underrepresented in the study sample. To account for this potential bias, we used case weights equal to the inverse of the probability of selection and adjusted to known crash totals to account for the oversampling of serious crashes in NASS-CDS. These weights are provided with the data. Thus for all the vari-ance estimates and tests of association reported below, Taylor series linearization estimates were used to account for the complex sample design of the NASS-CDS (Binder, 1983;Rao and Scott, 1981).

EDR data
Electronic data records (EDRs) are devices installed in motor-vehicle to continu-ously record technical vehicle and occupant information before, during, and after a crash. Examples of key variables recorded by all EDRs include vehicle longitu-dinal acceleration/deceleration, seat belt use, vehicle speed, brake use, and airbag deployment. Newer EDRs can record acceleration data at higher frequencies, ad-ditional events, acceleration and roll angle data from multiple crash directions, and activation of collision notification systems. Vehicle identification number (VIN) is also contained in EDRs. VINs can be decoded to identify standard and equipment (both safety and otherwise) in a vehicle.
A key feature of the EDR is that, in the event of a crash, airbag deployment timing and 300 ms during a crash event is locked in memory, for up to two events. These data, along with "baseline" measures of belt use and basic vehicle information obtainable from VINs, can be broadcast in real time to emergency responders.

Crash data eligible for analysis
Our analysis relies on the 2002-2012 NASS-CDS and associated EDR datasets. We first merged the NASS-CDS datasets and crash data extracted from EDR reports. Next, we removed crashes where the driver's age and sex were both missing, resulting in 5,736 crashes. We then restricted the crashes to 3,460 frontal crashes defined as those crashes where the principal direction of force was 0 • to 40• and 320• to 350•. We removed 3,198 unusable crashes that were missing crash pulse data or that recorded a constant crash pulse throughout the crash duration. We then removed one crash that involved a vehicle type that was neither a car, SUV, pickup, nor van. Finally, we removed 12 crashes that had missing injury outcomes (including unknown AIS scores) to arrive at 249 crashes eligible for analysis. Figure 1 shows the flowchart of our selection criteria for eligible crashes.  type, curb weight (kgs), EDR driver belt status (Yes/No/Not reported), EDR front-seat passenger belt status (Yes/No/No front-seat passenger/Not reported), NASS-CDS total delta-v, multiple event indication (Yes/No), and two crash pulse variables: time (milliseconds) and longitudinal (frontal) deceleration (m/s 2 ). The sampling weight and VMAIS were obtained from NASS-CDS datasets while the seat-belt statuses and crash pulse were obtained from the EDR datasets. Presence of multiple events was derived from EDR reports. Although curb weight and vehicle type were obtained from NASS-CDS datasets, these two variables could, in principle, be derived easily using the VIN.

Variables used in analysis
We modified two of the time-invariant variables as follows. For the vehicle type, we recategorized the 64 different vehicle types to three new categories: cars (category 0-9), SUVs (category [14][15][16][17][18][19], and pickups or vans (category 20-33). We grouped pickups and vans together because the number of severe injuries in the pickup category was small. VMAIS was recorded from 0 (no injury) to 6(maximum injury). We considered crashes with VMAIS ≥ 3 as having a severe injury in the vehicle. Finally, if two events were reported in the same crash report, the multiple crash indicator was coded as "Yes".

Comparing demographics of eligible and ineligible crashes
From Figure 1, we observed that nearly 92.8% of our frontal crashes were ineligi-ble. This ineligibility is largely due to the fact that early EDRs usually provided a cruder velocity change that is not easily converted to a continuous acceleration measure; nonetheless, this high percentage of ineligible crashes could create po-tential bias in our prediction model. To verify that our prediction model would still be valid, we compared six NASS-CDS demographic variables from the eli-gible and ineligible crashes. The six variables compared were: driver's age, sex, belt use (from crash investigation), intrusion, front-seat passenger belt use (again from crash investigation), and maximum vehicle injury severity.
We observed no significant difference between eligible and ineligible crashes for the weighted mean and proportion of driver's age, sex, and belt use, front-seat passenger belt use, maximum vehicle injury severity, and sampling weights. (Note that the analysis of the sampling weights was unweighted.) The weighted proportion of driver intrusion was significantly different between eligible and inel-igible crashes, with the excluded cases having somewhat greater rates of intrusion, suggesting that our analytic dataset may contain somewhat less severe crashes than the set of frontal crashes as a whole (Table 1). We would address this issue in the discussion section.

Functional Data Analysis
Figure 2 provides a typical example of the deceleration data available from an EDR report. To utilize these crash pulses in a regression model, we employed FDA. We considered the deceleration for the i th crash to be a continuous and smooth function of t, y i (t), observed with error E ij at time t ij , where j = 1, . . . , m i for the m i observations of the i th crash, i = 1, . . . , n. Let y ij be the observed deceleration. FDA attempts to convert y i1 , y i2 , . . . , y im i to the function y i (t). This means that for a given time t ij , y ij = y i (t ij ) + E ij where y i (t ij ) is the value of y i (t) evaluated at t ij and y ij is estimated using ŷ ij = y i (t ij ). This conversion can be accomplished using several methods (Ramsay et al., 2009). For this study, we considered four approaches: two methods using unpenalized splines, the "3-millisecond" (3mil) method and "combinatorial" (combi) method; and two methods using penalized splines, the "penalized natural cubic spline" (PCS) method, and "mixed-model" (MM) method.

3 Millisecond Method
The first unpenalized method we used was based on basis splines (de Boor, 2001). Let , ( ) ≡ , ( ) are basis splines matrices of degree d, and , ( ) is obtained by the recursion relation: are the internal knots of the basis splines, and and choice of the κ il . We chose a cubic spline d = 4 and defined κ il at equally spaced intervals of k. Our objective in determining k was to choose an interval so that the least squares estimator ŷ i was smooth and yet followed the original crash pulse closely. Based on a visual inspection of the resulting curves for a variety of crash pulses, we chose to define k = 3 (i.e., knots at every at every 3 milliseconds) to retain as much potential information from the original crash pulse as possible while still smoothing crash pulse to some degree.

Combinatorial Method
As an alternative to defining κ il at fixed intervals, we could fix the number of κ il minus the two boundaries as ≡ o, and allow κ il to be freely placed anywhere along the interval [ti2, ti,m i −1]. By defining o as a small value, we induce additional smoothing while hopefully still being able to find strong data signals by allowing the placement of κ il to vary. To maintain a reasonable number of submodels to consider, we defined the potential interval [t i2 , t i5 , . . . , t i,m i −1 ], corresponding to the knots for the 3mil model, for the o κ ij s to be placed. We set o = 5 to allow as much complexity in the fitted curve as possible while still keeping the computation manageable. We considered all possible combinations of κ ij on the reduced interval [ti2, ti5, . . . , ti,m 1], and calculated ∑ ( −̂) 2 =1 combination. We then chose the for the ℎ crash

Penalized Cubic Spline Method
To reduce the risk of overfitting, we considered a penalized spline approach that downweighted values of ĉ il that lead to highly variable values of ŷ ij by introducing a penalty in the objective function to be minimized and [ ( )] is the ℎ derivative of function . The first term in (3) is the objective function that minimizes the least squares estimation of c i in (1); the second term induces a penalty for large amounts of "acceleration" or "deceleration" when estimating the function y i (t). Least-squares estimation can be adapted to accommodate (3) by rewriting as  (Golub et al., 1979), an extension of standard cross-validation procedure that predicts the − ℎ observation after leaving out the − ℎ observation in the estimation of .
(Because ( ) = (̃), we do not need to convert ( )back to the ̃ scale.) We obtained by estimating the Cholesky decomposition of Ω using Smith's(1995) method. We then estimated as where, using standard mixed-model results (Harville, 1976), we have , where 2 ̂2 are estimated using restricted maximum likelihood (REML).

Slope of ( ) at
In addition to obtaining ( ), we were interested in calculating the slope of ( ) at . We hypothesized that vigorous fluctuations in the crash pulse may be an important predictor for severe injury outcomes. For the basis spline methods in 3.1.1-3.1.3, we computed the derivative of the least squares estimator as Φ̂ where For the slope of the MM method in 3.1.4, the first derivative of the equation 6 yields Where ̅ was the( , ) element of . To obtain ̅ , we took the first derivative of the Ω matrix with respect to . Then using the fact that

Extracting summary measures from the crash pulse
Typically in FDA, once y i (t ij ) has been estimated, researchers would use the coefficients c il as the variables in the second-stage regression model. However, in this study we faced two difficulties. First, the duration of each crash was different. Second, the time stamp interval was different for each crash. These two difficulties implied that defining a single set of linear combination of functions for all crashes and subsequently using their coefficients would be inappropriate since each crash should have its own linear combination of basis functions ( ) to define y i (t ij ). To overcome these difficulties, we extracted four summary measures from the y i (t ij ) of the four methods separately.
The first summary we extracted was ̂ defined as ̂= ∫ | ( )| absolute integral of yi(ti) over the duration ti. We took the absolutes because we did not want y i (t i ) to cancel each other during integration. Because closed form solutions were not available, we used numerical methods (Simpson's integral) to approximate ̂. ̂ sould approximate the "instan taneous deceleration" measure delta-v, a good predictor of severe injury risk. Indeed ̂ was correlated with the delta-v given in NASS-CDS (Kendall's correlation of 0.44, p-value<0.001) Thus ̂ should be an important predictor of injury risk. Figure 3 provides an illustration of the differences between crashes with a high ̂ and a low ̂.
The second summary measure that we extracted was ̂= ∫ | ′ ( )| the "absolute integral" of all the estimated slope of ( ) over the duration . This definition of ̂ was designed to capture the extent to which the crash pulse fluc-tuated. The absolute signs prevented positive slopes from canceling out negative slopes, ensuring that we did not underestimate the fluctuations in the crash pulse. Again, we used numerical methods (Simpson's integral) to approximate ̂. Fig-ure 4 illustrates the difference between a crash with higher ̂ vs. a crash with a lower ̂.
The third summary measure we extracted was 2 the variance of the residual. Residual variance provided an additional measure of the degree of variability in RECORDERS AND FUNCTIONAL DATA ANALYSIS the crash pulse not exp1ained by ̂ and ̂. Figure 5 illustrates this clearly. We estimated ̂2; for 3mil , combi , and PCS method by taking the samp1e variance of the residuals , whi1e for the MM method , we obtained the REML estimate of ̂2.
Lastly , we extracted a summary measure called time to 0.25Gs (tt025). Tt025 was defined as the time the crash pulse took to eventually return to ± 0.25 G s . This gives a measure of how long the force from a crash affects the vehicle. We computed tt025 by starting from the last crash pulse recorded for a given crash and moving backwards to the first deceleration value recorded , searching for the deceleration value that first exceeded ±0.25G s , 0.25G s .
We obtained tt025 by taking the difference of 0.25G s and the time at the start of the crash pulse. Figure 6 shows graphs of a crash pulse with a high tt025 vs. low tt025. A crash pulse with high tt025 tended to have their crash activity "stretched out" while crash pulses with low tt025 had their crash activity "compressed" . Hence , we expect a crash pulse that extends over a longer time (i.e., high tt025) should result in less severe injury than one that occurs in a shorter time frame.

Predicting risk of severe injury in motor vehicle crashes
We computed our prediction model by first merging each set of summary measures with the outcome, VMAIS and the other time-invariant variables (see subsection 2.4) to form four datasets. Then, we expanded these four sets to a further sixteen by taking different combinations of the following: 1. logarithm transform or no logarithm transform on the summary measures and 2. or no scaling of ̂. This was done to investigate whether logarithm transformation of the summary measures or whether the scaling of ̂ improves prediction performance. As presentation of coefficients and results from sixteen different models would be cumbersome, we used the highest leave-one-out cross-validation (CV) weighted area under the curve (AUC) of the receiver operating characteristic (ROC) curve to determine which model to present. The weighted AUC was computed based on the results of a weighted logistic regression with VMAIS as the outcome and summary measures and time-invariant (baseline) variables as the covariates.
To investigate the performance of our chosen model, henceforth called the `deceleration' model, we first compared the weighted AUC of the `deceleration' model with a `baseline' model. The `baseline' model was one which contained only time-invariant EDR variables. This was to determine if adding deceleration information to our model improved prediction performance. We also measured the variability of this difference by running 1,000 bootstrap samples of the eligible crashes and calculating the weighted AUC for each of them. We then took the difference of their weighted AUC, sorted this calculated difference, and took the 25th and 975th values as the lower and upper bounds for the 95% bootstrap confidence intervals of the weighted AUC difference.
To complement the weighted AUC results, we plotted the weighted ROC and precision-recall (PR) curves (Davis and Goadrich, 2006) and constructed a false discovery rate (FDR) table for both deceleration and baseline models. x-axis). More speciflcally, precision (yaxis) is defined as + , while recall (x-axis) is de_ned as + . PR curves share many properties with the ROC, but provide a focus on the trade-o_ between the fraction of AIS 3+ injuries captured at a given prediction level cutpoint versus the fraction of positive cases that will actually have an AIS 3+ injury at that level. These are particularly important measures in this setting, where the cost of missing severe injuries must be balanced against resources assigned unnecessarily to non-severe injuries. PR curves lying in the top-right imply good performance. We plotted weighted PR curves by first using the sampling weight of each crash as the number of either severe injury outcome or non-severe injury outcome a crash represents. Then, we smoothed these points using the locally weighted scatter-plot smoothing (loess) to curves the PR curves.
The FDR table we computed provides similar information as the PR curves.The difference was, instead of looking at the proportion of true-positives out of all positive results, our FDR table gave the proportion of false-positives out of all positive results. We presented this proportion in the right column of our FDR table with the predicted probability cutoff points in the left column, and the proportion of severe injuries categorized as positives in the middle column. For this study, we used predicted probability cutoff points of 1-10%, 15%, 20%, 25%, Finally, we ran a leave-one-out CV on the deceleration and baseline model to investigate how well they will perform given new data. Similarly, we looked at their weighted ROC and PR curves, and the FDR tables using the same predicted probability cutoff points mentioned in the previous paragraph. We completed all analyses using Rx64 3.0.3.

Results
To get a visual sense of how the four different FDA methods estimated the original crash pulse, we present four selected crashes with their estimated trajectories, ̂ in Figure 7. We observed that the 3mil and PCS methods tend to follow the original crash pulse closely while the combi method smooths almost all crash pulses. The MM method alternates between following the crash pulse when there are strong signals in the data (see the top two plots in Figure 7), or smoothing the pulses when they are \noisier" (see the bottom two plots in Figure 7).
From the results of the leave-one-out CV, we found that the model with summary measures estimated using penalized method, no logarithm transform on the summary measures, and no scaling produced the highest weight leave-one-out CV AUC, 0.75. Using this model, we next focus on the descriptive statistics of our variables ( Table 2). We found that the (weighted) mean of ̂, ̂ , ̂2, and curb weight in severe injury crashes were higher compared to non-severe injury crashes (847 vs. 499, 87 vs. 69, 59 vs. 12, and 1,694 vs. 1,681 respectively).
However, the mean of tt025 for severe injury crashes was slightly lower compared to non-severe injury crashes (110 vs. 112). The proportion of cars and pickups or vans were higher in severe injury crashes compared to non-severe injury crashes (45% vs. 42% and 26% vs. 24% respectively) while SUVs had a lower proportion involved in severe injury crashes (29% vs. 34%). More than half (59%) of the eligible crashes did not have seat-belt statuses of the driver and front-passenger reported in their EDR reports. We observed a lower proportion of the known-belted occupants involved in severe injury crashes (12% vs. 38% and 6% vs. 7% for driver and front-seat passenger belt status respectively). The proportion of unbelted drivers involved in severe injury crashes was higher than in non-severe injury crashes (21% vs. 3%), while we found a lower rate on unbelted front-seat passengers in severe injury crashes vs. non-severe injury crashes (4% vs. 17%). Finally, we observed a higher proportion multiple crash recorded of severe injury crashes (18% vs. 6%).  i Although estimated delta-v from databases like NASS-CDS underestimate actual delta-v, we still wanted to check whether our summary measures had good linear correlation to NASS-CDS delta-v because many studies have shown that NASS-CDS delta-v is a good predictor of severe injury risk. We found that all summary measures had a weak to moderately positive correlation with NASS-CDS delta-v with Ĝ having the highest linear correlation 0.44 while tt025 had the weakest, 0.05. Ĝ was moderately but significantly correlated with ĝ, σ 2 , and tt025 respectively (0.31, p-value< 0.001; 0.34, p-value< 0.001; 0.18, p-value< 0.001).

Prediction of Injury Risk
The coefficients of our deceleration model are given in Table 4. Ĝ , ̂2 , tt0.25 belted driver, and multiple crash indicator were significant, with a 10 m/s2 in-crease in Ĝ associated with a 4% increase in odds of AIS 3+ injury (95% CI 2%,6%), a 1-g increase in the error of the estimation method associated with a 2% increase in the odds of injury (95% CI 1%, 3%), a 1 millisecond increase in the duration of the crash associated with a 3% decrease in odds of AIS 3+ injury (95% CI 1%, 5%), presence of a belted driver associated with a 97% decrease in odds of AIS 3+ injury (95% CI 82%, 99%), and presence of multiple crashes associated with a 42fold increase in the odds of AIS 3+ injury (95% CI 2-fold,780-fold). Surprisingly, the coefficient for belted front-seat passenger was positive (1.49) although not significant (p-value=0.43). The weighted AUC for our decel-eration model was 0.93. We also computed the weighted AUC of the deceleration model with the summary measures centered and the square of these centered summary measures added. The weighted AUC was similarly 0.93 suggesting no improvement when we include squared terms of our summary measures. For the baseline model, only the coefficients of the belted driver and multiple crashes were significant (-3.23 and 2.04 respectively with p-values of < 0.001 and 0.04). The belted front-seat passenger coefficient was still positive, not significant (p-value=0.46), and the magnitude of the estimate was smaller than the deceler-ation model. The weighted AUC of our baseline model was 0.78, 0.15 less than the deceleration model. The 95% bootstrap confidence interval of the weighted AUC difference between the deceleration model and the baseline model was (0.04,

0.23).
We next checked each observed crash for influential points using the adjusted Cook's D for survey data. We found two influential observations where their covariate pattern suggested a non-severe injury crash but were in fact severe crashes. We removed them as part of a sensitivity analysis and found the resulting coefficients, 95% confidence intervals, and AUC similar to those presented in Table 4. Hence, we choose to present the model without the removal of the two influential points.
The ROC and PR curves suggested that our deceleration model performed better than the baseline model ( Figure 8). We observed that the weighted ROC for the baseline model was to the right of the ROC of the deceleration model for false positive rates between 0 and 0.8 implying better performance for the deceleration model. Similarly, we found that the loess PR curve of the baseline model was on the left of the deceleration model for all recall values indicating lower prediction performance of the baseline model. Similar results were found in the FDR table (Table 5). The deceleration model had a large decrease in capture rate at a predicted probability cutoff of 6%, at that cutoff, we would capture 7 in 10 of the AIS 3+ injuries, with about one in four crashes actually having an AIS 3+ injury. To capture a similar fraction of AIS 3+ crashes in the baseline model, we would have to use a cutpoint of 2%, at which point fewer than 1 in 20 crashes would have an AIS 3+ injury. Figure 9 provides the "leave-one-out" cross-validation results for the AUC and PR plots, and Table 6 provides the equivalent for the FDR table. These results suggest less gain for the deceleration model relative to the baseline model, although a cutoff of 5% risk still captures one-third of the AIS 3+ injuries with one in ten being AIS 3+; capturing an equivalent fraction of AIS 3+ injuries using baseline data would yield a "hit rate" of less than one AIS 3+ crash in 50.

Discussion
In this study, using only information available in EDR reports, we successfully developed a new severe injury risk prediction model able to estimate the proba-bility of a severe injury in motor vehicle crashes near real-time. Although such real-time prediction models are not new in trauma literature, our sole use of data from EDR guaranteed that our model reflects the actual situation when EDR data is received by EMS from a dispatcher. Another strength of our model lay in the use of a novel variablethe crash pulseas the main variable in our regression model. To use this crash pulse, we employed a two-stage approach, with 656 IMPROVING TRAUMA TRIAGE MODELS FOR MOTOR VEHICLE CRASHES USING EVENT DATA RECORDERS AND FUNCTIONAL DATA ANALYSIS summary measures from the first-stage FDA model, as inputs for the second-stage weighted logistic regression model. We were not aware of any such applications of FDA to crash pulses for predicting severe injury risk in trauma literature. By extracting summary measures from the estimated deceleration, we were able to avoid the problem of trying to find a common set of linear combination of functions for trajectories which have different duration and time stamp interval. With the use of FDA to extract information from trajectories gaining popularity, we believe more researchers will face similar situations where application of usual FDA procedures would be inappropriate.
Our deceleration model performed well when predicting the population of eligible crashes compared to the baseline model. In addition, we found that our deceleration model was able to give predicted probability cut-offs where the capture rate of true severe injury crashes was high, while incurring a reasonable amount of false discoveries. This could be seen from Figures  8 and 9 and Tables 5 and 6. We found less optimistic results using a leave-one-out cross-validation procedure. While suggesting a similar optimal cut-off point of 6-10% predicted risk for classifying crashes, our cross-validation analysis suggested that such a cutpoint would capture slightly less than 3 in 10 of AIS 3+ crashes, but yield a false discovery rate of approximately 9 in 10.
Not surprisingly, we found that Ĝ , driver belt use, and presence of multiple crashes increased the risk of a severe injury crash significantly. However, our finding that σ 2 and tt025 were significantly associated with severe injury in motor vehicle crashes was not entirely expected and could contribute to the further understanding of injury in motor vehicle crashes. A significant and positive σ 2 value implies that differences between the estimated deceleration and the original crash pulse are associated with higher risk of severe injury. This means that rapid changes in deceleration unexplained by the estimated pulse signal causes additional injury above and beyond what we would expect if there were no rapid changes.
The second significant coefficient that we felt could further the understand-ing of severe injury crashes was tt025. A negative and significantly coefficient implies that as the duration of the force of the crash decreases, the risk of severe injury increases. This may seem counterintuitive but with careful thought will be clear. Observe Figure 6. With similar force (Ĝ of 448 [left] and 478 [right]), the deceleration duration on the left took 149 milliseconds while the deceleration on the right took 65 milliseconds. As a result, the crash pulse on the right was 'squashed' and it had to 'deliver' the similar amount of force but within a shorter duration. This most likely leads to increased risk in having a severe injury. With the understanding that longer effective duration of the force of the crash would reduce the risk of severe injury in a crash, motor vehicles which prolong the effective duration of the crash could be developed.
Our success in developing a fairly good near real-time prediction model was limited by three issues: a high proportion of unusable crashes, non-reporting of seat-belt statuses, and inflated sampling weights. In our study, we found that 3,198 out of 3,460 (about 92.4%) crashes were unusable. This was because different manufacturers present different variables in their EDR report and most reports contained velocity change instead of deceleration because earlier EDRs usually provide cruder velocity change that is not easily converted to a continuous acceleration measure. As newer vehicles equipped with newer EDRs enter the car fleet and EDR penetration increasing, we expect this ineligibility issue to subside in the future.
While we found only one statistically significant difference between the usable and unusable crashes (a tendency for usable crashes to be somewhat less severe), we did not present a post-stratified model (a model that uses weights such that the rates of intrusion between eligible and ineligible crashes would be similar) because first, although post-stratification balanced the rates of intrusion between eligible and ineligible crashes, the resulting coefficients, 95% confidence intervals, and AUC were similar to the results reported in Table 4. Second, we noted that the intrusion rate of 46 cm ≤ x < 61 cm in ineligible crashes was 3%, an exceedingly high proportion for such severe injury crashes in the population. We suspect this is an artificial inflation of the rates due to highly variable survey weights being allocated to these ineligible crashes, an issue we will return to in subsequent paragraphs.
The second issue we faced was that many EDRs did not report driver and front-seat passenger belt status. From Table 2, we can see that nearly 59% (weighted) of our eligible crashes did not have their seat-belt status reported. Unlike deceleration, we were unable to observe any pattern or determine the rea-son why such a large proportion of EDR reports did not report seat-belt statuses. The consequence of this was likely the biasing of seat-belt status coefficients to-wards the null. The assumption we made here was that EDR reports were equally likely to not record seat-belt statuses of belted or not belted occupants. This im-plied that the belted effect for drivers could have been more protective but for a belted front-seat passenger, the contradictory effect of a belted front-seat pas-senger being a risk for severe injury during a crash would have been increased. A second consequence of EDR reports not recording seat-belt status was likely the reduction in prediction performance of cross-validation.
Finally, the third issue we faced was highly variable sampling weights. Since the NASS-CDS sampling weights reflected how likely a crash was to be selected, we would expect crashes resulting in a severe injury or driver intrusion > 3cm to have smaller sampling weights since severe crashes were oversampled. How-ever, we found that four of our eligible severe injury crashes (about 15%) had large weights of more than 100. Similarly, we found large sampling weights for driver intrusions of > 3cm for ineligible crashes. These inflated sampling weights affected the precision of our model's coefficients by increasing their variability.
The possible remedies for these three issues we faced were limited. We could have used some imputation method to recover the crash pulses in unusable crashes or perform sensitivity analyses to determine the effects of non-reporting seat-belt statuses and inflated sampling weights on our model. However, imputation of crash pulses in unusable crashes did not seem appropriate because the proportion of crashes that did not report crash pulse was high, about 92%, and crash pulses were usually very variable, making imputation difficult and likely inaccurate. Sensitivity analysis of seat-belt status may have given us a better sense of the variability in our model but we felt that sensitivity analysis would not be helpful because we did not have a good idea of the mechanism that generated the non-reporting of seat-belt statuses in EDR reports.
From the issues presented, we can see prediction of severe injury crashes in motor vehicles is still a challenging field. The model we have developed here should not be considered as the final model to be used for all motor vehicle crashes, but we believe it can be considered a promising seminal model for real-time severe injury crash prediction. To develop a comprehensive model we suggest the following steps: 1. Develop a separate prediction model that uses information from the ad-justed velocity change instead of deceleration. This step will be essential if some EDRs do not report deceleration.
2. Add the lateral, vertical, and rollover of crash pulse and adjusted velocity components into the model described in 1.
3. Develop a joint model to compute the FDA estimates and logistic regression model in a single step. Two-stage regression models are known to suffer from bias and loss of efficiency, due to the error in the estimation procedure not be accounting for in the second stage of estimation. However, joint model estimation is currently a time-demanding procedure, one that is un-likely to be carried out in practice given the time constraints involved in responding to a crash event. Methods that can utilize aspects of a joint estimation procedure outside of a full model-fitting would have great potential value in this setting and are the target of future research.
In conclusion, we believe that improve and reliably report crash pules and seat-belt statuses, the goal of developing a good and accurate near real-time prediction model of severe injury in a motor vehicle crash will be "just around the corner."