R/myStepFlexmix.R

In schiffner/locClass: Collection of Local Classification Methods

#' Runs \code{\link[flexmix]{flexmix}} repeatedly for different cluster values and returns the result with highest likelihood value.
#'
#' @title Run \code{\link[flexmix]{flexmix}} Repeatedly
#' 
#' @param \dots Passed to \code{\link[flexmix]{flexmix}}.
#' @param cluster List of initial cluster assignments of observations at the start of the EM algorithm.
#' @param verbose If \code{TRUE}, show progress information during computations.
#' @param drop If \code{TRUE} and k is of length 1, then a single flexmix object is returned instead of a "stepFlexmix" object.
#' @param unique If \code{TRUE}, then \code{unique()} is called on the result, see below.
#'
#' @return
#' An object of class \code{"stepFlexmix"} containing the best models with respect to the log likelihood for the different number of components in a slot if \code{length(k)>1}, else directly an object of class \code{"flexmix"}.
#' If \code{unique=FALSE}, then the resulting object contains one model per element of \code{k} (which is the number of clusters the EM algorithm started with). 
#' If \code{unique=TRUE}, then the result is resorted according to the number of clusters contained in the fitted models (which may be less than the number with which 
#' the EM algorithm started), and only the maximum likelihood solution for each number of fitted clusters is kept. This operation can also be done manually by 
#' calling \code{unique()} on objects of class \code{"stepFlexmix"}.
#' 
#' @import flexmix
#' @export
#' 
#' @rdname myStepFlexmix
#' @aliases myStepFlexmix 
#' 
#' @family mixtures
#'
#' @examples
#' library(benchData)
#' data <- xor3Data(500)
#' model <- FLXMCLmultinom(trace = FALSE, decay = 0.1)
#' cluster <- lapply(1:5, kmeans(data$x, centers = 3)$cluster)
#' fit <- myStepFlexmix(y ~ ., data = as.data.frame(data), concomitant = FLXPmultinom(~ x.1 + x.2), model = model, 
#'     cluster = cluster, control = list(verb = 1, tolerance = 10^-1))
#' fit2 <- flexmix(y ~ ., data = as.data.frame(data), concomitant = FLXPmultinom(~ x.1 + x.2), model = model, 
#'     cluster = posterior(fit), control = list(verb = 1))
#' pred <- mypredict(fit2, aggregate = TRUE)
#' mean(max.col(pred[[1]]) != data$y)


myStepFlexmix <- function (..., cluster, verbose = FALSE, drop = TRUE, 
    unique = FALSE) 
{
    MYCALL <- match.call()
    MYCALL1 <- MYCALL
    nrep <- length(cluster)
    bestFlexmix <- function(..., cluster) {
        z = new("flexmix", logLik = -Inf)
        logLiks = rep(NA, length.out = nrep)
        for (m in seq_len(nrep)) {
            if (verbose) 
                cat(" *")
            x = try(flexmix(..., cluster = cluster[[m]]))
            if (!is(x, "try-error")) {
                logLiks[m] <- logLik(x)
                if (logLik(x) > logLik(z)) 
                  z = x
            }
        }
        return(list(z = z, logLiks = logLiks))
    }
    z = list()
    # if (is.null(k)) {
        RET = bestFlexmix(..., cluster = cluster)
        z[[1]] <- RET$z
        logLiks <- as.matrix(RET$logLiks)
        z[[1]]@call <- MYCALL
        z[[1]]@control@nrep <- nrep
        names(z) <- as.character(z[[1]]@k)
        if (verbose) 
            cat("\n")
    # }
    # else {
        # k = as.integer(k)
        # logLiks <- matrix(nrow = length(k), ncol = nrep)
        # for (n in seq_along(k)) {
            # ns <- as.character(k[n])
            # if (verbose) 
                # cat(k[n], ":")
            # RET <- bestFlexmix(..., k = k[n])
            # z[[ns]] = RET$z
            # logLiks[n, ] <- RET$logLiks
            # MYCALL1[["k"]] <- as.numeric(k[n])
            # z[[ns]]@call <- MYCALL1
            # z[[ns]]@control@nrep <- nrep
            # if (verbose) 
                # cat("\n")
        # }
    # }
    logLiks <- logLiks[is.finite(sapply(z, logLik)), , drop = FALSE]
    z <- z[is.finite(sapply(z, logLik))]
    rownames(logLiks) <- lapply(cluster, function(x) {
    	ifelse(is.matrix(x), ncol(x), length(unique(x)))
    })##
# print(logLiks)
    if (!length(z)) 
        stop("no convergence to a suitable mixture")
    if (drop & (length(z) == 1)) {
        return(z[[1]])
    }
    else {
        z <- return(new("stepFlexmix", models = z, k = as.integer(names(z)), 
            nrep = as.integer(nrep), logLiks = logLiks, call = MYCALL))
        if (unique) 
            z <- unique(z)
        return(z)
    }
}