Comparing assault rates at one facility with two other facilities

Question

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?

Answer 1

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.