Levey-Jennings chart with NIST MAVRO.DAT transmittance data

Python Code

"""Clinical Lab QC — Levey-Jennings chart with NIST MAVRO.DAT transmittance data."""


import numpy as np
import pandas as pd
import requests
import warnings
import plotly.graph_objects as go


MAVRO_URL = "https://www.itl.nist.gov/div898/handbook/datasets/MAVRO.DAT"

N_HEADER = 25   # header lines to skip
N_VALUES = 50   # data points to read


def fetch_mavro():
    """Fetch NIST MAVRO.DAT — 25 header lines, then 50 float values."""
    warnings.filterwarnings("ignore")
    r = requests.get(MAVRO_URL, verify=False, timeout=30)
    r.raise_for_status()
    lines = r.text.strip().split("\n")
    data_lines = lines[N_HEADER:]
    values = []
    for line in data_lines:
        for token in line.strip().split():
            try:
                values.append(float(token))
            except ValueError:
                pass
            if len(values) >= N_VALUES:
                break
        if len(values) >= N_VALUES:
            break
    return np.array(values[:N_VALUES])


def westgard_violations(values, mean, sd):
    """Return dict of Westgard rule violations."""
    n = len(values)
    violations = {}
    z = (values - mean) / sd

    # 1_3s: any point outside ±3SD
    rule_1_3s = np.where(np.abs(z) > 3)[0]
    if len(rule_1_3s):
        violations["1₃s"] = rule_1_3s.tolist()

    # 2_2s: two consecutive points on same side beyond ±2SD
    rule_2_2s = []
    for i in range(n - 1):
        if z[i] > 2 and z[i + 1] > 2:
            rule_2_2s.extend([i, i + 1])
        elif z[i] < -2 and z[i + 1] < -2:
            rule_2_2s.extend([i, i + 1])
    if rule_2_2s:
        violations["2₂s"] = list(set(rule_2_2s))

    return violations


def generate():
    print("fetching NIST MAVRO transmittance data …")
    values = None
    try:
        values = fetch_mavro()
        print(f"  fetched {len(values)} values: min={values.min():.4f} max={values.max():.4f}")
    except Exception as e:
        print(f"  fetch failed ({e}), using synthetic QC data")

    if values is None or len(values) < 10:
        # Synthetic QC data: mostly in-control with a few violations
        rng = np.random.default_rng(41)
        values = rng.normal(0, 1, N_VALUES)
        values[20] = 3.5    # 1_3s violation
        values[35] = 2.2    # start of 2_2s
        values[36] = 2.3
        mean = 0.0
        sd   = 1.0
    else:
        mean = float(np.mean(values))
        sd   = float(np.std(values, ddof=1))

    runs = np.arange(1, len(values) + 1)
    z    = (values - mean) / sd
    violations = westgard_violations(values, mean, sd)

    print(f"  mean={mean:.5f}, SD={sd:.5f}")
    print(f"  violations: {violations}")

    fig = go.Figure()

    # ±3SD zone (red)
    fig.add_hrect(y0=mean + 2 * sd, y1=mean + 3 * sd,
                  fillcolor="rgba(218, 85, 151, 0.12)", line_width=0, layer="below")
    fig.add_hrect(y0=mean - 3 * sd, y1=mean - 2 * sd,
                  fillcolor="rgba(218, 85, 151, 0.12)", line_width=0, layer="below")

    # ±2SD zone (yellow/orange)
    fig.add_hrect(y0=mean + sd, y1=mean + 2 * sd,
                  fillcolor="rgba(233, 162, 59, 0.13)", line_width=0, layer="below")
    fig.add_hrect(y0=mean - 2 * sd, y1=mean - sd,
                  fillcolor="rgba(233, 162, 59, 0.13)", line_width=0, layer="below")

    # ±1SD zone (green)
    fig.add_hrect(y0=mean - sd, y1=mean + sd,
                  fillcolor="rgba(85, 182, 133, 0.12)", line_width=0, layer="below")

    # Reference lines
    for k, color, label in [
        (3,  PINK,   "+3SD"), (2, ORANGE, "+2SD"), (1, GREEN, "+1SD"),
        (-1, GREEN,  "−1SD"), (-2, ORANGE, "−2SD"), (-3, PINK, "−3SD"),
    ]:
        fig.add_hline(
            y=mean + k * sd,
            line=dict(color=color, width=1.2, dash="dot"),
            annotation_text=label,
            annotation_position="right",
            annotation_font=dict(size=9, color=color),
        )

    # Mean line
    fig.add_hline(
        y=mean,
        line=dict(color=VIOLET, width=2),
        annotation_text="Mean",
        annotation_position="right",
        annotation_font=dict(size=10, color=VIOLET),
    )

    # Main data line
    fig.add_trace(go.Scatter(
        x=runs, y=values,
        mode="lines+markers",
        name="QC measurement",
        line=dict(color="#1c2024", width=1.5),
        marker=dict(color="#1c2024", size=6),
        hovertemplate="Run %{x}: %{y:.5f}<extra></extra>",
    ))

    # Flag violations
    all_violation_indices = set()
    for rule, idxs in violations.items():
        all_violation_indices.update(idxs)

    if all_violation_indices:
        viol_runs = runs[list(all_violation_indices)]
        viol_vals = values[list(all_violation_indices)]
        fig.add_trace(go.Scatter(
            x=viol_runs, y=viol_vals,
            mode="markers",
            name="Westgard violation",
            marker=dict(color=PINK, size=12, symbol="circle-open",
                        line=dict(width=2.5, color=PINK)),
            hovertemplate="Run %{x}: %{y:.5f} — VIOLATION<extra></extra>",
        ))

    # Violation annotation
    rule_strs = ", ".join(violations.keys()) if violations else "none"
    fig.add_annotation(
        x=0.02, y=0.97, xref="paper", yref="paper",
        text=f"Violations: {rule_strs}",
        showarrow=False,
        font=dict(size=11, color=PINK if violations else MUTED),
        bgcolor="rgba(255,255,255,0.85)",
        bordercolor=GRID,
        borderwidth=1,
        align="left",
    )

    fig.update_layout(
        xaxis=dict(title="Run", dtick=5),
        yaxis=dict(title="Transmittance"),
        legend=dict(orientation="h", y=1.08),
        margin=dict(t=50, b=50, l=70, r=80),
        height=460,
    )
    return fig


if __name__ == "__main__":
    generate()

Made with Plotly