1. Executive Summary

This study develops and validates a predictive model for workers' compensation injury frequency using 692,412 establishment-year observations from OSHA's Injury Tracking Application (ITA) spanning 2016-2024. The model uses year Y features (prior injury rates, enforcement history, industry classification, employer size) to predict year Y+1 recordable injury counts.

The model architecture combines a Negative Binomial GLM with exposure offset (providing the interpretable base rate), Buhlmann-Straub credibility weighting (blending individual experience with industry priors), and a LightGBM model trained on GLM residuals (capturing nonlinear interactions, improving validation MAE by 19.1%). The GLM layer is fully transparent and standard in actuarial practice. The gradient boosting layer improves discrimination but introduces model governance considerations and explainability requirements discussed in Section 6. Model deployment is contingent on calendar-year adjustment and governance sign-off.

2. Data Description

2.1 Primary Data Source

The OSHA Injury Tracking Application (ITA) collects annual injury and illness data from establishments meeting reporting thresholds (generally 250+ employees, or 20+ employees in designated high-hazard industries). Each record represents one establishment-year and includes:

• Exposure: Total hours worked, annual average employee count
• Outcomes: Total recordable cases, days-away-from-work cases, job transfer/restriction cases, fatalities
• Derived rates: TRIR (Total Recordable Incident Rate = cases × 200,000 / hours), DART rate
• Classification: 6-digit NAICS code, EIN, establishment identifier

2.2 Supplementary Data Sources

• OSHA Enforcement: 5.1M inspections and 13.2M violations with severity classifications and penalty amounts
• OSHA Severe Injury Reports (SIR): 103K events including hospitalizations, amputations, fatalities

2.3 Panel Construction

The model uses a lagged panel design: for each establishment with consecutive-year ITA records, we pair year Y features with year Y+1 outcomes. This produces 692,412 observation pairs across 267,091 unique establishments.

Exhibit 1: Panel size by feature year (left) and observed TRIR trends (right). The COVID-19 dip in 2020-2021 is visible. 2017 has lower volume due to ITA reporting changes.

2.4 Temporal Split (No Data Leakage)

2.5 Data Quality Filters

• Minimum 2,000 annual hours worked (~1 FTE) — below this threshold, rates are unreliable
• Maximum 10M hours — 942 records with implausible values (up to 16 trillion) were capped
• Requires valid NAICS code (4+ digits) and non-null establishment linkage
• TRIR capped at 99th percentile (27.86) to limit outlier influence on the GLM

3. Methodology

3.1 Model Architecture

The model uses a three-stage architecture standard in actuarial pricing:

1. Stage 1 — Negative Binomial GLM: Provides the interpretable base prediction with industry fixed effects and an exposure offset. This is the actuarially transparent component.
2. Stage 2 — Buhlmann-Straub Credibility: Blends individual entity predictions with NAICS-4 group rates based on exposure volume, properly handling the small-entity problem.
3. Stage 3 — LightGBM Residual Model: A gradient boosting model trained on GLM Pearson residuals to capture nonlinear interactions and feature interactions not explained by the base GLM, improving MAE by 19.1% on validation.

3.2 Target Variable and Exposure Treatment

Total recordable cases (integer count) for year Y+1, with log(total hours worked) as the exposure offset. This is the standard actuarial approach: modeling claim frequency per unit of exposure. The Negative Binomial distribution accommodates overdispersion (injuries tend to cluster more than a Poisson process would predict).

3.3 Feature Set

4. Stage 1: Negative Binomial GLM

4.1 Specification

Where log(hours) is the offset (coefficient fixed at 1.0), ensuring the model predicts a rate per unit of exposure. The NAICS-2 sector dummies capture the industry base rate. The model was fit using IRLS (Iteratively Reweighted Least Squares) and converged in ~110 seconds on 346,191 training observations.

4.2 Coefficient Estimates

Table 1: GLM coefficients. Positive coefficients increase predicted injury counts. All features except serious_plus_count and sir_count are statistically significant. These two are captured via the GBM stage.

4.3 Key Coefficient Interpretations

• prior_trir (+0.098): A 1-point increase in prior TRIR predicts a ~10.2% increase in next-year injury count (exp(0.0975) = 1.102). This is the dominant predictor.
• log_prior_employees (-0.050): Holding hours constant, larger employers have slightly lower per-hour injury rates — consistent with the well-known large-employer safety advantage.
• prior_severity_ratio (+0.045): When a higher fraction of injuries are lost-time cases (vs. medical-only), the entity is predicted to have more injuries the following year.
• log_penalty (+0.018): Higher OSHA penalty amounts, controlling for industry and violations, predict more future injuries — penalties are a lagging indicator of hazardous conditions.

4.4 GLM Performance

5. Stage 2: Buhlmann-Straub Credibility

5.1 Motivation

An entity's observed injury rate is a noisy estimate of its true underlying risk. For a large employer with 500,000+ annual hours, the observed rate is highly credible. For a small employer with 10,000 hours, year-to-year variation can swamp the signal. Buhlmann-Straub credibility provides the actuarially standard solution: blend individual experience with a group prior, weighted by exposure volume.

5.2 Formula

5.3 Variance Component Estimation

Using multi-year entities (establishments with 2+ consecutive years of data), we estimated:

Exhibit 5: Left — distribution of credibility weights across the test set. Right — the credibility function showing Z vs. hours worked. At 50,000 hours (~25 FTE) the entity gets 50% weight on its own experience. At 200,000 hours (~100 FTE), Z ≈ 0.80.

5.4 Credibility Validation

To verify the credibility mechanism works correctly, we compared individual, group, and blended MAE within Z-bands on the test set:

Table 3: The blended prediction equals or beats both individual and group predictions in every Z-band, confirming the credibility mechanism is correctly calibrated.

6. Stage 3: LightGBM Residual Model

6.1 Approach

After the GLM produces a base prediction, we compute Pearson residuals: (observed - predicted) / sqrt(predicted). A LightGBM gradient boosting model is trained on these residuals using the full feature set (including nonlinear features the GLM cannot capture). The final prediction is: GLM_prediction + GBM_residual × sqrt(GLM_prediction).

6.2 Hyperparameter Tuning

Optuna Bayesian optimization with 5 trials selected: learning rate = 0.020, 87 leaves, max depth = 5. Early stopping at 30 rounds of no improvement.

6.3 Top Features by Gain

6.4 Stage 2 Improvement

The GBM stage provides meaningful improvement by capturing nonlinear interactions between features that the linear GLM cannot model (e.g., industry-specific penalty effects, size-dependent trends).

6.5 Governance Considerations

The GLM base layer is fully transparent: each coefficient has a clear directional interpretation, and the model can be expressed as a simple formula. The GBM residual layer improves discrimination but reduces pure interpretability — individual predictions cannot be decomposed into additive factor contributions as cleanly as the GLM alone.

For model governance purposes:

• The GLM layer can be deployed independently (with ~19% higher MAE) if full transparency is required
• The GBM layer should be monitored for feature drift (distributional shifts in input features over time), prediction stability (large changes in output for small input changes), and segment-level performance (degradation in specific NAICS or size segments)
• Periodic retraining (recommended annually as new ITA data is released) should include validation against the prior model version to detect performance changes
• SHAP values are available for individual predictions to provide post-hoc explainability of the GBM contribution, but these are approximations, not exact decompositions

7. Model Performance

7.1 Out-of-Time Test Set Results

All results below are on the held-out test set (feature year 2022, predicting 2023 outcomes, N = 116,957 establishments).

7.2 Calibration

Exhibit 8: Left — calibration plot showing predicted vs observed TRIR by decile. The model is well-calibrated with monotonically increasing actual TRIR across deciles. Right — Lorenz curve showing the model's ability to separate low-risk from high-risk establishments. The area between the model curve and the diagonal represents discriminatory power.

7.3 Lift Analysis

8. Size Bias Analysis

8.1 Is This Just Measuring Entity Size?

A critical question for any injury prediction model: is the predictive power simply "bigger entities have more injuries"? The evidence suggests the primary signal is persistent establishment-level risk differences (as reflected in prior observed experience), not headcount.

Exhibit 6: Size predicts case counts (trivially, ρ = 0.54) but barely predicts rates (ρ = 0.12). The real predictive signal is prior rate → future rate (ρ = 0.48).

8.2 Discrimination Within Size Bands

To confirm the model works beyond size, we tested predictive power within size bands — comparing same-size entities against each other:

Exhibit 4: Left — Spearman ρ between prior TRIR and next-year TRIR within each size band. All are statistically significant (p < 10^-64). Right — Lift (top vs bottom quartile of prior TRIR) within each size band. Even micro entities show 1.7x lift.

9. Small Entity Considerations

Exhibit 7: Left — Small entities (<25 FTE) with prior injuries have a next-year TRIR of 7.14 vs 2.82 for those without (2.5x signal). Right — distribution of non-zero injury rates for small entities.

9.1 The Zero-Inflation Problem

For small entities (<25 FTE, 19.8% of the test set):

• 59.8% have zero injuries in any given target year
• Of those with zero prior-year injuries, 71.5% remain at zero the following year
• Median target TRIR = 0.00 (more than half have no injuries)

9.2 What Works for Small Entities

• Binary signal is strong: Whether an entity had ANY prior injuries is a powerful discriminator (next-year TRIR of 7.14 vs 2.82, p < 0.001)
• NAICS classification is critical: When individual data is sparse, the industry base rate carries almost all the weight — NAICS-4 selection matters enormously
• Enforcement data helps: Even without ITA history, violation counts and penalties from OSHA inspections provide a regulatory history signal correlated with future injury experience

9.3 Recommended Approach by Entity Size

10. Comparison to Current Peer Grading

The current system ranks establishments by percentile within their NAICS-4 cohort using a weighted composite of 7 factors (violation severity 0.26, TRIR 0.20, SIR events 0.16, penalties 0.08, financial health 0.10, litigation 0.09, news sentiment 0.11). This is a point-in-time relative ranking, not a predictive model.

Exhibit 10: The new predictive model captures nearly twice the injury volume in the top risk decile.

11. Industry Rate Tables

11.1 Major Sector (NAICS-2) Base Rates

Exhibit 2: Exposure-weighted TRIR by major NAICS sector. Healthcare, Warehousing, and Public Administration have the highest observed rates. Construction's rate appears low because the construction ITA reporters tend to be larger, better-organized firms. Note on Mining (NAICS 21): The low TRIR shown (0.74) reflects only OSHA-regulated mining operations, primarily oil & gas extraction and surface mining support services. Underground coal mines, metal/nonmetal mines, and most quarries are regulated by MSHA (Mine Safety and Health Administration), not OSHA, and report injuries through a separate system. MSHA-regulated operations — which account for the majority of mining fatalities and severe injuries — are excluded from this study. The MSHA injury data (msha_accident) is available in a parallel dataset and could be integrated in a future extension of this model.

11.2 Minor Industry (NAICS-4) — Top 25 by Volume

11.3 COVID-19 Impact

Exhibit 9: Exposure-weighted TRIR by year. The 2020-2021 dip reflects reduced workplace exposure during the pandemic. The model includes a COVID binary indicator (coefficient = -0.088) to account for this structural break.

12. Limitations & Caveats

1. ITA Selection Bias (Critical): The OSHA ITA survey targets establishments with 250+ employees or in designated high-hazard industries. The model is trained on a non-representative sample biased toward larger, riskier employers. Predictions for employers outside this profile should be treated as class-rate priors with limited individual calibration (Z → 0).
2. Frequency Only, Not Severity: This model predicts recordable case counts, not dollar losses. Workers' compensation pricing requires both frequency and severity. Severity modeling would require carrier loss data not available in OSHA datasets.
3. Aggregate Over-Prediction: The model over-predicts total cases by ~25% on the test set (695,877 predicted vs 556,251 actual). This reflects the secular improvement in workplace safety — the 2016-2020 training data has higher baseline rates than 2023 outcomes. A calendar-year trend factor should be applied in production.
4. NAICS-6 Sparsity: Many 6-digit NAICS codes have fewer than 10 ITA reporters. The model uses NAICS-4 (415 codes) and NAICS-2 (50 sectors) to ensure adequate cohort sizes, but cannot differentiate within narrow sub-industries.
5. State-Plan Variation: Approximately half of US states operate their own OSHA programs. ITA reporting completeness may vary by jurisdiction, and state-plan inspection/violation data may have different patterns than federal OSHA. The model includes state as a feature but cannot fully control for reporting differences.
6. MSHA Exclusion: Mining operations regulated by the Mine Safety and Health Administration (MSHA) — including underground coal, metal/nonmetal mines, and most quarries — do not report to OSHA's ITA system. The "Mining" sector rates in this study reflect only OSHA-regulated operations (primarily oil & gas extraction and surface mining support). MSHA-regulated mining has substantially higher injury and fatality rates that are not captured here.
7. Not a Premium Rate: This model produces an expected injury frequency index, not an insurance premium. Converting to premium requires loss development factors, trend factors, expense loading, and profit provisions that are outside the scope of this study.
8. Stationarity Assumption: The model assumes that the relationship between year Y features and year Y+1 outcomes is approximately stationary. Structural changes (new OSHA regulations, industry shifts, economic cycles) may require periodic retraining.

13. Deployment & Governance Recommendations

1. Deploy as primary frequency predictor: The model outperforms the current peer grading system on every discrimination metric and should be used for forward-looking injury frequency estimation in underwriting.
2. Apply calendar-year trend factor: Fit a simple multiplicative trend to the TRIR time series to adjust for the secular improvement in workplace safety (~2-3% per year).
3. Retrain annually: As new ITA data becomes available, retrain the model to keep it current. The holdout set (2023→2024) is available for the first retrain validation.
4. Use credibility Z as a data-quality flag: Expose the credibility weight to underwriters so they know how much trust to place in individual experience vs. class rate.
5. Consider severity modeling: If carrier loss data becomes available, a parallel severity model (average cost per claim by injury type) would complete the actuarial pricing picture.
6. Extend to non-ITA entities: For the ~3.6M establishments without ITA records, the NAICS-4 group rate from this model can serve as a class-rated prior. Enforcement and SIR data (which exists for many non-ITA entities) can further adjust via the GBM stage.

Exhibit 3: Left — predicted risk decile (based on prior TRIR) vs actual next-year TRIR. The monotonic increase confirms that the model's risk ordering translates directly into observed outcomes. Right — percentage of all injuries captured by each decile. The top decile (D10) captures a disproportionate share of injuries while the bottom decile (D1) captures very few, demonstrating strong discrimination. Individual injury outcomes are inherently noisy (rare events follow a Negative Binomial process), but aggregate decile-level patterns are highly stable and predictable.

14. What This Model Can and Cannot Be Used For

15. Recommended Next Steps Before Production Use

1. Calendar-year trend adjustment: Fit a multiplicative trend factor to eliminate the ~25% aggregate over-prediction. Re-estimate this factor each year as new ITA data becomes available.
2. Parallel run: Deploy alongside the current peer grading system for 6-12 months. Compare predictions against realized outcomes for the parallel period before retiring peer grading for frequency estimation on ITA-like accounts.
3. Exposure projection protocol: Define how submitted or projected hours/payroll will be used at quote time in place of the realized hours used in backtesting. Validate that the rate prediction remains calibrated when paired with exposure estimates.
4. Segment-level validation: Evaluate discrimination and calibration separately by NAICS sector, state, and size band. Identify segments where the model underperforms and consider segment-specific adjustments or exclusions.
5. Holdout validation: Run the model on the reserved 2023→2024 holdout set to confirm out-of-time stability before production deployment.
6. Model governance documentation: Prepare model risk management documentation including ongoing monitoring plan, retraining schedule, performance thresholds that would trigger model review, and escalation procedures.
7. Credibility and applicability flags: Expose the credibility weight Z and an ITA-applicability flag in the model output so underwriters can see how much of the prediction is individually experience-rated vs. class-rated, and whether the account falls within the model's training population.