R/bootclustrange.R

In WeightedCluster: Clustering of Weighted Data

# mysample <- strataSampling(precluststrata, sample.size=sample.size, prob=ac$probs, randomize=0)






bootclustrange <- function(object, seqdata, seqdist.args = list(method = "LCS"),
                           R = 100, sample.size = 1000,
                           parallel = FALSE, progressbar = FALSE,
                           sampling = "clustering", strata = NULL) {
  ## Preparing
  start.time <- Sys.time() # debut du processus
  if (inherits(object, "seqclararange")) {
    medoids <- object$clara[[length(object$clara)]]$medoids
    clustering <- as.data.frame(object$clustering)
  } else {
    clustering <- as.data.frame(object)
    medoids <- NULL
  }
  if (is.null(strata) && sampling == "clustering") {
    strata <- clustering[, ncol(clustering)]
    tt <- prop.table(table(strata))
    if (any(round(tt * sample.size) < 2)) {
      stop("[!] sample size is to small for the stratified sampling of clustering. Consider a minimum value of ", 2 / min(tt))
    }
    sampling <- "strata"
  }
  if (parallel) {
    oplan <- plan(multisession)
    on.exit(plan(oplan), add = TRUE)
  }
  if (progressbar) {
    if (requireNamespace("progress", quietly = TRUE)) {
      old_handlers <- handlers(handler_progress(format = "(:spin) [:bar] :percent | Elapsed: :elapsed | ETA: :eta | :message"))
      if (!is.null(old_handlers)) {
        on.exit(handlers(old_handlers), add = TRUE)
      }
    } else {
      message(" [>] Install the progress package to see estimated remaining time.")
    }
    oldglobal <- handlers(global = TRUE)
    if (!is.null(oldglobal)) {
      on.exit(handlers(global = oldglobal), add = TRUE)
    }
  }
  p <- progressor(R)

  ## Parallel loop
  calc_cqi <- foreach(loop = 1:R, .options.future = list(packages = c("TraMineR", "WeightedCluster"), seed = TRUE, globals = structure(TRUE, add = c("sample.size", "seqdist.args")))) %dofuture% { # on stocke chaque bootstrap

    ## function for stratified subsampling
    stratsampling <- function() {
      tt <- table(strata)
      tosample <- round(prop.table(tt) * sample.size)
      correction <- sample.size - sum(tosample)
      if (correction != 0) {
        correct <- sample(1:length(tosample), abs(correction))
        tosample[correct] <- tosample[correct] + sign(correction)
      }
      mysample <- NULL
      for (u in unique(strata)) {
        mysample <- c(mysample, sample(which(strata == u), tosample[as.character(u)]))
      }
      return(mysample)
    }
    ltime <- proc.time()

    ## Ensure that we have at least one observation per cluster in our bootstrap
    done <- FALSE
    while (!done) {
      if (sampling == "strata") {
        mysample <- stratsampling()
      } else if (sampling == "medoids" && !is.null(medoids)) {
        mysample <- c(medoids, sample.int(nrow(clustering), sample.size - length(medoids)))
      } else {
        mysample <- sample.int(nrow(clustering), sample.size)
      }
      clust_sample <- clustering[mysample, ]
      done <- all(sapply(1:ncol(clustering), function(x) length(unique(clustering[, x])) == length(unique(clust_sample[, x]))))
    }

    ## Subsample
    seqdist.args$seqdata <- seqdata[mysample, ]
    suppressMessages(diss <- do.call(seqdist, seqdist.args))
    cqi <- as.clustrange(clustering[mysample, ], diss = diss)
    rm(diss)
    gc()
    cqi$stats
  }

  reframeData <- function(ll, clust, matrix = FALSE) {
    xx <- lapply(ll, FUN = function(x) x[clust, ])
    if (matrix) {
      return(do.call(rbind, xx))
    }
    return(do.call(c, xx))
  }

  ## Build the corresponding clustrange object.
  ret <- list(clustering = clustering)
  numclust <- ncol(ret$clustering)
  ret$kvals <- numeric(numclust)

  ret$meant <- matrix(-1, nrow = numclust, ncol = 10)
  ret$stderr <- matrix(-1, nrow = numclust, ncol = 10)
  ret$boot <- list()
  for (i in 1:numclust) {
    ret$boot[[i]] <- reframeData(calc_cqi, i, TRUE)
    ret$kvals[i] <- length(unique(ret$clustering[, i]))
    ret$meant[i, ] <- colMeans(ret$boot[[i]])
    ret$stderr[i, ] <- apply(ret$boot[[i]], 2L, function(x) sqrt(var(x)))
  }
  clnames <- paste0("cluster", ret$kvals)
  colnames(ret$clustering) <- clnames
  ret$meant <- as.data.frame(ret$meant)
  ret$stderr <- as.data.frame(ret$stderr)
  colnames(ret$meant) <- colnames(ret$boot[[i]])
  colnames(ret$stderr) <- colnames(ret$boot[[i]])
  rownames(ret$meant) <- clnames
  rownames(ret$stderr) <- clnames
  ret$stats <- ret$meant
  class(ret) <- c("bootclustrange", "clustrange", class(ret))
  return(ret)
}

plot.bootclustrange <- function(x, stat = "noCH", legendpos = "bottomright", norm = "none",
                                withlegend = TRUE, lwd = 1, col = NULL, ylab = "Indicators",
                                xlab = "N clusters", conf.int = 0.95, ci.method = "perc", ci.alpha = .3, line = "median", ...) {
  # 	NextMethod("plot")
  getS3method("plot", "clustrange")(x, stat = stat, legendpos = legendpos, norm = norm,
    withlegend = withlegend, lwd = lwd, col = col, ylab = ylab,
    xlab = xlab, conf.int = conf.int, ci.method = ci.method, ci.alpha = ci.alpha, line = line, ...)
}

print.bootclustrange <- function(x, digits = 2, bootstat = c("mean"), ...) {
  getS3method("print", "clustrange")(x, digits = digits, bootstat = bootstat, ...)
}