R/scoreyj.R

In fsdaR: Robust Data Analysis Through Monitoring and Dynamic Visualization

######
##  VT::11.05.2026
##
##
##  roxygen2::roxygenise("C:/users/valen/onedrive/myrepo/R/fsdaR", load_code=roxygen2:::load_installed)
##
#'  Computes the score test for Yeo and Johnson transformation
#'
#' @description Computes the score test for Yeo and Johnson transformation
#'
#' @param y Response variable. A vector with \code{n} elements that
#'  contains the response variable.
#'
#' @param x An \code{n x p} data matrix (\code{n} observations and \code{p} variables).
#'  Rows of \code{x} represent observations, and columns represent variables.
#'
#'  Missing values (NA's) and infinite values (Inf's) are allowed,
#'  since observations (rows) with missing or infinite values will
#'  automatically be excluded from the computations.
#'
#' @param intercept wheather to use constant term (default is \code{intercept=TRUE}
#'
#' @param la values of the transformation parameter for which it is necessary
#'  to compute the score test. Default value of lambda is
#'  \code{la=c(-1, -0.5, 0, 0.5, 1)}, i.e., the five most common values of lambda.
#'
#' @param lik likelihood for the augmented model. If true the value of the likelihood
#'  for the augmented model will be calculated and returend otherwise (default) only
#'  the value of the score test will be given
#'
#' @param nocheck Whether to check input arguments. If \code{nocheck=TRUE} no check is performed
#'  on matrix \code{y} and matrix \code{X}. Notice that \code{y} and \code{X}
#'  are left unchanged. In other words the additional column of ones for the
#'  intercept is not added. The default is \code{nocheck=FALSE}.
#'
#' @param tukey1df Tukey's one degree of freedome test for non-additivity.
#'  The constructed variable is given by
#'  \deqn{w_T(\lambda)= (\hat z(\lambda) - \overline  z(\lambda))^2 / 2 \overline  z(\lambda)}{ascii}
#'  where \eqn{z(\lambda)} is the transformed response, and
#'  \eqn{\hat z(\lambda)} are the fitted values on the transformed response. 
#'  The t test on the constructed variable above provides a test from departures from the
#'  assumed linear model and is known in the literature as Tukey's one degree of 
#'  freedome test for non-additivity. If \code{tukey1df=TRUE} the test is computed 
#'  and returned in \code{ScoreT} else (default) the value of the test is not computed.
#' @param trace Whether to print intermediate results. Default is \code{trace=FALSE}.
#' @param \dots potential further arguments passed to lower level functions.
#'
#' @return  An S3 object of class \code{\link{scoreYJ.object}} will be returned which is basically a list
#'  containing the following elements:
#'  \itemize{
#'  \item \code{Score}: score test. A vector of length \code{length(lambda)} which 
#'      contains the value of the score test for each value of lambda specified
#'      in the optional input parameter \code{la}. If \code{la} is not specified, 
#'      the vector will be of length 5 and contains the values of the score test for the
#'      5 most common values of \code{lambda}.
#'  \item \code{ScoreT}: value of the Tukey's one degree of freedome test for 
#'      non-additivity. This output is produced only if \code{tukey1df=TRUE}.
#'  \item \code{Lik}: value of the likelihood. This output is produced only if \code{lik=TRUE}.
#'  }
#'
#' @references
#'  Yeo, I.K. and Johnson, R. (2000), A new family of power
#'   transformations to improve normality or symmetry, "Biometrika", Vol. 87,  pp. 954-959.
#'
#' @examples
#'
#'  \dontrun{
#'  ##  ScoreYJ with all default options for the wool data.
#'  ##  Load the wool data.
#'  
#'    data(wool)
#'    XX <- wool
#'    y <- XX[, ncol(XX)]
#'    X <- XX[, 1:(ncol(XX)-1), drop=FALSE]
#'
#'    (out <- scoreYJ(X, y))                    # call 'scoreYJ' with all default parameters
#'    (out <- scoreYJ(X, y, lik=TRUE))          # return the likelihood
#'
#'    data(loyalty)
#'    head(loyalty)
#'    y=loyalty[, 4]
#'    X=loyalty[, 1:3]
#'
#'    ##    la is a vector containing the values of \lambda which have to be tested
#'    (out <- scoreYJ(X, y, la=c(0.25, 1/3, 0.4, 0.5)))
#'    (out <- scoreYJ(X, y, la=c(0.25, 1/3, 0.4, 0.5), lik=TRUE))
#'  }
#'
#' @export
#' @author FSDA team, \email{valentin.todorov@@chello.at}

scoreYJ <- function(x, y, intercept=TRUE, la=c(-1, -0.5, 0, 0.5, 1), lik=FALSE, nocheck=FALSE, tukey1df=FALSE, trace=FALSE, ...)
{
    if(is.data.frame(x))
      x <- data.matrix(x)
    else if(!is.matrix(x))
      x <- matrix(x, length(x), 1,
                  dimnames = list(names(x), deparse(substitute(x))))
    if(!is.numeric(x)) stop("x is not a numeric")

    if(is.data.frame(y))
      y <- data.matrix(y)
    else if(!is.matrix(y))
      y <- matrix(y, length(y), 1,
                  dimnames = list(names(y), deparse(substitute(y))))
    if(!is.numeric(y)) stop("y is not a numeric")

    storage.mode(x) <- "double"
    storage.mode(y) <- "double"

    dx <- dim(x)
    xn <- (dnx <- dimnames(x))[[2]]
    xn <- if (!is.null(xn))
        xn
    else if (dx[2] > 1)
        paste("X", 1:dx[2], sep = "")
    else if(dx[2])
        "X"
    dimnames(x) <- list(dnx[[1]], xn)

    n <- nrow(x)
    p <- ncol(x)

    control <- list()
    control$intercept <- ifelse(intercept, 1, 0)
    control$la <- la
    if(lik)
        control$Lik <- 1

    if(!is.numeric(nocheck) && !is.logical(nocheck) || length(nocheck) != 1)
        stop("'nocheck' must be logical or numeric of length 1!")
    control$nocheck <- ifelse(nocheck, 1, 0)
    if(!is.numeric(tukey1df) && !is.logical(tukey1df) || length(tukey1df) != 1)
        stop("'tukey1df' must be logical or numeric of length 1!")
    else if(tukey1df)
        control$tukey1df <- 1    

    outclass <- "scoreYJ"

    parlist = c(.jarray(y, dispatch=TRUE), .jarray(x, dispatch=TRUE))
    paramNames = names(control)
    if(trace)
        print(control)

    if(length(paramNames) > 0)
    {
        for (i in 1:length(paramNames)) {
            paramName = paramNames[i]
            paramValue = control[[i]]

            matlabValue = rType2MatlabType(paramName, paramValue)
            parlist = c(parlist, .jnew("java/lang/String", paramName), matlabValue)
        }
    }

    out <- callFsdaFunction("ScoreYJ", "[Ljava/lang/Object;", 1, parlist)
    if(is.null(out))
        return(NULL)

    arr1 = .jcast(out[[1]], "com/mathworks/toolbox/javabuilder/MWStructArray")
    arr = .jnew("org/jrc/ipsc/globesec/sitaf/fsda/FsdaMWStructArray", arr1)

    if(trace) {
        cat("\nReturning from MATLAB ScoreYJ().  Fields returned by MATLAB: \n")
        print(arr$fieldNames())
    }

    Score <- if(as.integer(arr$hasField("Score", as.integer(1))) != 1) NULL
             else as.matrix(.jevalArray(arr$get("Score", as.integer(1)), "[[D", simplify = TRUE))

    ScoreT <- if(as.integer(arr$hasField("ScoreT", as.integer(1))) != 1) NULL
             else as.matrix(.jevalArray(arr$get("ScoreT", as.integer(1)), "[[D", simplify = TRUE))

    Lik <- if(as.integer(arr$hasField("Lik", as.integer(1))) != 1) NULL
             else as.matrix(.jevalArray(arr$get("Lik", as.integer(1)), "[[D", simplify = TRUE))

    la_names <- c(paste0("la=", format(la, digits=2, nsmall=2)))
    if(!is.null(Score)) {
        Score <- Score[, 1]
        names(Score) <- la_names

    }

    ans <- list(call=match.call(), la=la, Score=Score)
    if(lik && !is.null(Lik)) {
        Lik <- Lik[, 1]
        names(Lik) <- la_names
        ans$Lik <- Lik
    }
    if(tukey1df && !is.null(ScoreT)) {
        ans$ScoreT <- ScoreT
    }

    freeMatlabResources(out)

    class(ans) <- outclass
    return (ans)
}