Introduction

A deterministic risk assessment gives you one number: THQ = 0.83. Below 1, so "safe" — case closed. But that number was built from a mean concentration, a mean ingestion rate, and a mean body weight. Real consumers don't eat the mean fish with the mean cadmium content at the mean body weight. Some eat more fish, some fish are more contaminated, some consumers weigh less. The honest question is not "is the mean consumer safe?" but "what fraction of the exposed population exceeds the threshold?"

That is exactly what Monte Carlo probabilistic risk assessment answers — and it is built into Labreador's Risk Assessment module. This post explains the method, the distribution choices, and how to read the output without a statistics degree.

The problem with point estimates

The deterministic chain is familiar:

EDI = (C × IR) / BW
THQ = EDI / RfD

Every input on the right side is uncertain:

C (concentration) — varies between individual fish, sampling sites, and seasons. A mean of 0.3 mg/kg can hide individuals at 1.2 mg/kg.
IR (ingestion rate) — consumption surveys show heavy right tails: most people eat little fish, a few eat a lot.
BW (body weight) — varies across the population, and children are a distinct scenario entirely.

A point estimate uses one value for each and produces one THQ. It tells you nothing about the spread — and regulators increasingly ask for exactly that spread.

How Monte Carlo works

The idea is simple: instead of one value per input, define a probability distribution per input. Then repeat the calculation thousands of times, each time drawing a random value from each distribution:

for each iteration i = 1..N:
    draw C_i  from the concentration distribution
    draw IR_i from the ingestion-rate distribution
    draw BW_i from the body-weight distribution
    THQ_i = (C_i × IR_i) / (BW_i × RfD)
collect all THQ_i → distribution of outcomes

With N = 10,000 iterations you get a full picture: the median outcome, the 95th percentile, and — most importantly — the probability of exceeding the threshold:

P(THQ > 1) = (number of iterations with THQ > 1) / N

Which distributions, and why

The distribution choices follow the US EPA (2001) Risk Assessment Guidance for Superfund, Volume III — Process for Conducting Probabilistic Risk Assessment:

Input	Distribution	Reason
Concentration (C)	Lognormal	Concentrations are strictly positive and right-skewed; a few highly contaminated samples drive the tail
Ingestion rate (IR)	Lognormal	Consumption is positive and right-skewed; high-end consumers matter most
Body weight (BW)	Normal, truncated	Roughly symmetric, but truncated at a lower bound so impossible values (negative or near-zero weights) never occur

The spread of each distribution is controlled by a coefficient of variation (CV = SD / mean). A CV of 0.30 for concentration means the standard deviation is 30% of the mean — a reasonable starting point for within-site variability when you don't have the raw replicate data. If you do have replicates, compute the CV from them.

One detail worth knowing: the body-weight truncation bound scales with the scenario mean. A fixed floor would silently turn a child scenario into an adult one and underestimate the dose — the bound adapts instead.

A worked example

Take cadmium in fish: mean C = 0.3 mg/kg (CV 0.30), IR = 0.0546 kg/day (CV 0.30), BW = 70 kg (CV 0.15), RfD = 0.001 mg/kg/day (US EPA).

Deterministic result:

EDI = (0.3 × 0.0546) / 70 = 2.34e-4 mg/kg/day
THQ = 2.34e-4 / 0.001 = 0.234

Comfortably below 1. But run 10,000 Monte Carlo iterations and the picture gains a dimension: the distribution of THQ is right-skewed, and the 95th percentile sits well above the mean. The key output is:

P(THQ > 1) ≈ 0.4 %

Interpretation: under these variability assumptions, roughly 4 in 1,000 random consumer-fish combinations exceed the reference dose. Whether that is acceptable is a risk-management decision — but now it is an informed one, not hidden inside a single average.

Reading the output

A Monte Carlo risk report typically gives you, per metal and for the hazard index (HI):

Mean and median — if they diverge strongly, the distribution is heavily skewed and the mean overstates the "typical" case.
P5–P95 interval — the 90% credible range of outcomes.
P(THQ > 1) / P(HI > 1) — the headline probability of exceeding the non-carcinogenic threshold.
P(CR > 1e-6) and P(CR > 1e-4) — for carcinogens, the probability of exceeding the de minimis (10⁻⁶) and upper-bound (10⁻⁴) cancer risk levels.

The HI histogram is usually plotted on a logarithmic x-axis, because the distribution spans orders of magnitude.

Reproducibility matters

A simulation that gives a different answer every time you run it is hard to defend in a paper or a report. Labreador's Monte Carlo engine uses a seeded pseudo-random generator: with the same seed and the same inputs, you get bit-identical results. Cite the seed, the iteration count, and the CV values in your methods section, and anyone can reproduce your numbers exactly.

Common mistakes to avoid

Using Monte Carlo to dress up bad data. Simulation propagates uncertainty — it does not reduce it. If your mean concentration is based on three samples, the CV assumption dominates the result.
Reporting only the mean of the distribution. That throws away the point of the exercise. Report the percentiles and the exceedance probability.
Ignoring non-detects. If half your samples are below the LOD, the concentration distribution is censored, and a naive lognormal fit is wrong. Resolve non-detects explicitly first (see our post on non-detects and LOD/2 substitution).
Mixing bases. Concentrations on a dry-weight basis must be converted to wet weight before they meet a wet-weight ingestion rate — see dry weight vs wet weight.

Try it yourself

The Risk Assessment module runs the full Monte Carlo analysis in your browser — per-metal THQ distributions, the HI histogram, exceedance probabilities, and carcinogenic risk percentiles — with a fixed seed for reproducibility. Your data never leaves your device.

References

US EPA (2001). Risk Assessment Guidance for Superfund: Volume III — Part A, Process for Conducting Probabilistic Risk Assessment. EPA 540-R-02-002.
US EPA IRIS — Cadmium oral RfD.