R/internal.R

In growthcharts/centile: Translate Measurements, Z-Scores and Centiles with the RIF format

validate_distribution <- function(names, dist) {
  if (length(dist) != 1L) return(FALSE)
  switch(dist,
         NO = all(c("mean", "sd") %in% names),
         LMS = all(c("L", "M", "S") %in% names),
         BCCG = all(c("mu", "sigma", "nu") %in% names),
         BCT = all(c("mu", "sigma", "nu", "tau") %in% names),
         BCPE = all(c("mu", "sigma", "nu", "tau") %in% names),
         MP = any(substr(names, 1L, 1L) == "p") && "mean" %in% names,
         MEA = !any(substr(names, 1L, 1L) == "p") && "mean" %in% names,
         PCT = any(substr(names, 1L, 1L) == "p"),
         ORD = any(substr(names, 1L, 3L) == "cat"),
         FALSE)
}

check_names <- function(df, needed) {
  if (missing(df)) stop("Required argument 'df' not found")
  if (missing(needed)) stop("Required argument 'needed' not found")

  notfound <- is.na(match(needed, names(df)))
  if (any(notfound)) stop("Not found: ", paste(needed[notfound], collapse = ", "))
}

# COPIED FROM dscore Package dd. 12apr2023
# Two functions adapted from gamlss.dist to work better with out-of-range data
#
# To cite package ‘gamlss.dist’ in publications use:
# Stasinopoulos M, Rigby R (2022). _gamlss.dist: Distributions for
# Generalized Additive Models for Location Scale and Shape_. R package version 6.0-3,
# <https://CRAN.R-project.org/package=gamlss.dist>.
#
# @Manual{,
#   title = {gamlss.dist: Distributions for Generalized Additive Models for Location Scale
#     and Shape},
#   author = {Mikis Stasinopoulos and Robert Rigby},
#   year = {2022},
#   note = {R package version 6.0-3},
#   url = {https://CRAN.R-project.org/package=gamlss.dist}}

qBCT <- function (p, mu = 5, sigma = 0.1, nu = 1, tau = 2, lower.tail = TRUE,
                  log.p = FALSE, na.rm = TRUE)
{
  if (any(mu < 0, na.rm = na.rm))
    stop(paste("mu must be positive", "\n", ""))
  if (any(sigma < 0, na.rm = na.rm))
    stop(paste("sigma must be positive", "\n", ""))
  if (any(tau < 0, na.rm = na.rm))
    stop(paste("tau must be positive", "\n", ""))
  if (log.p == TRUE)
    p <- exp(p)
  else p <- p
  if (any(p <= 0, na.rm = na.rm) | any(p >= 1, na.rm = na.rm))
    stop(paste("p must be between 0 and 1", "\n", ""))
  if (lower.tail == TRUE)
    p <- p
  else p <- 1 - p
  if (length(nu) > 1) {
    z <- ifelse((nu <= 0), qt(p * pt(1/(sigma * abs(nu)),
                                     tau), tau), qt(1 - (1 - p) * pt(1/(sigma * abs(nu)),
                                                                     tau), tau))
  }
  else {
    z <- if (nu <= 0)
      qt(p * pt(1/(sigma * abs(nu)), tau), tau)
    else qt(1 - (1 - p) * pt(1/(sigma * abs(nu)), tau), tau)
  }
  if (length(nu) > 1)
    ya <- ifelse(nu != 0, mu * (nu * sigma * z + 1)^(1/nu),
                 mu * exp(sigma * z))
  else if (nu != 0)
    ya <- mu * (nu * sigma * z + 1)^(1/nu)
  else ya <- mu * exp(sigma * z)
  ya[ya < 1.0 & !is.na(ya)] <- 1.0
  ya
}


pBCT <- function (q, mu = 5, sigma = 0.1, nu = 1, tau = 2, lower.tail = TRUE,
                  log.p = FALSE, na.rm = TRUE)
{
  if (any(mu < 0, na.rm = na.rm))
    stop(paste("mu must be positive", "\n", ""))
  if (any(sigma < 0, na.rm = na.rm))
    stop(paste("sigma must be positive", "\n", ""))
  if (any(tau < 0, na.rm = na.rm))
    stop(paste("tau must be positive", "\n", ""))
  # recode D-score < 1 to 1.
  q[q < 1.0 & !is.na(q)] <- 1.0
  # added check
  if (length(nu) == 1L && is.na(nu))
    return(NA)
  if (length(nu) > 1)
    z <- ifelse(nu != 0, (((q/mu)^nu - 1)/(nu * sigma)),
                log(q/mu)/sigma)
  else if (nu != 0)
    z <- (((q/mu)^nu - 1)/(nu * sigma))
  else if (is.na(nu))
    z <- NA
  else z <- log(q/mu)/sigma
  FYy1 <- pt(z, tau)
  if (length(nu) > 1)
    FYy2 <- ifelse(nu > 0, pt(-1/(sigma * abs(nu)), df = tau),
                   0)
  else if (nu > 0)
    FYy2 <- pt(-1/(sigma * abs(nu)), df = tau)
  else FYy2 <- 0
  FYy3 <- pt(1/(sigma * abs(nu)), df = tau)
  FYy <- (FYy1 - FYy2)/FYy3
  if (lower.tail == TRUE)
    FYy <- FYy
  else FYy <- 1 - FYy
  if (log.p == FALSE)
    FYy <- FYy
  else FYy <- log(FYy)
  FYy
}