- Email: [email protected]
Accounting for sampling variability, injury under-reporting, and sensor error in concussion injury risk curves
Journal of Biomechanics 48 (2015) 3059–3065
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Accounting for sampling variability, injury under-reporting, and sensor error in concussion injury risk curves Michael R. Elliott a,b, Susan S. Margulies c, Matthew R. Maltese d, Kristy B. Arbogast e,n a
Department of Biostatistics, School of Public Health, University of Michigan, Ann Arbor, MI 48109, United States Survey Methodology Program, Institute for Social Research, University of Michigan, Ann Arbor, MI 48109, United States c Department of Bioengineering, University of Pennsylvania, Philadelphia, PA 19104, United States d Department of Anesthesiology and Critical Care Medicine, Children's Hospital of Philadelphia, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104, United States e Center for Injury Research and Prevention, Children’s Hospital of Philadelphia, Department of Pediatrics, University of Pennsylvania, 34th and Civic, Center Blvd, Suite 1150, Philadelphia, PA 19104, United States b
art ic l e i nf o
a b s t r a c t
Article history: Accepted 23 July 2015
There has been recent dramatic increase in the use of sensors affixed to the heads or helmets of athletes to measure the biomechanics of head impacts that lead to concussion. The relationship between injury and linear or rotational head acceleration measured by such sensors can be quantified with an injury risk curve. The utility of the injury risk curve relies on the accuracy of both the clinical diagnosis and the biomechanical measure. The focus of our analysis was to demonstrate the influence of three sources of error on the shape and interpretation of concussion injury risk curves: sampling variability associated with a rare event, concussion under-reporting, and sensor measurement error. We utilized Bayesian statistical methods to generate synthetic data from previously published concussion injury risk curves developed using data from helmet-based sensors on collegiate football players and assessed the effect of the three sources of error on the risk relationship. Accounting for sampling variability adds uncertainty or width to the injury risk curve. Assuming a variety of rates of unreported concussions in the non-concussed group, we found that accounting for under-reporting lowers the rotational acceleration required for a given concussion risk. Lastly, after accounting for sensor error, we find strengthened relationships between rotational acceleration and injury risk, further lowering the magnitude of rotational acceleration needed for a given risk of concussion. As more accurate sensors are designed and more sensitive and specific clinical diagnostic tools are introduced, our analysis provides guidance for the future development of comprehensive concussion risk curves. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Sensor error Injury risk Traumatic brain injury
1. Introduction Helmet-based head impact sensors have been widely utilized to understand the link between the biomechanics of head impact and clinical outcomes of concussion in athletes (e.g. Rowson et al., 2009; Brainard et al., 2012; Crisco et al., 2010; Mihalik et al., 2007; Wilcox et al., 2015). Such sensors have accelerometers and in some cases, gyroscopes embedded in the padding or attached to the helmet shell to estimate the magnitude of linear and rotational acceleration experienced by the athlete during head impact. Similar systems have been developed for non-helmeted sports and are integrated into headbands, skull-caps or mouthguards or directly attached to the athlete near the mastoid process (King n
Corresponding author. Tel.: þ 1 215 590 6075; fax: þ1 215 590 5425. E-mail address: [email protected] (K.B. Arbogast).
http://dx.doi.org/10.1016/j.jbiomech.2015.07.026 0021-9290/& 2015 Elsevier Ltd. All rights reserved.
et al., 2015; Hernandez et al., 2014; Bartsch et al., 2014). Quantifying the relationship between biomechanical input and clinical outcome is critical to the advancement of concussion prevention principles including the assessment of injury risk, the design of protective equipment such as helmets, and the development of training and policies intended to limit exposure to head impacts and injury risk. The most common approach to quantifying the link between biomechanical input and concussion is through injury risk curves. Injury risk curves describe the injury probability given a specific mechanical input – for example, concussion risk given a particular head acceleration. This relationship is not linear but rather sigmoid, often in the form of a Weibull distribution or some other cumulative distribution function (CDF). Previous research has constructed injury risk curves for concussion. Pellman et al. (2003) used head impact data from the NFL and Rowson and Duma (2011) used collegiate football head impact data to describe the relationship
3060
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
between linear acceleration and concussion risk. Recognizing that rotational acceleration or velocity is an important part of the causal pathway for concussion, Rowson et al. (2012) expanded this work and quantified a relationship between rotational acceleration and concussion risk. Specifically, they used data from 300,977 head impacts among US football players (57 of which resulted in diagnosed concussions) to estimate a logistic risk function for concussion as a function of rotational acceleration. The same authors proposed a combined injury risk function that incorporated both linear and rotational components (Rowson and Duma, 2013). The utility of injury risk curves relies on accuracy of the metrics that lie on both the abscissa and ordinate axes: accuracy of clinical diagnosis (dichotomous variable on the ordinate) and accuracy of biomechanical measure (continuous variable on the abscissa). There are three sources of uncertainty or error that influence concussion injury risk curve development: sampling variability of the clinical outcome, concussion under-reporting, and head impact sensor error. The effect of these errors on injury risk curves has not been rigorously explored in the literature and therefore detailed technical treatment of these three issues is the focus of this analysis. For exemplar purposes, we utilize the Rowson et al. (2012) data as the basis for our analyses however these approaches can be applied to other injury risk curves. Each source of error is discussed below. First, injury risk curves do not typically account for uncertainty in the finite sample of head impacts. This might seem minor, given that typically large numbers of impacts are collected over a given period of games and practices. For example, 300,977 observations were made in the Rowson et al. (2012) analysis. However, concussion remains a rare event (e.g. Rowson et al. report 57 diagnosed concussions), thus there is substantial sampling variability that influences the relationship between concussion risk and rotational acceleration. Standard logistic regression methods (Hosmer et al., 2013) will account for this sampling variability; which is non-trivial when the outcome is rare. A second source of uncertainty is misclassification of the concussion outcome itself, with under-reporting estimated to be as high as 53% (McCrea et al., 2004). The exact under-reporting proportion is not known so parametric analysis using a rational range of estimates is important. A third issue is that the measure of acceleration obtained from helmet-based head impact systems is subject to non-trivial error levels between what is measured by the sensor and the true head kinematics (Allison et al., 2014, 2015; Jadischke et al., 2013; Funk et al., 2012). If this error is non-differential (that is, not correlated with the true acceleration), as suggested by Allison et al. (2014), then accounting for it in the creation of injury risk curves will strengthen the estimated relationship between rotational acceleration and concussion risk, and lead to more accurate concussion thresholds. Therefore, the objective of this work is to demonstrate, using accepted analytical methods, the magnitude of the effect of these three sources of error on the shape and interpretation of concussion injury risk curves. This analysis provides a roadmap for future research efforts as new sensors are designed and clinical diagnostic tools become more accurate. We estimate concussion risk curves as a function of rotational acceleration and identify acceleration thresholds associated with a specific injury risk that can be derived from these relationships. 2. Methods We use summative published data as the basis for our assumptions about true concussion rates and information about rotational acceleration measurement error (e.g. Rowson et al., 2012; Allison et al., 2014) and generate synthetic datasets to provide a platform for statistical analysis to assess the sources of uncertainty related to sampling, under-reporting, and sensor error. The synthetic dataset was developed using the Bayesian statistical principles, which allow generation of posterior predictive
distributions that can be used to obtain the injury risk curve measures of interest. The synthetic origin of the dataset is of minor importance; rather, the major objective of this communication was to demonstrate the importance of accounting for uncertainty in conclusions drawn from the dataset.
2.1. Dataset generation Using statistics about the distributions among the concussed and non-concussed individuals in the Rowson et al. data (Rowson et al., 2012, Table 1), we can use Bayesian methods to estimate the relationship between concussion risk and rotational acceleration described by their data. We use their reported distributions to generate posterior predictive distributions (synthetic data) of the rotational acceleration by observed concussion status for analysis. Specifically, we construct a posterior predictive distribution of the rotational acceleration of non-concussive injuries from the HIT system and 6DOF measurement device,1 and of the concussive injuries, also measured through the HIT system. A posterior predictive distribution is a Bayesian concept, and can be thought of as a draw of a future set of observations, yrep , conditional on observed data y, and assuming a model for the data f (y|θ ) and a prior distribution p(θ ) for the model parameters (Gelman et al., 2014). The mathematical formalism is as follows:
p(y rep |y ) =
∫ f (yrep |θ )p(θ |y)dθ
(1)
where p(θ|y ) = f (y|θ )p(θ )/ ∫ f (y|θ )p(θ )dθ is the posterior distribution of θ . The integral in (1) is usually not computable analytically, but p(yrep |y ) can be obtained by simulation by drawing θ rep~p(θ|y ), and then yrep ~f (y|θ rep). This approach provides a principled way to accommodate uncertainty in the data distributions of head rotational acceleration measures in football players. Although we do not actually have the data from Rowson et al., we can still approximate the posterior distribution of θ as N (θ^, J −1(θ^ )) , where θ^ is the maximum likelihood estimator of θ , and J (θ ) is the Fisher information matrix associated with the distribution of y , computed as J (θ ) = E (−
∂2 ∂θ ∂θ T
log f (y|θ )) for the assumed
distribution f (y|θ ) . 2.2. Assessing sampling variability Because concussion is a rare event when considering all sensor-detected head movements on the sports field, conclusions on the relationship between injury and head acceleration may be subject to substantial uncertainty due to the variability in a finite dataset attributable to the statistical likelihood associated with an underlying rate of concussion. We can assess this uncertainty by creating many examples of our synthetic population and observing x concussions in one example and y (slightly more or less than x) concussions in another. The non-concussed data is assumed to follow a Weibull distribution; the concussed data is assumed to follow a Rician distribution (Rice, 1945). The large sample size of the non-concussed data means that the normal approximation will be highly accurate; for the concussed data, the smaller sample size is at least partially offset by the fact that the Rician distribution will be approximately normal for the parameter values reported. To compute an estimated concussion risk as a function of rotational acceleration, we follow Rowson et al. and assume a logistic regression model, so that
⎛ e β0 + β1y ⎞ ⎟⎟ P (Z = 1|Y = y ) = ⎜⎜ ⎝ 1 + e β0 + β1y ⎠
(2)
where Z is an indicator for concussion, Y is the rotational acceleration, and Z is set deterministically to 1 if Y is drawn from the Rician distribution and 0 if drawn from the Weibull distribution. We can obtain inference about the parameters β that govern the risk function (2) by extending the posterior predictive distribution in (1):
p(β|y ) =
∫ p(β|yrep )p(yrep |θ )p(θ |y)dθdyrep .
(3)
We can obtain draws from the posterior distribution in (3) under the assumption of a flat prior by fitting a logistic regression model to a draw of yrep rep rep rep rep from (1), and then drawing β ~N (β^ , I −1(β^ )), where β^ and I −1(β^ ) are
obtained from logistic regression using yrep . In this analysis we obtain 1000 separate samples or “draws” of parameters of interest to make inference (compute posterior means, credible intervals [the Bayesian analog of confidence intervals]). Assuming that the distributional assumptions are correct, we are able to obtain estimates of the concussion/rotational acceleration risk that account for the influence of sampling variability without any loss of information. We use these estimates to develop concussion risk curves as a function of rotational acceleration that
1 Two different variants of a similar head impact sensor were used in the Rowson study. Most data were collected by the commercially available HIT system (Simbex, Lebanon, NH) however some data were collected using a version modified for research use that implements additional sensors towards improving accuracy.
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065 incorporate uncertainty due to sampling variability, and compare them to the risk curves that do not account for this uncertainty.
The Bayesian paradigm is also particularly convenient for handling random under-reporting of concussion in our subjects in the non-concussed group, because it is the basis behind accepted missing data methodological approaches (Little and Rubin, 2002). We thus extend our work accounting for uncertainty in the association between rotational acceleration and concussion risk to account for underreporting of concussions by identifying missing concussions in the non-concussed group and “converting” them to the concussed group. We assume that missing concussions are missing at random (Little and Rubin, 2002) conditional on the rotational acceleration distribution. We require an assumption about the fraction of concussions that are missed. We consider three fractions: (1) a fraction based on the overall concussion rate assumed in Rowson et al. (2012) – 0.0726%; (2) a fraction based on one-half of that assumed in Rowson et al. (2012), and (3) a fraction twice that assumed in Rowson et al. (2012). If the true concussion rate in the sample is 0.0726%, then there are a total of 219 concussed players, and 162 ( ¼ 219–57) of the 300,920 “non-concussed” players are actually concussed. We assume that there are no misclassified observed concussions, that is, all rotational acceleration observations Y generated from the Rician distribution are associated with concussions (Z ¼1). Next, we assume that the unobserved concussions are drawn from the same distribution of Y as the observed concussions. The posterior distribution of the logistic parameters governing the relationship between rotational acceleration and concussion risk is given by
∫ p(β|yrep , z rep)p(yrep |θ )p(θ |y)dθdyrep z rep
(4)
Draws from p(β|yrep , z rep) are obtained using a Monte Carlo Markov Chain (MCMC) (Gelman et al., 2014). The full algorithm proceeds as follows: 1. Draw θ rep~p(θ|y ) from N (θ^, J −1(θ^ )) , where 286,579 draws are from the Weibull distribution estimated using non-concussive rotational accelerations obtained from the HIT system, 14,341 draws are from the Weibull distribution estimated using non-concussive rotational accelerations obtained from the 6DOF system, and 57 draws are from the Rician distribution estimated using concussive rotational accelerations. 2. Draw yrep ~f (y|θ rep) under the appropriate distributions. 3. Obtain a draw from p(β|y
rep
,z
rep
) using MCMC by iterating as follows:
a. Impute the sth z rep, s|yrep , β (s − 1) , where β (s − 1) is the draw of β from the previous iteration of the MCMC ( β (0) is set to be previous saved draw in the overall algorithm, i.e., the final drawn from the previous MCMC). z rep, s is set to 1 if
yrep P (Z
is
rep
= 1|Y
drawn rep
)=
from
the
Rician
exp(β (s − 1) + β (s − 1)Y rep) 0 1 , 1 + exp(β (s − 1) + β (s − 1)Y rep) 0 1
distribution;
otherwise
subject to the constraints that
∑i zirep, s = 219. b. Draw the sth β (s)|yrep , z rep, s by first fitting a logistic regression model to draw rep rep rep of z rep, s and yrep , and then drawing β ~N (β^ , I −1(β^ )), where β^ and
I −1(β^
rep
) are obtained from this logistic regression.
c. Repeat (a) and (b) for S times, until the autocorrelation between the draws of early values of β (s) and later values of β (s) is trivial. We chose S = 100 in the analysis below. To account for the uncertainty in the actual concussion rate, we repeat both of the above analyses assuming a true concussion rate twice that (0.1452%) and onehalf that (0.0363%) assumed in Rowson et al. (2012). This yields an estimate of an additional 380 and 52 missed concussions in the “non-concussed” group, respectively, above the 57 reported concussions. 2.4. Accounting for sensor error We use the SIMEX method to account for the measurement error (Cook and Stefanski, 1994). SIMEX, short for Simulation Extraction, begins from the observation that the ordinary least squares estimator of a regression slope Z = β0 + βxX + ε estimates
βx σ 2 σ 2 + σu2
do here, estimate an approximate quadratic regression model β^(λ ) = γ1 + γ2λ + γ3λ2, computing β^x = γ^1 − γ^2 + γ3 (Carrol et al., 1995). Variance estimates of β^x are computed via a bootstrap (Stefanski and Cook, 1995). We assume a multiplicative error form for the rotational acceleration data – specifically, the observed rotational acceleration Y equals XU , where U follows a log-normal distribution with expec-
2.3. Accounting for under-reporting
p(β|y ) =
3061
when X is replaced with Y = X + U , where σu2 is the variance of the
(mean zero) measurement error U and σ 2 is the variance of the true residual error ε . By adding error with variance λσ 2 to W, one can compute a series of values β^(λ ) u
from different values of λ as a function of
βx σ 2 σ 2 + (1 + λ )σu2
and extrapolate to β^( − 1) ,
which estimates βx . One can work directly with this non-linear equation, or, as we
2
tation 1 and variance e σ − 1. Thus log(U ) is generated from a normal distribution with mean −σ /2 and variance σ 2 , and Y is replaced with exp(log(Y ) + λ log(U )) to 2
implement SIMEX. We estimate σ 2 as 0.2706 by method of moments using the alldirection impact data from Allison et al. (2014), where instrumented anthropomorphic test devices (ATDs) fitted with HIT sensor-outfitted helmets were impacted, and helmet rotational accelerations were compared to values measured within the ATD head. We embed this helmet sensor measurement error in our multiple imputation method by replacing the draw the β (s)|yrep , z rep, s at step (3b) (s ) (s ) (s ) (s ) with a draw of βx(s)|yrep , z rep, s from N (β^x , cov(β^x )) , where β^x and cov(β^x ) are
obtained from the SIMEX algorithm utilizing the draws yrep , z rep, s .
3. Results The consolidated results of all simulations with and without accounting for the uncertainties due to sample variability, underreporting and sensor error are provided in Tables 1 and 2, with the mean logistic regression parameters in Table 1 and injury risk curve data in Table 2. Importantly, Table 2 presents the median values for rotational acceleration thresholds (and confidence intervals) for a range of injury risk levels, from 10 to 90%. To orient the reader, column A in Tables 1 and 2 is the synthetic dataset (based on Rowson et al. (2012) published observations of 57 concussions of 300,977 impacts), without accounting for any sample variability, under-reporting, or sensor measurement error. 3.1. Sampling variability The empirical CDFs of 20 draws from the posterior predictive distributions of the 286,579 non-concussive rotational accelerations obtained from the HIT system (Fig. 1) and the 14,341 non-concussive rotational accelerations obtained from the 6DOF system (Fig. 2) show that there is essentially trivial sampling variability in the estimates of the distribution of non-concussive rotational accelerations due to the very large sample sizes. In contrast, the CDF of rotational acceleration at respective maximum likelihood estimates of the 57 concussive rotational accelerations exhibit substantial sampling variability in the estimates of the concussive rotational accelerations obtained from the HIT system (Fig. 3). Assuming that the only concussions among the 300,977 impacts are the 57 confirmed, diagnosed concussions as reported in Rowson et al. (2012), the posterior means and 95% credible intervals for the logistic regression intercept and slope relative risk of concussion to rotational acceleration (rad/s2) are 13.7 ( 15.2, 12.4) and 0.00193 (0.00166,0.00223) (see Table 1, column A). To illustrate the degree of uncertainty associated with this sampling variability on the estimation of the risk function, 50 draws of the risk function associated with β are shown in Fig. 4, along with the overall risk curve function estimated from the median value of β from 1000 draws. Column A of Table 2 provides the associated median and 95% credible intervals (in parentheses) for the injury risk curve for our synthetic dataset, with the rotational acceleration threshold associated with concussion risk p, given by
log(p / (1 − p)) − β0 β1
. Note that the
acceleration thresholds required to sustain a given injury risk are higher here than in the Rowson paper (Table 2, column H), because Rowson et al. (2012) weighted up the impacts resulting in concussions in their analysis so that they represented 0.0726% of all hits, rather than the 0.0189% of hits that were actual clinically diagnosed concussions. Weighting the data has the effect of shifting the risk
3062
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
Table 1 Posterior mean and 95% credible interval of slopes and intercepts for logistic regression model (2): (A) using only 57 observed concussions; (B) through (G): imputing additional concussions using three fractions of underreporting, with and without accounting for sensor error. Using 57 observed concussions (A)
Imputing 162 additional concussions
Imputing 52 additional concussions
Imputing 380 additional concussions
Without accounting for sensor error (B)
Accounting for sensor error (C)
Without accounting for sensor error (D)
Accounting for sensor error (E)
Accounting for senWithout accounting for sensor error sor error (G) (F)
Intercept ( β0 )
13.7 ( 15.2, 12.4)
11.7 ( 12.8, 10.6)
17.5 ( 21.5, 14.6)
12.9 ( 14.4, 11.7)
20.2 ( 26.0, 16.0)
10.2 ( 11.0, 9.5)
14.2 ( 16.8, 12.3)
Slope ( β1)
1.93 (1.66,2.23)
1.82 (1.56,2.08)
3.45 (2.69,4.37)
1.92 (1.64,2.24)
3.83 (2.86,5.07)
1.64 (1.42,1.85)
2.81 (2.24,3.43)
( ×10−3 )
Table 2 Posterior median and 95% credible interval of rotational acceleration (rad/s2) associated with a risk p of sustaining a concussion: (A) using only 57 observed concussions; (B) through (G): imputing additional concussions using three fractions of underreporting, with and without accounting for sensor error; (H) original estimates by Rowson et al. (2012) (note they were provided without confidence intervals). All values in rad/s2. P
0.1 0.25 0.5 0.75 0.9
Using 57 observed concussions (A)
5951 6525 7093 7658 8229
(5715,6220) (6226,6882) (6726,7542) (7230,8208) (7733,8860)
Imputing 162 additional concussions
Imputing 52 additional concussions
Imputing 380 additional concussions
Accounting for Without accounting for sensor error sensor error (C) (B)
Accounting for Without accounting for sensor error sensor error (E) (D)
Accounting for Without accounting for sensor error sensor error (G) (F)
5190 5794 6395 6994 7599
5560 6132 6700 7273 7846
4883 5553 6228 6899 7571
(5015,5430) (5532,6124) (6057,6837) (6588,7542) (7119,8249)
4437 4757 5077 5404 5725
(4322,4600) (4602,5018) (4861,5416) (5126,5817) (5380,6215)
(5359,5838) (5863,6526) (6360,7176) (6858,7821) (7340,8498)
curve to the left, so that lower-intensity impacts are predicted to have a greater injury risk than if the data were left unweighted. We address this issue of concussion under-reporting in the next section.
4680 4960 5247 5533 5816
(4529,4873) (4792,5231) (5030,5604) (5263,5974) (5475,6343)
(4719,5105) (5316,5869) (5912,6637) (6504,7407) (7107,8173)
4254 4646 5039 5438 5828
(4133,4490) (4460,4983) (4772,5470) (5094,5957) (5419,6445)
Original Rowson (2012) estimates (H)
5260 5821 6383 6945 7483
concussion risk of 90% is essentially unchanged compared to the assumption of a missing concussion rate of 0.0726%. 3.3. Accounting for sensor error
3.2. Accounting for unreported concussions The estimates for the posterior distributions of β obtained under (4) multiply imputing 162 additional concussions (to match the overall concussion rate assumed in Rowson et al. (2012)) are presented in Table 1 column B. The presence of under-reporting lowers the estimates of rotational acceleration required for a given concussion risk (compare columns A and B in Table 2), because the overall concussion rate in the imputed model is greater than for the un-imputed results; see also Fig. 5. Also note that the point estimates under this assumption are similar to those of Rowson et al. (compare Table 2 column B with Table 2 column H), although we now provide intervals that account for the uncertainty due to both sampling variability and the misclassification of concussions at the same rate as reported by Rowson et al. (0.0726%). Finally, to evaluate the sensitivity of the risk curve shape and thresholds to decreases or increases in the underlying actual concussion rate, we considered analyses in which the true rate of concussions was considered to be either one half (0.0363%, column D in Tables 1 and 2) or twice (0.1452%, column F in Tables 1 and 2) the rate reported by Rowson et al. (2012). Compared to the assumption of concussion rate of 0.0726% (column B, Tables 1 and 2), reducing the true rate of concussions to 0.0363% increases the rotational acceleration required for a 10% concussion risk by about 7%, and increases the rotational acceleration required associated with a 90% concussion risk by about 3%. When the true rate of concussions is increased to 0.1452%, the rotational acceleration associated with a 10% concussion risk is about 6% lower after accounting for underreporting, and the rotational acceleration associated with a
After accounting for the sensor measurement error as estimated by Allison et al. (2014) for underlying concussion rates of 0.0726% (162 missing concussions), 0.0363% (52 missing concussions), and 0.1452% (380 missing concussions), presented in columns C, E and G, respectively in Tables 1 and 2, we find strengthened relationships between rotational acceleration and injury risk, with steeper slopes in the risk curves. Specifically, accounting for the measurement error decreases the intercept and increases the slope, as it removes “bias toward the null” that the measurement error induces; at the same time the uncertainty in the association is increased, as evidenced by the widened confidence intervals, as we are now treating the rotational acceleration as having uncertainty. Fig. 6(a) illustrates this concept, with the steeper slope with wider confidence intervals after accounting for measurement error. Accounting for measurement error also reduces the rotational acceleration estimated to induce a 10%, 25%, 50%, 75%, and 90% risk of concussion at the three levels of assumed true concussion rates (see Fig. 7, and compare columns C, E and G in Table 2 which include measurement error correction to columns B, D, and F, respectively, which do not include sensor measurement error correction).
4. Discussion There are three sources of error to consider when assessing the relationship between rotational acceleration and concussion risk. The first is sampling variability. While the large number of impacts analyzed might seem to make this trivial, the fact that concussion is a rare outcome means that there remains substantial uncertainty in
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
Fig. 1. 20 empirical CDFs of posterior predictive distributions of rotation acceleration based on non-concussive HIT data. Red line gives CDF of rotational acceleration at maximum likelihood estimates of non-concussive HIT data.
Fig. 2. 20 empirical CDFs of posterior predictive distributions of rotation acceleration based on non-concussive 6DOF data. Red line gives CDF of rotational acceleration at maximum likelihood estimates of non-concussive 6DOF data.
the association between rotational acceleration and concussion risk even if there is no under-reporting or sensor error. The second source of error is under-reporting of the injury, resulting in concussion misclassification – that is, a number of concussions that occur are recorded as non-concussive impacts. The rate at which this occurs is unknown so a sensitivity analysis under differing assumptions of this rate must be made. Finally, all sensors measure rotational acceleration with some error, often due to limited precision or accuracy as well as non-rigid attachment to the head. If one assumes that this error is “non-differential” or not associated with concussion risk or
3063
Fig. 3. 20 empirical CDFs of posterior predictive distributions of rotation acceleration based on concussive HIT data. Red line gives CDF of rotational acceleration at maximum likelihood estimates of concussive HIT data.
Fig. 4. Posterior median of risk of concussion based on 1000 draws of risk function parameters β (red), together with 50 draws of the risk function based on 50 draws of β (black).
acceleration level, presence of this error in the data will tend to blur any relationship between load and injury, and bias associations between acceleration and concussion toward the null. Based on initial point estimates for concussion risk calculated by Rowson et al. (2012), accounting only for sampling error suggests that rotational acceleration of 5715–6220 rad/s2 would be sufficient to provide a 10% concussion risk, 6726–7542 rad/s2 would be sufficient to provide a 50% concussion risk, and 7733– 8860 rad/s2 would be sufficient to provide a 90% concussion risk. Accounting for injury under-reporting shifts and widens these ranges to 4719–5838, 5912–7176, and 7107–8498 rad/s2
3064
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
respectively, depending on the degree of mis-classification. The decrease in acceleration threshold is a function of the fact that more concussions are assumed to occur, and the increase in width of the estimated interval associated with a given risk is due to the increase in uncertainty due to the misclassification as well as the range of misclassification rate explored. Finally, accounting for measurement error in the rotational acceleration measures yields estimates of 4133–4873, 4772–5604, and 5380–6445 rad/s2 respectively for 10%, 50%, and 90% concussion risk, respectively, again depending on the degree of misclassification assumed. We have a further decline in rotational acceleration required for a given injury risk due to the fact that the assumed skewness of the error terms means that observed high accelerations are more likely to actually represent smaller true accelerations than vice-
Fig. 5. Posterior median of risk of concussion based on 1000 draws of risk function parameters β : (A) using only 57 observed concussions; (B) multiply imputing 162 additional concussions under missing at random assumption.
versa. The cumulative effect of these three forms of error is demonstrated in Fig. 8 for both 10% and 50% concussion risk. By accounting for the cumulative effect of these three forms of error or uncertainty, the final model describes the fundamental biomechanical relationship between injury and acceleration if a perfect sensor – i.e. one with no error – was utilized. This hypothetical scenario suggests how associations might appear as sensors improve and measurement error is reduced. Until a sensor with
Fig. 7. Posterior median of risk of concussion based on 200 draws of risk function parameters β : (B) multiply imputed 162 unobserved concussions under missing at random assumption without measurement error, (D) multiply imputed 52 unobserved concussions under missing at random assumption without measurement error, (F) multiply imputed 380 unobserved concussions under missing at random assumption without measurement error, (C) multiply imputed 162 unobserved concussions under missing at random assumption with measurement error, (E) multiply imputed 52 unobserved concussions under missing at random assumption with measurement error, and (G) multiply imputed 380 unobserved concussions under missing at random assumption with measurement error.
Fig. 6. Posterior mean and associated 95% credible intervals for risk of concussion, multiply imputing 162 unobserved concussions under missing at random assumption, without accounting for measurement error (black) and removing measurement error (red), (a) on log-odds scale, (b) on probability scale.
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
3065
Appendix A. Supplementary Information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jbiomech.2015.07.026
References
Fig. 8. Cumulative effect of three types of error on the 10% risk (solid bars) and 50% (grey bars) of concussion. The left most data points represent the point estimates from Rowson et al. (2012) and moving to the right, represent the progressive effect of adding sampling variability, under-reporting, and accounting for sensor error. The ranges shown represent the smallest lower end of a 95% confidence interval with the largest upper end of a 95% confidence interval across the three levels of under-reporting considered.
negligible error exists, we recommend that future reports should account for both sampling variability and concussion underreporting in their measured data to predict injury risk from observed rotational acceleration and should incorporate sensor error to define fundamental biomechanical relationships from such data. 4.1. Limitations There are a number of assumptions in this analysis. First, for concussion under-reporting, we have assumed (a) that misclassification is monotonic – that is, there are true concussions classed as non-concussions but no true non-concussions classed as concussions, (b) the misclassification rate is known (although allowed to vary over three rates as a sensitivity analysis), and (c) random, conditional on rotational acceleration. Second, we assume that rotational acceleration sensor measurement error is non-differential (again, random conditional on true rotational acceleration), multiplicative and log-normally distributed, with known variance. This assumption differs from a previous report by Funk et al. (2012), in which sensor error was presumed to occur predominantly at higher accelerations. Third, we utilized sensor error calculated on a hockey–helmet based system on this analysis of football data. While the magnitude of the sensor error for the football and hockey based systems is similar, the specific value of the error for football may be slightly different. However, for illustrative purposes, the effect of sensor error on the injury risk curves is adequately demonstrated by the hockey-based helmet system sensor error. Finally, because we are generating hypothetical posterior population predictive distributions using results of Rowson et al. (2012), we assumed a linear relationship between the log-odds of concussion and rotational acceleration.
Conflict of interest statement Kristy B. Arbogast serves as a consultant for the National Football League Players Association on topics related to head injury and player safety.
Allison, M.A., Kang, Y.S., Bolte 4th, J.H., Maltese, M.R., Arbogast, K.B., 2014. Validation of a helmet-based system to measure head impact biomechanics in ice hockey. Med. Sci. Sports Exerc. 46, 115–123. Allison, M.A., Kang, Y.S., Bolte 4th, J.H., Maltese, M.R., Arbogast, K.B., 2015. Measurement of Hybrid III impact kinematics using an accelerometer and gyroscope system in ice hockey helmets. Ann. Biomed. Eng. 43, 1896–1906. Bartsch, A., Samorezov, S., Benzel, E., Miele, V., Brett, D., 2014. Validation of an “intelligent mouthguard” single event head impact dosimeter. Stapp Car Crash J. 58, 1–28. Brainard, L.L., Beckwith, J.G., Chu, J.J., Crisco, J.J., McAllister, T.W., Duhaime, A.C., Maerlender, A.C., Greenwald, R.M., 2012. Gender differences in head impacts sustained by collegiate ice hockey players. Med. Sci. Sports Exerc. 44, 297–304. Carrol, R.J., Ruppert, D., Stefanski, L.A., 1995. Measurement Error in Nonlinear Models: A Modern Perspective. Chapman and Hall/CRC Press, Boca Raton, FL. Cook, J.R., Stefanski, L.A., 1994. Simulation-extrapolation estimation in parametric measurement error models. J. Am. Stat. Assoc. 89, 1314–1328. Crisco, J.J., Fiore, R., Beckwith, J.G., Chu, J.J., Brolinson, P.G., Duma, S., McAllister, T. W., Duhaime, A.C., Greenwald, R.M., 2010. Frequency and location of head impact exposures in individual collegiate football players. J. Athl. Train. 45, 549–559. Funk, J.R., Rowson, S., Daniel, R.W., Duma, S.M., 2012. Validation of concussion risk curves for collegiate football players derived from HITS data. Ann. Biomed. Eng. 40, 79–89. Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B., 2014. Bayesian Data Analysis, 3rd edition. CDC Press, Boca Raton, FL. Hernandez, F., Wu, L.C., Yip, M.C., Laksari, K., Hoffman, A.R., Lopez, J.R., Grant, G.A., Kleivem, S., Camarillo, D.B., 2014. Six degree-of-freedom measurements of human mild traumatic brain injury. Ann. Biomed. Eng. 43, 1918–1934. Hosmer, D.W., Lemeshow, S., Sturdivant, R.X., 2013. Applied Logistic Regression, 3rd Edition Wiley Press, Hoboken, NJ. Jadischke, R., Viano, D.C., Dau, N., King, A.I., McCarthy, J., 2013. On the accuracy of the head impact telemetry (HIT) system used in football helmets. J. Biomech. 46, 2310–2315. King, D., Hume, P.A., Brughelli, M., Gisanne, C., 2015. Instrumented mouthguard acceleration analyses for head impacts in amateur rugby union players over a season of matches. Am. J. Sports Med. 43, 614–624. Little, R.J.A., Rubin, D.B., 2002. Statistical Analysis with Missing Data, 2nd edition. Wiley, New York. McCrea, M., Hammeke, T., Olsen, G., Leo, P., Guskiewicz, K., 2004. Unreported concussion in high school football players: implications for prevention. Clin. J. Sport Med. 14, 13–17. Mihalik, J., Bell, D., Marshall, S., Guskiewicz, K., 2007. Measurement of head impacts in collegiate football players: an investigation of positional and event-type differences. Neurosurgery 6, 1229–1235. Pellman, E.J., Viano, D.C., Tucker, A.M., Casson, I.R., Waeckerle, J.F., 2003. Concussion in professional football: reconstruction of game impacts and injuries. Neurosurgery 53, 799–812. Rice, S.O., 1945. Mathematical analysis of random noise. Bell Syst. Tech. J. 24, 46–156. Rowson, S., Brolinson, G., Goforth, M., Dietter, D., Duma, S., 2009. Linear and angular head acceleration measurements in collegiate football. J. Biomech. Eng. 131, 061016. Rowson, S., Duma, S.M., 2011. Development of the STAR evaluation system for football helmets: integrating player head impact exposure and risk of concussion. Ann. Biomed. Eng. 39, 2130–2140. Rowson, S., Duma, S.M., 2013. Brain injury prediction: assessing the combined probability of concussion using linear and rotational head acceleration. Ann. Biomed. Eng. 41, 873–882. Rowson, S., Duma, S.M., Beckwith, J.G., Chu, J.J., Greenwald, R.M., Crisco, J.J., Brolinson, P.G., Duhaime, A.C., McAllister, T.W., Maerlender, A.C., 2012. Rotational head kinematics in football impacts: an injury risk function for concussion. Ann. Biomed. Eng. 40, 1–13. Stefanski, L.A., Cook, J.R., 1995. Simulation-extrapolation: the measurement error jackknife. J. Am. Stat. Assoc. 90, 1247–1256. Wilcox, B.J., Beckwith, J.G., Greenwald, R.M., Raukar, N.P., Chu, J.J., McAllister, T.W., Flashman, L.A., Maerlender, A.C., Duhaime, A.C., Crisco, J.J., 2015. Biomechanics of head impacts associated with diagnosed concussion in female collegiate ice hockey players. J. Biomech. 48, 2201–2204.
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Accounting for sampling variability, injury under-reporting, and sensor error in concussion injury risk curves Michael R. Elliott a,b, Susan S. Margulies c, Matthew R. Maltese d, Kristy B. Arbogast e,n a
Department of Biostatistics, School of Public Health, University of Michigan, Ann Arbor, MI 48109, United States Survey Methodology Program, Institute for Social Research, University of Michigan, Ann Arbor, MI 48109, United States c Department of Bioengineering, University of Pennsylvania, Philadelphia, PA 19104, United States d Department of Anesthesiology and Critical Care Medicine, Children's Hospital of Philadelphia, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104, United States e Center for Injury Research and Prevention, Children’s Hospital of Philadelphia, Department of Pediatrics, University of Pennsylvania, 34th and Civic, Center Blvd, Suite 1150, Philadelphia, PA 19104, United States b
art ic l e i nf o
a b s t r a c t
Article history: Accepted 23 July 2015
There has been recent dramatic increase in the use of sensors affixed to the heads or helmets of athletes to measure the biomechanics of head impacts that lead to concussion. The relationship between injury and linear or rotational head acceleration measured by such sensors can be quantified with an injury risk curve. The utility of the injury risk curve relies on the accuracy of both the clinical diagnosis and the biomechanical measure. The focus of our analysis was to demonstrate the influence of three sources of error on the shape and interpretation of concussion injury risk curves: sampling variability associated with a rare event, concussion under-reporting, and sensor measurement error. We utilized Bayesian statistical methods to generate synthetic data from previously published concussion injury risk curves developed using data from helmet-based sensors on collegiate football players and assessed the effect of the three sources of error on the risk relationship. Accounting for sampling variability adds uncertainty or width to the injury risk curve. Assuming a variety of rates of unreported concussions in the non-concussed group, we found that accounting for under-reporting lowers the rotational acceleration required for a given concussion risk. Lastly, after accounting for sensor error, we find strengthened relationships between rotational acceleration and injury risk, further lowering the magnitude of rotational acceleration needed for a given risk of concussion. As more accurate sensors are designed and more sensitive and specific clinical diagnostic tools are introduced, our analysis provides guidance for the future development of comprehensive concussion risk curves. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Sensor error Injury risk Traumatic brain injury
1. Introduction Helmet-based head impact sensors have been widely utilized to understand the link between the biomechanics of head impact and clinical outcomes of concussion in athletes (e.g. Rowson et al., 2009; Brainard et al., 2012; Crisco et al., 2010; Mihalik et al., 2007; Wilcox et al., 2015). Such sensors have accelerometers and in some cases, gyroscopes embedded in the padding or attached to the helmet shell to estimate the magnitude of linear and rotational acceleration experienced by the athlete during head impact. Similar systems have been developed for non-helmeted sports and are integrated into headbands, skull-caps or mouthguards or directly attached to the athlete near the mastoid process (King n
Corresponding author. Tel.: þ 1 215 590 6075; fax: þ1 215 590 5425. E-mail address: [email protected] (K.B. Arbogast).
http://dx.doi.org/10.1016/j.jbiomech.2015.07.026 0021-9290/& 2015 Elsevier Ltd. All rights reserved.
et al., 2015; Hernandez et al., 2014; Bartsch et al., 2014). Quantifying the relationship between biomechanical input and clinical outcome is critical to the advancement of concussion prevention principles including the assessment of injury risk, the design of protective equipment such as helmets, and the development of training and policies intended to limit exposure to head impacts and injury risk. The most common approach to quantifying the link between biomechanical input and concussion is through injury risk curves. Injury risk curves describe the injury probability given a specific mechanical input – for example, concussion risk given a particular head acceleration. This relationship is not linear but rather sigmoid, often in the form of a Weibull distribution or some other cumulative distribution function (CDF). Previous research has constructed injury risk curves for concussion. Pellman et al. (2003) used head impact data from the NFL and Rowson and Duma (2011) used collegiate football head impact data to describe the relationship
3060
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
between linear acceleration and concussion risk. Recognizing that rotational acceleration or velocity is an important part of the causal pathway for concussion, Rowson et al. (2012) expanded this work and quantified a relationship between rotational acceleration and concussion risk. Specifically, they used data from 300,977 head impacts among US football players (57 of which resulted in diagnosed concussions) to estimate a logistic risk function for concussion as a function of rotational acceleration. The same authors proposed a combined injury risk function that incorporated both linear and rotational components (Rowson and Duma, 2013). The utility of injury risk curves relies on accuracy of the metrics that lie on both the abscissa and ordinate axes: accuracy of clinical diagnosis (dichotomous variable on the ordinate) and accuracy of biomechanical measure (continuous variable on the abscissa). There are three sources of uncertainty or error that influence concussion injury risk curve development: sampling variability of the clinical outcome, concussion under-reporting, and head impact sensor error. The effect of these errors on injury risk curves has not been rigorously explored in the literature and therefore detailed technical treatment of these three issues is the focus of this analysis. For exemplar purposes, we utilize the Rowson et al. (2012) data as the basis for our analyses however these approaches can be applied to other injury risk curves. Each source of error is discussed below. First, injury risk curves do not typically account for uncertainty in the finite sample of head impacts. This might seem minor, given that typically large numbers of impacts are collected over a given period of games and practices. For example, 300,977 observations were made in the Rowson et al. (2012) analysis. However, concussion remains a rare event (e.g. Rowson et al. report 57 diagnosed concussions), thus there is substantial sampling variability that influences the relationship between concussion risk and rotational acceleration. Standard logistic regression methods (Hosmer et al., 2013) will account for this sampling variability; which is non-trivial when the outcome is rare. A second source of uncertainty is misclassification of the concussion outcome itself, with under-reporting estimated to be as high as 53% (McCrea et al., 2004). The exact under-reporting proportion is not known so parametric analysis using a rational range of estimates is important. A third issue is that the measure of acceleration obtained from helmet-based head impact systems is subject to non-trivial error levels between what is measured by the sensor and the true head kinematics (Allison et al., 2014, 2015; Jadischke et al., 2013; Funk et al., 2012). If this error is non-differential (that is, not correlated with the true acceleration), as suggested by Allison et al. (2014), then accounting for it in the creation of injury risk curves will strengthen the estimated relationship between rotational acceleration and concussion risk, and lead to more accurate concussion thresholds. Therefore, the objective of this work is to demonstrate, using accepted analytical methods, the magnitude of the effect of these three sources of error on the shape and interpretation of concussion injury risk curves. This analysis provides a roadmap for future research efforts as new sensors are designed and clinical diagnostic tools become more accurate. We estimate concussion risk curves as a function of rotational acceleration and identify acceleration thresholds associated with a specific injury risk that can be derived from these relationships. 2. Methods We use summative published data as the basis for our assumptions about true concussion rates and information about rotational acceleration measurement error (e.g. Rowson et al., 2012; Allison et al., 2014) and generate synthetic datasets to provide a platform for statistical analysis to assess the sources of uncertainty related to sampling, under-reporting, and sensor error. The synthetic dataset was developed using the Bayesian statistical principles, which allow generation of posterior predictive
distributions that can be used to obtain the injury risk curve measures of interest. The synthetic origin of the dataset is of minor importance; rather, the major objective of this communication was to demonstrate the importance of accounting for uncertainty in conclusions drawn from the dataset.
2.1. Dataset generation Using statistics about the distributions among the concussed and non-concussed individuals in the Rowson et al. data (Rowson et al., 2012, Table 1), we can use Bayesian methods to estimate the relationship between concussion risk and rotational acceleration described by their data. We use their reported distributions to generate posterior predictive distributions (synthetic data) of the rotational acceleration by observed concussion status for analysis. Specifically, we construct a posterior predictive distribution of the rotational acceleration of non-concussive injuries from the HIT system and 6DOF measurement device,1 and of the concussive injuries, also measured through the HIT system. A posterior predictive distribution is a Bayesian concept, and can be thought of as a draw of a future set of observations, yrep , conditional on observed data y, and assuming a model for the data f (y|θ ) and a prior distribution p(θ ) for the model parameters (Gelman et al., 2014). The mathematical formalism is as follows:
p(y rep |y ) =
∫ f (yrep |θ )p(θ |y)dθ
(1)
where p(θ|y ) = f (y|θ )p(θ )/ ∫ f (y|θ )p(θ )dθ is the posterior distribution of θ . The integral in (1) is usually not computable analytically, but p(yrep |y ) can be obtained by simulation by drawing θ rep~p(θ|y ), and then yrep ~f (y|θ rep). This approach provides a principled way to accommodate uncertainty in the data distributions of head rotational acceleration measures in football players. Although we do not actually have the data from Rowson et al., we can still approximate the posterior distribution of θ as N (θ^, J −1(θ^ )) , where θ^ is the maximum likelihood estimator of θ , and J (θ ) is the Fisher information matrix associated with the distribution of y , computed as J (θ ) = E (−
∂2 ∂θ ∂θ T
log f (y|θ )) for the assumed
distribution f (y|θ ) . 2.2. Assessing sampling variability Because concussion is a rare event when considering all sensor-detected head movements on the sports field, conclusions on the relationship between injury and head acceleration may be subject to substantial uncertainty due to the variability in a finite dataset attributable to the statistical likelihood associated with an underlying rate of concussion. We can assess this uncertainty by creating many examples of our synthetic population and observing x concussions in one example and y (slightly more or less than x) concussions in another. The non-concussed data is assumed to follow a Weibull distribution; the concussed data is assumed to follow a Rician distribution (Rice, 1945). The large sample size of the non-concussed data means that the normal approximation will be highly accurate; for the concussed data, the smaller sample size is at least partially offset by the fact that the Rician distribution will be approximately normal for the parameter values reported. To compute an estimated concussion risk as a function of rotational acceleration, we follow Rowson et al. and assume a logistic regression model, so that
⎛ e β0 + β1y ⎞ ⎟⎟ P (Z = 1|Y = y ) = ⎜⎜ ⎝ 1 + e β0 + β1y ⎠
(2)
where Z is an indicator for concussion, Y is the rotational acceleration, and Z is set deterministically to 1 if Y is drawn from the Rician distribution and 0 if drawn from the Weibull distribution. We can obtain inference about the parameters β that govern the risk function (2) by extending the posterior predictive distribution in (1):
p(β|y ) =
∫ p(β|yrep )p(yrep |θ )p(θ |y)dθdyrep .
(3)
We can obtain draws from the posterior distribution in (3) under the assumption of a flat prior by fitting a logistic regression model to a draw of yrep rep rep rep rep from (1), and then drawing β ~N (β^ , I −1(β^ )), where β^ and I −1(β^ ) are
obtained from logistic regression using yrep . In this analysis we obtain 1000 separate samples or “draws” of parameters of interest to make inference (compute posterior means, credible intervals [the Bayesian analog of confidence intervals]). Assuming that the distributional assumptions are correct, we are able to obtain estimates of the concussion/rotational acceleration risk that account for the influence of sampling variability without any loss of information. We use these estimates to develop concussion risk curves as a function of rotational acceleration that
1 Two different variants of a similar head impact sensor were used in the Rowson study. Most data were collected by the commercially available HIT system (Simbex, Lebanon, NH) however some data were collected using a version modified for research use that implements additional sensors towards improving accuracy.
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065 incorporate uncertainty due to sampling variability, and compare them to the risk curves that do not account for this uncertainty.
The Bayesian paradigm is also particularly convenient for handling random under-reporting of concussion in our subjects in the non-concussed group, because it is the basis behind accepted missing data methodological approaches (Little and Rubin, 2002). We thus extend our work accounting for uncertainty in the association between rotational acceleration and concussion risk to account for underreporting of concussions by identifying missing concussions in the non-concussed group and “converting” them to the concussed group. We assume that missing concussions are missing at random (Little and Rubin, 2002) conditional on the rotational acceleration distribution. We require an assumption about the fraction of concussions that are missed. We consider three fractions: (1) a fraction based on the overall concussion rate assumed in Rowson et al. (2012) – 0.0726%; (2) a fraction based on one-half of that assumed in Rowson et al. (2012), and (3) a fraction twice that assumed in Rowson et al. (2012). If the true concussion rate in the sample is 0.0726%, then there are a total of 219 concussed players, and 162 ( ¼ 219–57) of the 300,920 “non-concussed” players are actually concussed. We assume that there are no misclassified observed concussions, that is, all rotational acceleration observations Y generated from the Rician distribution are associated with concussions (Z ¼1). Next, we assume that the unobserved concussions are drawn from the same distribution of Y as the observed concussions. The posterior distribution of the logistic parameters governing the relationship between rotational acceleration and concussion risk is given by
∫ p(β|yrep , z rep)p(yrep |θ )p(θ |y)dθdyrep z rep
(4)
Draws from p(β|yrep , z rep) are obtained using a Monte Carlo Markov Chain (MCMC) (Gelman et al., 2014). The full algorithm proceeds as follows: 1. Draw θ rep~p(θ|y ) from N (θ^, J −1(θ^ )) , where 286,579 draws are from the Weibull distribution estimated using non-concussive rotational accelerations obtained from the HIT system, 14,341 draws are from the Weibull distribution estimated using non-concussive rotational accelerations obtained from the 6DOF system, and 57 draws are from the Rician distribution estimated using concussive rotational accelerations. 2. Draw yrep ~f (y|θ rep) under the appropriate distributions. 3. Obtain a draw from p(β|y
rep
,z
rep
) using MCMC by iterating as follows:
a. Impute the sth z rep, s|yrep , β (s − 1) , where β (s − 1) is the draw of β from the previous iteration of the MCMC ( β (0) is set to be previous saved draw in the overall algorithm, i.e., the final drawn from the previous MCMC). z rep, s is set to 1 if
yrep P (Z
is
rep
= 1|Y
drawn rep
)=
from
the
Rician
exp(β (s − 1) + β (s − 1)Y rep) 0 1 , 1 + exp(β (s − 1) + β (s − 1)Y rep) 0 1
distribution;
otherwise
subject to the constraints that
∑i zirep, s = 219. b. Draw the sth β (s)|yrep , z rep, s by first fitting a logistic regression model to draw rep rep rep of z rep, s and yrep , and then drawing β ~N (β^ , I −1(β^ )), where β^ and
I −1(β^
rep
) are obtained from this logistic regression.
c. Repeat (a) and (b) for S times, until the autocorrelation between the draws of early values of β (s) and later values of β (s) is trivial. We chose S = 100 in the analysis below. To account for the uncertainty in the actual concussion rate, we repeat both of the above analyses assuming a true concussion rate twice that (0.1452%) and onehalf that (0.0363%) assumed in Rowson et al. (2012). This yields an estimate of an additional 380 and 52 missed concussions in the “non-concussed” group, respectively, above the 57 reported concussions. 2.4. Accounting for sensor error We use the SIMEX method to account for the measurement error (Cook and Stefanski, 1994). SIMEX, short for Simulation Extraction, begins from the observation that the ordinary least squares estimator of a regression slope Z = β0 + βxX + ε estimates
βx σ 2 σ 2 + σu2
do here, estimate an approximate quadratic regression model β^(λ ) = γ1 + γ2λ + γ3λ2, computing β^x = γ^1 − γ^2 + γ3 (Carrol et al., 1995). Variance estimates of β^x are computed via a bootstrap (Stefanski and Cook, 1995). We assume a multiplicative error form for the rotational acceleration data – specifically, the observed rotational acceleration Y equals XU , where U follows a log-normal distribution with expec-
2.3. Accounting for under-reporting
p(β|y ) =
3061
when X is replaced with Y = X + U , where σu2 is the variance of the
(mean zero) measurement error U and σ 2 is the variance of the true residual error ε . By adding error with variance λσ 2 to W, one can compute a series of values β^(λ ) u
from different values of λ as a function of
βx σ 2 σ 2 + (1 + λ )σu2
and extrapolate to β^( − 1) ,
which estimates βx . One can work directly with this non-linear equation, or, as we
2
tation 1 and variance e σ − 1. Thus log(U ) is generated from a normal distribution with mean −σ /2 and variance σ 2 , and Y is replaced with exp(log(Y ) + λ log(U )) to 2
implement SIMEX. We estimate σ 2 as 0.2706 by method of moments using the alldirection impact data from Allison et al. (2014), where instrumented anthropomorphic test devices (ATDs) fitted with HIT sensor-outfitted helmets were impacted, and helmet rotational accelerations were compared to values measured within the ATD head. We embed this helmet sensor measurement error in our multiple imputation method by replacing the draw the β (s)|yrep , z rep, s at step (3b) (s ) (s ) (s ) (s ) with a draw of βx(s)|yrep , z rep, s from N (β^x , cov(β^x )) , where β^x and cov(β^x ) are
obtained from the SIMEX algorithm utilizing the draws yrep , z rep, s .
3. Results The consolidated results of all simulations with and without accounting for the uncertainties due to sample variability, underreporting and sensor error are provided in Tables 1 and 2, with the mean logistic regression parameters in Table 1 and injury risk curve data in Table 2. Importantly, Table 2 presents the median values for rotational acceleration thresholds (and confidence intervals) for a range of injury risk levels, from 10 to 90%. To orient the reader, column A in Tables 1 and 2 is the synthetic dataset (based on Rowson et al. (2012) published observations of 57 concussions of 300,977 impacts), without accounting for any sample variability, under-reporting, or sensor measurement error. 3.1. Sampling variability The empirical CDFs of 20 draws from the posterior predictive distributions of the 286,579 non-concussive rotational accelerations obtained from the HIT system (Fig. 1) and the 14,341 non-concussive rotational accelerations obtained from the 6DOF system (Fig. 2) show that there is essentially trivial sampling variability in the estimates of the distribution of non-concussive rotational accelerations due to the very large sample sizes. In contrast, the CDF of rotational acceleration at respective maximum likelihood estimates of the 57 concussive rotational accelerations exhibit substantial sampling variability in the estimates of the concussive rotational accelerations obtained from the HIT system (Fig. 3). Assuming that the only concussions among the 300,977 impacts are the 57 confirmed, diagnosed concussions as reported in Rowson et al. (2012), the posterior means and 95% credible intervals for the logistic regression intercept and slope relative risk of concussion to rotational acceleration (rad/s2) are 13.7 ( 15.2, 12.4) and 0.00193 (0.00166,0.00223) (see Table 1, column A). To illustrate the degree of uncertainty associated with this sampling variability on the estimation of the risk function, 50 draws of the risk function associated with β are shown in Fig. 4, along with the overall risk curve function estimated from the median value of β from 1000 draws. Column A of Table 2 provides the associated median and 95% credible intervals (in parentheses) for the injury risk curve for our synthetic dataset, with the rotational acceleration threshold associated with concussion risk p, given by
log(p / (1 − p)) − β0 β1
. Note that the
acceleration thresholds required to sustain a given injury risk are higher here than in the Rowson paper (Table 2, column H), because Rowson et al. (2012) weighted up the impacts resulting in concussions in their analysis so that they represented 0.0726% of all hits, rather than the 0.0189% of hits that were actual clinically diagnosed concussions. Weighting the data has the effect of shifting the risk
3062
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
Table 1 Posterior mean and 95% credible interval of slopes and intercepts for logistic regression model (2): (A) using only 57 observed concussions; (B) through (G): imputing additional concussions using three fractions of underreporting, with and without accounting for sensor error. Using 57 observed concussions (A)
Imputing 162 additional concussions
Imputing 52 additional concussions
Imputing 380 additional concussions
Without accounting for sensor error (B)
Accounting for sensor error (C)
Without accounting for sensor error (D)
Accounting for sensor error (E)
Accounting for senWithout accounting for sensor error sor error (G) (F)
Intercept ( β0 )
13.7 ( 15.2, 12.4)
11.7 ( 12.8, 10.6)
17.5 ( 21.5, 14.6)
12.9 ( 14.4, 11.7)
20.2 ( 26.0, 16.0)
10.2 ( 11.0, 9.5)
14.2 ( 16.8, 12.3)
Slope ( β1)
1.93 (1.66,2.23)
1.82 (1.56,2.08)
3.45 (2.69,4.37)
1.92 (1.64,2.24)
3.83 (2.86,5.07)
1.64 (1.42,1.85)
2.81 (2.24,3.43)
( ×10−3 )
Table 2 Posterior median and 95% credible interval of rotational acceleration (rad/s2) associated with a risk p of sustaining a concussion: (A) using only 57 observed concussions; (B) through (G): imputing additional concussions using three fractions of underreporting, with and without accounting for sensor error; (H) original estimates by Rowson et al. (2012) (note they were provided without confidence intervals). All values in rad/s2. P
0.1 0.25 0.5 0.75 0.9
Using 57 observed concussions (A)
5951 6525 7093 7658 8229
(5715,6220) (6226,6882) (6726,7542) (7230,8208) (7733,8860)
Imputing 162 additional concussions
Imputing 52 additional concussions
Imputing 380 additional concussions
Accounting for Without accounting for sensor error sensor error (C) (B)
Accounting for Without accounting for sensor error sensor error (E) (D)
Accounting for Without accounting for sensor error sensor error (G) (F)
5190 5794 6395 6994 7599
5560 6132 6700 7273 7846
4883 5553 6228 6899 7571
(5015,5430) (5532,6124) (6057,6837) (6588,7542) (7119,8249)
4437 4757 5077 5404 5725
(4322,4600) (4602,5018) (4861,5416) (5126,5817) (5380,6215)
(5359,5838) (5863,6526) (6360,7176) (6858,7821) (7340,8498)
curve to the left, so that lower-intensity impacts are predicted to have a greater injury risk than if the data were left unweighted. We address this issue of concussion under-reporting in the next section.
4680 4960 5247 5533 5816
(4529,4873) (4792,5231) (5030,5604) (5263,5974) (5475,6343)
(4719,5105) (5316,5869) (5912,6637) (6504,7407) (7107,8173)
4254 4646 5039 5438 5828
(4133,4490) (4460,4983) (4772,5470) (5094,5957) (5419,6445)
Original Rowson (2012) estimates (H)
5260 5821 6383 6945 7483
concussion risk of 90% is essentially unchanged compared to the assumption of a missing concussion rate of 0.0726%. 3.3. Accounting for sensor error
3.2. Accounting for unreported concussions The estimates for the posterior distributions of β obtained under (4) multiply imputing 162 additional concussions (to match the overall concussion rate assumed in Rowson et al. (2012)) are presented in Table 1 column B. The presence of under-reporting lowers the estimates of rotational acceleration required for a given concussion risk (compare columns A and B in Table 2), because the overall concussion rate in the imputed model is greater than for the un-imputed results; see also Fig. 5. Also note that the point estimates under this assumption are similar to those of Rowson et al. (compare Table 2 column B with Table 2 column H), although we now provide intervals that account for the uncertainty due to both sampling variability and the misclassification of concussions at the same rate as reported by Rowson et al. (0.0726%). Finally, to evaluate the sensitivity of the risk curve shape and thresholds to decreases or increases in the underlying actual concussion rate, we considered analyses in which the true rate of concussions was considered to be either one half (0.0363%, column D in Tables 1 and 2) or twice (0.1452%, column F in Tables 1 and 2) the rate reported by Rowson et al. (2012). Compared to the assumption of concussion rate of 0.0726% (column B, Tables 1 and 2), reducing the true rate of concussions to 0.0363% increases the rotational acceleration required for a 10% concussion risk by about 7%, and increases the rotational acceleration required associated with a 90% concussion risk by about 3%. When the true rate of concussions is increased to 0.1452%, the rotational acceleration associated with a 10% concussion risk is about 6% lower after accounting for underreporting, and the rotational acceleration associated with a
After accounting for the sensor measurement error as estimated by Allison et al. (2014) for underlying concussion rates of 0.0726% (162 missing concussions), 0.0363% (52 missing concussions), and 0.1452% (380 missing concussions), presented in columns C, E and G, respectively in Tables 1 and 2, we find strengthened relationships between rotational acceleration and injury risk, with steeper slopes in the risk curves. Specifically, accounting for the measurement error decreases the intercept and increases the slope, as it removes “bias toward the null” that the measurement error induces; at the same time the uncertainty in the association is increased, as evidenced by the widened confidence intervals, as we are now treating the rotational acceleration as having uncertainty. Fig. 6(a) illustrates this concept, with the steeper slope with wider confidence intervals after accounting for measurement error. Accounting for measurement error also reduces the rotational acceleration estimated to induce a 10%, 25%, 50%, 75%, and 90% risk of concussion at the three levels of assumed true concussion rates (see Fig. 7, and compare columns C, E and G in Table 2 which include measurement error correction to columns B, D, and F, respectively, which do not include sensor measurement error correction).
4. Discussion There are three sources of error to consider when assessing the relationship between rotational acceleration and concussion risk. The first is sampling variability. While the large number of impacts analyzed might seem to make this trivial, the fact that concussion is a rare outcome means that there remains substantial uncertainty in
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
Fig. 1. 20 empirical CDFs of posterior predictive distributions of rotation acceleration based on non-concussive HIT data. Red line gives CDF of rotational acceleration at maximum likelihood estimates of non-concussive HIT data.
Fig. 2. 20 empirical CDFs of posterior predictive distributions of rotation acceleration based on non-concussive 6DOF data. Red line gives CDF of rotational acceleration at maximum likelihood estimates of non-concussive 6DOF data.
the association between rotational acceleration and concussion risk even if there is no under-reporting or sensor error. The second source of error is under-reporting of the injury, resulting in concussion misclassification – that is, a number of concussions that occur are recorded as non-concussive impacts. The rate at which this occurs is unknown so a sensitivity analysis under differing assumptions of this rate must be made. Finally, all sensors measure rotational acceleration with some error, often due to limited precision or accuracy as well as non-rigid attachment to the head. If one assumes that this error is “non-differential” or not associated with concussion risk or
3063
Fig. 3. 20 empirical CDFs of posterior predictive distributions of rotation acceleration based on concussive HIT data. Red line gives CDF of rotational acceleration at maximum likelihood estimates of concussive HIT data.
Fig. 4. Posterior median of risk of concussion based on 1000 draws of risk function parameters β (red), together with 50 draws of the risk function based on 50 draws of β (black).
acceleration level, presence of this error in the data will tend to blur any relationship between load and injury, and bias associations between acceleration and concussion toward the null. Based on initial point estimates for concussion risk calculated by Rowson et al. (2012), accounting only for sampling error suggests that rotational acceleration of 5715–6220 rad/s2 would be sufficient to provide a 10% concussion risk, 6726–7542 rad/s2 would be sufficient to provide a 50% concussion risk, and 7733– 8860 rad/s2 would be sufficient to provide a 90% concussion risk. Accounting for injury under-reporting shifts and widens these ranges to 4719–5838, 5912–7176, and 7107–8498 rad/s2
3064
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
respectively, depending on the degree of mis-classification. The decrease in acceleration threshold is a function of the fact that more concussions are assumed to occur, and the increase in width of the estimated interval associated with a given risk is due to the increase in uncertainty due to the misclassification as well as the range of misclassification rate explored. Finally, accounting for measurement error in the rotational acceleration measures yields estimates of 4133–4873, 4772–5604, and 5380–6445 rad/s2 respectively for 10%, 50%, and 90% concussion risk, respectively, again depending on the degree of misclassification assumed. We have a further decline in rotational acceleration required for a given injury risk due to the fact that the assumed skewness of the error terms means that observed high accelerations are more likely to actually represent smaller true accelerations than vice-
Fig. 5. Posterior median of risk of concussion based on 1000 draws of risk function parameters β : (A) using only 57 observed concussions; (B) multiply imputing 162 additional concussions under missing at random assumption.
versa. The cumulative effect of these three forms of error is demonstrated in Fig. 8 for both 10% and 50% concussion risk. By accounting for the cumulative effect of these three forms of error or uncertainty, the final model describes the fundamental biomechanical relationship between injury and acceleration if a perfect sensor – i.e. one with no error – was utilized. This hypothetical scenario suggests how associations might appear as sensors improve and measurement error is reduced. Until a sensor with
Fig. 7. Posterior median of risk of concussion based on 200 draws of risk function parameters β : (B) multiply imputed 162 unobserved concussions under missing at random assumption without measurement error, (D) multiply imputed 52 unobserved concussions under missing at random assumption without measurement error, (F) multiply imputed 380 unobserved concussions under missing at random assumption without measurement error, (C) multiply imputed 162 unobserved concussions under missing at random assumption with measurement error, (E) multiply imputed 52 unobserved concussions under missing at random assumption with measurement error, and (G) multiply imputed 380 unobserved concussions under missing at random assumption with measurement error.
Fig. 6. Posterior mean and associated 95% credible intervals for risk of concussion, multiply imputing 162 unobserved concussions under missing at random assumption, without accounting for measurement error (black) and removing measurement error (red), (a) on log-odds scale, (b) on probability scale.
M.R. Elliott et al. / Journal of Biomechanics 48 (2015) 3059–3065
3065
Appendix A. Supplementary Information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jbiomech.2015.07.026
References
Fig. 8. Cumulative effect of three types of error on the 10% risk (solid bars) and 50% (grey bars) of concussion. The left most data points represent the point estimates from Rowson et al. (2012) and moving to the right, represent the progressive effect of adding sampling variability, under-reporting, and accounting for sensor error. The ranges shown represent the smallest lower end of a 95% confidence interval with the largest upper end of a 95% confidence interval across the three levels of under-reporting considered.
negligible error exists, we recommend that future reports should account for both sampling variability and concussion underreporting in their measured data to predict injury risk from observed rotational acceleration and should incorporate sensor error to define fundamental biomechanical relationships from such data. 4.1. Limitations There are a number of assumptions in this analysis. First, for concussion under-reporting, we have assumed (a) that misclassification is monotonic – that is, there are true concussions classed as non-concussions but no true non-concussions classed as concussions, (b) the misclassification rate is known (although allowed to vary over three rates as a sensitivity analysis), and (c) random, conditional on rotational acceleration. Second, we assume that rotational acceleration sensor measurement error is non-differential (again, random conditional on true rotational acceleration), multiplicative and log-normally distributed, with known variance. This assumption differs from a previous report by Funk et al. (2012), in which sensor error was presumed to occur predominantly at higher accelerations. Third, we utilized sensor error calculated on a hockey–helmet based system on this analysis of football data. While the magnitude of the sensor error for the football and hockey based systems is similar, the specific value of the error for football may be slightly different. However, for illustrative purposes, the effect of sensor error on the injury risk curves is adequately demonstrated by the hockey-based helmet system sensor error. Finally, because we are generating hypothetical posterior population predictive distributions using results of Rowson et al. (2012), we assumed a linear relationship between the log-odds of concussion and rotational acceleration.
Conflict of interest statement Kristy B. Arbogast serves as a consultant for the National Football League Players Association on topics related to head injury and player safety.
Allison, M.A., Kang, Y.S., Bolte 4th, J.H., Maltese, M.R., Arbogast, K.B., 2014. Validation of a helmet-based system to measure head impact biomechanics in ice hockey. Med. Sci. Sports Exerc. 46, 115–123. Allison, M.A., Kang, Y.S., Bolte 4th, J.H., Maltese, M.R., Arbogast, K.B., 2015. Measurement of Hybrid III impact kinematics using an accelerometer and gyroscope system in ice hockey helmets. Ann. Biomed. Eng. 43, 1896–1906. Bartsch, A., Samorezov, S., Benzel, E., Miele, V., Brett, D., 2014. Validation of an “intelligent mouthguard” single event head impact dosimeter. Stapp Car Crash J. 58, 1–28. Brainard, L.L., Beckwith, J.G., Chu, J.J., Crisco, J.J., McAllister, T.W., Duhaime, A.C., Maerlender, A.C., Greenwald, R.M., 2012. Gender differences in head impacts sustained by collegiate ice hockey players. Med. Sci. Sports Exerc. 44, 297–304. Carrol, R.J., Ruppert, D., Stefanski, L.A., 1995. Measurement Error in Nonlinear Models: A Modern Perspective. Chapman and Hall/CRC Press, Boca Raton, FL. Cook, J.R., Stefanski, L.A., 1994. Simulation-extrapolation estimation in parametric measurement error models. J. Am. Stat. Assoc. 89, 1314–1328. Crisco, J.J., Fiore, R., Beckwith, J.G., Chu, J.J., Brolinson, P.G., Duma, S., McAllister, T. W., Duhaime, A.C., Greenwald, R.M., 2010. Frequency and location of head impact exposures in individual collegiate football players. J. Athl. Train. 45, 549–559. Funk, J.R., Rowson, S., Daniel, R.W., Duma, S.M., 2012. Validation of concussion risk curves for collegiate football players derived from HITS data. Ann. Biomed. Eng. 40, 79–89. Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B., 2014. Bayesian Data Analysis, 3rd edition. CDC Press, Boca Raton, FL. Hernandez, F., Wu, L.C., Yip, M.C., Laksari, K., Hoffman, A.R., Lopez, J.R., Grant, G.A., Kleivem, S., Camarillo, D.B., 2014. Six degree-of-freedom measurements of human mild traumatic brain injury. Ann. Biomed. Eng. 43, 1918–1934. Hosmer, D.W., Lemeshow, S., Sturdivant, R.X., 2013. Applied Logistic Regression, 3rd Edition Wiley Press, Hoboken, NJ. Jadischke, R., Viano, D.C., Dau, N., King, A.I., McCarthy, J., 2013. On the accuracy of the head impact telemetry (HIT) system used in football helmets. J. Biomech. 46, 2310–2315. King, D., Hume, P.A., Brughelli, M., Gisanne, C., 2015. Instrumented mouthguard acceleration analyses for head impacts in amateur rugby union players over a season of matches. Am. J. Sports Med. 43, 614–624. Little, R.J.A., Rubin, D.B., 2002. Statistical Analysis with Missing Data, 2nd edition. Wiley, New York. McCrea, M., Hammeke, T., Olsen, G., Leo, P., Guskiewicz, K., 2004. Unreported concussion in high school football players: implications for prevention. Clin. J. Sport Med. 14, 13–17. Mihalik, J., Bell, D., Marshall, S., Guskiewicz, K., 2007. Measurement of head impacts in collegiate football players: an investigation of positional and event-type differences. Neurosurgery 6, 1229–1235. Pellman, E.J., Viano, D.C., Tucker, A.M., Casson, I.R., Waeckerle, J.F., 2003. Concussion in professional football: reconstruction of game impacts and injuries. Neurosurgery 53, 799–812. Rice, S.O., 1945. Mathematical analysis of random noise. Bell Syst. Tech. J. 24, 46–156. Rowson, S., Brolinson, G., Goforth, M., Dietter, D., Duma, S., 2009. Linear and angular head acceleration measurements in collegiate football. J. Biomech. Eng. 131, 061016. Rowson, S., Duma, S.M., 2011. Development of the STAR evaluation system for football helmets: integrating player head impact exposure and risk of concussion. Ann. Biomed. Eng. 39, 2130–2140. Rowson, S., Duma, S.M., 2013. Brain injury prediction: assessing the combined probability of concussion using linear and rotational head acceleration. Ann. Biomed. Eng. 41, 873–882. Rowson, S., Duma, S.M., Beckwith, J.G., Chu, J.J., Greenwald, R.M., Crisco, J.J., Brolinson, P.G., Duhaime, A.C., McAllister, T.W., Maerlender, A.C., 2012. Rotational head kinematics in football impacts: an injury risk function for concussion. Ann. Biomed. Eng. 40, 1–13. Stefanski, L.A., Cook, J.R., 1995. Simulation-extrapolation: the measurement error jackknife. J. Am. Stat. Assoc. 90, 1247–1256. Wilcox, B.J., Beckwith, J.G., Greenwald, R.M., Raukar, N.P., Chu, J.J., McAllister, T.W., Flashman, L.A., Maerlender, A.C., Duhaime, A.C., Crisco, J.J., 2015. Biomechanics of head impacts associated with diagnosed concussion in female collegiate ice hockey players. J. Biomech. 48, 2201–2204.