R/sir_utils.R

In popEpi: Functions for Epidemiological Analysis using Population Data

#' @title Confidence intervals for the ratio of two SIRs/SMRs
#' @author Matti Rantanen
#' @description Calculate ratio of two SIRs/SMRs and the confidence intervals of the ratio.
#' 
#' @details Function works with pooled sir-objects i.e. the \code{print} argument in \code{sir} is ignored.
#' Also \code{x} and \code{y} can be a vector of two where first index is the
#' observed cases and second is expected cases (see examples).
#' Note that the ratio of two SIR's is only applicable when the age distributions are similar
#' in both populations.
#'
#' \strong{Formula}
#' 
#' The observed number of first sir \code{O1} is considered as a Binomial variable with sample 
#' size of \code{O1+O2}. The confidence intervals for Binomial proportion \code{A} 
#' is solved using \code{exact} or \code{asymptotic} 
#' method. Now the CI for ratio \code{O1/O2} is \code{B = A/(1 - A)}. And further the CI for SIR/SMR 
#' is B*E2/E1. (Ederer and Mantel)
#' 
#' @param x a sir-object or a vector of two; observed and expected cases.
#' @param y a sir-object or a vector of two; observed and expected cases.
#' @param conf.level the type-I error in confidence intervals, default 0.95 for 95\% CI.
#' @param type How the binomial confidence intervals are calculated (default:) \code{exact} or \code{asymptotic}.
#' @param alternative The null-hypothesis test: (default:) \code{two.sided}, \code{less}, \code{greater}
#' @param digits number of digits in the output
#' 
#' @note
#' Parameter \code{alternative} is always \code{two.sided} when parameter 
#' \code{type} is set to \code{asymptotic}.
#' 
#' @examples 
#' ## Ratio for sir-object and the same values given manually:
#' 
#' 
#' ## create example dataset
#' dt1 <- data.frame(obs = rep(c(5,7), 10),
#'                   pyrs = rep(c(250,300,350,400), 5),
#'                   var = 1:20)
#' Ref <- data.frame(obs = rep(c(50,70,80,100), 5),
#'                  pyrs = rep(c(2500,3000,3500,4000), 5),
#'                  var = 1:20)
#' ## sir using the function
#' s1 <- sir(coh.data = dt1, coh.obs = obs, coh.pyrs = pyrs, 
#'           ref.data = Ref, ref.obs = obs, ref.pyrs = pyrs,
#'           adjust = var)
#'
#' ## Ratio is simply 1:
#' sir_ratio(s1, c(120, 150))
#' 
#' @seealso \code{\link{sir}}
#' \href{../doc/sir.html}{A SIR calculation vignette}
#' 
#' @references 
#' Statistics with Confidence: Confidence Intervals and Statistical Guidelines, 
#' Douglas Altman, 2000. ISBN: 978-0-727-91375-3
#' 
#' @family sir functions
#' 
#' @return A vector length of three: sir_ratio, and lower and upper confidence intervals.
#' 
#' @export sir_ratio
#' 
#' @import data.table
#' @import stats


sir_ratio <- function(x, y, digits = 3, alternative = 'two.sided', 
                      conf.level = 0.95, type = 'exact') {
  # prepare input values: x
  # Tests are located in test_sir script.
  if(inherits(x = x, what = 'sir')){
    O1 <- sum(x$observed)
    E1 <- sum(x$expected)
  }
  else if(is.vector(x) && length(x) == 2) {
    O1 <- x[1]
    E1 <- x[2]
  }
  else{
    stop('Input x is not correct: x is neighter a vector of 2 nor sir-object')
  }
  # prepare y:
  if(inherits(y,'sir')){
    O2 <- sum(y$observed)
    E2 <- sum(y$expected)
  }
  else if(is.vector(y) && length(y) == 2) {
    O2 <- y[1]
    E2 <- y[2]
  }
  else{
    stop('Input y is not correct: y is neighter a vector of 2 nor sir-object')
  }
  
  type <- match.arg(type, c('asymptotic', 'exact'), several.ok = FALSE)
  alternative <- match.arg(alternative, c('two.sided','less', 'greater'), several.ok = FALSE)
  # conf.level
  
  p <- O1/(O1+O2)
  if(type == 'asymptotic') {
    alpha <- (1 - conf.level)/2
    Ex <- p + c(-qnorm(1-alpha),qnorm(1-alpha)) * sqrt((1/(O1+O2))*p*(1-p))
    if( alternative != 'two.sided') {
      message('Test changed to two.sided when asymptotic.')
      alternative <- 'two.sided'
    }
  }
  if(type == 'exact') {
    Ex <- binom.test(c(O1,O2), p = 0.5, alternative = alternative, conf.level = conf.level)$conf.int
  }
  B = Ex/(1-Ex)
  
  res <- round(c(sir_ratio = (O1/E1)/(O2/E2), lower=(B*(E2/E1))[1], upper = (B*(E2/E1))[2]), digits = digits)
  return(res)
}