qikit

SPC charts and funnel plots for healthcare quality, in Python.

Two questions, two functions, one vocabulary:

• qic(): is one measure changing over time? Control charts and run charts, with qicharts2's rules or NHS England's Making Data Count rules.
• funnel(): which providers differ from their peers? Spiegelhalter funnel plots with exact limits, over-dispersion, p- and q-values, and power.

Both return a result whose .data is a table ready to write back to a database, whose .summary says what was found, and whose .plot() draws a Plotly figure.

Install

pip install clinical-qikit            # tables only: numpy and pandas
pip install "clinical-qikit[plot]"    # adds plotly for .plot()

The core needs only numpy and pandas, so a Snowpark stored procedure or UDTF can compute limits without plotly. import qikit never loads plotly; .plot() raises an ImportError naming the extra when it is missing.

Quickstart

from qikit import qic, funnel

r = qic(df, x="month", y="deaths", n="admissions", chart="p")
r.data          # one row per month: cl, ucl, lcl, signal, runs_signal
r.summary       # which rules fired, and how the limits were set
r.plot()        # a Plotly figure

f = funnel(df, y="deaths", n="expected", measure="oe", unit="surgeon", worse="higher")
print(f)        # what was compared, what flagged, what to know
f.data[f.data.signal != "none"]

One vocabulary

The two functions answer different questions, so they stay separate, but every shared idea has one name.

Arguments removed or renamed in 1.0 raise a TypeError naming the replacement, so an old call fails loudly rather than drawing something different. See CHANGELOG.md.

Over time: qic()

Rows that share an x are aggregated into one subgroup. Numeric and date x are sorted. x_period="month" buckets dates.

Rules (method=):

• "anhoej" (default): qicharts2's run rules. A run on one side of the centre line longer than round(log₂ n) + 3, or fewer crossings than the 5th percentile of Binomial(n − 1, ½).
• "nhs": NHS England's Making Data Count, for XmR boards. A point beyond the limits, 7 on one side, 7 rising or falling, or 2 of 3 near a limit. With worse=, each point is classed improving, concerning or common cause; with target=, the summary gives the assurance verdict (pass, fail or variable).
• "ihi": Provost and Murray. "weco", "nelson": the classic zone rules.

Baseline and phases. freeze=k sets the limits from points 1..k and judges later points against them. part= splits the chart into phases, each with its own centre line, limits and rules. exclude= drops points from the limits and the rules.

qic() reproduces qicharts2 to 1e-6 on its own data. Where it departs, and why, is in docs/qic-methods.md.

Across providers: funnel()

Each provider is plotted at its own volume, against Spiegelhalter's exact 99.8% and 95% probability contours (Statistics in Medicine 2005;24:1185-1202). Every method, reference and design choice is recorded in docs/methods.md, including where qikit departs from the paper.

Choosing the settings

Rows that share a unit are summed, so stay-level or year-level data aggregates itself.

Reading the result

r.data has one row per provider:

r.summary holds the model and how to read the counts:

• n_outlier, n_alert, and expected_false_flags.
• tau: the estimated spread of true values between providers; tau_method, tau_scale, tau_source ("estimated" or "pre-specified"), tau_estimated, and tau_ci under tau_method="ml".
• phi: the over-dispersion ratio as the chosen method defines it.
• count_model, limit_method, within_dispersion, dispersion_shape, underpowered, n_filtered, n_fragile.
• With benchmark=: benchmark, a dict of method, share and value. Under "abc" it also holds:
  • providers and group_denominator: who sets it and their volumes.
  • coverage: their share of patients.
  • msd and n_below_msd: Kiefe's minimum sufficient denominator.
  • without_small: the ABC without providers below it.
• notes: sentences ready to show a reader.

On a faceted funnel the same keys are under by_facet.

What decides a flag

Each contour is the boundary of an exact test. A provider is beyond the 99.8% line exactly when its one-sided p-value is below 0.001, and beyond the 95% line exactly when it is below 0.025. For the unadjusted O/E chart that is the same as "the provider's own exact 95% interval excludes the center line". The line, the p-value, the colour and the interval always agree.

signal= picks one test for both tiers:

• "limits" (default): outlier beyond the outer contour, alert beyond the inner one. About 0.1% of average providers land beyond each 99.8% line by chance, and 2.5% beyond each 95% line. summary["expected_false_flags"] counts that for your chart. With 1,000 providers on a one-sided chart it is about 1 outlier and 24 alerts, so "35 alerts" is read against that.
• "fdr": outlier when q_value < 0.05, alert when q_value < 0.10 (fdr_levels=). At most that share of each tier is expected to be false, however many providers there are. The flags then no longer match the drawn lines, and the caption says so.

There is deliberately no "limit OR q" mode. That union only ever adds flags, so it never removes the chance alerts that make FDR worth using, and neither guarantee holds for the combined tier (methods §4.11).

limits=(0.998, 0.95) moves the contours, for example to (0.997, 0.95).

Over-dispersion: real differences between providers

No risk model captures everything, so providers differ by more than chance allows. With overdispersion=True (the default) each provider is judged against that spread as well as chance, using Spiegelhalter's additive random-effects model. The smallest providers' limits barely move. The busiest level off at about the centre ± z·τ instead of closing onto the line. When τ̂ is 0, which is common at physician scale, the chart is identical to overdispersion=False.

How τ is estimated (tau_method=, methods §4.16):

The published estimator understates a small real spread, so average providers flag somewhat more often than the stated rate. That is the right error for a screen whose flags are reviewed. funnel_sensitivity() runs each option and marks the providers whose tier depends on the choice:

from qikit import funnel_sensitivity

s = funnel_sensitivity(df, y="deaths", n="expected", measure="oe", unit="surgeon")
s[~s.consistent]      # flagged under one model and not another: a judgement call

Treat a non-zero τ̂ as a finding about the system: unmeasured case mix, data differences, or real variation in practice.

Length of stay and other clustered counts

A funnel of observed/expected patient-days treats every day as an independent event. Days come in stays, so a provider's day total varies 5–10 times more than that assumption allows. Every small provider's limits come out far too narrow, and every p-value and q-value too small. τ cannot absorb this, because the missing variance shrinks with volume while τ does not.

Pass stay-level rows and let qikit measure it:

stays = session.sql("""
    SELECT provider_id, provider_name, los_days, expected_los_days
    FROM inpatient_stays WHERE discharge_date >= '2026-01-01'
""").to_pandas()

r = funnel(stays, y="LOS_DAYS", n="EXPECTED_LOS_DAYS", measure="oe",
           unit="PROVIDER_ID", label="PROVIDER_NAME",
           within_dispersion="estimate", worse="higher")
r.summary["within_dispersion"]     # φ̂: how much more variable each count is than Poisson

• The estimate: φ̂ comes from each stay's residual around its own provider's rate, so real differences between providers stay with τ.
• Already-aggregated data: pass a number, computed once in SQL: SUM(POWER(los - exp_los, 2) / exp_los) / (COUNT(*) - COUNT(DISTINCT provider)).
• Event counts (deaths, readmissions) need nothing; leave the default of 1.
• The guard: if the spread looks clustered and within_dispersion was left at 1, qikit raises a UserWarning and says so in summary["notes"] and on the figure.
• Upstream choices: cap extreme stays at a stated percentile before aggregating. If φ̂ differs a lot between service lines, facet by service line; "estimate" is then taken per facet.
• Multi-item surveys are the same problem on a proportion: one patient answers the three doctor questions alike, so item counts look like more data than they are. Pass one row per survey with within_dispersion="estimate"; the model then sees y/φ̂ out of n/φ̂ (Rao–Scott effective counts), interpolated rather than rounded (methods §4.13).

In one simulated length-of-stay cohort of 1,000 providers with no real differences between them, the chart flags 24 outliers and 65 alerts without the fix. With within_dispersion="estimate" it flags 3 and 20, against 1 and 24 expected by chance (methods §5 has the full simulation).

Physician-level comparison: what the chart could not find

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 absence of a signal is not evidence of acceptable performance. That is exactly how a blank funnel tends to get read. So every provider carries min_detectable and an exact interval. summary["underpowered"] says whether the typical provider could detect even a doubling, and an underpowered figure captions itself so the caveat survives a slide deck. A physician at O/E 2.5 whose min_detectable is 4.6 has not been cleared; the chart cannot see that far at their volume.

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

import pandas as pd
from qikit import funnel

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

# Any single year: no outlier, and underpowered is True.
funnel(df[df.year == 2025], y="observed_deaths", n="expected_deaths",
       measure="oe", unit="physician").summary["underpowered"]       # True

# All three years (rows per physician are summed): three real signals, SURG-07 and
# SURG-15 above, and a high performer, SURG-21, below the lower limit.
r = funnel(df, y="observed_deaths", n="expected_deaths", measure="oe", unit="physician")
r.data[r.data.signal != "none"]

An achievable goal: benchmark=

The centre line is the average, and the average is not a goal. benchmark= adds one the best providers in the same data already reach:

• "abc", the achievable benchmark of care (Kiefe et al. 1998), for process measures. It is the pooled rate of the best providers who together hold at least 10% of patients. Providers are ranked by (x + 1)/(n + 2) so that 3 of 3 cannot top the list on luck.
• "eb", for outcomes. It is the 90th percentile, toward better, of the shrunk estimates eb_rate (O'Brien et al. 2008).

r = funnel(df, y="hp_within_24h", n="admissions", measure="proportion", unit="physician",
           worse="lower", benchmark="abc", multiply=100)
b = r.summary["benchmark"]
b["value"]                           # e.g. 97.0: the % the best physicians achieve
r.data[r.data.benchmark_group]       # who set it
b["without_small"]                   # does it lean on low-volume physicians?
r.plot()                             # the benchmark as a dashed line

funnel(df, y="deaths", n="expected", measure="oe", unit="surgeon", worse="higher",
       benchmark="eb").summary["benchmark"]["value"]   # an achievable O/E

The ABC's authors designed it for process measures. They urge extreme caution with outcomes, and discourage its use for pay, credentialing or public report cards (methods §4.18). Both are goals for confidential feedback, not tests, so signal is unchanged.

Figures

HOUSE = {"font": {"family": "Segoe UI", "size": 13}, "limits": {"color": "#555"}}

r.plot(
    styles=[HOUSE, {"lcl": {"visible": False}, "lcl_95": {"visible": False}}],
    x_range=(2, None),          # start at 2 expected deaths; the y axis refits
    label_signals="outlier",
    width=1000, height=560,
)

• Defaults for a projector: contours at 3.5:1 contrast, the 95% line dashed, points at 0.85 opacity, and the legend on. Outliers are red circles and alerts amber diamonds, so the tier never depends on colour alone.
• Axes and caption: the y axis is fitted to the providers and the narrow end of the funnel, not its low-volume flare. An unset caption states the flag rule.
• styles= restyles by role: data, outlier, alert, runs, improving, concerning, neutral, cl, ucl, lcl, ucl_95, lcl_95, target, the groups limits and limits_95, end_labels and font. A list of dicts is merged in order, so a house style is one constant.
• r.limits(n=None) returns the contours at any denominators, for drawing elsewhere.
• Unknown names raise. An unknown .plot() option raises TypeError.

Options at a glance

In Snowflake

qic(...).data and funnel(...).data are the tables a pipeline writes back. With facets= a list of columns, one call covers every measure and service line, and the facet columns stay in the output, so the result joins back to the source table. docs/snowflake.md has a stored procedure and a UDTF that add limits to a pre-aggregated metrics table, packaging steps, and Streamlit in Snowflake notes.

Also included

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

For AI agents

• AGENTS.md ships inside the package. It maps each module, says which function answers which question, and lists the conventions an agent must not break. llms.txt is the same map at the repository root.
• Docstrings are the reference. help(qikit.qic) and help(qikit.funnel) open with a decision table and document every result column.
• Skills for Snowflake Cortex Code or any agent that loads them: skills/spc-monitoring/SKILL.md and skills/provider-funnels/SKILL.md.
• Errors name the fix, and removed arguments name their replacement.

Excel add-in

A Microsoft Excel task-pane add-in built on the same engine, ported to TypeScript, lives in excel-addin/. See its README and docs/excel-addin-roadmap.md.

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 .
uv run mypy
uv run pytest tests/

The Python engine is authoritative. The TypeScript port in excel-addin/packages/engine conforms to it through the shared JSON fixtures in fixtures/spc. After an intentional engine change, regenerate the snapshots with uv run python scripts/update_snapshots.py. fixtures/reference holds output from R's qicharts2, FunnelPlotR and NHSRplotthedots, regenerated by scripts/make_reference_fixtures.R; tests/test_reference.py checks qikit against it to 1e-6, with each accepted departure listed and explained.

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.
3. Anhøj J, Olesen AV. Run charts revisited: a simulation study of run chart rules for detection of non-random variation in health care processes. PLoS ONE 2014;9(11):e113825.
4. Spiegelhalter DJ. Funnel plots for comparing institutional performance. Statistics in Medicine 2005;24(8):1185–1202.
5. NHS England. Making Data Count. Strengthening your decisions, 2019.

License

MIT