# Confidence interval for difference between predicted vs actual rates

I would like to generate a confidence interval for predicted vs actual rates.

I am auditing my group of anaesthetists (aka anesthesiologists) to see how we compare on a number of potentially preventable complications (eg post-operative nausea, severe pain, hypothermia).

I have 20000 surgical operation records and I can make a GLM to make a "case-mix adjusted risk" (using age, gender, type of surgery, duration of surgery as risk factors) and thus generate a risk for each patient.

I can then aggregate the risk and actual per clinician I can generate an actual and predicted rate for each complication. I can make a confidence interval for my actual - but it seems a bit simplistic to just test to see if the confidence interval on the actual rate includes the rate generated from summing the glm-predicted risks.

This question has some pointers to a package but I am hoping for some more specific suggestions.

To clarify what I have already (using "requirement for pain protocol" as an example):

# make model (dependent variable has values 1/0)
model.pp = glm(
pain_protocol1 ~
age + log_age + age2 + inv_age
+ op_time + log_op_time + op_time2
+ gender
+ category
+ thimble,
family = "binomial",
data=d4)
# calculate predicted PACU time and then difference between predicted and actual:
d4$pred_pp = predict(model.pp, newdata=d4, type="response", na.action="na.pass") d4$extra_pp = d4$pain_protocol1 - d4$pred_pp
# aggregate deviation from predicted rate
ppr_pa <- aggregate(extra_pp ~ adult_anaesthetist, data=d4, FUN=mean)
barplot(ppr_pa$extra_pp, name=ppr_pa$adult_anaesthetist,las=2)


So I can make this plot for my colleagues, showing the variation we have in how much pain our patients experience in the "post anaesthesia care unit". These variations are great enough so that they definitely represent a material difference in patient experience, and most of the difference will also be unlikely to be variation due to chance (ie "statistically significant"). However, as I examine smaller subgroups and other complications that are less frequent it would be good to be able to calculate confidence intervals.

Note that each clinician has a different number of cases, and each clinician is given a bird code-name for anonymity.

• I'm not sure I follow your scenario. What is your "actual" data? Is it yes/no of something? What? – gung - Reinstate Monica Nov 11 '14 at 2:16
• @gung - sorry yes it is binomial. I want to produce a series of reports/tables/plots for each complication. Each complication will have its own model. – ErichBSchulz Nov 11 '14 at 2:19

mmm... so this is the best answer I have currently - I'll post it here because it gives others a basis to improve on.

My concerns with this approach is that I'm assuming a normal distribution (which it clearly isn't when i plot(density(na.omit(d4$extra_pp)))) and failing to account for the idiosyncrasies of the GLM. That said, it feels like the intervals provide some guidance to my colleagues on how seriously they should take the variation from their colleagues... require(ggplot2) alpha<-0.05 dsbc <- aggregate(d4$extra_pp, by=list(d4$adult_anaesthetist), FUN=mean) dsbc$mean <- dsbc$x dsbc$n <- aggregate(d4$extra_pp, by=list(d4$adult_anaesthetist), FUN=length)$x dsbc$sd <- aggregate(d4$extra_pp, by=list(d4$adult_anaesthetist), FUN=sd)$x dsbc$sem <- dsbc$sd/sqrt(dsbc$n) # standard error of the mean (SEM)
dsbc$me <- qt(1-alpha/2, df=dsbc$n)*dsbc\$sem # margin of error for ci
print(
ggplot(dsbc, aes(x = Group.1, y = mean)) +
geom_bar(position = position_dodge(), stat="identity", fill="blue") +
geom_errorbar(aes(ymin=mean-me, ymax=mean+me))
+ ggtitle("Risk of patients requiring PACU pain-protocols above expected")
+ ylab("Additional rate (with 95% confidence intervals)")
+ xlab("Anaesthetist")
+ theme(axis.text.x  = element_text(angle=90, vjust=0.5, size=8))
)