R/tclustregIC.R

In fsdaR: Robust Data Analysis Through Monitoring and Dynamic Visualization

######
##  VT::22.11.2019
##
##
##  roxygen2::roxygenise("C:/users/valen/onedrive/myrepo/R/fsdaR", load_code=roxygen2:::load_installed)
##
#'  Computes \code{tclustreg} for different number of groups \code{k}
#'  and restriction factors \code{c}.
#'
#' @description The last two letters stand for 'Information Criterion'. This function computes
#'  the values of BIC (MIXMIX), ICL (MIXCLA) or CLA (CLACLA), for different values
#'  of \code{k} (number of groups) and different values of \code{c}
#'  (restriction factor for the variances of the residuals), for
#'  a prespecified level of trimming. In order to minimize randomness, given \code{k},
#'  the same subsets are used for each value of \code{c}.
#'
#' @param y Response variable. A vector with \code{n} elements that
#'  contains the response variable.
#'
#' @param x An n x p data matrix (n observations and p variables).
#'  Rows of 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 alphaLik Trimming level, a number between 0 and 0.5 or an
#'  integer number specifying the number of observations which have to be trimmed.
#'  If \code{alphaLik=0}, there is no trimming.  More in detail, if \code{0 < alphaLik < 1}
#'  clustering is based on \code{h = floor(n * (1 - alphaLik))} observations.
#'  If \code{alphaLik} is an integer greater than 1 clustering is
#'  based on \code{h = n - floor(alphaLik)}. The likelihood
#'  contributions are sorted and the units associated with the smallest \code{n - h}
#'  contributions are trimmed.
#'
#' @param alphaX Second-level trimming or constrained weighted model for \code{x}.
#   \code{alphaX} is a value in the interval \code{[0 1]}.
#'  \itemize{
#'      \item If \code{alphaX=0} there is no second-level trimming.
#'      \item If \code{alphaX} is in the interval [0, 0.5] it indicates the
#'              fixed proportion of units subject to second level trimming.
#'              In this case alphaX is usually smaller than alphaLik.
#'              For further details see Garcia-Escudero et. al. (2010).
#'      \item If \code{alphaX} is in the interval [0.5, 1], it indicates a
#'              Bonferronized confidence level to be used to identify the
#'              units subject to second level trimming. In this case the
#               proportion of units subject to second level trimming is not
#'              fixed a priori, but is determined adaptively.
#'              For further details see Torti et al. (2018).
#'      \item If \code{alphaX=1}, constrained weighted model for \code{X} is assumed
#'              (Gershenfeld, 1997). The CWM estimator is able to
#'              take into account different distributions for the explanatory
#'              variables across groups, so overcoming an intrinsic limitation
#'              of mixtures of regression, because they are implicitly
#'              assumed equally distributed. Note that if \code{alphaX=1} it is
#'              also possible to apply using restrfactor(2) the constraints
#'              on the cov matrices of the explanatory variables.
#'              For further details about CWM see Garcia-Escudero et al.
#'              (2017) or Torti et al. (2018).
#'  }
#'
#' @param intercept wheather to use constant term (default is \code{intercept=TRUE}
#'
#' @param cc a vector specifying the values of the restriction factor which have to 
#'  be considered for the variances of the residuals of the regression lines. 
#'  By default \code{cc=c(1, 2, 4, 8, 16, 32, 64, 128)}.
#' @param ccSigmaX A number specifying the value of the restriction factor which has to be
#'  considered for the covariance matrices of the explanatory variables. The default value 
#'  is \code{ccsigmaX=12}. Note that this option is used only if \code{alphaX=1}, that is 
#'  if constrained weighted model (CWM) for \code{x} is assumed.
#' @param kk an integer vector specifying the number of mixture components (clusters) 
#'  for which the information criteria are be calculated. By default \code{kk=1:5}.
#' @param whichIC A character value which specifies which information criteria must be computed
#'  for each \code{k} (number of groups) and each value of the restriction factor 
#'  \code{c}. Possible values for \code{whichIC} are:
#'  \itemize{
#'   \item "MIXMIX": a mixture model is fitted and for computing the information criterion
#'      the mixture likelihood is used. This option corresponds to the use of the Bayesian
#'      Information criterion (BIC). In output just the matrix \code{MIXMIX} is given.
#'  \item "MIXCLA": a mixture model is fitted but to compute the information criterion
#'      the classification likelihood is used. This option corresponds to the use of the
#'      Integrated Complete Likelihood (ICL). In the output just the matrix \code{MIXCLA} is given.
#'  \item "CLACLA": everything is based on the classification likelihood. This information
#'      criterion will be called CLA. In the output just the matrix \code{CLACLA} is given.
#'  \item "ALL": both classification and mixture likelihood are used. In this case all
#'      three information criteria CLA, ICL and BIC are computed. In the output all
#'      three matrices \code{MIXMIX}, \code{MIXCLA} and \code{CLACLA} are given.
#'  }
#'
#' @param nsamp If a scalar, it contains the number of subsamples which will be extracted.
#'  If \code{nsamp = 0} all subsets will be extracted. Remark - if the number of all possible
#'  subset is greater than 300 the default is to extract all subsets, otherwise just 300.
#'  If \code{nsamp} is a matrix it contains in the rows the indexes of the subsets which
#'  have to be extracted. \code{nsamp} in this case can be conveniently generated by
#'  function \code{subsets()}. \code{nsamp} must have \code{k * p} columns. The first \code{p}
#'  columns are used to estimate the regression coefficient of group 1, ..., the last \code{p}
#'  columns are used to estimate the regression coefficient of group \code{k}.
#'
#' @param RandNumbForNini pre-extracted random numbers to initialize proportions.
#'  Matrix of size k-by-nrow(nsamp) containing the random numbers which
#'  are used to initialize the proportions of the groups. This option is effective only if
#'  \code{nsamp} is a matrix which contains pre-extracted subsamples. The purpose of this
#'  option is to enable the user to replicate the results when the function \code{tclustreg()}
#'  is called using a parfor instruction (as it happens for example in routine IC, where
#'  \code{tclustreg()} is called through a parfor for different values of the restriction factor).
#'  The default is that \code{RandNumbForNini} is empty - then uniform random numbers are used.
#'
#' @param refsteps Number of refining iterations in each subsample. Default is \code{refsteps=10}.
#'  \code{refsteps = 0} means "raw-subsampling" without iterations.
#'
#' @param reftol Tolerance of the refining steps. The default value is 1e-14
#'
#' @param equalweights A logical specifying wheather cluster weights in the concentration
#'  and assignment steps shall be considered. If \code{equalweights=TRUE} we are (ideally)
#'  assuming equally sized groups, else if \code{equalweights = false} (default) we allow for
#'  different group weights. Please, check in the given references which functions
#' are maximized in both cases.
#'
#' @param startv1 How to initialize centroids and covariance matrices. Scalar.
#'  If \code{startv1=1} then initial centroids and covariance matrices are based
#'  on \code{(p+1)} observations randomly chosen, else each centroid is initialized
#'  taking a random row of input data matrix and covariance matrices are initialized
#'  with identity matrices. The default value is\code{startv1=1}.
#'
#'  Remark 1: in order to start with a routine which is in the required parameter space,
#'  eigenvalue restrictions are immediately applied.
#'
#' Remark 2 - option \code{startv1} is used only if \code{nsamp} is a scalar
#'  for more details see the help associated with \code{nsamp}).
#'
#' @param plot If \code{plot=FALSE} (default) or \code{plot=0}  no plot is produced.
#'  If \code{plot=TRUE} a plot with the final allocation is shown (using the spmplot function).
#'  If \code{X} is 2-dimensional, the lines associated to the groups are shown too.
#' @param msg  Controls whether to display or not messages on the screen If \code{msg==TRUE} (default)
#'  messages are displayed on the screen. If \code{msg=2}, detailed messages are displayed,
#'  for example the information at iteration level.
#' @param nocheck Check input arguments. If \code{nocheck=TRUE} no check is performed
#'  on matrix \code{X}. The default is \code{nocheck=FALSE}.
#' @param UnitsSameGroup List of the units which must (whenever possible) have
#'  a particular label. For example \code{UnitsSameGroup=c(20, 26)}, means that
#'  group which contains unit 20 is always labelled with number 1. Similarly,
#'  the group which contains unit 26 is always labelled with number 2, (unless
#'  it is found that unit 26 already belongs to group 1).
#'  In general, group which contains unit \code{UnitsSameGroup(r)} where \code{r=2, ...length(kk)-1}
#'  is labelled with number \code{r} (unless it is found that unit \code{UnitsSameGroup(r)}
#'  has already been assigned to groups \code{1, 2, ..., r-1}. The default value of \code{UnitsSameGroup} is
#'  an empty list, that is consistent labels are not imposed.
#' @param we Weights. A vector of size n-by-1 containing application-specific weights
#'    Default is a vector of ones.
#' @param commonslope wheather to impose a constraint of common slope on the regression coefficients.
#'  If \code{commonslope=TRUE}, the groups are forced to have the same regression coefficients 
#'  (apart from the intercepts). The default value of commonslope is \code{commonslope=FALSE}. 
#' @param Ysave weather to save on output the unput response variable \code{y} and matrix of predictors \code{x}.
#' @param trace Whether to print intermediate results. Default is \code{trace=FALSE}.
#'
#' @param ... potential further arguments passed to lower level functions.
#'
#' @return  An S3 object of class \code{tclustregic} which is basically a list with the following componnts
#'  \itemize{
#'      \item call the matched call
#'      \item CLACLA A matrix of size 5-times-8 if \code{kk} and \code{cc} are not
#'          specififed else it is a matrix of size \code{length(kk)}-times-\code{length(cc)}
#'          containinig the value of the penalized classification likelihood. 
#'          This output is present only if \code{whichIC="CLACLA")} or \code{whichIC="ALL")}.
#'      \item IDXCLA array of size 5-times-8 if \code{kk} and \code{cc} are not
#'          specififed else it is an array of size \code{length(kk)}-times-\code{length(cc)}.
#'          Each element of the array is a list with one element which is a vector of length \code{n} containinig the assignment 
#'          of each unit using the classification model. This output is present only
#'          if \code{whichIC="CLACLA")} or \code{whichIC="ALL")}. 
#'      \item MIXMIX A matrix of size 5-times-8 if \code{kk} and \code{cc} are not
#'          specififed else it is a matrix of size \code{length(kk)}-times-\code{length(cc)}
#'          containinig the value of the penalized mixture likelihood. 
#'          This output is present only if \code{whichIC="MIXMIX")} or \code{whichIC="ALL")}.
#'      \item MIXCLA A matrix of size 5-times-8 if \code{kk} and \code{cc} are not
#'          specififed else it is a matrix of size \code{length(kk)}-times-\code{length(cc)}
#'          containinig the value of the ICL. 
#'          This output is present only if \code{whichIC="MIXCLA")} or \code{whichIC="ALL")}.
#'      \item IDXMIX array of size 5-times-8 if \code{kk} and \code{cc} are not
#'          specififed else it is an array of size \code{length(kk)}-times-\code{length(cc)}.
#'          Each element of the array is a list with one element which is a vector of length \code{n} containinig the assignment 
#'          of each unit using the mixture model. This output is present only
#'          if \code{whichIC="MIXMIX")}, \code{whichIC="MIXCLA")} or \code{whichIC="ALL")}. 
#'      \item kk a vector containing the values of \code{k} (number of components) which have been considered.
#'          This vector is identical to the argument \code{kk} (default is \code{kk=1:5}.
#'      \item cc a vector containing the values of \code{c} (values of the restriction factor) which
#'          have been considered for the variance of the residuals. This vector is identical 
#'          to the argument \code{cc} (defalt is \code{cc=c(1, 2, 4, 8, 16, 32, 64, 128)}.
#'      \item ccSigmaX values of the restriction factor which
#'          have been considered for the covariance matrices of the esplnatory variables. 
#'          This vector is identical the argument \code{ccsigmaX}.
#'      \item alpha the trimming level which has been used in the likelidood (it stores the values of input alphaLik).
#'      \item alphaX second-level trimming or constrained weighted model for \code{X}.
#'      \item X original data matrix of explanatory variables. Present if \code{Ysave=TRUE}.
#'      \item y original vector containing the response. Present if \code{Ysave=TRUE}.
#'  }
#'
#' @references
#'
#'      Torti F., Perrotta D., Riani, M. and Cerioli A. (2019). Assessing Robust Methodologies for Clustering Linear Regression Data,
#'      Advances in Data Analysis and Classification, Vol. 13, pp 227-257.
#'
#' @examples
#'  \dontrun{

#'  ## The X data have been introduced by Gordaliza, Garcia-Escudero & Mayo-Iscar (2013).
#'  ## The dataset presents two parallel components without contamination.
#'
#'  data(X)
#'  y1 = X[, ncol(X)]
#'  X1 = X[,-ncol(X), drop=FALSE]
#'
#'  (out <- tclustregIC(y1, X1, plot=TRUE))
#'
#'  tclustICplot(out, whichIC="MIXMIX")
#'
#'  }
#' @export
#' @author FSDA team, \email{valentin.todorov@@chello.at}

tclustregIC <- function(y, x, alphaLik=0, alphaX=1, intercept=TRUE, 
        whichIC=c("ALL", "MIXMIX", "MIXCLA", "CLACLA"),
        kk=1:5, cc=c(1, 2, 4, 8, 16, 32, 64, 128), ccSigmaX=12,
        plot=FALSE,
        nsamp, refsteps=10, reftol=10e-14, equalweights=FALSE, we, msg=TRUE, nocheck=FALSE, 
        RandNumbForNini, startv1=1, UnitsSameGroup, commonslope=FALSE, Ysave=TRUE,
        trace=FALSE, ...)
{

    whichIC <- match.arg(whichIC)

    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")

    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)

    if(length(alphaLik) != 1 || as.double(alphaLik) < 0 || as.double(alphaLik) > 1)
        stop("'alphaLik must be a number in the interval [0,1]")

    if(length(alphaX) != 1  || as.double(alphaX) < 0 || as.double(alphaX) > 1)
        stop("'alphaX must be a number in the interval [0,1]")
        
    control <- list(...)

    control$kk <- kk
    control$cc <- cc
    control$ccSigmaX <- ccSigmaX
    control$alphaLik <- as.double(alphaLik)
    control$alphaX <- as.double(alphaX)
    control$whichIC <- whichIC
    control$intercept <- ifelse(intercept, 1, 0)
    control$plots <- ifelse(plot, 1, 0)

    if(!missing(nsamp)){
        control$nsamp <- nsamp
        if(is.matrix(nsamp) && !missing(RandNumbForNini))
            control$RandNumbForNini <- RandNumbForNini
    }
    control$refsteps <- refsteps
    control$reftol <- reftol
    control$equalweights <- equalweights
    if(!missing(we))
    {
        if(!is.numeric(we) || length(we) != n)
            stop("Parameter 'we' must be a numeric vector of length ", n)
        control$we <- we
    }

    if(!missing(UnitsSameGroup))
        control$UnitsSameGroup <- UnitsSameGroup
    xmsg <- 0
    if(is.logical(msg))
        xmsg <- ifelse(msg, 1, 0)
    else  if(is.numeric(msg) && msg >= 0 && msg <= 2)
        xmsg <- msg
    else
        stop("Invalid parameter 'msg'. Should be TRUE/FALSE or 0, 1, 2.")
    control$msg <- xmsg
    xtemp <- 0
    if(is.logical(nocheck))
        xtemp <- ifelse(nocheck, 1, 0)
    else  if(is.numeric(nocheck) && nocheck >= 0 && nocheck <= 1)
        xtemp <- nocheck
    else
        stop("Invalid parameter 'nocheck'. Should be TRUE/FALSE or 0, 1.")
    control$nocheck <- xtemp
    control$startv1 <- startv1
    control$commonslope <- commonslope
    control$Ysave <- Ysave
    outclass <- "tclustregic"

    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("tclustregIC", "[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 tclustregIC().  Fields returned by MATLAB: \n")
        print(arr$fieldNames())
    }

    kk_ret <- as.vector(as.matrix(.jevalArray(arr$get("kk", as.integer(1)), "[[D", simplify = TRUE)))
    cc_ret <- as.vector(as.matrix(.jevalArray(arr$get("cc", as.integer(1)), "[[D", simplify = TRUE)))
    alpha_ret = as.vector(as.matrix(.jevalArray(arr$get("alpha", as.integer(1)), "[[D", simplify = TRUE)))[1]
    alphaX_ret <- as.vector(as.matrix(.jevalArray(arr$get("alphaX", as.integer(1)), "[[D", simplify = TRUE)))
    ccSigmaX_ret <- as.vector(as.matrix(.jevalArray(arr$get("ccSigmaX", as.integer(1)), "[[D", simplify = TRUE)))

    MIXMIX <- if(as.integer(arr$hasField("MIXMIX", as.integer(1))) != 1) NULL
                 else as.matrix(.jevalArray(arr$get("MIXMIX", as.integer(1)), "[[D", simplify = TRUE))
    MIXCLA <- if(as.integer(arr$hasField("MIXCLA", as.integer(1))) != 1) NULL
                 else as.matrix(.jevalArray(arr$get("MIXCLA", as.integer(1)), "[[D", simplify = TRUE))
    IDXMIX <- if(as.integer(arr$hasField("IDXMIX", as.integer(1))) != 1) NULL
              else unwrapComplexNumericCellArray(as.matrix(.jevalArray(arr$get("IDXMIX", as.integer(1)))))
    CLACLA <- if(as.integer(arr$hasField("CLACLA", as.integer(1))) != 1) NULL
                 else as.matrix(.jevalArray(arr$get("CLACLA", as.integer(1)), "[[D", simplify = TRUE))
    IDXCLA <- if(as.integer(arr$hasField("IDXCLA", as.integer(1))) != 1) NULL
              else unwrapComplexNumericCellArray(as.matrix(.jevalArray(arr$get("IDXCLA", as.integer(1)))))

    xkk <- paste0("k=", kk_ret)
    xcc <- paste0("c=", cc_ret)
    
    if(!is.null(MIXMIX))
        dimnames(MIXMIX) <- list(xkk, xcc)
    if(!is.null(MIXCLA))
        dimnames(MIXCLA) <- list(xkk, xcc)
    if(!is.null(CLACLA))
        dimnames(CLACLA) <- list(xkk, xcc)
    if(!is.null(IDXMIX))
        dimnames(IDXMIX) <- list(xkk, xcc)
    if(!is.null(IDXCLA))
        dimnames(IDXCLA) <- list(xkk, xcc)

    ans <- list(call=match.call(), CLACLA=CLACLA, IDXCLA=IDXCLA, MIXMIX=MIXMIX, MIXCLA=MIXCLA, IDXMIX=IDXMIX,
            kk=kk_ret, cc=cc_ret,
            alpha=alpha_ret, alphaX=alphaX_ret, ccSigmaX=ccSigmaX_ret, whichIC=whichIC)

    if(Ysave) {
        ans$y <- y
        ans$X <- x
    }
    
    freeMatlabResources(out)

    ## Remove any NULL elements (this happens if whichIC != ALL)
    if(length(del <- which(unlist(lapply(ans, FUN=is.null)))) > 0)
        ans <- ans[-del]

    if(trace) {
        cat("\ntclustIC(): object 'out' after removing the NULL objects:")
        print(names(ans))
    }

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