Structural equation modelling for safety risk

Drag the sample size

Every estimate below was produced by fitting data simulated from known path weights, 300 replications per cell. Drag the sample size. The wide band is where estimates actually land; the thin bar is the interval the estimator reports. When the thin bar is shorter, the reported confidence is false. API reference →

Sample size n = 100

Where estimates land (mean ± 1.96 empirical SD)Interval the estimator reports (± 1.96 analytic SE)True weight

All four paths at this n select a row

-

True weight: -
Mean estimate: -
Empirical SD: -
Mean analytic SE: -
SE ratio: -
95% coverage: -

Tick marks the nominal 95%.

The selected path, at every sample size tested

Source artifact: results/study01_recovery.csv

Why this exists

A formula that fits perfectly, recovers its own coefficients, and passes every conventional cutoff can still point the wrong way down its own arrows.

A particular formula circulates in practitioner writing on safety analytics: risk as a weighted sum of unsafe acts, operational stress, system condition and safety response capability, with weights around 0.45, 0.30, -0.25 and -0.20, described as a structural equation model with machine-constructed latent variables.

Read on

This repository builds the estimator that would be needed to produce that formula, from scratch in numpy, and then tests what such an estimator can support. Most of the answer is negative, and the negative results are the contribution.

The clearest one is here. Fitting the same data with the causal arrow reversed returns an identical coefficient (0.493 both ways), an identical chi-square, an identical CFI and an identical SRMR. Both models meet every conventional cutoff. Goodness of fit cannot tell you which way causation runs, because it never could.

Everything here is simulated

Read on

Every number on this page comes from data generated by a model whose parameters were fixed in advance. That is what a simulation study is for: ground truth is known, so estimator behaviour can be measured against it rather than assumed.

No number on this page is a claim about any real workplace, any real programme, or any real injury. Nothing in this repository touches real injury data at any point.

What this does and does not show

What this shows

Point estimates are nearly unbiased at every sample size tested, from n=100 to n=5,000.
Interval coverage sits at 76% to 86% against a nominal 95%, and does not improve with n.
Telling the two negative coefficients apart at 80% power needs a sample of 6,884.
Fitting the same data with the causal arrow reversed returns an identical coefficient (0.493 both ways) and identical fit statistics.

What this does not show

Simulation establishes estimator properties, not empirical facts. These studies say nothing about whether any particular safety programme works, and they are not evidence about any real site.
The generating models are the author's, chosen to resemble the structure the circulating formula implies. A different generating model would give different numbers, though the identification results do not depend on the parameterisation.
Replication counts are modest in places and each table records its own. The bootstrap study runs 40 replications, which is enough to show the direction of the coverage gap and not enough to pin its size.
Nothing here has been peer reviewed.

The practical conclusion is narrow and worth stating plainly: illustrative path weights of the kind that circulate in practitioner writing cannot be read causally, cannot be transferred between sites, and cannot be validated by goodness of fit.

Recovery, and the intervals around it

Start with the friendliest possible test.

Why

Generate data from a known model, then estimate it. If the estimator cannot recover coefficients it generated itself, nothing downstream is worth reading.

Coverage sits at 76% to 86% against a nominal 95%, at every sample size.

Source artifact: results/study01_recovery.csv

AUC holds near 0.69 at every base rate; alerts per true event climb to 385 at the rarest.

Source artifact: results/study02_discrimination.csv

Interval coverage against sample size

Coverage of the nominal 95% confidence interval, by sample size, for each of the four path coefficients.

Why

Coverage should sit on the dashed line. It sits between 76% and 86% instead, and it does not improve with sample size: the analytic standard errors are 63% to 76% of the empirical ones at every n tested. The point estimates are nearly unbiased. It is the uncertainty around them that is wrong, which is the more dangerous failure, because an interval that is too narrow reads as confidence.

Source artifact: results/study01_recovery.csv

n	predictor	true beta	mean estimate	bias	empirical sd	coverage 95	SE ratio
100	UnsafeActs	0.45	0.4562	0.006205	0.1519	0.78	0.6311
100	OperationalStress	0.3	0.3074	0.007413	0.1431	0.7933	0.6696
100	SystemCondition	-0.25	-0.2571	-0.007051	0.1466	0.7967	0.6579
100	SafetyResponseCapability	-0.2	-0.2021	-0.002123	0.1424	0.8233	0.6746
200	UnsafeActs	0.45	0.4537	0.00375	0.09747	0.8233	0.6981
200	OperationalStress	0.3	0.3021	0.00206	0.09338	0.85	0.7279
200	SystemCondition	-0.25	-0.2569	-0.006931	0.101	0.7933	0.6741
200	SafetyResponseCapability	-0.2	-0.1951	0.004912	0.09568	0.8467	0.7118
500	UnsafeActs	0.45	0.4522	0.002231	0.0565	0.8633	0.7572
500	OperationalStress	0.3	0.3015	0.001504	0.06038	0.8333	0.7081
500	SystemCondition	-0.25	-0.2532	-0.003214	0.06663	0.7567	0.6424
500	SafetyResponseCapability	-0.2	-0.2018	-0.001839	0.0583	0.8567	0.7331
1000	UnsafeActs	0.45	0.4509	0.0008792	0.04103	0.86	0.7357
1000	OperationalStress	0.3	0.3007	0.0006662	0.04383	0.8067	0.6885
1000	SystemCondition	-0.25	-0.2484	0.001632	0.04073	0.8433	0.7421
1000	SafetyResponseCapability	-0.2	-0.2016	-0.00155	0.04136	0.84	0.7297
2000	UnsafeActs	0.45	0.4527	0.002746	0.02892	0.8433	0.7366
2000	OperationalStress	0.3	0.2984	-0.001575	0.0305	0.8233	0.6988
2000	SystemCondition	-0.25	-0.2503	-0.0002661	0.03195	0.8	0.6667
2000	SafetyResponseCapability	-0.2	-0.2013	-0.001287	0.03029	0.8333	0.703
5000	UnsafeActs	0.45	0.4514	0.001357	0.01862	0.86	0.7241
5000	OperationalStress	0.3	0.3002	0.000187	0.01942	0.8367	0.6939
5000	SystemCondition	-0.25	-0.2514	-0.00143	0.01937	0.83	0.6964
5000	SafetyResponseCapability	-0.2	-0.2	-3.099e-05	0.01955	0.82	0.6894

Bias, spread and interval coverage for each coefficient at each sample size, over 300 replications per cell.

Why

Bias is negligible everywhere. Coverage never reaches its nominal level anywhere.

Source artifact: results/study01_recovery.csv

How much data the formula would need

Estimating a coefficient and being able to tell it apart from the coefficient next to it are different tasks with very different sample size requirements.

requirement	basis	n required
global fit: 80% power for the test of close fit	MacCallum-style RMSEA power, df = 80	153
path precision: SE = 0.05 on a single coefficient	standardized regression SE, R^2 = 0.288, VIF = 1.28	369
path precision: SE = 0.025 on a single coefficient	standardized regression SE, R^2 = 0.288, VIF = 1.28	1459
distinguish beta = 0.45 from beta = 0.30 at 80% power	z test on the coefficient difference, correlated predictors	770
distinguish beta = -0.25 from beta = -0.20 at 80% power	z test on the coefficient difference, correlated predictors	6884

Sample size required for each analytic goal, computed from the same generating model.

Why

Distinguishing the two negative coefficients from one another needs a sample of 6,884, which is larger than most sites will ever have. Reporting them as separate drivers at any realistic n is reporting noise.

Source artifact: results/study01_requirements.csv

coefficient a	beta a	coefficient b	beta b	absolute difference	n for 80pct power	se of difference at n 250
UnsafeActs	0.45	OperationalStress	0.3	0.15	770	0.09456
UnsafeActs	0.45	SystemCondition	-0.25	0.7	41	0.09456
UnsafeActs	0.45	SafetyResponseCapability	-0.2	0.65	46	0.09456
OperationalStress	0.3	SystemCondition	-0.25	0.55	62	0.09456
OperationalStress	0.3	SafetyResponseCapability	-0.2	0.5	74	0.09456
SystemCondition	-0.25	SafetyResponseCapability	-0.2	0.05	6884	0.09456

Every pairwise comparison between coefficients, with the sample size that would be needed to separate the pair at 80% power.

Why

Ranking the drivers is a much harder task than estimating them.

Source artifact: results/study01_distinguishability.csv

What rare events do

Safety incidents are rare.

Why

Rare outcomes break naive estimation in a specific way: discrimination survives, calibration does not, and calibration is what an operational alert depends on.

training sample	train n	train event fraction	auc test	mean predicted test	observed rate test	overprediction ratio	calibration slope	ppv at top 5pct
none (as collected)	100,000	0.00506	0.7111	0.005072	0.00516	0.983	1.215	0.0192
random undersampling of majority	1012	0.5	0.7075	0.454	0.00516	87.98	1.1	0.0188
random oversampling of minority	198,988	0.5	0.7109	0.4521	0.00516	87.62	1.187	0.0198
SMOTE-style interpolation	198,988	0.5	0.7108	0.4485	0.00516	86.92	1.102	0.0198

What resampling does to a rare-outcome model.

Why

Discrimination is essentially unchanged by rebalancing, while the predicted probabilities become meaningless: undersampling the majority class overpredicts by 88x if the intercept is not corrected.

Source artifact: results/study02_imbalance.csv

shift hours	base rate per shift	worker shifts per event	sensitivity	specificity	ppv	alerts per true event
8	9.2e-05	10,870	0.8	0.95	0.00147	680.3
8	9.2e-05	10,870	0.8	0.99	0.007307	136.9
8	9.2e-05	10,870	0.9	0.999	0.07647	13.08
8	9.2e-05	10,870	0.95	0.9999	0.4664	2.144
12	0.000138	7246	0.8	0.95	0.002203	453.8
12	0.000138	7246	0.8	0.99	0.01092	91.57
12	0.000138	7246	0.9	0.999	0.1105	9.05
12	0.000138	7246	0.95	0.9999	0.5673	1.763

Alerts per true event at realistic base rates.

Why

This is the operational meaning of a rare outcome: at a base rate this low, even a highly specific model buries the true events under false ones, and the people receiving the alerts stop reading them.

Source artifact: results/study02_alert_burden.csv

unit of analysis	base rate	relative precision target	events required	units required
worker-shift (8 h)	9.2e-05	0.1	384.1	4,175,499
worker-shift (8 h)	9.2e-05	0.25	61.46	668,080
worker-shift (12 h)	0.000138	0.1	384.1	2,783,666
worker-shift (12 h)	0.000138	0.25	61.46	445,387
worker-year	0.023	0.1	384.1	16,702
worker-year	0.023	0.25	61.46	2672

Events and units of observation required to estimate a rare-event rate to a given relative precision.Source artifact: results/study02_sizing.csv

Zero inflation

Most establishments record no injuries in a period.

Why

That is not a small deviation from a count model, it is a different data-generating process.

structural zero fraction	observed zero fraction	poisson expected zero fraction	excess zeros	mean count	variance count	variance to mean ratio	ols slope on latent 1
0	0.2367	0.2357	0.000997	1.625	2.098	1.291	0.5693
0.2	0.3884	0.2357	0.1527	1.31	2.109	1.61	0.5325
0.38	0.5204	0.2357	0.2847	1.042	2.03	1.948	0.4772

Zero-inflated outcomes are the norm in injury counts.

Why

As the structural zero fraction rises, the observed zeros pull away from what a Poisson would produce and the variance-to-mean ratio climbs, while the slope on the latent variable attenuates.

Source artifact: results/study02_zero_inflation.csv

Misspecification and equivalence

This is the section that matters.

Why

A model can be wrong in the specific way that inflates the number you care about, and still fit better than the correct one.

fitted model	beta predictor on outcome	true direct effect	bias	chi square	df	cfi	rmsea	meets conventional cutoffs
confounder included	0.173	0.2	-0.02705	28.46	24	0.9997	0.00682	True
confounder omitted	0.4955	0.2	0.2955	5.702	8	1	0	True

The same data, fitted with and without a confounder.

Why

Omitting it inflates the coefficient from 0.173 to 0.495 against a true direct effect of 0.20. Both models meet every conventional fit cutoff, and the misspecified one fits better.

Source artifact: results/study03_confounder.csv

fitted model	beta	se	p value	chi square	df	cfi	rmsea	srmr	meets conventional cutoffs
Cause -> Effect (correct direction)	0.4927	0.01589	0	7.696	8	1	0	0.007859	True
Effect -> Cause (reversed)	0.4927	0.01589	0	7.696	8	1	0	0.007859	True

Two models with opposite causal directions, fitted to the same data.

Why

Every fit statistic is identical to the last decimal place. No amount of model comparison separates them; only an assumption from outside the data can.

Source artifact: results/study03_equivalent.csv

generating structure	fitted structure	estimated beta	se	p value	chi square	df	cfi	rmsea	srmr	meets conventional cutoffs
Injuries -> Climate, beta = -0.55	Climate -> Injuries	-0.5698	0.01501	0	8.442	8	0.9999	0.004294	0.006726	True

Data generated with injuries driving climate, then fitted as climate driving injuries.

Why

The reversed model recovers a large, highly significant coefficient and passes every cutoff.

Source artifact: results/study03_reverse.csv

What a coefficient means

Even a correctly estimated coefficient does not tell you what happens if you intervene, and does not survive a change of scale or a change of site.

world	fitted coefficients	effect of intervening on UnsafeActs	note
1. four direct causes, correlated exogenously	[0.448, 0.306, -0.254, -0.2]	0.45	the only world in which the coefficient is an intervention effect
2. SystemCondition also acts through UnsafeActs	[0.453, 0.307, -0.26, -0.208]	0.45	total effect of SystemCondition is -0.070, not -0.250; the reported coefficient is the direct effect only
3. an unmeasured common cause of UnsafeActs and Risk	[0.448, 0.303, -0.246, -0.203]	0.281	generating direct effects are [ 0.281, 0.252, -0.298, -0.248]; the population regression on the measured predictors returns the headline vector exactly

Four causal worlds that produce nearly identical fitted coefficients, and what intervening on a variable would actually do in each.

Why

The coefficient is the same. What it licenses you to do is not.

Source artifact: results/study04_ambiguity.csv

site	crew	score raw inputs	score z inputs	ranking on z inputs
Site 1 (overtime varies widely)	A	2.55	0.05625	crew A scores higher
Site 1 (overtime varies widely)	B	3.15	-0.5812	crew A scores higher
Site 2 (overtime tightly controlled)	A	2.55	3.45	crew B scores higher
Site 2 (overtime tightly controlled)	B	3.15	6.45	crew B scores higher

Two crews carrying identical raw scores at two sites.

Why

Standardising the inputs flips the ranking at the site where overtime varies widely and leaves it intact where overtime is tightly controlled. Nothing about the crews changed between the rows; only the scaling did. A weighted-sum score is not scale-free, so the choice of scaling is a decision about who gets flagged.

Source artifact: results/study04_raw_scale.csv

quantity	value
Site A loadings	[0.806, 0.764, 0.769]
Site B loadings	[0.774, 0.793, 0.301]
Site A composite reliability	0.823
Site B composite reliability	0.675
correlation of Site B scores under its own vs Site A's measurement model	0.9463
share of Site B units whose top-decile flag flips when Site A's measurement model is imposed	0.046
correlation of Site A score with the generating latent	0.9103
correlation of Site B score with the generating latent	0.8759

Measurement loadings at two sites.

Why

When an indicator loads differently between sites, the latent variable does not mean the same thing at both, and comparing their scores compares two different constructs.

Source artifact: results/study04_invariance.csv

Run it yourself

ehs-risk-sem

git clone https://github.com/priyatham9/ehs-risk-sem
cd ehs-risk-sem
python3 -m unittest discover -s tests -v
python3 simulations/run_all.py --quick       # writes results_quick/
python3 simulations/run_all.py               # full run, writes results/

--quick writes to results_quick/, so a reduced-replication run cannot silently replace the checked-in tables.

What SEM cannot tell you about safety risk

Drag the sample size

All four paths at this n select a row

-

The selected path, at every sample size tested

Why this exists

Everything here is simulated

What this does and does not show

+What this shows

−What this does not show

Recovery, and the intervals around it

Interval coverage against sample size

How much data the formula would need

What rare events do

Zero inflation

Misspecification and equivalence

What a coefficient means

Run it yourself

What this shows

What this does not show