R/JohnsonFit.R

In Haplin: Analyzing Case-Parent Triad and/or Case-Control Data with SNP Haplotypes

#' The Johnson distributions
#'
#' Density of the Johnson distribution; adapted from the orphaned SuppDists package.
#' 
#' 	The Johnson system (Johnson 1949) is a very flexible system for describing statistical
#' 	 distributions. It is defined by
#' 	 \deqn{z=\gamma+\delta \log{f(u)}, u=(x-\xi)/\lambda}{z=gamma+delta log(f(u)), with u=(x-xi)/lambda}
#'   and where \eqn{f( )} has four possible forms:
#'   \tabular{ll}{
#'   	SL:\tab \eqn{f(u)=u} the log normal \cr
#'   	SU:\tab \eqn{f(u)=u+\sqrt{1+u^2}}{f(u)=u+sqrt(1+u^2)} an unbounded distribution\cr
#'   	SB:\tab \eqn{f(u)=u/(1-u)} a bounded distribution\cr
#'   	SN:\tab \eqn{\exp(u)} the normal
#'   }
#'   Estimation of the Johnson parameters may be done from quantiles. The procedure of
#'    Wheeler (1980) is used. They may also be estimated from the moments.
#'    Applied Statistics algorithm 99, due to Hill, Hill, and Holder (1976) has been
#'    translated into C for this implementation.  
#' 
#' @aliases pJohnson JohnsonFit
#' 
#' @param q vector of quantities.
#' @param t observation vector, t=x.
#' @param parms list or list of lists each containing output of \code{JohnsonFit}.
#' @param moment character scalar specifying t: for now only "quant".
#' @param lower.tail logical vector; if TRUE (default), probabilities are \eqn{P[X <= x]},
#'   otherwise, \eqn{P[X > x]}.
#' @param log.p logical vector; if TRUE, probabilities p are given as log(p).
#'
#' @return \code{pJohnson()} gives the distribution function.
#'   \code{JohnsonFit()} outputs a list containing the Johnson parameters
#'   (gamma, delta, xi, lambda, type), where type is one of the Johnson types: "SN", "SL",
#'   "SB", or "SU". \code{JohnsonFit()} does this using 5 order statistics when
#'   moment="quant".
#'
#' @references Hill, I.D., Hill, R., and Holder, R.L. (1976). Fitting Johnson curves
#'    by moments. \emph{Applied Statistics.} AS99;
#'  Johnson, N.L. (1949). Systems of frequency curves generated by methods of translation.
#'    \emph{Biometrika,} \bold{36.} 149-176;
#'  Wheeler, R.E. (1980). Quantile estimators of Johnson curve parameters. \emph{Biometrika.}
#'    \bold{67-3} 725-728
#'
#' @author Bob Wheeler


# THESE FUNCTIONS ARE TAKEN FROM SuppDists (https://cran.r-project.org/package=SuppDists)

pJohnson <-
function (q, parms, lower.tail = TRUE, log.p = FALSE) 
{
    tfun <- function(x) if (x == "SN") 
        1
    else if (x == "SL") 
        2
    else if (x == "SU") 
        3
    else 4
    vecFromList <- function(item, aList) {
        if (!is.list(aList[[1]])) 
            return(aList[[item]])
        else {
            tVec <- vector(length = 0)
            for (i in 1:length(aList)) {
                tVec <- append(tVec, (aList[[i]])[[item]])
            }
        }
        tVec
    }
    gamma <- vecFromList(1, parms)
    delta <- vecFromList(2, parms)
    xi <- vecFromList(3, parms)
    lambda <- vecFromList(4, parms)
    type <- vecFromList(5, parms)
    type <- sapply(type, tfun)
    N <- max(length(gamma), length(q))
    q <- rep(q, length.out = N)
    gamma <- rep(gamma, length.out = N)
    delta <- rep(delta, length.out = N)
    xi <- rep(xi, length.out = N)
    lambda <- rep(lambda, length.out = N)
    type <- rep(type, length.out = N)
    if (lower.tail == TRUE) {
        value <- .C("pJohnsonR", as.double(q), as.double(gamma), 
            as.double(delta), as.double(xi), as.double(lambda), 
            as.integer(type), as.integer(N), val = double(N),PACKAGE="Haplin")$val
    }
    else {
        value <- .C("uJohnsonR", as.double(q), as.double(gamma), 
            as.double(delta), as.double(xi), as.double(lambda), 
            as.integer(type), as.integer(N), val = double(N),PACKAGE="Haplin")$val
    }
    if (log.p == TRUE) 
        value <- log(value)
    value
}

#' @rdname pJohnson
JohnsonFit <-
function (t, moment = "quant")
{
#     firstChar=substring(moment,1,1)
#     if (firstChar=="f") {
#         mom <- moments(t)
#         mu <- mom[[1]]
#         sigma <- mom[[2]]
#         skew <- mom[[3]]
#         kurt <- mom[[4]]
#         value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma), 
#             as.double(skew), as.double(kurt), gamma = double(1), 
#             delta = double(1), xi = double(1), lambda = double(1), 
#             type = integer(1),PACKAGE="Haplin")
#     }
#     else if (firstChar=="u") {
#       mu<-t[1]
#       sigma<-sqrt(t[2])
#       skew<-t[3]/sigma^3
#       kurt<-(t[4]/t[2]^2)-3
#       value <- .C("JohnsonMomentFitR", as.double(mu), as.double(sigma), 
#           as.double(skew), as.double(kurt), gamma = double(1), 
#           delta = double(1), xi = double(1), lambda = double(1), 
#           type = integer(1),PACKAGE="Haplin")
#      }
#     else if (firstChar=="q") {
        input <- quantile(t, probs = c(0.05, 0.206, 0.5, 0.794, 
            0.95), names = FALSE)
        x5 <- input[[1]]
        x20.6 <- input[[2]]
        x50 <- input[[3]]
        x79.4 <- input[[4]]
        x95 <- input[[5]]
        value <- .C("JohnsonFitR", as.double(x95), as.double(x79.4), 
            as.double(x50), as.double(x20.6), as.double(x5), 
            gamma = double(1), delta = double(1), xi = double(1), 
            lambda = double(1), type = integer(1),PACKAGE="Haplin")
#     }
#     else return(NA)
    types <- c("SN", "SL", "SU", "SB")
    list(gamma = value$gamma, delta = value$delta, xi = value$xi, 
        lambda = value$lambda, type = types[value$type])
}