R/fit_functions.R

In wangyuchen/jtrans: Johnson Transformation for Normality

#' Fit Functions for Johnson Curves
#' 
#' Fit functions for three Johnson Curves.
#' 
#' Three types of transformations are SB, SL and SU. Their forms are described 
#' below: \deqn{S_B:   Z = \gamma + \eta * ln((X-\epsilon) / (\lambda + \epsilon
#' - X))} \deqn{S_L:   Z = \gamma + \eta * ln(X - \epsilon)} \deqn{S_U:   Z = 
#' \gamma + \eta * asinh((X-\epsilon) / \lambda)} in whihc Z is the standard 
#' normal varible, and X is the non-normal original data, all the necessary 
#' parameters will be returned. Before fitting these curves, sample quantiles 
#' should be calculated according to z values. the \code{qtls} function here is 
#' to provide every useful parameters for Johnson curve fitting.
#' 
#' These functions could also be used for predicting new values when you have 
#' already fitted a model and obtained a \code{jtrans} object. This could be 
#' done by set the \code{newx} parameter. See examples for details.
#' 
#' Note that when predicting new data, the new data should be from the same 
#' distribution as the original data used for fitting the model. All fits have 
#' certain restrictions on data range, if the new data is outside the range of 
#' the model, the prediction will return NA for all the values. Try to exclude 
#' some out-of-range values and predict again may fix this problem.
#' 
#' @param x the non-normal numerical data.
#' @param q the quantiles and some statistics generated by quantiles, it must be
#'   the return value of \code{qtls} function.
#' @param z a single z value for model fitting. It's returned by 
#'   \code{\link{jtrans}}
#'   
#' 
#' @return return NA when the prediction failed, return a list with 2 component 
#'   when fit succeeded. The first component \code{trans} is the transformed 
#'   value and the second component \code{params} is the parameters used in the 
#'   transformation.
#'   
#' 
#' 
#' @rdname fit_functions
fit_sb <- function(x, q) {
  
  eta <- q$z / acosh(.5 * sqrt((1 + q$xm / q$xu) * (1 + q$xm / q$xl)))  
  gamma <- eta * asinh((q$xm / q$xl - q$xm / q$xu) * 
                         sqrt((1 + q$xm / q$xu) * 
                                (1 + q$xm / q$xl) - 4) / 
                         (2 * (q$xm^2 / q$xl / q$xu - 1))) 
  lambda <- (q$xm * sqrt(((1 + q$xm / q$xu) * 
                            (1 + q$xm / q$xl) - 2)^2 - 4) / 
               (q$xm^2 / q$xl / q$xu - 1))
  epsilon <- .5 * (q$x2 + q$x3 - lambda + 
                     q$xm * (q$xm / q$xl - q$xm / q$xu) / 
                     (q$xm^2 / q$xl / q$xu - 1))
  
  if (is.nan(gamma) | is.nan(epsilon) | eta <= 0 | lambda <= 0) 
    return(NA)
  
  if (all(x > epsilon) & all(x < epsilon + lambda)) {
    return(list(trans  = gamma + eta *  log((x - epsilon) / 
                                              (lambda + epsilon - x)),
                params = c(eta, gamma, lambda, epsilon, q$z)))
  } else return(NA)
}


#' @rdname fit_functions
fit_su <- function (x, q) {
  
  eta <- 2 * q$z / acosh(.5 * (q$xu / q$xm + q$xl / q$xm))
  gamma <- eta * asinh((q$xl / q$xm - q$xu / q$xm) / 
                         (2 * sqrt(q$xu * q$xl / q$xm^2 - 1)))
  lambda <- (2 * q$xm * sqrt(q$xu * q$xl / q$xm^2 - 1) / 
               (q$xu / q$xm + q$xl / q$xm - 2) / 
               sqrt(q$xu / q$xm + q$xl / q$xm + 2))
  epsilon <- .5 * (q$x2 + q$x3 + q$xm * (q$xl / q$xm - q$xu / q$xm) / 
                     (q$xu / q$xm + q$xl / q$xm - 2))

  if (is.nan(gamma) | is.nan(epsilon) | eta <= 0 | lambda <= 0) 
    return(NA)
  
  return(list(trans  = gamma + eta * asinh((x - epsilon) / lambda),
              params = c(eta, gamma, lambda, epsilon, q$z)))
}

#' @rdname fit_functions
fit_sl <- function(x, q) {
  
  if (q$xu / q$xm <= 1) return(NA)
  
  eta <- 2 * q$z / log(q$xu / q$xm)
  gamma <- eta * log((q$xu / q$xm - 1) / sqrt(q$xu * q$xm))
  epsilon <- .5 * (q$x2 + q$x3 - q$xm * (q$xu / q$xm + 1) / 
                     (q$xu / q$xm - 1))
  
  if (is.nan(gamma) | is.nan(epsilon) | eta <= 0) return(NA)
  
  if (all(x > epsilon)) {
    return(list(trans  = gamma + eta * log(x - epsilon), 
                params = c(eta, gamma, NA, epsilon, q$z)))
  } else return(NA)
}