Computing Diversion Z-Score Thresholds

Compute robust per-entity diversion baselines with median and MAD, flag events above a modified z-score, and add hysteresis to stop alert flapping

The naive way to flag a diverting prescriber is to compute a mean and standard deviation of their dispensing history and alert when a new event lands more than three standard deviations out. It looks statistically respectable and it fails in the field for one reason: the outliers you are hunting are the same outliers that inflate the standard deviation you measure them against. Two or three large removals widen the band until the next large removal scores an unremarkable z = 1.8, and the detector quietly stops detecting. This recipe implements the estimator that survives contact with real diversion data — a median and median absolute deviation (MAD) baseline scored with a modified z-score — and adds the hysteresis that stops a borderline entity from flapping in and out of alert on every event. It is the statistical engine behind Diversion Alert Thresholds & Suspicious-Order Monitoring, which defines how the score becomes a hold and a review.

The problem to solve precisely: given a stream of per-entity events (a dispense quantity, an order size, an override count), decide for each new event whether it deviates enough from that entity’s own history to warrant an alert — robustly, on small samples, without heavy dependencies, and without oscillating around the threshold. The output is a score, a band, and an alert state that a downstream monitor can act on and an auditor can reproduce.

Prerequisites & environment

  • Python 3.10+ — the code uses statistics from the standard library, dataclasses, and modern type syntax.
  • Standard library only for the core path: statistics, collections.deque, math, hashlib, logging. numpy is optional and only worthwhile at very high event volume; the pure-stdlib path is intentionally the default so a compliance node needs no scientific stack. Add pytest for the test block.
  • Regulatory grounding you should already hold: 21 CFR § 1301.74(b) requires detecting orders of unusual size, frequency, or pattern deviation — the “pattern deviation” clause is exactly what a z-score operationalizes. Retention and completeness of the resulting records fall under 21 CFR § 1304.21, and integrity of the audit trail under 45 CFR § 164.312(b).
  • Conceptual dependency: this recipe assumes the per-entity, per-schedule framing established in the parent topic. It computes the score; the parent decides the band-to-action mapping and the escalation path.

Why plain mean/standard deviation is the wrong tool

A quick numeric example makes the failure concrete. Suppose an entity’s true routine dispensing is around 5 units, and a diverter begins removing 60 units at a time.

  • Mean/std: after three 60-unit events mixed into a 5-unit history, the mean climbs and the standard deviation balloons. A fourth 60-unit event may score below 3σ and pass. The estimator has been poisoned by the very behavior it should catch.
  • Median/MAD: the median stays near 5 and the MAD stays small until more than half the window is corrupted. Each 60-unit event scores a large modified z-score every time, because the baseline refuses to move for a minority of anomalies.

The robustness difference is not cosmetic — it is the 50% breakdown point of the median versus the 0% breakdown point of the mean. In diversion monitoring, where the adversary controls some of the data, only a high-breakdown estimator is defensible.

Robust baseline band and a modified z-score threshold breach A horizontal axis of event magnitudes shows a grouping of routine events near a median marked M, surrounded by a shaded acceptance band whose half-width is a multiple of the median absolute deviation. Two later events sit far to the right beyond the band edge, each labelled with a modified z-score above the alert threshold, illustrating that a robust median and MAD band flags anomalies that a mean and standard deviation would absorb. event magnitude acceptance band: median ± k·MAD M threshold z′ = 6.2 z′ = 4.1 flagged: beyond k·MAD

Implementation: robust baseline + modified z-score with hysteresis

The function computes the modified z-score z' = 0.6745 · (x − median) / MAD. The 0.6745 constant makes MAD a consistent estimator of the standard deviation for normally distributed data, so a threshold like 3.5 carries roughly its familiar “3.5 sigma” meaning. Hysteresis is implemented with two thresholds: an entity enters the alert state at enter_z and only leaves it when a later event falls back below the lower exit_z. A single band would let an entity hovering at the boundary toggle alert/clear on every event, burying the real signal in noise.

python
from __future__ import annotations

import hashlib
import logging
import statistics
from collections import defaultdict, deque
from dataclasses import dataclass, field
from datetime import datetime, timezone
from typing import Deque, Optional

logger = logging.getLogger("cs.diversion.zscore")

_MAD_SCALE = 0.6745          # consistency constant: MAD -> sigma for normal data
_EPS = 1e-9                  # guards a zero MAD (all-identical window)


@dataclass(frozen=True)
class ZConfig:
    window: int = 60         # trailing observations retained per entity
    min_obs: int = 12        # below this the baseline is unstable -> WARMING
    enter_z: float = 3.5     # modified z at/above which alert engages
    exit_z: float = 2.5      # modified z below which alert releases (hysteresis)
    magnitude_floor: float = 0.0  # optional absolute backstop; 0 disables


@dataclass(frozen=True)
class Score:
    entity_id: str
    ndc: str
    quantity: float
    modified_z: Optional[float]   # None while WARMING
    state: str                    # WARMING | NORMAL | ALERT
    ts: str = field(default_factory=lambda: datetime.now(timezone.utc).isoformat())

    def audit_hash(self) -> str:
        raw = f"{self.entity_id}|{self.ndc}|{self.quantity}|{self.modified_z}|{self.state}|{self.ts}"
        return hashlib.sha256(raw.encode("utf-8")).hexdigest()


class RobustBaseline:
    """Per-(entity, ndc) median/MAD baseline with hysteretic alert state."""

    def __init__(self, cfg: ZConfig = ZConfig()) -> None:
        self.cfg = cfg
        self._hist: dict[tuple[str, str], Deque[float]] = defaultdict(
            lambda: deque(maxlen=cfg.window)
        )
        self._alerting: set[tuple[str, str]] = set()

    def _modified_z(self, key: tuple[str, str], x: float) -> Optional[float]:
        hist = self._hist[key]
        if len(hist) < self.cfg.min_obs:
            return None
        med = statistics.median(hist)
        mad = statistics.median([abs(v - med) for v in hist])
        return _MAD_SCALE * (x - med) / (mad if mad > _EPS else _EPS)

    def score(self, entity_id: str, ndc: str, quantity: float) -> Score:
        """Score one event, apply hysteresis, then update the baseline.

        Scoring BEFORE update means an event is never judged against itself,
        preserving a clean pattern-deviation test per 21 CFR § 1301.74(b).
        """
        key = (entity_id, ndc)
        z = self._modified_z(key, abs(quantity))

        if z is None:
            state = "WARMING"
        else:
            engaged = key in self._alerting
            floor_hit = (self.cfg.magnitude_floor > 0
                         and abs(quantity) >= self.cfg.magnitude_floor)
            if not engaged and (abs(z) >= self.cfg.enter_z or floor_hit):
                self._alerting.add(key)
                state = "ALERT"
            elif engaged and abs(z) < self.cfg.exit_z and not floor_hit:
                self._alerting.discard(key)
                state = "NORMAL"
            else:
                state = "ALERT" if engaged else "NORMAL"

        self._hist[key].append(abs(quantity))  # update AFTER scoring

        result = Score(entity_id, ndc, quantity,
                       None if z is None else round(z, 3), state)
        # Structured, PHI-free audit line: entity is a role/location key.
        logger.info(
            "zscore_evaluated",
            extra={"entity": entity_id, "ndc": ndc, "state": state,
                   "modified_z": result.modified_z,
                   "audit_hash": result.audit_hash()},
        )
        return result

The numpy optional path only changes _modified_z: numpy.median over a window is marginally faster once windows exceed a few thousand points, but the semantics are identical, so keeping stdlib as the default avoids shipping a scientific stack to every compliance node. The scoring-before-update ordering is the one invariant you must not “optimize” away — folding the current event into the baseline first is the classic bug that lets a large event hide inside the statistics meant to catch it.

Verification & testing

The tests prove three properties: MAD resists poisoning where std does not, hysteresis actually suppresses flapping, and the audit hash binds the score to its inputs.

python
import pytest


def _feed(bl, entity, ndc, values):
    return [bl.score(entity, ndc, v) for v in values]


def test_mad_flags_what_mean_std_would_absorb():
    bl = RobustBaseline(ZConfig(min_obs=8, window=40))
    # Routine ~5-unit history, then repeated 60-unit "diversion" events.
    _feed(bl, "op-7", "00093-0058-01", [5, 4, 6, 5, 5, 4, 6, 5])
    hits = _feed(bl, "op-7", "00093-0058-01", [60, 60, 60, 60])
    # Every large event stays flagged; a mean/std band would fade after a few.
    assert all(s.state == "ALERT" for s in hits)
    assert all(abs(s.modified_z) > 3.5 for s in hits)


def test_hysteresis_prevents_flapping():
    bl = RobustBaseline(ZConfig(min_obs=6, enter_z=3.5, exit_z=2.5, window=30))
    _feed(bl, "presc-1", "11111-1111-11", [10, 11, 9, 10, 10, 11])
    # A borderline event engages the alert...
    (s_enter,) = _feed(bl, "presc-1", "11111-1111-11", [40])
    assert s_enter.state == "ALERT"
    # ...a mid-band event does NOT immediately clear it (must fall below exit_z).
    mid = bl.score("presc-1", "11111-1111-11", 20)
    assert mid.state == "ALERT"          # still engaged: hysteresis holds
    low = bl.score("presc-1", "11111-1111-11", 10)
    assert low.state == "NORMAL"         # now safely below exit band


def test_warming_does_not_score_small_samples():
    bl = RobustBaseline(ZConfig(min_obs=12))
    scores = _feed(bl, "loc-3", "22222-2222-22", [5, 200, 5])
    assert all(s.state == "WARMING" and s.modified_z is None for s in scores)


def test_audit_hash_is_bound_to_inputs():
    bl = RobustBaseline(ZConfig(min_obs=4))
    _feed(bl, "op-9", "33333-3333-33", [3, 3, 3, 3])
    (s,) = _feed(bl, "op-9", "33333-3333-33", [90])
    assert s.audit_hash() == Score(
        s.entity_id, s.ndc, s.quantity, s.modified_z, s.state, s.ts
    ).audit_hash()

A representative emitted log line — note it carries only operational identifiers, never patient data:

text
INFO cs.diversion.zscore zscore_evaluated entity=op-7 ndc=00093-0058-01 \
     state=ALERT modified_z=8.094 audit_hash=b41d...7ce2

Re-deriving the hash from the persisted entity_id, ndc, quantity, modified_z, state, and ts and comparing it to the stored value is the tamper check: any post-hoc edit to a score breaks the digest, which is the integrity signal 45 CFR § 164.312(b) expects.

Gotchas & compliance pitfalls

  • Small-sample instability. With fewer than ~12 observations the MAD is a coarse estimate and a single unusual-but-legitimate event can score enormously. Enforce a min_obs warm-up during which events accumulate the baseline but are not scored; scoring against three points invents both false alerts and false clears.
  • Seasonality and dosing cycles. A chemo clinic’s controlled-substance use spikes predictably on infusion days; a flat rolling baseline will flag every cycle. Either bucket the baseline by day-of-week/shift or widen the window to span several full cycles so the periodic peak is inside the normal band rather than compared against the trough.
  • Gaming the baseline. A patient diverter can slowly inflate their own routine so a later bulk removal scores as normal. Median/MAD blunts this because more than half the window must move first, but pair it with an absolute magnitude_floor so no amount of baseline drift can lower the hard ceiling on a single Schedule II removal.
  • A zero MAD is real. If every event in the window is identical (common for unit-dose items), MAD is exactly zero and the score would divide by zero. Floor the denominator with a small epsilon so an identical-history entity does not silently throw, and treat the first genuine deviation as a strong signal.
  • Do not delete cleared alerts. Hysteresis will move an entity from ALERT back to NORMAL; that transition is a record, not an erasure. Append the state change to the ledger so the history shows the system detected, tracked, and released — the evidence that the § 1301.74(b) program is live.
  • Confirmed anomalies flow onward. A score is not a finding. When review upgrades an alert to confirmed loss, the numbers feed Automating Form 106 from Ledger Snapshots, which renders the filing from an immutable snapshot so the reported figure cannot drift after discovery.

Frequently Asked Questions

What z-score threshold should I start with?

A modified z-score of 3.5 to enter and 2.5 to exit is a defensible starting band for Schedule II, because the 0.6745 scaling makes the modified z comparable to standard deviations for well-behaved data. Tighten enter_z for the most tightly controlled drugs and loosen it for high-volume, low-risk items, but always keep exit_z below enter_z so hysteresis has room to work.

Do I need numpy?

No. The statistics module in the standard library computes the median and MAD correctly, and for the window sizes a per-entity diversion baseline uses — tens to low hundreds of points — pure Python is fast enough. Reach for numpy only if you are scoring millions of events per batch, and even then the semantics are unchanged, so the stdlib path remains the correctness reference.

Why score before updating the baseline?

Because folding the current event into the median and MAD first lets a large diversion event partially normalize its own score, weakening exactly the detection you want. Scoring the event against the prior history and only then appending it keeps the pattern-deviation test honest, which is what 21 CFR § 1301.74(b) is asking for.

How is this different from a simple high-water limit?

A high-water limit (“alert above 100 units”) is scale-fixed: it is deaf to a low-volume drug being drained in small amounts and noisy for a high-volume clinic. A per-entity modified z-score is scale-free — it measures deviation relative to that entity’s own established pattern — so it catches the relative anomaly a fixed limit misses. In practice you run both: the z-score for pattern deviation and a magnitude floor as an absolute backstop.