qikit

SPC charts and quality improvement tools for healthcare, in Python.

One function. One result object. One .plot() call.

Install

pip install clinical-qikit

Quickstart

from qikit import qic

result = qic(y=values, chart="i")
result.plot()   # returns a Plotly figure

qic() mirrors R's qicharts2, the package most healthcare QI practitioners already know — same chart types, same run-chart signal detection, same mental model, now in Python.

Chart types

All chart types are computed by the same qic() entry point via the chart= argument:

qic(..., funnel=True) renders a funnel plot instead, for cross-sectional comparison across p/pp/u/up/oe/oep charts sorted by denominator, with Spiegelhalter's exact 95% and 99.8% probability contours in place of 3σ. Each provider is plotted at its own volume (sample size, exposure or expected events), with the contours drawn as smooth curves across that range. Every method, reference and design choice behind the funnels is recorded in docs/methods.md, including where qikit follows Spiegelhalter (2005) and where it departs from him.

Risk-adjusted O/E funnel plots

Which chart for a cross-sectional comparison:

chart="oe" compares observed events against model-derived expected events — the form used for provider-level benchmarking of mortality, readmissions, and complications. Pass observed counts as y= and expected events as n=:

qic(x=physicians, y=observed_deaths, n=expected_deaths,
    chart="oe", funnel=True, show_95=True)

Limits are exact Poisson probability contours at 95% and 99.8% rather than 3σ (Spiegelhalter, Statistics in Medicine 2005;24:1185-1202), which matters at the volumes individual physicians actually have: the limits are asymmetric about the center line, and the lower limit stays above zero where a normal approximation would clip it and lose better-than-expected outliers. limit_method="byar" selects the closed-form Wilson-Hilferty approximation instead; it drifts from exact below ~20 expected events.

The center line defaults to the pooled ΣO/ΣE, so points are judged against their peers and any overall miscalibration of the risk model is absorbed. Pass cl=1.0 to judge them against the risk model itself.

Over-dispersion

Providers often differ from one another by more than chance allows, because no risk model captures everything. Across a few hundred hospitals, exact limits can flag far more than the nominal 0.2%. Every funnel, and every oe chart, estimates that excess spread the way Spiegelhalter published it (Stat Med 2005, §5). It uses Pearson residuals, winsorized at the 10th/90th percentiles so the outliers being screened for cannot inflate it, and a method-of-moments between-provider variance τ²:

r = qic(y=observed, n=expected, chart="oe", funnel=True)
r.summary["dispersion_tau"]    # spread of true O/E between providers; 0.0 when none is detected
r.summary["overdispersed"]     # True when τ̂ > 0
r.summary["dispersion_note"]   # one quotable sentence, present only when overdispersed

An unadjusted funnel (p, u, oe) reports the estimate but does not use it. It answers is this provider consistent with the target, given chance? The adjusted funnels (pp and up with funnel=True, and oep) answer is this provider unusual, given the spread between providers? They draw the additive random-effects model of his eq. 11: each provider's true rate varies around the target with SD τ̂. The limits are the exact Poisson–gamma or beta-binomial quantiles of that model rather than its normal approximation:

r = qic(y=observed, n=expected, chart="oep", funnel=True)
r.summary["dispersion_adjusted"]  # False when τ̂ = 0, and then the limits are oe's exactly
r.summary["dispersion_method"]    # "additive"

Under the additive model the smallest providers' limits barely move, because their own noise dominates. The busiest providers' limits level off at roughly θ₀ ± z·τ instead of closing onto the center line. At physician volumes τ̂ is usually 0, and the two charts are then identical. On the bundled 24-surgeon data it is 0.

Choose which question you are asking before looking at the results. Treat a non-zero τ̂ as a finding about the system: unmeasured case mix, data differences or real variation in practice. Spiegelhalter calls the adjustment a temporary measure, not a fix. Outside a funnel, pp/up remain Laney's p′/u′ charts, whose moving-range σ_z is designed for time order.

On a small cohort, one extreme provider drags the pooled center line toward itself and inflates τ̂ even after its own residual is clipped. Pass cl=1.0 to anchor the O/E center line if that matters. docs/methods.md covers this and the estimator's known downward bias.

Physician-level comparison

At individual-physician volumes a funnel is mostly a test with no power. The chance of flagging a physician whose true O/E is 2.0, at the 99.8% limit:

Most physicians sit at 5–20 expected events, so a genuinely poor performer is more likely to be missed than caught. That is a property of the sample size, not of the method — but it means absence of a signal is not evidence of acceptable performance, and that is exactly how a blank funnel tends to get read.

O/E charts therefore report what they could not have found. Every point carries its own detectability and interval, and the summary says whether the chart had any reach:

r = qic(x=physicians, y=observed_deaths, n=expected_deaths, chart="oe", funnel=True)

r.data[["observed", "expected", "y", "min_detectable_oe",
        "ci_95_lower", "ci_95_upper", "sigma_signal"]]

r.summary["underpowered"]               # True when the median point can't detect a doubling
r.summary["min_detectable_oe_median"]
r.summary["n_underpowered"]
r.summary["power_note"]                 # one quotable sentence, present only when underpowered

• min_detectable_oe — the smallest true O/E this point would flag 80% of the time. A physician at O/E 2.5 who does not signal, but whose min_detectable_oe is 4.3, has not been cleared; the chart simply cannot see that far at their volume.
• ci_95_lower / ci_95_upper — an exact (Garwood) interval on this point's own O/E. This is a different object from the funnel band: the band asks where the point would fall if the risk model were right and the provider average, while the interval asks what range of true performance is consistent with what the provider actually did. The band moves with the center line; the interval does not.
• An underpowered chart also captions itself, so the caveat survives being pasted into a slide deck. Supplying your own caption= replaces it.

Note too that at 99.8% about 0.2% of points flag by chance alone, so in a 500-physician funnel roughly one signal is expected even when every physician is average.

p and u funnels report the same way. A complication rate over a surgeon's 30 cases is no better powered than an O/E over their 5 expected deaths, so funnel=True on p/pp/u/up adds the same columns and summary keys:

r = qic(x=surgeons, y=complications, n=cases, chart="p", funnel=True)

r.data[["y", "ucl", "min_detectable", "ci_95_lower", "ci_95_upper", "sigma_signal"]]
r.summary["underpowered"]                # True when the median point can't detect a doubling
r.summary["min_detectable_ratio_median"] # median of min_detectable / cl
r.summary["n_underpowered"], r.summary["power_note"]

• min_detectable is a proportion (p, pp) or rate (u, up) on the chart's own y scale, so it scales with multiply= like ucl. A surgeon at 4 of 12 (33%) against an 11.6% center line does not signal, and their min_detectable of 0.58 says why.
• The interval is exact: Clopper–Pearson for proportions, Garwood for rates.
• "Underpowered" is relative to the center line, because a proportion has no natural unit the way O/E does: a point is underpowered when min_detectable / cl > 2. On a p funnel whose center line is 50% or more, a doubling is not a proportion at all, so only points that cannot flag at any rate count there.

The workflow that actually works is to aggregate until the chart has power. Using the bundled example (24 surgeons, 3 years):

import pandas as pd
from qikit import qic

df = pd.read_csv("data/examples/physician_mortality_oe.csv")

# Any single year: not one surgeon flags, and underpowered is True.
y2025 = df[df.year == 2025]
qic(x="physician", y="observed_deaths", n="expected_deaths", data=y2025,
    chart="oe", funnel=True).summary["underpowered"]        # True

# Three years pooled: the chart gains resolution and finds three real signals,
# including a high performer below the lower limit.
agg = df.groupby("physician", as_index=False)[["observed_deaths", "expected_deaths"]].sum()
qic(x="physician", y="observed_deaths", n="expected_deaths", data=agg,
    chart="oe", funnel=True).plot()

The volume axis

A funnel plots each point against its denominator, as Spiegelhalter draws them: sample size on p/pp, exposure on u/up, expected events on oe/oep. The limits are drawn as smooth curves on a grid across the volume range, and each curve passes exactly through every point's own limit, so a flagged point never appears inside the funnel.

r = qic(x="physician", y="observed_deaths", n="expected_deaths", data=agg,
        chart="oe", funnel=True, show_95=True,
        x_log=True)                  # log-scale the volume axis
r.plot()                             # x-axis titled "expected_deaths"; xlab= to rename
r.funnel_limits()                    # the drawn contours: n, cl, ucl, lcl, ucl_95, lcl_95
r.data[["x", "n", "y", "ucl"]]       # "n" is the denominator on every chart

• xlab titles whatever x-axis is drawn. Left unset, a volume axis takes the n= column name, or "Expected events", "Sample size" or "Exposure".
• x_spacing="x" plots by x value instead: the categorical funnel of releases before 0.6.0. x_order= orders categories, so on a funnel it needs x_spacing="x". A funnel whose denominators are all equal has no volume to spread, and falls back to "x" on its own.
• funnel_limits(n=...) evaluates the contours at any denominators, for a figure of your own. Use it rather than reading x and y off the traces: the trace names are stable, but not how many points each trace has.

Labeling flagged providers

When x is an ID (an NPI, say) and the reader needs names, keep the ID as the key and pass the display text separately:

r = qic(x="npi", y="observed_deaths", n="expected_deaths", data=agg,
        chart="oep", funnel=True, show_95=True,
        point_labels="physician_name",          # hover, tick text, signal labels
        hover_extra=["encounters", "practice"], # extra tooltip lines
        label_signals="both",                   # label 99.8% and 95% points
        x_spacing="rank")                       # even spacing in funnel order
r.data[["x", "label", "sigma_signal", "warn_signal"]]

• warn_signal marks points past the 95% contour that are not already sigma_signal, so the two tiers never overlap. It is present on every oe/oep chart and on p/pp/u/up funnels, counted in summary["n_warn"], and kept out of signals: at one chance alarm in twenty, a large funnel would always report one. With show_95=True these points are drawn amber.
• label_signals= takes "outlier", "warn", "both" or None. Labels above the center line point up and those below point down, and neighbours are staggered.
• x_spacing="rank" plots positions 1…k in funnel order with the labels as tick text, for when the names matter more than the volumes. Every trace and annotation uses the same positions, so numeric-looking IDs are never coerced into a stretched axis (the volume axis avoids that too, since IDs never reach it).
• Trace names ("data", "CL", "99.8%+", "95%-", …) are a stable contract for restyling, unless you rename a trace with styles={role: {"name": ...}} (below). The tooltip's customdata layout is not, so use the parameters above rather than editing it.

A last caution for anyone reading these in a credentialing context: at 5 expected deaths, a confidence interval that contains 1.0 means insufficient information, not exoneration. The honest report of a low-volume physician is the interval, not a verdict.

Flagging one direction, or by p-value and FDR

r = qic(x="physician", y="observed_deaths", n="expected_deaths", data=agg,
        chart="oe", funnel=True, cl=1.0,
        signal_side="above",     # only excess deaths flag
        signal_method="fdr")     # flag by Benjamini–Hochberg q-value
r.data[["x", "y", "p_value", "q_value", "sigma_signal"]]
r.summary["count_model"]         # "poisson": the model behind the limits and p-values

• signal_side= is "both" (default), "above" or "below". It filters sigma_signal and warn_signal, and with them the colors, signals and labels. The other side's limits are still drawn; hide them with styles (below). Pooled over 2023–2025, the bundled surgeons flag SURG-07 and SURG-15 above, and SURG-21 below. It is valid only on charts with funnel contours.

• p_value and q_value are on every funnel and O/E chart:

  • p_value is the exact p-value against the same Poisson, Poisson–gamma, binomial or beta-binomial model, and the same τ, that the limits use: P(X ≥ O) above, P(X ≤ O) below, doubled for two-sided.
  • q_value is its Benjamini–Hochberg adjustment within the chart (each facet separately).

• signal_method= chooses what decides the flags:

  • "limits" (default): the drawn contours.
  • "exact": p < 0.001 for red and p < 0.025 for amber, doubled for two-sided.
  • "fdr": the same levels, applied to q.

  At small counts "exact" is stricter than the line above the center: 13 deaths against 5 expected sits above the 99.8% line, but P(X ≥ 13) = 0.2%. docs/methods.md §4.11 explains why. When the curves do not decide the colors, an uncaptioned figure says what did.

Slide-quality figures

r.plot(
    show_95=True,
    x_range=(2, None),          # start at 2 expected deaths; the y axis refits
    styles={
        "limits_95": {"color": "#777", "width": 1.5, "dash": "solid"},
        "lcl": {"visible": False}, "lcl_95": {"visible": False},
        "ucl": {"name": "99.8% limit"}, "ucl_95": {"name": "95% limit"},
    },
    show_legend=True,
    width=1000, height=560,
)

• styles= restyles by role.
  • Roles: data, sigma, warn, runs, cl, ucl, lcl, ucl_95, lcl_95, target and end_labels, plus the groups limits and limits_95.
  • Attributes: color, width, dash, visible, name, and size on end_labels.
  • A specific role overrides its group.
  • name renames the trace, its legend entry and its end label.
  • The defaults draw exactly the figure you get without styles=.
• x_range= / y_range= take (low, high), and either end may be None.
  • x_range is in the units drawn: the denominator on a volume axis, positions under "rank", numbers or dates otherwise. It is rejected on a categorical axis.
  • With x_range set, the y axis fits only the points inside it, the center line, and each limit's narrowest in-range value. Trimmed points and the funnel's low-volume flare then stop stretching it.
• show_legend=True groups each limit pair into one entry and adds an entry per signal color: "Beyond 99.8%" and "Beyond 95%", or "p < …"/"q < …" under signal_method.
• show_labels=False hides the end-of-line labels ("CL=…", "99.8%+=…").
• width/height default to None, which keeps the figure responsive in notebooks and Streamlit. Plotly's static export then draws it at 700×500, so set both for slides.
• Unknown names raise. Every option here can be passed to qic() or overridden per call with r.plot(...). An unknown name raises TypeError rather than being ignored.

Funnel options at a glance

Beyond qic(), the package also provides:

• paretochart(x, data=...) — Pareto chart of categorical frequencies.
• bchart(...) — Bernoulli CUSUM chart for early detection of shifts in binary (pass/fail) data.
• design(factors, ...) / analyze(design_obj, response, ...) — full-factorial design of experiments (DOE) and effect analysis.

Excel add-in

A Microsoft Excel task-pane add-in built on the same charting engine (ported to TypeScript) lives in excel-addin/ — see its README for dev-harness and sideload instructions, and docs/excel-addin-roadmap.md for the roadmap.

Web app

app.py is a Streamlit front end over the same engine:

pip install -e .[app]
streamlit run app.py

Development

uv sync --extra dev
uv run ruff check .        # lint (CI-enforced)
uv run pytest tests/       # 176 tests incl. cross-language fixture conformance

The Python engine is authoritative; the TypeScript port in excel-addin/packages/engine conforms to it via the shared JSON fixtures in fixtures/. After intentional engine changes, regenerate the fixture snapshots with uv run python scripts/update_snapshots.py.

References

1. Montgomery DC. Introduction to Statistical Quality Control, 8th ed. Wiley, 2019.
2. Provost LP, Murray SK. The Health Care Data Guide, 2nd ed. Jossey-Bass, 2022.

License

MIT