4. Feature engineering for crash injury prediction

Raw crash records rarely arrive in a form a random forest or XGBoost model can consume directly. A typical police report or event data recorder (EDR) export contains free-text narratives, categorical codes, timestamps, and a handful of continuous measurements whose physical meaning depends on the collision configuration. Feature engineering is the step that converts these heterogeneous fields into a compact, physically interpretable, and model-ready matrix. In injury prediction, this step matters more than model selection: a well-engineered feature set with a modest gradient-boosted tree will usually outperform a poorly engineered set with a heavily tuned one, because the signal that separates a minor injury from a serious one is concentrated in a small number of kinematic and occupant-related quantities.

This chapter works through three feature families — vehicle kinematics, occupant characteristics, and environmental risk factors — and shows how each is derived, validated, and encoded. Throughout, we assume a tabular dataset with one row per occupant (not per vehicle), since injury outcomes are occupant-level labels. This choice has a direct consequence: vehicle-level features must be broadcast to every occupant in that vehicle, and occupant-level features must never leak across rows.

Vehicle kinematics: from Δv and PDOF to occupant-level exposure

The single most predictive kinematic quantity in crash injury literature is the change in velocity, Δv, experienced by the occupant compartment. If the dataset provides pre- and post-impact speed vectors, Δv is the vector difference; if it provides only a scalar collision severity code, you must reconstruct an estimate. A practical derivation uses the principle of conservation of momentum for a two-vehicle collision:

# Example: reconstructing Δv for vehicle 1 in a two-vehicle collision
# Inputs (example values, illustrative only)
m1, m2 = 1500.0, 1800.0          # kg, vehicle masses
v1_pre = 50.0 / 3.6              # m/s, vehicle 1 pre-impact speed (50 km/h)
v2_pre = 0.0                     # m/s, vehicle 2 stationary
restitution = 0.1                # coefficient of restitution, typical for car-to-car

# Common post-impact velocity (perfectly inelastic baseline, then corrected)
v_common = (m1 * v1_pre + m2 * v2_pre) / (m1 + m2)
dv1 = (1 + restitution) * (v_common - v1_pre)   # negative: deceleration
dv1_kmh = abs(dv1) * 3.6
print(f"Estimated Δv for vehicle 1: {dv1_kmh:.1f} km/h")

With the example inputs above, the estimate is approximately 27.5 km/h. Treat this as a teaching value, not a measured result. The point is the workflow: mass ratio, closing speed, and restitution jointly determine Δv, and the model should receive Δv as a continuous feature rather than a coarse severity bucket. Where EDR data is available, prefer the recorded longitudinal Δv and use the reconstruction only as a fallback, flagging reconstructed rows with a boolean dv_is_estimated feature so the tree can learn to discount them.

Δv alone ignores direction. The principal direction of force (PDOF) determines which body regions are loaded: a frontal 0° impact loads the head and chest through the restraint system, while a 60° near-side impact loads the thorax and pelvis through the door. Encode PDOF as a continuous angle in degrees plus a small set of derived indicators — is_frontal (|PDOF| ≤ 30°), is_side (60°–120°), is_rear (≥ 150°) — and, critically, an interaction term dv_x_pdof or separate dv_frontal and dv_side columns. Trees can split on interactions, but supplying the interaction explicitly reduces the depth needed and improves stability on small datasets.

Two further kinematic features are worth deriving when the data supports them. First, delta_v_per_mass = Δv / vehicle mass, a proxy for the momentum transfer per unit structure. Second, time_to_collision or pre-impact braking duration, which captures whether the occupant was braced. If only a binary braking flag exists, keep it as a categorical feature and let the model decide its value.

Occupant characteristics: restraint use, seating, and vulnerability

Occupant features are where most leakage errors occur. Restraint status (belted, airbag_deployed) must be recorded at the occupant level, not inferred from the vehicle. A common mistake is to fill missing belt status with the vehicle's majority value; this silently imputes the outcome-correlated variable and inflates validation scores. Instead, encode missingness explicitly with a belt_status_missing indicator and impute with a neutral category such as unknown.

Age and sex are standard, but their effect is non-linear and interacts with crash direction. Rather than feeding raw age, derive age_band (e.g., 16–24, 25–44, 45–64, 65+) and a is_pediatric / is_elderly pair. Seating position should be one-hot encoded (seat_driver, seat_front_passenger, seat_rear_left, …) because the same Δv produces very different injury risk for a driver versus a rear-seat passenger. Height and weight, when present, yield BMI, which is a recognized modifier of restraint fit; keep BMI continuous and add a bmi_obese flag at the 30 kg/m² threshold.

Environmental risk factors: road, light, and time

Environmental features describe the conditions under which the kinematics occurred. Road type (highway, urban arterial, rural), surface condition (dry, wet, snow/ice), speed limit, and lighting (daylight, dark-lit, dark-unlit) are the core set. Speed limit deserves special treatment: it is a strong prior for impact speed and should be kept continuous, with a derived speed_over_limit = estimated impact speed − posted limit when both are available. Lighting and surface condition interact with Δv — a wet rural road at night raises the probability of a higher-energy impact — so consider a small number of explicit interaction flags such as dark_and_wet.

Time-of-day and day-of-week are cyclical and should be encoded with sine/cosine pairs rather than raw integers, otherwise the model treats 23:00 and 00:00 as maximally distant. Weather codes from police reports are notoriously inconsistent; map them to a controlled vocabulary of five to seven categories and record the mapping in a versioned lookup table so the transformation is reproducible.

Assembling and validating the feature matrix

Once the three families are derived, join them on crash_id and occupant_id, then verify the matrix before training. A short validation script catches the most common defects: duplicate occupant rows, features with zero variance, and columns whose correlation with the label exceeds 0.95 (a strong sign of leakage). The example below is a reproducible check you can adapt:

import pandas as pd

df = pd.read_parquet("features_occupant_level.parquet")

# 1. Uniqueness
assert not df.duplicated(subset=["crash_id", "occupant_id"]).any()

# 2. Zero-variance columns
nunique = df.nunique()
print("Constant columns:", nunique[nunique <= 1].index.tolist())

# 3. Suspiciously high label correlation (numeric features only)
num = df.select_dtypes("number")
corr = num.corrwith(df["mais_3plus"]).abs().sort_values(ascending=False)
print(corr.head(10))

Expected output is a clean uniqueness assertion, an empty constant-column list, and a correlation ranking in which no single feature exceeds roughly 0.6. If a feature such as hospital_admitted appears at the top, it is a post-outcome variable and must be removed — it is not available at prediction time. This is the most frequent cause of a model that scores 0.95 AUC in cross-validation and collapses in deployment.

Finally, document the feature set as a versioned schema: column name, dtype, derivation rule, source field, and missing-value policy. Random forests and XGBoost are forgiving of scale but not of silent schema drift; a feature that changes meaning between data releases will quietly degrade the model. Keeping the derivation rules in code, with the schema checked at load time, is what makes the feature engineering reproducible rather than a one-off notebook exercise.