Monitoring Surgical Outcomes Using a Risk-Adjusted CUSUM Chart
In this example, we demonstrate the use of a risk-adjusted CUSUM chart. The data comes from a center for cardiac surgery in the UK and is available in the cardiacsurgery.csv
file at this Zenodo link.
First, the required packages are loaded
using StatisticalProcessMonitoring, Distributions, Random, Parameters, CSV, DataFrames, CategoricalArrays, MixedModels, Plots
The dataset is then loaded as a DataFrame
and the surgeon
variable is categorized.
julia> dat = CSV.read("cardiacsurgery.csv", DataFrame)
julia> dat.surgeon = categorical(dat.surgeon)
We divide the data into Phase I (first two years) for model estimation and Phase II (following year) for prospective monitoring.
julia> dat_ic = dat[dat.date .<= 730, :]
We estimate the post-operative mortality rate on the Phase I data using a logistic regression model with the Parsonnet score as a covariate and the surgeon as a random intercept.
julia> using MixedModels
julia> mod = fit(MixedModel, @formula(status ~ Parsonnet + (1|surgeon)),
dat_ic, Bernoulli())
julia> print(mod)
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
status ~ 1 + Parsonnet + (1 | surgeon)
Distribution: Bernoulli{Float64}
Link: LogitLink()
logLik deviance AIC AICc BIC
-388.8235 777.6471 783.6471 783.6607 800.0816
Variance components:
Column Variance Std.Dev.
surgeon (Intercept) 0.037837 0.194518
Number of obs: 1769; levels of grouping factors: 6
Fixed-effects parameters:
─────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
─────────────────────────────────────────────────────
(Intercept) -3.65655 0.17509 -20.88 <1e-96
Parsonnet 0.0818093 0.00723527 11.31 <1e-28
─────────────────────────────────────────────────────
For Phase II monitoring, we consider data from the following year.
julia> dat_oc = dat[730 .< dat.date .<= 1095, :]
We use a risk-adjusted CUSUM control chart to monitor potential increases in post-operative patient mortality.
julia> Random.seed!(239184367)
julia> STAT = RiskAdjustedCUSUM(Δ = 0.75, model = mod, response = :status)
We set the IC average run length to 1000 and initialize the control limit.
julia> NOM = ARL(1000)
julia> LIM = OneSidedFixedLimit(1.0, true)
In-control run lengths for estimating the control limit are simulated by resampling the Phase I data using bootstrap.
julia> PH2 = Phase2(Bootstrap(), dat_ic)
After creating the ControlChart
object, we find the control limit using stochastic approximations.
julia> CH = ControlChart(STAT, LIM, NOM, PH2)
julia> h = saCL!(CH, verbose = true, gamma = 0.05)
Running SA ...
Running adaptive gain ...
Estimated gain D = 0.4117860274331409
Running optimization ...
i: 0/50000 h: 2.77856 hm: 0.0 stop: 0
i: 1000/50000 h: 2.86375 hm: 2.9541 stop: 3457
i: 2000/50000 h: 2.94813 hm: 2.95045 stop: 3435
i: 3000/50000 h: 2.92253 hm: 2.95877 stop: 3420
i: 3359/50000 Convergence!
(h = 2.9568877844997226, iter = 3359, status = "Convergence")
We apply the control chart to the Phase II data and plot the results of the monitoring.
julia> proc = apply_chart(CH, dat_oc)
julia> plt = plot_series(proc, dpi = 300, label = "", xlab = L"t",
ylab = L"C_t")