R/CVfromCI_BL.R

In PowerTOST: Power and Sample Size for (Bio)Equivalence Studies

#------------------------------------------------------------------------------
# function to calculate the CV from a given 90% confidence interval
# 
# Author: dlabes, hschuetz
# adapted to unbalanced studies by Benjamin Lang
#------------------------------------------------------------------------------

CVfromCI <- function(pe, lower, upper, n, design = "2x2", alpha = 0.05, 
                     robust = FALSE) 
{
  if (missing(pe)) pe <- NULL
  if (missing(n)) stop("Sample size n must be given!", call. = FALSE)
  # according to Helmut's suggestion
  if ((missing(lower) && missing(upper)) || ((missing(lower) && 
       missing(pe)) || (missing(upper) && is.null(pe)))) {
    stop("At least both CLs or PE and one CL must be given!", call. = FALSE)
  }
  # calculate missing lower or upper
  if (!is.null(pe) && (missing(lower) || missing(upper))) {
    ifelse(missing(lower), lower <- pe^2/upper, upper <- pe^2/lower)
  }
  # calculate missing PE
  if (is.null(pe)) pe <- sqrt(lower * upper)
  # should we really do that?
  if (length(pe) > 1 || length(lower) > 1 || length(upper) > 1)
    stop("pe, lower, upper must be scalars.", call. = FALSE)
  # check design
  d.no <- .design.no(design)
  if (is.na(d.no)) stop("Design ", design, " unknown!", call. = FALSE)
  ades <- .design.props(d.no)
  if (length(n) == 1) {
    # n given as ntotal
    n <- nvec(n = n, grps = ades$steps)
    if (n[1] != n[length(n)]) {
      message("Unbalanced ", design, " design. n(i)= ", paste(n, collapse = "/"),
              " assumed.")
    }
  }
  else {
    # n given as vector of # of subjects in (sequence) groups    
    if (length(n) != ades$steps) stop("Length of n vector must be ", ades$steps, "!")
  }
  nc <- sum(1/n)
  n <- sum(n)
  se.fac <- sqrt(ades$bkni * nc)
  df <- eval(.design.df(ades, robust = robust))
  tval <- qt(1 - alpha, df)
  s1 <- (log(pe) - log(lower))/se.fac/tval
  s2 <- (log(upper) - log(pe))/se.fac/tval
  sw <- 0.5 * (s1 + s2)
  if (abs(s1 - s2)/sw > 0.1) {
    warning(paste("sigma based on pe & lower CL more than 10% different than\n", 
            "sigma based on pe & upper CL. Check input."), call. = FALSE)
  }
  return(se2CV(sw))
}
# ---------------------------------------------------------------------------
# alias to CVfromCI
CI2CV <- function(pe, lower, upper, n, design="2x2", alpha=0.05, robust=FALSE)
{
  CVfromCI(pe, lower, upper, n, design=design, alpha=alpha, robust=robust)
}