R/qq.R

In lessR: Less Code with More Comprehensive Results

# Quantile-quantile chart support for XY()
#
# A q-q chart is an ordinary XY() scatterplot in which the x-variable is
#   not a column of data but the quantiles of a theoretical distribution,
#   requested with a keyword in the  x  role, as  XY(.normal, Salary).
#   The construction follows the run chart, for which the keyword .index
#   generates the consecutive integers 1 to the number of rows of data.
#
# Each distribution is fit to the data values of  y  by the method of
#   moments, so the generated quantiles are on the scale of  y. Both axes
#   then carry the same units and the reference is the 45-degree line
#   through the origin. Every distribution here is fit in closed form, with
#   no shape parameter left to estimate or to choose.
#
# Called from: XY.R


# the keywords recognized in the  x  role, one per distribution. A keyword
#   names no variable, so XY() matches it without regard to case and passes
#   it on here already lowered, which is the form kept in this file
# Called from: XY.R
.qq.dists <- function()
  c(".normal", ".lognormal", ".exponential", ".uniform")


# the x-axis label for a keyword
# Called from: XY.R
.qq.lab <- function(keyword) {
  switch(keyword,
    ".normal"      = "Normal Quantiles",
    ".lognormal"   = "Lognormal Quantiles",
    ".exponential" = "Exponential Quantiles",
    ".uniform"     = "Uniform Quantiles",
    "Theoretical Quantiles")
}


# Theoretical quantiles, one for each value of  y, returned in the row
#   order of  y  rather than sorted, so that the pairing of a data value
#   with its theoretical quantile survives the deletion of missing data,
#   the grouping by  by, and the paneling by  facet  that XY() applies
#   after the x-variable is established.
#   The i-th smallest value of  y  pairs with the i-th smallest quantile,
#   which is the sorted pairing that defines a q-q chart.
# Called from: XY.R
.qq.quantiles <- function(y, keyword=".normal") {

  if (!is.numeric(y)) {
    cat("\n"); stop(call.=FALSE, "\n------\n",
      "A q-q chart displays the distribution of a continuous variable,\n",
      "  but the  y  variable is not numeric.\n\n")
  }

  ok <- !is.na(y)
  n <- sum(ok)
  if (n < 3) {
    cat("\n"); stop(call.=FALSE, "\n------\n",
      "A q-q chart requires at least 3 non-missing values of  y,\n",
      "  but only ", n, " found.\n\n")
  }
  yv <- y[ok]

  m <- mean(yv)
  s <- stats::sd(yv)
  if (s == 0) {
    cat("\n"); stop(call.=FALSE, "\n------\n",
      "All values of  y  are the same, so there is no distribution\n",
      "  to compare against a theoretical distribution.\n\n")
  }

  # the support of the distribution bounds the data it can be fit to
  if (keyword == ".lognormal"  &&  any(yv <= 0)) {
    cat("\n"); stop(call.=FALSE, "\n------\n",
      "A lognormal distribution is defined only for positive values,\n",
      "  but  y  has ", sum(yv <= 0), " value(s) of zero or less.\n\n",
      "For data that include zero, consider  .exponential\n\n")
  }
  if (keyword == ".exponential"  &&  any(yv < 0)) {
    cat("\n"); stop(call.=FALSE, "\n------\n",
      "An exponential distribution is defined only for non-negative\n",
      "  values, but  y  has ", sum(yv < 0), " negative value(s).\n\n")
  }

  p <- stats::ppoints(n)
  q <- switch(keyword,
    ".normal"      = stats::qnorm(p, m, s),
    ".lognormal"   = stats::qlnorm(p, mean(log(yv)), stats::sd(log(yv))),
    ".exponential" = stats::qexp(p, rate=1/m),
    ".uniform"     = stats::qunif(p, m - sqrt(3)*s, m + sqrt(3)*s))

  # rank of each value of  y  indexes its theoretical quantile
  r <- rank(y, na.last="keep", ties.method="first")
  q[r]
}