# Comparing assault rates at one facility with two other facilities

Assaults on staff are fairly uncommon in our facilities with between 70-140 in a year at a given facility. Time, similar to 'person years' in epidemiology, is tracked as 'number of bed days' (i.e., number of days a bed is occupied). Bed days per facility range from 75,000 to 135,000.
For the comparisons, I'm thinking to calculate Confidence Intervals modeled with Poisson distribution with a Bonferonni correction on the confidence level for multiple comparisons (comparing x facility with both y and z facilities; hence confidence level set to 0.975[see below]).

1. facility x assaults=71; days stay=93,516
2. facility y assaults=61; days stay=74,272
3. facility z assaults=142; days stay=133,699

R procedures and output (EpiTools package):

1. > pois.exact(71, 93516, conf.level = 0.975)
x    pt         rate        lower        upper conf.level
1    71 93516 0.0007592284 0.0005717496 0.0009874694      0.975

2. > pois.exact(61, 74272, conf.level = 0.975)
x    pt         rate        lower       upper conf.level
1 61 74272 0.0008213055 0.0006038488 0.001090219      0.975

3. > pois.exact(142, 133699, conf.level = 0.975)
x     pt        rate        lower       upper conf.level
1 142 133699 0.001062087 0.0008724176 0.001279991      0.975


As there is overlap of CIs between facility x and both y and z, there is no significant difference in assault rates.

Does this analysis look alright?

• But you haven't really told us what you are trying to do, so I'm not sure if we can tell you if your analysis looks right. May 14, 2019 at 18:21
• One issue you might want to deal with is the people in the facility. Presumably some people have a higher propensity to assault (perhaps on their first day). If stays are shorter, you have more people passing through. if that's the case, you might expect more assaults. This violation of independence means your CIs are not correct. May 14, 2019 at 18:22
• Trying to answer the question: Does facility x have lower staff assault rates than facility y or facility z? May 14, 2019 at 18:22
• Jeremy, thank you for fixing the code formatting! We don't have good reason to believe that people at facility x are at higher risk of committing assaults due to personal factors. We have a reasonably predictive actuarial risk assessment that suggests (in aggregate) there are not demonstrable differences between clients at any of these facilities. Of course, there may be some unknown issues that violate independence assumption, but we don't have any evidence to suggest this is currently the case. May 14, 2019 at 18:31
• I think you are moving toward the right conclusion, but something I'd caution about is that CI overlap doesn't exclude statistical significance a priori. Though if there is significant overlap as I believe is the case here, then it sorta points that way. May 14, 2019 at 20:22

As said in comments, judging statistical significance by overlap of confidence intervals is problematic. So I would prefer to build a model, in this case a Poisson regression model for rate, that is, using the logarithm of beddays as an offset. Then we can test the hypothesis of some differences between facilities with a 2 df test, as one hypothesis, avoiding multiplicity concerns. In R I do that as

    fac_df  <- data.frame(fac=c("X", "Y", "Z"), bed_days=c(93516, 74272, 133699), assaults=c(71, 61, 142))
fac_mod <- glm( assaults ~ offset(log(bed_days))+fac, family=poisson(), data=fac_df)
fac_mod0 <- glm( assaults ~ offset(log(bed_days)), family=poisson(), data=fac_df)
anova(fac_mod0, fac_mod, test="Chi")
Analysis of Deviance Table

Model 1: assaults ~ offset(log(bed_days))
Model 2: assaults ~ offset(log(bed_days)) + fac
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1         2      6.364
2         0      0.000  2    6.364   0.0415 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


so this indicates some differences between facilities. Then you could continue with a post-hoc analysis, looking at individual contrasts. See Post-hoc after GLM.

But see

summary(fac_mod)

Call:
glm(formula = assaults ~ offset(log(bed_days)) + fac, family = poisson(),
data = fac_df)

Deviance Residuals:
[1]  0  0  0

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.18321    0.11868  -60.53   <2e-16 ***
facY         0.07859    0.17458    0.45   0.6526
facZ         0.33569    0.14535    2.31   0.0209 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 6.3640e+00  on 2  degrees of freedom
Residual deviance: 2.6645e-15  on 0  degrees of freedom
AIC: 24.849

Number of Fisher Scoring iterations: 2


this is a saturated model, as evidenced by the deviance residuals of zero. That means that model criticism based on residuals is impossible, which always is a concern. In particular, with Poisson regression overdispersion is always a possibility. To look into that you would need to use data with more details, maybe per week or per month.

• This is awesome Kjetil! I will take some time to digest, but it looks very promising. Thank you! May 15, 2019 at 18:55