vectorized_confidence_intervals.R

# vectorised confidence intervals for poisson count in R
# wrapper for base::poisson.test
# returns rate and 95% confidence interval
# for use with vectors of counts of person-times
# x = vector of counts
# t = vector of corresponding person-times
# form = output format. T = formatted; F = matrix of rate, lower, upper
# digs = number of digits if form = T
# uses 'mapply' so not actually vectorized for performance
# ... additional arguments passed to poisson.test (e.g. 'confidence.level')
# returns 0, 0, 0 for input values where x or t is NA
# returns errors where a value in x is negative or not an integer

vpt <- function(x, t, form = F, digs = 2, ...) {
  f <- function(xl, tl) {
    if (is.na(xl) | is.na(tl)) return(c(0, 0, 0))
    y <- poisson.test(xl, tl, ...)
    c(xl/tl, y$conf.int[1:2])
  }
  res <- mapply(f, xl = x, tl = t)
  return(if (form) {
    res <- format(round(res, digs), digits = digs, nsmall = digs)
    res <- apply(res, 2, function(x) paste0(x[1], '(', x[2], '-', x[3], ')'))
    res <- gsub(' ', '', res)
    gsub('\\(', ' (', res)
  } else {
    `colnames<-`(t(res), c('rate', 'lower', 'upper'))
  })
}

# for binomial proportions:

propCI <- function(X, N, ..., form = T, digs = 0, pc = T) {
  r <- t(rbind(X/N, mapply(function(x, n) prop.test(x, n, ...)$conf.int[1:2], x = X, n = N)))
  colnames(r) <- c('prop', 'lower', 'upper')
  if (form == T) {
    if (pc == T) {r <- r * 100}
    r <- format(round(r, digs), nsmall = digs, digits = digs)
    r <- paste0(r[,1], '(', r[,2], '-', r[,3], ')')
    gsub('\\(', ' (', gsub(' ', '', r))
  } else {
    r
  }
}

# for rate ratios                           
                           
vrr <- function (c1, c2, t1, t2, form = T, digs = 2) { # regression-based vectorised confidence intervals for rate ratio
  f <- function (c1, c2, t1, t2) {
    m <- glm(c(c1, c2) ~ c('a', 'b') + offset(log(c(t1, t2))), family = 'poisson')
    m <- `dimnames<-`(exp(cbind(coef(m), confint(m))), list(c('intercept', 'rate_ratio'), c('est', 'lower', 'upper')))
    m[2,]
  }
  r <- mapply(f, c1 = c1, c2 = c2, t1 = t1, t2 = t2)
  if (form == T) {
    r <- format(round(r, digs), nsmall = digs, digits = digs)
    r <- paste0(r[1,], '(', r[2,], '-', r[3,], ')')
    gsub('\\(', ' (', gsub(' ', '', r))
  } else {
    t(r)
  }
}