All data were stored on a secure drive accessible only to the study team. Statistical analyses were performed in Rstudio version 4.0.1 (Rstudio Team 2022 ). Using the NEISS-Work dataset, national ED-treated occupational injury count estimates were produced using the R packages “survey” and “srvyr” (Ellis et al. 2021 ; Lumley 2021 ) using the aforementioned NEISS-Work survey weights. ED-treated occupational injury count estimates were generated for all injuries and by injury event type, a categorical variable denoting the way an injury was incurred and is based on the aforementioned BLS OIICS v 2.01 classification system (National Institute for Occupational Safety and Health (NIOSH) Division of Safety Research 2021b ); all analyses were conducted both for total injury rate estimates and stratified by injury event type. ED-treated occupational injury rates were calculated per 10,000 FTE using Current Population Survey (CPS) estimates which were generated using NIOSH’s Employed Labor Force (ELF) query system; as NEISS-Work includes all work-related ED-treated injuries, FTE estimates were generated for all jobs (as opposed to “primary” or “secondary” jobs only) (National Institute for Occupational Safety and Health (NIOSH) Division of Safety Research 2021c ). Standard errors (SE) for FTE estimates were generated using generalized variance functions provided by BLS; standard errors were used to calculate monthly FTE variances by multiplying the square of the SE by corresponding ELF-generated monthly FTE estimates (i.e., the corresponding monthly sample size) (National Institute for Occupational Safety and Health (NIOSH) Division of Safety Research 2021c ). Variances of both numerator (injury count estimates) and denominator (FTE) data were used to calculate 95% confidence intervals (CI) for ED-treated occupational injury rate estimates based on Taylor series expansion (National Institute for Occupational Safety and Health (NIOSH) Division of Safety Research 2021d ) and were reported as injury rate estimates ± margin of error.
Seasonality of injury rate estimates was assessed by calculating seasonality indices per month. Seasonality indices were calculated by dividing the mean rate for each month by the mean monthly occupational injury rate for the entire dataset; seasonality indices of greater and less than one indicate higher than and lower than expected injury rates for a given month, respectively (Zhang et al. 2014 ).
To assess linear trends in injury rates over time, we fit a linear regression model to monthly injury rate estimates and adjusted for autocorrelation and serially correlated error terms using autoregressive integrated moving average (ARIMA) modeling. This analysis was conducted using both monthly total injury rate estimates and monthly estimates stratified by injury event type. In data violating the linear regression assumption of no autocorrelation, ARIMA models are used to control for serial correlation (e.g., seasonality) by including lagged dependent variable values and errors, including in studies of injury data (Box et al. 2016 ; Zhu et al. 2015 (link)). An ARIMA model takes the form ARIMA(p,d,q)(P,D,Q)m, where p is the order of autocorrelation, d is the number of differences applied to the data, q is the order of moving average terms, P, D, and Q are the seasonal versions of these terms, and m is the order of seasonality (e.g., 12 for annually seasonality in monthly data) (Hyndman and Athanasopoulos 2018a ). ARIMA models were fit to monthly injury rates by examining autocorrelation and partial autocorrelation plots. A lagged regression estimate was included if it showed statistical significance (p < 0.05) and was necessary to control for serial correlation. Finally, significance of each model’s Ljung-Box Q statistic was observed to ensure proper model fit, with a non-significant value considered a properly fit model (Ljung and Box 1978 ). The conditional sum of squares method was used to estimate all models. To assess temporal trends, a trend regressor with slope of one was included in each ARIMA model as a covariate and reported with 95% CIs (Hyndman and Athanasopoulos 2018b ). A total percent decrease in injury rates throughout the study period was estimated by multiplying this term by 96 (i.e., the total number of months in the study period) and calculating the percent difference from the model’s intercept; an analogous calculation using each trend parameter’s 95% CI was performed to determine each percent decrease’s 95% CI.
Seasonality of injury rate estimates was assessed by calculating seasonality indices per month. Seasonality indices were calculated by dividing the mean rate for each month by the mean monthly occupational injury rate for the entire dataset; seasonality indices of greater and less than one indicate higher than and lower than expected injury rates for a given month, respectively (Zhang et al. 2014 ).
To assess linear trends in injury rates over time, we fit a linear regression model to monthly injury rate estimates and adjusted for autocorrelation and serially correlated error terms using autoregressive integrated moving average (ARIMA) modeling. This analysis was conducted using both monthly total injury rate estimates and monthly estimates stratified by injury event type. In data violating the linear regression assumption of no autocorrelation, ARIMA models are used to control for serial correlation (e.g., seasonality) by including lagged dependent variable values and errors, including in studies of injury data (Box et al. 2016 ; Zhu et al. 2015 (link)). An ARIMA model takes the form ARIMA(p,d,q)(P,D,Q)m, where p is the order of autocorrelation, d is the number of differences applied to the data, q is the order of moving average terms, P, D, and Q are the seasonal versions of these terms, and m is the order of seasonality (e.g., 12 for annually seasonality in monthly data) (Hyndman and Athanasopoulos 2018a ). ARIMA models were fit to monthly injury rates by examining autocorrelation and partial autocorrelation plots. A lagged regression estimate was included if it showed statistical significance (p < 0.05) and was necessary to control for serial correlation. Finally, significance of each model’s Ljung-Box Q statistic was observed to ensure proper model fit, with a non-significant value considered a properly fit model (Ljung and Box 1978 ). The conditional sum of squares method was used to estimate all models. To assess temporal trends, a trend regressor with slope of one was included in each ARIMA model as a covariate and reported with 95% CIs (Hyndman and Athanasopoulos 2018b ). A total percent decrease in injury rates throughout the study period was estimated by multiplying this term by 96 (i.e., the total number of months in the study period) and calculating the percent difference from the model’s intercept; an analogous calculation using each trend parameter’s 95% CI was performed to determine each percent decrease’s 95% CI.
Free full text: Click here
