## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----echo=FALSE, warning=FALSE, message=FALSE, include = FALSE----------------
library(dplyr)
library(boot)
library(SciViews)
x <- c(
  0.2, 0.5, 1.1, 1.4, 1.8, 2.3, 2.5, 2.7, 3.5, 4.4,
  4.6, 5.4, 5.4, 5.7, 5.8, 5.9, 6.0, 6.6, 7.1, 7.9
)
cqv_versatile <- function(
  x,
  na.rm = FALSE,
  digits = 1,
  method = NULL,
  R = NULL,
  ...
) {
  # require(dplyr)
  # require(SciViews)
  # require(boot)
  if (missing(x) || is.null(x)) {
    stop("object 'x' not found")
  } else if (!missing(x)) {
    x <- x
  }
  if (!is.numeric(x)) {
    stop("argument is not a numeric vector: returning NA")
    return(NA_real_)
  }
  na.rm <- na.rm # removes NAs if TRUE
  if (na.rm == TRUE) {
    x <- x[!is.na(x)]
  } else if (anyNA(x)) {
    stop(
      "missing values and NaN's not allowed if 'na.rm' is FALSE"
    )
  }
  if (is.null(digits)) {
    digits <- 1
  }
  digits <- digits # digits required for rounding
  method <- method # returns 95% confidence interval
  if (is.null(R)) {
    R <- 1000
  }
  q3 <- unname(
    quantile(
      x,
      probs = 0.75, # third quartile (0.75 percentile)
      na.rm = na.rm
    )
  )
  q1 <- unname(
    quantile(
      x,
      probs = 0.25, # first quartile (0.25 percentile)
      na.rm = na.rm
    )
  )
  if (q3 == 0) { # to avoid NaNs when q3 and q1 are zero
    warning(
      "cqv is NaN because q3 and q1 are 0, max was used instead of q3"
    )
    q3 <- max(x, na.rm = na.rm)
  }
  a <- ceiling(
    (length(x) / 4) - (1.96 * (((3 * length(x)) / 16)^(0.5)))
  )
  b <- round(
    (length(x) / 4) + (1.96 * (((3 * length(x)) / 16)^(0.5))),
    digits = 0
  )
  c <- length(x) + 1 - b
  d <- length(x) + 1 - a
  Ya <- dplyr::nth(x, a, order_by = x)
  Yb <- dplyr::nth(x, b, order_by = x)
  Yc <- dplyr::nth(x, c, order_by = x)
  Yd <- dplyr::nth(x, d, order_by = x)
  star <- 0
  for (i in a:(b - 1)) {
    star[i] <- (
      (choose(length(x), i)) * (0.25^(i)) * (0.75^(length(x) - i))
    )
    alphastar <- 1 - sum(star[i], na.rm = na.rm)
  }
  zzz <- qnorm((1 - ((1 - alphastar) / 2)))
  f1square <- (3 * (zzz)^2) / (4 * length(x) * ((Yb - Ya)^2))
  f3square <- (3 * (zzz)^2) / (4 * length(x) * ((Yd - Yc)^2))
  D <- q3 - q1
  S <- q3 + q1
  v <- (
    (1 / (16 * length(x))) * (
      (((3 / f1square) + (3 / f3square) - (2 / sqrt(f1square * f3square))) / D^2) +
        (((3 / f1square) + (3 / f3square) + (2 / sqrt(f1square * f3square))) / S^2) -
        ((2 * ((3 / f3square) - (3 / f1square))) / (D * S))
    )
  )
  ccc <- length(x) / (length(x) - 1)
  upper.tile <- exp(((SciViews::ln((D / S)) * ccc)) + (zzz * (v^(0.5))))
  lower.tile <- exp(((SciViews::ln((D / S)) * ccc)) - (zzz * (v^(0.5))))
  if (
    unname(quantile(x, probs = 0.75, na.rm = na.rm)) != 0
  ) {
    boot.cqv <- boot::boot(
      x,
      function(x, i) {
        round(((
          unname(quantile(x[i], probs = 0.75, na.rm = na.rm)) -
            unname(quantile(x[i], probs = 0.25, na.rm = na.rm))
        ) / (
          unname(quantile(x[i], probs = 0.75, na.rm = na.rm)) +
            unname(quantile(x[i], probs = 0.25, na.rm = na.rm))
        )) * 100, digits = digits)
      },
      R = R
    )
  } else if (
    unname(quantile(x, probs = 0.75, na.rm = na.rm)) == 0
  ) {
    boot.cqv <- boot::boot(
      x,
      function(x, i) {
        round(((
          max(x[i], na.rm = na.rm) -
            unname(quantile(x[i], probs = 0.25, na.rm = na.rm))
        ) / (
          max(x[i], na.rm = na.rm) +
            unname(quantile(x[i], probs = 0.25, na.rm = na.rm))
        )) * 100, digits = digits)
      },
      R = R
    )
  }

  if (is.null(method)) {
    boot.cqv.ci <- NA
  } else if (method == "bonett") {
    boot.cqv.ci <- NA
  } else if (method == "norm") {
    boot.norm.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "norm")
  } else if (method == "basic") {
    boot.basic.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "basic")
  } else if (method == "perc") {
    boot.perc.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "perc")
  } else if (method == "bca") {
    boot.bca.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "bca")
  } else if (method == "all") {
    boot.norm.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "norm")
    boot.basic.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "basic")
    boot.perc.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "perc")
    boot.bca.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "bca")
  }

  if (is.null(method)) {
    cqv <- round(
      100 * ((q3 - q1) / (q3 + q1)),
      digits = digits
    )
  } else if (method == "bonett") {
    cqv <- round(
      100 * ((q3 - q1) / (q3 + q1)),
      digits = digits
    )
  } else if (method == "norm") {
    cqv <- round(
      100 * ((q3 - q1) / (q3 + q1)),
      digits = digits
    )
  } else if (method == "basic") {
    cqv <- round(
      100 * ((q3 - q1) / (q3 + q1)),
      digits = digits
    )
  } else if (method == "perc") {
    cqv <- round(
      100 * ((q3 - q1) / (q3 + q1)),
      digits = digits
    )
  } else if (method == "bca") {
    cqv <- round(
      100 * ((q3 - q1) / (q3 + q1)),
      digits = digits
    )
  } else if (method == "all") {
    cqv <- round(
      100 * ((q3 - q1) / (q3 + q1)),
      digits = digits
    )
  }

  if (is.null(method)) {
    lower <- NA
    upper <- NA
  } else if (method == "bonett") {
    lower <- round(lower.tile * 100, digits = digits)
    upper <- round(upper.tile * 100, digits = digits)
  } else if (method == "norm") {
    lower <- round(boot.norm.ci$normal[2], digits = digits)
    upper <- round(boot.norm.ci$normal[3], digits = digits)
  } else if (method == "basic") {
    lower <- round(boot.basic.ci$basic[4], digits = digits)
    upper <- round(boot.basic.ci$basic[5], digits = digits)
  } else if (method == "perc") {
    lower <- round(boot.perc.ci$percent[4], digits = digits)
    upper <- round(boot.perc.ci$percent[5], digits = digits)
  } else if (method == "bca") {
    lower <- round(boot.bca.ci$bca[4], digits = digits)
    upper <- round(boot.bca.ci$bca[5], digits = digits)
  }

  if (is.null(method)) {
    return(
      list(
        method = "cqv = (q3-q1)/(q3+q1)",
        statistics = data.frame(
          est = cqv,
          row.names = c(" ")
        )
      )
    )
  } else if (method == "bonett" && cqv != 100) {
    return(
      list(
        method = "cqv with Bonett 95% CI",
        statistics = data.frame(
          est = cqv,
          lower = lower,
          upper = upper,
          row.names = c(" ")
        )
      )
    )
  } else if (method == "norm" && cqv != 100) {
    return(
      list(
        method = "cqv with normal approximation 95% CI",
        statistics = data.frame(
          est = cqv,
          lower = lower,
          upper = upper,
          row.names = c(" ")
        )
      )
    )
  } else if (method == "basic" && cqv != 100) {
    return(
      list(
        method = "cqv with basic bootstrap 95% CI",
        statistics = data.frame(
          est = cqv,
          lower = lower,
          upper = upper,
          row.names = c(" ")
        )
      )
    )
  } else if (method == "perc" && cqv != 100) {
    return(
      list(
        method = "cqv with bootstrap percentile 95% CI",
        statistics = data.frame(
          est = cqv,
          lower = lower,
          upper = upper,
          row.names = c(" ")
        )
      )
    )
  } else if (method == "bca" && cqv != 100) {
    return(
      list(
        method = "cqv with adjusted bootstrap percentile (BCa) 95% CI",
        statistics = data.frame(
          est = cqv,
          lower = lower,
          upper = upper,
          row.names = c(" ")
        )
      )
    )
  } else if (
    (
      method == "norm" | method == "bonett" | method == "basic" | method == "perc" |
        method == "bca" | method == "all"
    ) && cqv == 100
  ) {
    warning(
      "All values of t are equal to  100 \n Cannot calculate confidence intervals \n"
    )
    return(
      list(
        method = "cqv with Bonett 95% CI",
        statistics = data.frame(
          est = cqv,
          lower = round(lower.tile * 100, digits = digits),
          upper = round(upper.tile * 100, digits = digits),
          row.names = c(" ")
        )
      )
    )
  } else if (method == "all" && cqv != 100) {
    return(
      list(
        method = "All methods",
        statistics = data.frame(
          row.names = c(
            "bonett",
            "norm",
            "basic",
            "percent",
            "bca"
          ),
          est = c(cqv, cqv, cqv, cqv, cqv),
          lower = c(
            round(lower.tile * 100, digits = digits),
            round(boot.norm.ci$normal[2], digits = digits),
            round(boot.basic.ci$basic[4], digits = digits),
            round(boot.perc.ci$percent[4], digits = digits),
            round(boot.bca.ci$bca[4], digits = digits)
          ),
          upper = c(
            round(upper.tile * 100, digits = digits),
            round(boot.norm.ci$normal[3], digits = digits),
            round(boot.basic.ci$basic[5], digits = digits),
            round(boot.perc.ci$percent[5], digits = digits),
            round(boot.bca.ci$bca[5], digits = digits)
          ),
          description = c(
            "cqv with Bonett 95% CI",
            "cqv with normal approximation 95% CI",
            "cqv with basic bootstrap 95% CI",
            "cqv with bootstrap percentile 95% CI",
            "cqv with adjusted bootstrap percentile (BCa) 95% CI"
          )
        )
      )
    )
  } else {
    stop("method for confidence interval is not available")
    return(NA_real_)
  }
}

## ----eval = TRUE, warning=FALSE, message=FALSE--------------------------------
x <- c(
  0.2, 0.5, 1.1, 1.4, 1.8, 2.3, 2.5, 2.7, 3.5, 4.4,
  4.6, 5.4, 5.4, 5.7, 5.8, 5.9, 6.0, 6.6, 7.1, 7.9
)
cqv_versatile(
  x,
  na.rm = TRUE,
  digits = 3,
  method = "bonett"
)

## ----eval = TRUE, warning=FALSE, message=FALSE--------------------------------
cqv_versatile(
  x,
  na.rm = TRUE,
  digits = 3,
  method = "bca"
)

## ----eval = TRUE, warning=FALSE, message=FALSE--------------------------------
cqv_versatile(
  x,
  na.rm = TRUE,
  digits = 3,
  method = "all"
)