R/lmsUtils.R

In CrescNet/GrowthSDS: Utility Functions for Growth SDS Calculations

#' Build a list which contains the zscore, as well as user defined centiles.
#'
#' @param x        The x value (e.g., (gestational) age).
#' @param y        The y value (e.g., height).
#' @param ref      A data.frame, which contains x and LMS values.
#' @param centiles A vector with centiles to be calculated and added to the resulting list.
#' @param zscores  A vector with zscores to be calculated and added to the resulting list.
#' @return List with zscore and y values for the provided centiles and zscores.
#' @examples
#' evaluateCentiles(0, 42.1, data.frame(x = c(), l = c(), m = c(), s = c()), centiles = c(.03, .5, .9), zscores = c(1.881))
evaluateCentiles <- function(x, y, ref, centiles, zscores) {
  if (missing(ref) | nrow(ref) == 0) return(NA)

  if (nrow(ref) == 1) {
    lms <- ref
  } else {
    lms <- list(
      l = approx(ref$x, ref$l, x, rule = 2)$y,
      m = approx(ref$x, ref$m, x, rule = 2)$y,
      s = approx(ref$x, ref$s, x, rule = 2)$y
    )
  }

  result <- list()

  if (!missing(y))
    result$zscore <-
      qnorm(gamlss.dist::pBCCG(y, mu = lms$m, sigma = lms$s, nu = lms$l))

  if (!missing(centiles)) {
    for (centile in centiles) {
      result[[paste0('p', centile * 100)]] <-
        gamlss.dist::qBCCG(centile, mu = lms$m, sigma = lms$s, nu = lms$l)
    }
  }

  if (!missing(zscores)) {
    for (zscore in zscores) {
      result[[paste(zscore)]] <-
        gamlss.dist::qBCCG(pnorm(zscore), mu = lms$m, sigma = lms$s, nu = lms$l)
    }
  }

  return(result)
}

#' Calculate SDS values.
#' You can transform the resulting SDS values to a p-value by calling pnorm().
#' If you provide SDS values, y values will be returned instead.
#'
#' @param x           x value of the subject (e.g. (gestational) age)
#' @param y           y value of the subject (e.g. height)
#' @param sex         sex of the subject
#' @param measurement measurement from which y originates
#' @param ref         Reference object to be used to calculate the SDS, defaults to \code{\link{kromeyerHauschild}}
#'   either ref or refName must be given
#' @param refName     name of the reference to be used to calculate the SDS values
#' @param recodeSex   list containing a mapping for sex to available sexes of the reference
#' @return SDS value
#' @examples
#' sds(
#'   x   = 0,
#'   y   = 47,
#'   sex = 1,
#'   measurement = 'height',
#'   recodeSex   = list(`1` = 'male', `2` = 'female')
#' )
#' @export
sds <- function(x, y, sex = 'male', measurement = 'height', ref = GrowthSDS::kromeyerHauschild, refName, recodeSex) {
  if (missing(refName) & is.null(ref)) stop("no reference provided, please specify 'ref' or 'refName'")
  if (!missing(refName)) ref <- get(refName)
  if (class(ref) != 'Reference') stop("'ref' must be of class 'Reference'")

  if (!missing(recodeSex))
    sex <- factor(sex, names(recodeSex), unlist(recodeSex, use.names = FALSE))

  refSubset <- subsetReference(ref, measurement, unique(x), unique(sex))

  return(purrr::pmap_dbl(list(x, sex, y), function(x, sex, y) {
    if (is.na(x) | is.na(y)) return(NA)

    evaluateCentiles(
      x   = x,
      y   = y,
      ref = refSubset[[sex]][refSubset[[sex]]$x == x,]
    )$zscore
  }))
}

#' Calculate y values.
#'
#' @param x           x value of the subject (e.g. (gestational) age)
#' @param sds         sds value of the subject (e.g. height sds)
#' @param sex         sex of the subject
#' @param measurement measurement from which sds originates
#' @param ref         Reference object to be used to calculate the centile, defaults to \code{\link{kromeyerHauschild}}
#'   either ref or refName must be given
#' @param refName     name of the reference to be used to calculate the centile values
#' @param recodeSex   list containing a mapping for sex to available sexes of the reference
#' @return centile
#' @examples
#' sds(
#'   x   = 0,
#'   sds = 47,
#'   sex = 1,
#'   measurement = 'height',
#'   recodeSex   = list(`1` = 'male', `2` = 'female')
#' )
#' @export
centile <- function(x, sds, sex = 'male', measurement = 'height', ref = GrowthSDS::kromeyerHauschild, refName, recodeSex) {
  if (missing(refName) & is.null(ref)) stop("no reference provided, please specify 'ref' or 'refName'")
  if (!missing(refName)) ref <- get(refName)
  if (class(ref) != 'Reference') stop("'ref' must be of class 'Reference'")

  if (!missing(recodeSex))
    sex <- factor(sex, names(recodeSex), unlist(recodeSex, use.names = FALSE))

  refSubset <- subsetReference(ref, measurement, unique(x), unique(sex))

  return(purrr::pmap_dbl(list(x, sex, sds), function(x, sex, sds) {
    if (is.na(x) | is.na(sds)) return(NA)
    evaluateCentiles(
      x       = x,
      zscores = sds,
      ref     = refSubset[[sex]][refSubset[[sex]]$x == x,]
    )[[paste(sds)]]
  }))
}

#' Get the x value of the reference to be used for sds calculations.
#'
#' If the provided x value is outside the x range of the reference (see xRanges method),
#' the lowest or highest value of the range is returned respectively.
#'
#' @param x           x value of the subject (e.g. (gestational) age)
#' @param sex         Sex of the subject
#' @param measurement The measurement from which y originates
#' @param refName     The name of the reference to be used to calculate the SDS values
#' @param ref         Reference object to be used to calculate the SDS, defaults to \code{\link{kromeyerHauschild}}
#'   either ref or refName must be given
#' @param recodeSex   list containing a mapping for sex to available sexes of the reference
#' @return x value, adjusted to the selected reference
#' @examples
#'   nearestX(50, sex = 'male', measurement = 'height')
#' @export
nearestX <- function(x, sex = 'male', measurement = 'height', refName, ref = GrowthSDS::kromeyerHauschild, recodeSex = NULL) {
  if (missing(x)) stop("'x' missing with no default")
  if (missing(refName) & is.null(ref)) stop("no reference provided, please specify 'ref' or 'refName'")
  if (!missing(refName)) ref <- get(refName)
  if (class(ref) != 'Reference') stop("'ref' must be of class 'Reference'")

  if (!is.null(recodeSex))
    sex <- factor(sex, names(recodeSex), unlist(recodeSex, use.names = FALSE))

  purrr::pmap_dbl(list(x, sex), function(x, sex) {
    if (!sex %in% names(ref@values[[measurement]])) {
      warning(sprintf("measurement '%s' in reference '%s' has no sex '%s'", measurement, ref@title, sex), call. = FALSE)
      return(NA)
    }
    if (is.na(x)) return(NA)

    refX <- ref@values[[measurement]][[sex]]$x
    approx(refX, refX, xout = x, rule = 2)$y
  })
}