Nothing
#' An S4 class to represent the data to use in a DMBC model.
#'
#' @slot diss A list whose elements are the dissimilarity matrices corresponding
#' to the judgments expressed by the \emph{S} subjects/raters. These matrices
#' must be defined as a \code{dist} object.
#' @slot n A length-one character vector representing the number of objects
#' compared by each subject.
#' @slot S A length-one numeric vector representing the number of subjects.
#' @slot family A length-one character vector representing the type of data to
#' analyze. Currently, it accepts only the 'binomial' value, but future
#' developments will include the possibility to analyze continuous,
#' multinomial and count data.
#'
#' @name dmbc_data-class
#' @rdname dmbc_data-class
#' @aliases dmbc_data
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @references
#' Venturini, S., Piccarreta, R. (2021), "A Bayesian Approach for Model-Based
#' Clustering of Several Binary Dissimilarity Matrices: the \pkg{dmbc}
#' Package in \code{R}", Journal of Statistical Software, 100, 16, 1--35, <10.18637/jss.v100.i16>.
#'
#' @examples
#' showClass("dmbc_data")
#'
#' @exportClass dmbc_data
setClass(Class = "dmbc_data",
slots = c(
diss = "list",
n = "numeric",
S = "numeric",
family = "character"
)
)
#' Create an instance of the \code{dmbc_data} class using new/initialize.
#'
#' @param .Object Prototype object from the class \code{\link{dmbc_data}}.
#' @param diss A list whose elements are the dissimilarity matrices corresponding
#' to the judgments expressed by the \emph{S} subjects/raters. These matrices
#' must be defined as a \code{dist} object.
#' @param n A length-one character vector representing the number of objects
#' compared by each subject.
#' @param S A length-one numeric vector representing the number of subjects.
#' @param family A length-one character vector representing the type of data to
#' analyze. Currently, it accepts only the 'binomial' value, but future
#' developments will include the possibility to analyze continuous,
#' multinomial and count data.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases initialize,dmbc_data-method
#' @aliases dmbc_data-initialize
#'
#' @importFrom methods initialize
#' @exportMethod initialize
setMethod("initialize", "dmbc_data",
function(
.Object,
diss = list(),
n = numeric(),
S = numeric(),
family = character()
)
{
.Object@diss <- diss
.Object@n <- n
.Object@S <- S
# .Object@family <- family
.Object@family <- "binomial"
.Object
}
)
#' Show an instance of the \code{dmbc_data} class.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @param object An object of class \code{\link{dmbc_data}}.
#'
#' @aliases show,dmbc_data-method
#' @aliases dmbc_data-show
#'
#' @importFrom methods show
#' @exportMethod show
setMethod("show",
"dmbc_data",
function(object) {
cat("Observed dissimilarity matrices to use in a DMBC analysis\n")
cat("Number of objects (n):", object@n, "\n")
cat("Number of subjects (S):", object@S, "\n")
cat("Family:", object@family, "\n")
}
)
#' Provide a summary of a \code{dmbc_data} class instance.
#'
#' @param object An object of class \code{\link{dmbc_data}}.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases summary,dmbc_data-method
#' @aliases dmbc_data-summary
#'
#' @exportMethod summary
setMethod("summary",
"dmbc_data",
function(object) {
show(object)
cat("Observed dissimilarities:\n")
wd <- as.integer(options("width"))
for (s in 1:object@S) {
dash_string_L <- strrep("-", floor((wd - nchar(paste0("Subject ", object@S)) - 2)/2))
dash_string_R <- strrep("-", floor((wd - nchar(paste0("Subject ", object@S)) - 2)/2) - nchar(s) + 1)
cat(dash_string_L, " ", "Subject ", s, " ", dash_string_R, "\n", sep = "")
diss_nm <- attr(object@diss[[s]], "Labels")
print_matrix(as.matrix(object@diss[[s]]), rownm = diss_nm, colnm = diss_nm,
ndigits = ifelse(object@family != "normal", 0, 2),
colwidth = ifelse(is.null(diss_nm), 5, max(nchar(diss_nm))),
isint = (object@family != "normal"), between_cols = 1)
}
}
)
#' Provide a graphical summary of a \code{dmbc_data} class instance.
#'
#' @param x An object of class \code{\link{dmbc_data}}.
#' @param colors A character vector providing the colors to use in the plot.
#' @param font A length-one numeric vector for the font to use for text.
#' Can be a vector. \code{NA} values (the default) mean use \code{par("font")}.
#' @param cex.font A length-one numeric vector for the character expansion
#' factor. \code{NULL} and \code{NA} are equivalent to \code{1.0}. This is an
#' absolute measure, not scaled by \code{par("cex")} or by setting
#'' \code{par("mfrow")} or \code{par("mfcol")}. Can be a vector.
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases plot,dmbc_data-method
#' @aliases dmbc_data-plot
#'
#' @exportMethod plot
#'
#' @examples
#' data(simdiss)
#' library(bayesplot)
#' cols <- color_scheme_set("brightblue")
#' plot(simdiss, colors = unlist(cols)[c(1, 6)], font = 1, cex.font = 0.75)
setMethod("plot",
signature(x = "dmbc_data"),
function(x, colors = c("white", "black"), font = NA, cex.font = NA, ...) {
if (x@family == "binomial") {
plot_dmbc_data <- function(ncat, colors, crit, INPUT, MD_SOL, MD_DIM,
labxaxis = "", labyaxis = "") {
nsat <- length((1:nrow(INPUT))[crit])
if (nsat < 2) {
a <- matrix(0, ncol = ncol(INPUT), nrow = 1)
}
if (nsat > 1) {
a <- as.matrix(INPUT[crit, ][order(MD_SOL[crit, MD_DIM]), ])
}
graphics::image(x = 1:nrow(a), z = a, y = 1:ncol(a), axes = FALSE,
zlim = c(1, ncat), xlim = c(1, nsat), xlab = NA, ylab = NA, col = colors)
graphics::box()
graphics::axis(side = 1, tck = -.03, labels = NA)
graphics::axis(side = 2, tck = -.03, labels = NA)
graphics::axis(side = 1, lwd = 0, line = -.9, cex.axis = 0.6)
graphics::axis(side = 2, lwd = 0, line = -.6, las = 1, cex.axis = 0.6)
graphics::mtext(side = 1, labxaxis, line = 2, cex = 0.6)
graphics::mtext(side = 2, labyaxis, line = 2.5, cex = 0.6)
}
D <- x@diss
n <- x@n
S <- x@S
col_bn <- colors[1:2]
nr <- floor(sqrt(S))
nc <- S %/% nr
nr <- ifelse(S %% nr, nr + 1, nr)
opar <- graphics::par(c("mfrow", "mar", "oma"))
on.exit(par(opar))
graphics::par(mfrow = c(nr, nc), mar = c(1, .75, .75, .5) + 0.1, oma = c(1, 0, 0, 0))
for(s in 1:S) {
use <- as.matrix(D[[s]]) + 1
plot_dmbc_data(2, col_bn, use[, 1] > 0, use, as.matrix((1:nrow(use))), 1)
mtext(paste0("Subject ", s), side = 3, font = font, line = 0.05, cex = cex.font)
}
} else {
stop("no plot methods implemented yet for non-binary data.")
}
}
)
#' An S4 class to represent a DMBC model.
#'
#' @slot p A length-one character vector representing the number of dimensions
#' of the latent space to use in the MDS analysis.
#' @slot G A length-one numeric vector representing the number of clusters to
#' partition the subjects into.
#' @slot family A length-one character vector representing the type of data to
#' analyze. Currently, it accepts only the 'binomial' value, but future
#' developments will include the possibility to analyze continuous,
#' multinomial and count data.
#'
#' @name dmbc_model-class
#' @rdname dmbc_model-class
#' @aliases dmbc_model
#'
#' @references
#' Venturini, S., Piccarreta, R. (2021), "A Bayesian Approach for Model-Based
#' Clustering of Several Binary Dissimilarity Matrices: the \pkg{dmbc}
#' Package in \code{R}", Journal of Statistical Software, 100, 16, 1--35, <10.18637/jss.v100.i16>.
#'
#' @examples
#' showClass("dmbc_model")
#'
#' @exportClass dmbc_model
setClass(Class = "dmbc_model",
slots = c(
p = "numeric",
G = "numeric",
family = "character"
)
)
#' Create an instance of the \code{dmbc_model} class using new/initialize.
#'
#' @param .Object Prototype object from the class \code{\link{dmbc_model}}.
#' @param p A length-one character vector representing the number of dimensions
#' of the latent space to use in the MDS analysis.
#' @param G A length-one numeric vector representing the number of clusters to
#' partition the subjects into.
#' @param family A length-one character vector representing the type of data to
#' analyze. Currently, it accepts only the 'binomial' value, but future
#' developments will include the possibility to analyze continuous,
#' multinomial and count data.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases initialize,dmbc_model-method
#' @aliases dmbc_model-initialize
#'
#' @importFrom methods initialize
#' @exportMethod initialize
setMethod("initialize", "dmbc_model",
function(
.Object,
p = numeric(),
G = numeric(),
family = character()
)
{
.Object@p <- p
.Object@G <- G
.Object@family <- family
.Object
}
)
#' Show an instance of the \code{dmbc_model} class.
#'
#' @param object An object of class \code{\link{dmbc_model}}.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases show,dmbc_model-method
#' @aliases dmbc_model-show
#'
#' @importFrom methods show
#' @exportMethod show
setMethod("show",
"dmbc_model",
function(object) {
cat("Dissimilarity Model Based Clustering definition\n")
cat("Number of latent dimensions (p):", object@p, "\n")
cat("Number of clusters (G):", object@G, "\n")
cat("Family:", object@family, "\n")
}
)
#' An S4 class to represent the results of fitting DMBC model.
#'
#' @description
#' An S4 class to represent the results of fitting DMBC model using a single
#' Markov Chain Monte Carlo chain.
#'
#' @slot z.chain An object of class \code{array}; posterior draws from
#' the MCMC algorithm for the (untransformed) latent configuration \eqn{Z}.
#' @slot z.chain.p An object of class \code{array}; posterior draws from
#' the MCMC algorithm for the (Procrustes-transformed) latent configuration
#' \eqn{Z}.
#' @slot alpha.chain An object of class \code{matrix}; posterior draws
#' from the MCMC algorithm for the \eqn{\alpha} parameters.
#' @slot eta.chain An object of class \code{matrix}; posterior draws
#' from the MCMC algorithm for the \eqn{\eta} parameters.
#' @slot sigma2.chain An object of class \code{matrix}; posterior draws
#' from the MCMC algorithm for the \eqn{\sigma^2} parameters.
#' @slot lambda.chain An object of class \code{matrix}; posterior draws
#' from the MCMC algorithm for the \eqn{\lambda} parameters.
#' @slot prob.chain An object of class \code{array}; posterior draws
#' from the MCMC algorithm for the cluster membership probabilities.
#' @slot x.ind.chain An object of class \code{array}; posterior draws
#' from the MCMC algorithm for the cluster membership indicators.
#' @slot x.chain An object of class \code{matrix}; posterior draws from
#' the MCMC algorithm for the cluster membership labels.
#' @slot accept An object of class \code{matrix}; final acceptance rates
#' for the MCMC algorithm.
#' @slot diss An object of class \code{list}; list of observed
#' dissimilarity matrices.
#' @slot dens An object of class \code{list}; list of log-likelihood,
#' log-prior and log-posterior values at each iteration of the MCMC simulation.
#' @slot control An object of class \code{list}; list of the control
#' parameters (number of burnin and sample iterations, number of MCMC chains,
#' etc.). See \code{\link{dmbc_control}()} for more information.
#' @slot prior An object of class \code{list}; list of the prior
#' hyperparameters. See \code{\link{dmbc_prior}()} for more information.
#' @slot dim An object of class \code{list}; list of dimensions for
#' the estimated model, i.e. number of objects (\emph{n}), number of latent
#' dimensions (\emph{p}), number of clusters (\emph{G}), and number of
#' subjects (\emph{S}).
#' @slot model An object of class \code{\link{dmbc_model}}.
#'
#' @name dmbc_fit-class
#' @rdname dmbc_fit-class
#'
#' @references
#' Venturini, S., Piccarreta, R. (2021), "A Bayesian Approach for Model-Based
#' Clustering of Several Binary Dissimilarity Matrices: the \pkg{dmbc}
#' Package in \code{R}", Journal of Statistical Software, 100, 16, 1--35, <10.18637/jss.v100.i16>.
#'
#' @examples
#' showClass("dmbc_fit")
#'
#' @exportClass dmbc_fit
setClass(Class = "dmbc_fit",
slots = c(
z.chain = "array",
z.chain.p = "array",
alpha.chain = "matrix",
eta.chain = "matrix",
sigma2.chain = "matrix",
lambda.chain = "matrix",
prob.chain = "array",
x.ind.chain = "array",
x.chain = "matrix",
accept = "matrix",
diss = "list",
dens = "list",
control = "list",
prior = "list",
dim = "list",
model = "dmbc_model"
)
)
#' Create an instance of the \code{dmbc_fit} class using new/initialize.
#'
#' @param .Object Prototype object from the class \code{\link{dmbc_fit}}.
#' @param z.chain An object of class \code{array}; posterior draws from
#' the MCMC algorithm for the (untransformed) latent configuration \eqn{Z}.
#' @param z.chain.p An object of class \code{array}; posterior draws from
#' the MCMC algorithm for the (Procrustes-transformed) latent configuration
#' \eqn{Z}.
#' @param alpha.chain An object of class \code{matrix}; posterior draws
#' from the MCMC algorithm for the \eqn{\alpha} parameters.
#' @param eta.chain An object of class \code{matrix}; posterior draws
#' from the MCMC algorithm for the \eqn{\eta} parameters.
#' @param sigma2.chain An object of class \code{matrix}; posterior draws
#' from the MCMC algorithm for the \eqn{\sigma^2} parameters.
#' @param lambda.chain An object of class \code{matrix}; posterior draws
#' from the MCMC algorithm for the \eqn{\lambda} parameters.
#' @param prob.chain An object of class \code{array}; posterior draws
#' from the MCMC algorithm for the cluster membership probabilities.
#' @param x.ind.chain An object of class \code{array}; posterior draws
#' from the MCMC algorithm for the cluster membership indicators.
#' @param x.chain An object of class \code{matrix}; posterior draws from
#' the MCMC algorithm for the cluster membership labels.
#' @param accept An object of class \code{matrix}; final acceptance rates
#' for the MCMC algorithm.
#' @param diss An object of class \code{list}; list of observed
#' dissimilarity matrices.
#' @param dens An object of class \code{list}; list of log-likelihood,
#' log-prior and log-posterior values at each iteration of the MCMC simulation.
#' @param control An object of class \code{list}; list of the control
#' parameters (number of burnin and sample iterations, number of MCMC chains,
#' etc.). See \code{\link{dmbc_control}()} for more information.
#' @param prior An object of class \code{list}; list of the prior
#' hyperparameters. See \code{\link{dmbc_prior}()} for more information.
#' @param dim An object of class \code{list}; list of dimensions for
#' the estimated model, i.e. number of objects (\emph{n}), number of latent
#' dimensions (\emph{p}), number of clusters (\emph{G}), and number of
#' subjects (\emph{S}).
#' @param model An object of class \code{\link{dmbc_model}}.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases initialize,dmbc_fit-method
#' @aliases dmbc_fit-initialize
#'
#' @importFrom methods initialize
#' @exportMethod initialize
setMethod("initialize",
"dmbc_fit",
function(
.Object,
z.chain = array(),
z.chain.p = array(),
alpha.chain = matrix(),
eta.chain = matrix(),
sigma2.chain = matrix(),
lambda.chain = matrix(),
prob.chain = array(),
x.ind.chain = array(),
x.chain = matrix(),
accept = matrix(),
diss = list(),
dens = list(),
control = list(),
prior = list(),
dim = list(),
model = NA
)
{
.Object@z.chain <- z.chain
.Object@z.chain.p <- z.chain.p
.Object@alpha.chain <- alpha.chain
.Object@eta.chain <- eta.chain
.Object@sigma2.chain <- sigma2.chain
.Object@lambda.chain <- lambda.chain
.Object@prob.chain <- prob.chain
.Object@x.ind.chain <- x.ind.chain
.Object@x.chain <- x.chain
.Object@accept <- accept
.Object@diss <- diss
.Object@dens <- dens
.Object@control <- control
.Object@prior <- prior
.Object@dim <- dim
.Object@model <- model
.Object
}
)
#' Show an instance of the \code{dmbc_fit} class.
#'
#' @param object An object of class \code{\link{dmbc_fit}}.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases show,dmbc_fit-method
#' @aliases dmbc_fit-show
#'
#' @importFrom methods show
#' @exportMethod show
setMethod("show",
"dmbc_fit",
function(object) {
cat("Dissimilarity Model Based Clustering simulated chain\n")
cat("Number of latent dimensions (p):", object@model@p, "\n")
cat("Number of clusters (G):", object@model@G, "\n")
cat("Family:", object@model@family, "\n")
cat("\n")
cat("To get a summary of the object, use the 'summary()' function.")
}
)
#' Provide a summary of a \code{dmbc_fit} class instance.
#'
#' @param object An object of class \code{\link{dmbc_fit}}.
#' @param include.burnin A length-one logical vector. If \code{TRUE} the
#' burnin iterations (if available) are included in the summary.
#' @param summary.Z A length-one logical vector. If \code{TRUE} the summary
#' also includes the latent configuration coordinates.
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases summary,dmbc_fit-method
#' @aliases dmbc_fit-summary
#'
#' @exportMethod summary
setMethod("summary",
"dmbc_fit",
function(object, include.burnin = FALSE, summary.Z = FALSE, ...) {
control <- object@control
n <- object@dim[["n"]]
p <- object@dim[["p"]]
G <- object@dim[["G"]]
S <- object@dim[["S"]]
res.coda <- dmbc_fit_to_mcmc(object, include.burnin = include.burnin, verbose = FALSE)
if (summary.Z) {
res.coda <- res.coda[, (n*p*G + 1):(2*n*p*G + 4*G)]
} else {
res.coda <- res.coda[, (2*n*p*G + 1):(2*n*p*G + 4*G)]
}
out <- summary(res.coda)
return(out)
}
)
#' @export
setGeneric("subset", function(x) standardGeneric("subset"))
#' Subsetting a \code{dmbc_fit} object.
#'
#' @param x An object of class \code{\link{dmbc_fit}}.
#' @param pars An optional character vector of parameter names. If neither
#' \code{pars} nor \code{regex_pars} is specified, the default is to use all
#' parameters.
#' @param regex_pars An optional \code{\link[=grep]{regular expression}} to use
#' for parameter selection. Can be specified instead of \code{pars} or in addition
#' to \code{pars}.
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases subset,dmbc_fit-method
#' @aliases dmbc_fit-subset
#'
#' @export
setMethod("subset",
"dmbc_fit",
function(x, pars = character(), regex_pars = character(), ...) {
x_mcmc <- dmbc_fit_to_mcmc(x, include.burnin = TRUE, verbose = FALSE)
parnames <- colnames(x_mcmc)
pars <- select_pars(explicit = pars, patterns = regex_pars, complete = parnames)
out <- x_mcmc[, pars]
return(out)
}
)
#' Provide a graphical summary of a \code{dmbc_fit} class instance.
#'
#' @param x An object of class \code{\link{dmbc_fit}}.
#' @param what A length-one character vector providing the plot type to produce.
#' Admissible values are those provided by the \pkg{\link{bayesplot}} package,
#' that is: \code{acf}, \code{areas}, \code{dens}, \code{hex}, \code{hist},
#' \code{intervals}, \code{neff}, \code{pairs}, \code{parcoord}, \code{recover},
#' \code{rhat}, \code{scatter}, \code{trace}, \code{violin} or \code{combo}.
#' In particular, \code{combo} allows to mix different plot types. For more
#' details see the documentation of the \pkg{\link{bayesplot}} package,
#' starting from \code{\link[=MCMC-overview]{this overview page}}.
#' @param pars An optional character vector of parameter names. If neither
#' \code{pars} nor \code{regex_pars} is specified, the default is to use all parameters.
#' @param regex_pars An optional \code{\link[=grep]{regular expression}} to use for
#' parameter selection. Can be specified instead of \code{pars} or in addition to
#' \code{pars}.
#' @param include.burnin A length-one logical vector. If \code{TRUE} the
#' burnin iterations (if available) are included in the summary.
#' @param combo A character vector providing the plot types to combine (see
#' \code{\link[bayesplot]{mcmc_combo}}).
#' @param ... Further arguments to pass on.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases plot,dmbc_fit-method
#' @aliases dmbc_fit-plot
#'
#' @exportMethod plot
setMethod("plot",
signature(x = "dmbc_fit"),
function(x, what = "trace", pars = character(), regex_pars = "lambda", include.burnin = FALSE,
combo = NULL, ...) {
stopifnot(is.character(pars),
is.character(regex_pars),
is.character(what))
if (!(what %in% unlist(all_plots_list, use.names = FALSE)))
stop("the plot type specified is not available.")
x_mcmc <- dmbc_fit_to_mcmc(x, include.burnin = include.burnin, verbose = FALSE)
control <- x@control
if (what %in% acf_plot_list) {
if (what == "acf") {
p <- bayesplot::mcmc_acf(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "acf_bar") {
p <- bayesplot::mcmc_acf_bar(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% areas_plot_list) {
if (what == "areas") {
p <- bayesplot::mcmc_areas(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "areas_ridges") {
p <- bayesplot::mcmc_areas_ridges(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% dens_plot_list) {
if (what == "dens") {
p <- bayesplot::mcmc_dens(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "dens_overlay") {
p <- bayesplot::mcmc_dens_overlay(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "dens_chains") {
p <- bayesplot::mcmc_dens_chains(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% hex_plot_list) {
if (what == "hex") {
p <- bayesplot::mcmc_hex(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% hist_plot_list) {
if (what == "hist") {
p <- bayesplot::mcmc_hist(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "hist_by_chain") {
p <- bayesplot::mcmc_hist_by_chain(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% intervals_plot_list) {
if (what == "intervals") {
p <- bayesplot::mcmc_intervals(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% neff_plot_list) {
x_sub <- subset(x, pars = pars, regex_pars = regex_pars)
nsample <- floor((control[["burnin"]] + control[["nsim"]])/control[["thin"]])
neff <- coda::effectiveSize(x_sub)
ratio <- neff/nsample
if (what == "neff") {
p <- bayesplot::mcmc_neff(ratio = ratio, ...)
} else if (what == "neff_hist") {
p <- bayesplot::mcmc_neff_hist(ratio = ratio, ...)
}
}
if (what %in% pairs_plot_list) {
if (what == "pairs") {
p <- bayesplot::mcmc_pairs(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% parcoord_plot_list) {
if (what == "parcoord") {
p <- bayesplot::mcmc_parcoord(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% recover_plot_list) {
x_sub <- subset(x, pars = pars, regex_pars = regex_pars)
if (what == "recover_hist") {
p <- bayesplot::mcmc_recover_hist(x = x_sub, ...)
} else if (what == "recover_intervals") {
p <- bayesplot::mcmc_recover_intervals(x = x_sub, ...)
} else if (what == "recover_scatter") {
p <- bayesplot::mcmc_recover_scatter(x = x_sub, ...)
}
}
if (what %in% rhat_plot_list) {
x_sub <- subset(x, pars = pars, regex_pars = regex_pars)
rhat <- coda::gelman.diag(x_sub, multivariate = FALSE)$psrf[, 1]
if (what == "rhat") {
p <- bayesplot::mcmc_rhat(rhat = rhat, ...)
} else if (what == "rhat_hist") {
p <- bayesplot::mcmc_rhat_hist(rhat = rhat, ...)
}
}
if (what %in% scatter_plot_list) {
if (what == "scatter") {
p <- bayesplot::mcmc_scatter(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% trace_plot_list) {
if (what == "trace") {
p <- bayesplot::mcmc_trace(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "trace_highlight") {
p <- bayesplot::mcmc_trace_highlight(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% violin_plot_list) {
if (what == "violin") {
p <- bayesplot::mcmc_violin(x = x_mcmc, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what == "combo") {
if (!is.null(combo)) {
p <- bayesplot::mcmc_combo(x = x_mcmc, pars = pars, regex_pars = regex_pars, combo = combo, ...)
} else {
stop("to produce an 'mcmc_combo' plot, the 'combo' option must be specified.")
}
}
p
}
)
#' An S4 class to represent the results of fitting DMBC model.
#'
#' @description
#' An S4 class to represent the results of fitting DMBC model using multiple
#' Markov Chain Monte Carlo chains.
#'
#' @slot results An object of class \code{list}; list of \code{dmbc_fit}
#' objects corresponding to the parallel MCMC chains simulated during the
#' estimation.
#'
#' @name dmbc_fit_list-class
#' @rdname dmbc_fit_list-class
#' @aliases dmbc_fit_list
#'
#' @seealso
#' \code{\link{dmbc_fit}} for more details on the components of each element of
#' the list.
#'
#' @references
#' Venturini, S., Piccarreta, R. (2021), "A Bayesian Approach for Model-Based
#' Clustering of Several Binary Dissimilarity Matrices: the \pkg{dmbc}
#' Package in \code{R}", Journal of Statistical Software, 100, 16, 1--35, <10.18637/jss.v100.i16>.
#'
#' @examples
#' showClass("dmbc_fit_list")
#'
#' @exportClass dmbc_fit_list
setClass(Class = "dmbc_fit_list",
slots = c(
results = "list"
)
)
#' Create an instance of the \code{dmbc_fit_list} class using new/initialize.
#'
#' @param .Object Prototype object from the class \code{\link{dmbc_fit_list}}.
#' @param results A list whose elements are the \code{dmbc_fit} objects for
#' each simulated chain.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases initialize,dmbc_fit_list-method
#' @aliases dmbc_fit_list-initialize
#'
#' @importFrom methods initialize
#' @exportMethod initialize
setMethod("initialize", "dmbc_fit_list",
function(
.Object,
results = list()
)
{
.Object@results <- results
.Object
}
)
#' Show an instance of the \code{dmbc_fit_list} class.
#'
#' @param object An object of class \code{\link{dmbc_fit_list}}.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases show,dmbc_fit_list-method
#' @aliases dmbc_fit_list-show
#'
#' @importFrom methods show
#' @exportMethod show
setMethod("show",
"dmbc_fit_list",
function(object) {
cat("List of Dissimilarity Model Based Clustering simulated chains\n")
cat("Number of simulated chains:", length(object@results), "\n")
cat("Number of latent dimensions (p):", object@results[[1]]@model@p, "\n")
cat("Number of clusters (G):", object@results[[1]]@model@G, "\n")
cat("Family:", object@results[[1]]@model@family, "\n")
cat("\n")
cat("To get a summary of the object, use the 'summary()' function.")
}
)
#' Provide a summary of a \code{dmbc_fit_list} class instance.
#'
#' @param object An object of class \code{\link{dmbc_fit_list}}.
#' @param include.burnin A length-one logical vector. If \code{TRUE} the
#' burnin iterations (if available) are included in the summary.
#' @param summary.Z A length-one logical vector. If \code{TRUE} the summary
#' also includes the latent configuration coordinates.
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases summary,dmbc_fit_list-method
#' @aliases dmbc_fit_list-summary
#'
#' @exportMethod summary
setMethod("summary",
"dmbc_fit_list",
function(object, include.burnin = FALSE, summary.Z = FALSE, ...) {
control <- object@results[[1]]@control
nchains <- control[["nchains"]]
n <- object@results[[1]]@dim[["n"]]
p <- object@results[[1]]@dim[["p"]]
G <- object@results[[1]]@dim[["G"]]
S <- object@results[[1]]@dim[["S"]]
res.coda <- dmbc_fit_list_to_mcmc.list(object, include.burnin = include.burnin, verbose = FALSE)
for (c in 1:nchains) {
if (summary.Z) {
res.coda[[c]] <- res.coda[[c]][, (n*p*G + 1):(2*n*p*G + 4*G)]
} else {
res.coda[[c]] <- res.coda[[c]][, (2*n*p*G + 1):(2*n*p*G + 4*G)]
}
}
out <- summary(res.coda)
return(out)
}
)
#' Subsetting a \code{dmbc_fit_list} object.
#'
#' @param x An object of class \code{\link{dmbc_fit_list}}.
#' @param pars An optional character vector of parameter names. If neither
#' \code{pars} nor \code{regex_pars} is specified, the default is to use all
#' parameters.
#' @param regex_pars An optional \code{\link[=grep]{regular expression}} to use
#' for parameter selection. Can be specified instead of \code{pars} or in addition
#' to \code{pars}.
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases subset,dmbc_fit_list-method
#' @aliases dmbc_fit_list-subset
#'
#' @export
setMethod("subset",
"dmbc_fit_list",
function(x, pars = character(), regex_pars = character(), ...) {
x_mcmc.list <- dmbc_fit_list_to_mcmc.list(x, include.burnin = TRUE, verbose = FALSE)
parnames <- colnames(x_mcmc.list[[1]])
pars <- select_pars(explicit = pars, patterns = regex_pars, complete = parnames)
out <- coda::mcmc.list(lapply(x_mcmc.list, function(chain) chain[, pars]))
return(out)
}
)
#' Provide a graphical summary of a \code{dmbc_fit_list} class instance.
#'
#' @param x An object of class \code{\link{dmbc_fit_list}}.
#' @param what A length-one character vector providing the plot type to produce.
#' Admissible values are those provided by the \pkg{\link{bayesplot}} package,
#' that is: \code{acf}, \code{areas}, \code{dens}, \code{hex}, \code{hist},
#' \code{intervals}, \code{neff}, \code{pairs}, \code{parcoord}, \code{recover},
#' \code{rhat}, \code{scatter}, \code{trace}, \code{violin} or \code{combo}.
#' In particular, \code{combo} allows to mix different plot types. For more
#' details see the documentation of the \pkg{\link{bayesplot}} package,
#' starting from \code{\link[=MCMC-overview]{this overview page}}.
#' @param pars An optional character vector of parameter names. If neither
#' \code{pars} nor \code{regex_pars} is specified, the default is to use all parameters.
#' @param regex_pars An optional \code{\link[=grep]{regular expression}} to use for
#' parameter selection. Can be specified instead of \code{pars} or in addition to
#' \code{pars}.
#' @param include.burnin A length-one logical vector. If \code{TRUE} the
#' burnin iterations (if available) are included in the summary.
#' @param combo A character vector providing the plot types to combine (see
#' \code{\link[bayesplot]{mcmc_combo}}).
#' @param ... Further arguments to pass on.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases plot,dmbc_fit_list-method
#' @aliases dmbc_fit_list-plot
#'
#' @exportMethod plot
setMethod("plot",
signature(x = "dmbc_fit_list"),
function(x, what = "trace", pars = character(), regex_pars = "lambda", include.burnin = FALSE,
combo = NULL, ...) {
stopifnot(is.character(pars),
is.character(regex_pars),
is.character(what))
if (!(what %in% unlist(all_plots_list, use.names = FALSE)))
stop("the plot type specified is not available.")
x_mcmc.list <- dmbc_fit_list_to_mcmc.list(x, include.burnin = include.burnin, verbose = FALSE)
control <- x@results[[1]]@control
if (what %in% acf_plot_list) {
if (what == "acf") {
p <- bayesplot::mcmc_acf(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "acf_bar") {
p <- bayesplot::mcmc_acf_bar(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% areas_plot_list) {
if (what == "areas") {
p <- bayesplot::mcmc_areas(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "areas_ridges") {
p <- bayesplot::mcmc_areas_ridges(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% dens_plot_list) {
if (what == "dens") {
p <- bayesplot::mcmc_dens(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "dens_overlay") {
p <- bayesplot::mcmc_dens_overlay(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "dens_chains") {
p <- bayesplot::mcmc_dens_chains(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% hex_plot_list) {
if (what == "hex") {
p <- bayesplot::mcmc_hex(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% hist_plot_list) {
if (what == "hist") {
p <- bayesplot::mcmc_hist(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "hist_by_chain") {
p <- bayesplot::mcmc_hist_by_chain(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% intervals_plot_list) {
if (what == "intervals") {
p <- bayesplot::mcmc_intervals(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% neff_plot_list) {
x_sub <- subset(x, pars = pars, regex_pars = regex_pars)
nsample <- control[["nchains"]]*floor((control[["burnin"]] + control[["nsim"]])/control[["thin"]])
neff <- coda::effectiveSize(x_sub)
ratio <- neff/nsample
if (what == "neff") {
p <- bayesplot::mcmc_neff(ratio = ratio, ...)
} else if (what == "neff_hist") {
p <- bayesplot::mcmc_neff_hist(ratio = ratio, ...)
}
}
if (what %in% pairs_plot_list) {
if (what == "pairs") {
p <- bayesplot::mcmc_pairs(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% parcoord_plot_list) {
if (what == "parcoord") {
p <- bayesplot::mcmc_parcoord(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% recover_plot_list) {
x_sub <- subset(x, pars = pars, regex_pars = regex_pars)
if (what == "recover_hist") {
p <- bayesplot::mcmc_recover_hist(x = x_sub, ...)
} else if (what == "recover_intervals") {
p <- bayesplot::mcmc_recover_intervals(x = x_sub, ...)
} else if (what == "recover_scatter") {
p <- bayesplot::mcmc_recover_scatter(x = x_sub, ...)
}
}
if (what %in% rhat_plot_list) {
x_sub <- subset(x, pars = pars, regex_pars = regex_pars)
rhat <- coda::gelman.diag(x_sub, multivariate = FALSE)$psrf[, 1]
if (what == "rhat") {
p <- bayesplot::mcmc_rhat(rhat = rhat, ...)
} else if (what == "rhat_hist") {
p <- bayesplot::mcmc_rhat_hist(rhat = rhat, ...)
}
}
if (what %in% scatter_plot_list) {
if (what == "scatter") {
p <- bayesplot::mcmc_scatter(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% trace_plot_list) {
if (what == "trace") {
p <- bayesplot::mcmc_trace(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
} else if (what == "trace_highlight") {
p <- bayesplot::mcmc_trace_highlight(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what %in% violin_plot_list) {
if (what == "violin") {
p <- bayesplot::mcmc_violin(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, ...)
}
}
if (what == "combo") {
if (!is.null(combo)) {
p <- bayesplot::mcmc_combo(x = x_mcmc.list, pars = pars, regex_pars = regex_pars, combo = combo, ...)
} else {
stop("to produce an 'mcmc_combo' plot, the 'combo' option must be specified.")
}
}
p
}
)
#' An S4 class to represent the comparison of a set of DMBC models.
#'
#' @description
#' An S4 class to represent the comparison of a set of DMBC models through
#' the dissimilarity model-based clustering information criterion (DCIC).
#'
#' @slot logprior An object of class \code{matrix} providing the
#' log-prior values corresponding to different values of \emph{p} and
#' \emph{G}.
#' @slot logmlik An object of class \code{matrix} providing the
#' marginal log-likelihood values corresponding to different values of
#' \emph{p} and \emph{G}.
#' @slot logcorrfact An object of class \code{matrix} providing the
#' logarithm of the correction factors corresponding to different values of
#' \emph{p} and \emph{G}.
#' @slot DCIC An object of class \code{matrix} providing the values
#' of the DCIC index corresponding to different values of \emph{p} and
#' \emph{G}.
#' @slot post.est An object of class \code{list}; named list with
#' elements representing the parameter estimates corresponding to different
#' values of \emph{p} and \emph{G}.
#' @slot est A length-one character vector representing the estimate
#' type used in computing the DCIC index. Possible values are \code{mean},
#' \code{median}, \code{ml} and \code{map}. See \code{\link{dmbc_ic}()} for
#' more details about these values.
#' @slot res_last_p An object of class \code{list}; list of
#' \code{dmbc_fit_list} objects with the results of fitting the DMBC
#' models corresponding to the last value of \emph{p}. This is needed in case
#' of an update of the DCIC calculations using additional \emph{p} and/or
#' \emph{G} values.
#'
#' @name dmbc_ic-class
#' @rdname dmbc_ic-class
#' @aliases dmbc_ic
#'
#' @references
#' Venturini, S., Piccarreta, R. (2021), "A Bayesian Approach for Model-Based
#' Clustering of Several Binary Dissimilarity Matrices: the \pkg{dmbc}
#' Package in \code{R}", Journal of Statistical Software, 100, 16, 1--35, <10.18637/jss.v100.i16>.
#'
#' @examples
#' showClass("dmbc_ic")
#'
#' @exportClass dmbc_ic
setClass(Class = "dmbc_ic",
slots = c(
logprior = "matrix",
logmlik = "matrix",
logcorrfact = "matrix",
DCIC = "matrix",
post.est = "list",
est = "character",
res_last_p = "list"
)
)
#' Create an instance of the \code{dmbc_ic} class using new/initialize.
#'
#' @param .Object Prototype object from the class \code{\link{dmbc_ic}}.
#' @param logprior An object of class \code{matrix} providing the
#' log-prior values corresponding to different values of \emph{p} and
#' \emph{G}.
#' @param logmlik An object of class \code{matrix} providing the
#' marginal log-likelihood values corresponding to different values of
#' \emph{p} and \emph{G}.
#' @param logcorrfact An object of class \code{matrix} providing the
#' logarithm of the correction factors corresponding to different values of
#' \emph{p} and \emph{G}.
#' @param DCIC An object of class \code{matrix} providing the values
#' of the DCIC index corresponding to different values of \emph{p} and
#' \emph{G}.
#' @param post.est An object of class \code{list}; named list with
#' elements representing the parameter estimates corresponding to different
#' values of \emph{p} and \emph{G}.
#' @param est A length-one character vector representing the estimate
#' type used in computing the DCIC index. Possible values are \code{mean},
#' \code{median}, \code{ml} and \code{map}. See \code{\link{dmbc_ic}()} for
#' more details about these values.
#' @param res_last_p An object of class \code{list}; list of
#' \code{dmbc_fit_list} objects with the results of fitting the DMBC
#' models corresponding to the last value of \emph{p}. This is needed in case
#' of an update of the DCIC calculations using additional \emph{p} and/or
#' \emph{G} values.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases initialize,dmbc_ic-method
#' @aliases dmbc_ic-initialize
#'
#' @importFrom methods initialize
#' @exportMethod initialize
setMethod("initialize", "dmbc_ic",
function(
.Object,
logprior = matrix(),
logmlik = matrix(),
logcorrfact = matrix(),
DCIC = matrix(),
post.est = list(),
est = character(),
res_last_p = list()
)
{
.Object@logprior <- logprior
.Object@logmlik <- logmlik
.Object@logcorrfact <- logcorrfact
.Object@DCIC <- DCIC
.Object@post.est <- post.est
.Object@est <- est
.Object@res_last_p <- res_last_p
.Object
}
)
#' Show an instance of the \code{dmbc_ic} class.
#'
#' @param object An object of class \code{\link{dmbc_ic}}.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases show,dmbc_ic-method
#' @aliases dmbc_ic-show
#'
#' @importFrom methods show
#' @exportMethod show
setMethod("show",
"dmbc_ic",
function(object) {
pmax <- nrow(object@DCIC)
Gmax <- ncol(object@DCIC)
last <- length(object@res_last_p)
family <- object@res_last_p[[1]]@results[[1]]@model@family
cat("Dissimilarity Model Based Clustering -- Model choice\n")
cat("Latent dimensions (p) requested: from 1 to", pmax, "\n")
cat("Number of clusters (G) requested: from 1 to", Gmax, "\n")
cat("Family:", family, "\n")
cat("Dissimilarity Model Based Clustering Information Criterion (DCIC):\n")
dcic_rownm <- paste0("p = ", 1:pmax)
dcic_colnm <- paste0("G = ", 1:Gmax)
print_matrix(object@DCIC, dcic_rownm, dcic_colnm, colwidth = 8)
}
)
#' Provide a summary of a \code{dmbc_ic} class instance.
#'
#' @param object An object of class \code{\link{dmbc_ic}}.
#' @param p An optional length-one numeric vector providing the number of
#' latent space dimension to use in the summary.
#' @param G An optional length-one numeric vector providing the number of
#' clusters to use in the summary.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases summary,dmbc_ic-method
#' @aliases dmbc_ic-summary
#'
#' @exportMethod summary
setMethod("summary",
"dmbc_ic",
function(object, p = NULL, G = NULL) {
pmax <- nrow(object@DCIC)
Gmax <- ncol(object@DCIC)
dcic_rownm <- paste0("p = ", 1:pmax)
dcic_colnm <- paste0("G = ", 1:Gmax)
if (!is.null(p))
if (p > pmax) stop("the number of dimensions (p) specified is not available.")
if (!is.null(G))
if (G > Gmax) stop("the number of clusters (G) specified is not available.")
show(object)
cat("\n")
cat("Dissimilarity Model Based Clustering -- Logarithm of marginal likelihood values:\n")
print_matrix(object@logmlik, dcic_rownm, dcic_colnm)
cat("\n")
cat("Dissimilarity Model Based Clustering -- Logarithm of prior values:\n")
print_matrix(object@logprior, dcic_rownm, dcic_colnm)
cat("\n")
cat("Dissimilarity Model Based Clustering -- Logarithm of correction factors:\n")
print_matrix(object@logcorrfact, dcic_rownm, dcic_colnm)
if (is.null(p) || is.null(G)) {
G_best <- ifelse((which.min(t(object@DCIC)) %% Gmax) == 0, Gmax, (which.min(t(object@DCIC)) %% Gmax))
p_best <- which.min(object@DCIC[, G_best])
} else {
G_best <- G
p_best <- p
}
best_lab <- paste("p = ", p_best, " -- G = ", G_best, sep = "")
best_idx <- which(names(object@post.est) == best_lab)
best_res <- object@post.est[[best_idx]]
n <- dim(best_res$z.m)[1]
cat("\n")
est <- switch(object@est,
mean = "posterior mean",
median = "posterior median",
ml = "maximum likelihood",
map = "maximum a posteriori"
)
cat("Dissimilarity Model Based Clustering -- Parameter estimates (", est, ") -- p = ", p_best,
" - G = ", G_best, "\n", sep = "")
wd <- as.integer(options("width"))
for (g in 1:G_best) {
dash_string_L <- strrep("-", floor((wd - nchar(paste0("Cluster ", Gmax)) - 2)/2))
dash_string_R <- strrep("-", floor((wd - nchar(paste0("Cluster ", Gmax)) - 2)/2) - nchar(g) + 1)
cat(dash_string_L, " ", "Cluster ", g, " ", dash_string_R, "\n", sep = "")
cat("(1) Z:\n")
print_matrix(as.matrix(best_res$z.m[, , g]), paste0("n = ", 1:n), paste0("Dimension ", 1:p_best),
colwidth = 11, ndigits = 4, shift = 7)
cat("(2) alpha:", round(best_res$alpha.m[g], digits = 4), "\n")
cat("(3) eta:", round(best_res$eta.m[g], digits = 4), "\n")
cat("(4) sigma2:", round(best_res$sigma2.m[g], digits = 4), "\n")
cat("(5) lambda:", round(best_res$lambda.m[g], digits = 4), "\n")
}
}
)
#' Provide a graphical summary of a \code{dmbc_ic} class instance.
#'
#' @param x An object of class \code{\link{dmbc_ic}}.
#' @param size A length-two numeric vector providing the optional sizes of
#' points and lines in the plot.
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases plot,dmbc_ic-method
#' @aliases dmbc_ic-plot
#'
#' @exportMethod plot
setMethod("plot",
signature(x = "dmbc_ic"),
function(x, size = NULL, ...) {
data <- reshape2::melt(data = x@DCIC,
varnames = c("p", "G"),
value.name = "DCIC")
data$p <- factor(data$p)
n_p <- nlevels(data$p)
geom_args <- list()
if (is.null(size)) {
geom_args$size <- c(3, 1)
} else {
if (length(size) != 2)
stop("the 'size' argument must be a numeric vector with two elements.")
geom_args$size <- size
}
# mapping <- ggplot2::aes_(x = ~ G, y = ~ DCIC, color = ~ p, group = ~ p, shape = ~ p)
mapping <- ggplot2::aes(x = G, y = DCIC, color = p, group = p, shape = p)
graph <- ggplot2::ggplot(data, mapping) + scale_shape_manual(values = 1:n_p) +
bayesplot::bayesplot_theme_get()
graph <- graph + ggplot2::geom_point(size = geom_args$size[1])
if (n_p > 1) {
graph <- graph + ggplot2::geom_line(size = geom_args$size[2])
}
graph <- graph +
ggplot2::scale_color_manual("p", values = choose_colors(n_p))
graph <- graph +
ggplot2::scale_x_continuous(breaks = 1:ncol(x@DCIC)) +
bayesplot::legend_move(ifelse(n_p > 1, "right", "none")) +
bayesplot::xaxis_title(on = TRUE) +
bayesplot::yaxis_title(on = TRUE)
graph
}
)
#' Provide an update of a \code{dmbc_ic} class instance.
#'
#' @param object An object of class \code{\link{dmbc_ic}}.
#' @param pmax A length-one numeric vector indicating the maximum number of
#' dimensions of the latent space to consider.
#' @param Gmax A length-one numeric vector indicating the maximum number of
#' cluster to consider.
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases update,dmbc_ic-method
#' @aliases dmbc_ic-update
#'
#' @return A \code{dmbc_ic} object.
#'
#' @seealso \code{\link{dmbc}()} for fitting a DMBC model.
#' @seealso \code{\link{dmbc_ic}} for a description of the elements included
#' in the returned object.
#'
#' @references
#' Venturini, S., Piccarreta, R. (2021), "A Bayesian Approach for Model-Based
#' Clustering of Several Binary Dissimilarity Matrices: the \pkg{dmbc}
#' Package in \code{R}", Journal of Statistical Software, 100, 16, 1--35, <10.18637/jss.v100.i16>.
#'
#' @examples
#' \dontrun{
#' data(simdiss, package = "dmbc")
#'
#' pmax <- 2
#' Gmax <- 2
#' prm.prop <- list(z = 1.5, alpha = .75)
#' burnin <- 2000
#' nsim <- 1000
#' seed <- 1809
#'
#' set.seed(seed)
#'
#' control <- list(burnin = burnin, nsim = nsim, z.prop = prm.prop[["z"]],
#' alpha.prop = prm.prop[["alpha"]], random.start = TRUE, verbose = TRUE,
#' thin = 10, store.burnin = TRUE)
#' sim.ic <- dmbc_IC(data = simdiss, pmax = pmax, Gmax = Gmax, control = control,
#' est = "mean")
#'
#' pmax <- pmax + 1
#' Gmax <- Gmax + 2
#' new.ic <- update(sim.ic, pmax = pmax, Gmax = Gmax)
#' new.ic
#'
#' # plot the results
#' library(bayesplot)
#' library(ggplot2)
#' color_scheme_set("mix-yellow-blue")
#' p <- plot(new.ic, size = c(4, 1.5))
#' p + panel_bg(fill = "gray90", color = NA)
#' }
#'
#' @importFrom methods new
#' @importFrom stats4 update
#' @exportMethod update
setMethod("update", "dmbc_ic",
function(object, pmax = NULL, Gmax = NULL, ...) {
control <- object@res_last_p[[1]]@results[[1]]@control
verbose <- control[["verbose"]]
burnin <- control[["burnin"]]
nsim <- control[["nsim"]]
prior <- object@res_last_p[[1]]@results[[1]]@prior
pmax.old <- nrow(object@DCIC)
Gmax.old <- ncol(object@DCIC)
if (is.null(pmax) & is.null(Gmax))
stop("pmax and Gmax cannot be both null.")
if (is.null(pmax)) {
warning("you did not provide a pmax value, so it is set to its previous value.",
call. = FALSE, immediate. = TRUE)
pmax <- pmax.old
}
if (is.null(Gmax)) {
warning("you did not provide a Gmax value, so it is set to its previous value.",
call. = FALSE, immediate. = TRUE)
Gmax <- Gmax.old
}
if ((pmax <= pmax.old) & (Gmax <= Gmax.old))
stop("at least one of the new values for pmax and Gmax must be larger than the previous ones.")
if (pmax < pmax.old) {
warning("the pmax value (", pmax, ") is smaller than the previous one (", pmax.old, ") and hence it is set at that value.", call. = FALSE, immediate. = TRUE)
pmax <- pmax.old
}
if (Gmax < Gmax.old) {
warning("the Gmax value (", Gmax, ") is smaller than the previous one (", Gmax.old, ") and hence it is set at that value.", call. = FALSE, immediate. = TRUE)
Gmax <- Gmax.old
}
D <- object@res_last_p[[1]]@results[[1]]@diss
data <- new("dmbc_data", diss = D, n = object@res_last_p[[1]]@results[[1]]@dim[["n"]],
S = object@res_last_p[[1]]@results[[1]]@dim[["S"]])
est <- object@est
logprior <- logmlik <- logcorrfact <- dcic <- matrix(NA, nrow = pmax, ncol = Gmax)
logprior[1:pmax.old, 1:Gmax.old] <- object@logprior
logmlik[1:pmax.old, 1:Gmax.old] <- object@logmlik
logcorrfact[1:pmax.old, 1:Gmax.old] <- object@logcorrfact
dcic[1:pmax.old, 1:Gmax.old] <- object@DCIC
res_list <- list()
res_save <- list()
res_all <- object@post.est
res_last_p <- object@res_last_p
res.i <- pmax.old*Gmax.old + 1
p.seq <- if (pmax == pmax.old) pmax.old else seq(from = (pmax.old + 1), to = pmax, by = 1)
G.seq <- if (Gmax == Gmax.old) Gmax.old else seq(from = (Gmax.old + 1), to = Gmax, by = 1)
if (Gmax > Gmax.old) {
for (G.i in G.seq) {
for (p.i in 1:pmax.old) {
if (verbose)
message("--- p = ", p.i, " -- G = ", G.i, " ---")
prior.i <- update_prior(prior, p.i, G.i)
res <- dmbc(data = data, p = p.i, G = G.i, control = control, prior = prior.i)
res_list[[p.i]] <- res
if (p.i == pmax.old)
res_last_p[[G.i]] <- res
if (est == "mean") {
est.tmp <- dmbc_get_postmean(res)
} else if (est == "median") {
est.tmp <- dmbc_get_postmedian(res)
} else if (est == "ml") {
est.tmp <- dmbc_get_ml(res)
} else if (est == "map") {
est.tmp <- dmbc_get_map(res)
}
z.m <- est.tmp$z
alpha.m <- est.tmp$alpha
eta.m <- est.tmp$eta
sigma2.m <- est.tmp$sigma2
lambda.m <- est.tmp$lambda
class.m <- dmbc_get_postmean(res)$cluster
res_all[[res.i]] <- list(z.m = z.m, alpha.m = alpha.m, eta.m = eta.m, sigma2.m = sigma2.m,
lambda.m = lambda.m, class.m = class.m)
names(res_all)[res.i] <- paste("p = ", p.i, " -- G = ", G.i, sep = "")
res.i <- res.i + 1
logprior[p.i, G.i] <- log_marg_prior(res, z.m)
logmlik[p.i, G.i] <- log_marg_lik(res, z.m)
if (p.i > 1)
logcorrfact[p.i, G.i] <- log_corr_fact(res_list[[p.i - 1]], z.m)
dcic[p.i, G.i] <- -2*(logprior[p.i, G.i] + logmlik[p.i, G.i] +
ifelse(p.i > 1, sum(logcorrfact[2:p.i, G.i], na.rm = TRUE), 0))
if (verbose) {
message(" ")
}
}
res_list <- list()
}
}
if (pmax > pmax.old) {
for (G.i in 1:Gmax) {
for (p.i in p.seq) {
if (verbose)
message("--- p = ", p.i, " -- G = ", G.i, " ---")
prior.i <- update_prior(prior, p.i, G.i)
res <- dmbc(data = data, p = p.i, G = G.i, control = control, prior = prior.i)
res_list[[p.i]] <- res
if (p.i == p.seq[1])
res_list[[p.i - 1]] <- res_last_p[[G.i]]
if (p.i == pmax)
res_save[[G.i]] <- res
if (est == "mean") {
est.tmp <- dmbc_get_postmean(res)
} else if (est == "median") {
est.tmp <- dmbc_get_postmedian(res)
} else if (est == "ml") {
est.tmp <- dmbc_get_ml(res)
} else if (est == "map") {
est.tmp <- dmbc_get_map(res)
}
z.m <- est.tmp$z
alpha.m <- est.tmp$alpha
eta.m <- est.tmp$eta
sigma2.m <- est.tmp$sigma2
lambda.m <- est.tmp$lambda
class.m <- dmbc_get_postmean(res)$cluster
res_all[[res.i]] <- list(z.m = z.m, alpha.m = alpha.m, eta.m = eta.m, sigma2.m = sigma2.m,
lambda.m = lambda.m, class.m = class.m)
names(res_all)[res.i] <- paste("p = ", p.i, " -- G = ", G.i, sep = "")
res.i <- res.i + 1
logprior[p.i, G.i] <- log_marg_prior(res, z.m)
logmlik[p.i, G.i] <- log_marg_lik(res, z.m)
logcorrfact[p.i, G.i] <- log_corr_fact(res_list[[p.i - 1]], z.m)
dcic[p.i, G.i] <- -2*(logprior[p.i, G.i] + logmlik[p.i, G.i] +
ifelse(p.i > 1, sum(logcorrfact[2:p.i, G.i], na.rm = TRUE), 0))
if (verbose) {
message(" ")
}
}
res_list <- list()
}
}
out <- new("dmbc_ic",
logprior = logprior,
logmlik = logmlik,
logcorrfact = logcorrfact,
DCIC = dcic,
post.est = res_all,
est = est,
res_last_p = res_save
)
return(out)
}
)
#' An S4 class to represent the latent configuration estimate for a DMBC model.
#'
#' @description
#' An S4 class to represent the the latent configuration estimate for a DMBC
#' model.
#'
#' @slot Z.est An array containing the estimate of the latent
#' configuration for a DMBC model.
#' @slot Z.sd An array containing the standard deviation of the latent
#' configuration for a DMBC model.
#' @slot cluster A numeric vector providing the estimated group
#' membership for the \emph{S} subjects in the data.
#' @slot est A length-one character vector providing the estimate type
#' returned in \code{Z.est}. Possible values are \code{mean} (posterior
#' mean), \code{median} (posterior median), \code{ml} (maximum likelihood)
#' and \code{map} (maximum-a-posteriori).
#' @slot n A length-one numeric vector providing the number of objects.
#' @slot p A length-one numeric vector providing the number of latent
#' dimensions.
#' @slot S A length-one numeric vector providing the number of subjects.
#' @slot G A length-one numeric vector providing the number of clusters.
#' @slot family An object of class \code{list}; named list with
#' elements representing the parameter estimates corresponding to different
#' values of \emph{p} and \emph{G}.
#' @slot chain A length-one numeric vector representing the ID of
#' the MCMC chain used to compute the estimates.
#' @slot labels A character vector for the (optional) strings to use
#' in the plots for labeling the objects.
#'
#' @name dmbc_config-class
#' @rdname dmbc_config-class
#' @aliases dmbc_config
#'
#' @references
#' Venturini, S., Piccarreta, R. (2021), "A Bayesian Approach for Model-Based
#' Clustering of Several Binary Dissimilarity Matrices: the \pkg{dmbc}
#' Package in \code{R}", Journal of Statistical Software, 100, 16, 1--35, <10.18637/jss.v100.i16>.
#'
#' @examples
#' showClass("dmbc_config")
#'
#' @exportClass dmbc_config
setClass(Class = "dmbc_config",
slots = c(
Z.est = "array",
Z.sd = "array",
cluster = "numeric",
est = "character",
n = "numeric",
p = "numeric",
S = "numeric",
G = "numeric",
family = "character",
chain = "numeric",
labels = "character"
)
)
#' Create an instance of the \code{dmbc_config} class using new/initialize.
#'
#' @param .Object Prototype object from the class \code{\link{dmbc_config}}.
#' @param Z.est An array containing the estimate of the latent
#' configuration for a DMBC model.
#' @param Z.sd An array containing the standard deviation of the latent
#' configuration for a DMBC model.
#' @param cluster A numeric vector providing the estimated group
#' membership for the \emph{S} subjects in the data.
#' @param est A length-one character vector providing the estimate type
#' returned in \code{Z.est}. Possible values are \code{mean} (posterior
#' mean), \code{median} (posterior median), \code{ml} (maximum likelihood)
#' and \code{map} (maximum-a-posteriori).
#' @param n A length-one numeric vector providing the number of objects.
#' @param p A length-one numeric vector providing the number of latent
#' dimensions.
#' @param S A length-one numeric vector providing the number of subjects.
#' @param G A length-one numeric vector providing the number of clusters.
#' @param family An object of class \code{list}; named list with
#' elements representing the parameter estimates corresponding to different
#' values of \emph{p} and \emph{G}.
#' @param chain A length-one numeric vector representing the ID of
#' the MCMC chain used to compute the estimates.
#' @param labels A character vector for the (optional) strings to use
#' in the plots for labeling the objects.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases initialize,dmbc_config-method
#' @aliases dmbc_config-initialize
#'
#' @importFrom methods initialize
#' @exportMethod initialize
setMethod("initialize", "dmbc_config",
function(
.Object,
Z.est = array(),
Z.sd = array(),
cluster = numeric(),
est = character(),
n = numeric(),
S = numeric(),
p = numeric(),
G = numeric(),
family = character(),
chain = numeric(),
labels = character()
)
{
.Object@Z.est <- Z.est
.Object@Z.sd <- Z.sd
.Object@cluster <- cluster
.Object@est <- est
.Object@n <- n
.Object@S <- S
.Object@p <- p
.Object@G <- G
.Object@family <- family
.Object@chain <- chain
.Object@labels <- labels
.Object
}
)
#' Show an instance of the \code{dmbc_config} class.
#'
#' @param object An object of class \code{\link{dmbc_config}}.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases show,dmbc_config-method
#' @aliases dmbc_config-show
#'
#' @importFrom methods show
#' @exportMethod show
setMethod("show",
"dmbc_config",
function(object) {
est <- switch(object@est,
mean = "posterior mean",
median = "posterior median",
ml = "maximum likelihood",
map = "maximum a posteriori"
)
cat("Dissimilarity Model Based Clustering -- Latent configuration estimates\n")
cat("Number of objects (n):", object@n, "\n")
cat("Number of latent dimensions (p):", object@p, "\n")
cat("Number of subjects (S):", object@S, "\n")
cat("Number of clusters (G):", object@G, "\n")
cat("Family:", object@family, "\n")
cat("Chain ID:", object@chain, "\n")
cat("Estimate type:", est, "\n")
cat("Cluster sizes:\n")
cl <- factor(object@cluster, levels = 1:object@G)
cl_tbl <- table(cl)
cl_sizes <- as.matrix(trunc(cl_tbl))
rownames(cl_sizes) <- paste0("Cluster ", levels(cl))
colnames(cl_sizes) <- "Size"
print_matrix(cl_sizes, rownames(cl_sizes), colnames(cl_sizes), ndigits = 0, colwidth = 6, isint = TRUE)
cat("Subjects assigned to each cluster:\n")
for (g in 1:object@G) {
subj_clg <- which(object@cluster == g)
subj_str <- subj_clg[1]
if (length(subj_clg) > 1) {
for (h in 2:length(subj_clg)) {
subj_str <- paste0(subj_str, ", ", subj_clg[h])
}
}
cat("Cluster ", g, ": ", subj_str, "\n", sep = "")
}
}
)
#' Provide a summary of a \code{dmbc_config} class instance.
#'
#' @param object An object of class \code{\link{dmbc_config}}.
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases summary,dmbc_config-method
#' @aliases dmbc_config-summary
#'
#' @exportMethod summary
setMethod("summary",
"dmbc_config",
function(object) {
show(object)
cat("Latent configuration estimates (Z):\n")
wd <- as.integer(options("width"))
for (g in 1:object@G) {
dash_string_L <- strrep("-", floor((wd - nchar(paste0("Cluster ", object@G)) - 2)/2))
dash_string_R <- strrep("-", floor((wd - nchar(paste0("Cluster ", object@G)) - 2)/2) - nchar(g) + 1)
cat(dash_string_L, " ", "Cluster ", g, " ", dash_string_R, "\n", sep = "")
print_matrix(as.matrix(object@Z.est[, , g]),
if (length(object@labels)) object@labels else paste0("i = ", 1:object@n),
paste0("Dimension ", 1:object@p), colwidth = 11, ndigits = 4)
}
}
)
#' Provide a graphical summary of a \code{dmbc_config} class instance.
#'
#' @param x An object of class \code{\link{dmbc_config}}.
#' @param size A length-two numeric vector providing the optional sizes of
#' points and lines in the plot.
#' @param size_lbl A length-one numeric vector providing the size of labels.
#' @param nudge_x A length-one numeric vector providing the optional horizontal
#' adjustment to nudge labels by.
#' @param nudge_y A length-one numeric vector providing the optional vertical
#' adjustment to nudge labels by.
#' @param label_objects A length-one logical vector. If \code{TRUE}, labels are
#' added to the plot.
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases plot,dmbc_config-method
#' @aliases dmbc_config-plot
#'
#' @exportMethod plot
setMethod("plot",
signature(x = "dmbc_config"),
function(x, size = NULL, size_lbl = NULL, nudge_x = 0, nudge_y = 0, label_objects = TRUE, ...) {
n <- x@n
p <- x@p
G <- x@G
data <- prepare_data_to_plot(x)
n_G <- nlevels(data$G)
geom_args <- list()
if (is.null(size)) {
geom_args$size <- 2
} else {
geom_args$size <- size
}
if (is.null(size_lbl)) {
geom_args$size_lbl <- 3
} else {
geom_args$size_lbl <- size_lbl
}
mapping <- ggplot2::aes_(x = ~ p_1, y = ~ p_2, color = ~ G)
graph <- ggplot2::ggplot(data, mapping) + bayesplot::bayesplot_theme_get()
graph <- graph + ggplot2::geom_point(size = geom_args$size)
if (p <= 2) {
if (p == 1) {
min_val <- tapply(data[, 1], data$G, min, na.rm = TRUE)*1.15
max_val <- tapply(data[, 1], data$G, max, na.rm = TRUE)*1.15
lbl_cl <- unique(data$cl)
lbl_cl_rep <- margin.table(table(data$G, data$cl), margin = 2)/n
if (any(lbl_cl_rep > 1)) {
lbl_cl <- rep(lbl_cl, lbl_cl_rep[match(lbl_cl, names(lbl_cl_rep))])
}
lbl_data <- data.frame(min_val = min_val, max_val = max_val, G = factor(1:G, levels = 1:G),
cl = lbl_cl)
graph <- graph +
ggplot2::geom_segment(data = lbl_data,
mapping = ggplot2::aes_(x = ~ min_val, y = ~ min_val, xend = ~ max_val, yend = ~ max_val),
arrow = ggplot2::arrow(angle = 45, length = ggplot2::unit(0.05, "npc"), type = "open"), size = .5)
}
graph <- graph + ggplot2::facet_wrap(~ G + cl,
labeller = ggplot2::label_bquote(cols = paste("Cluster ", .(G), ", ", italic(S)[.(G)], " = ", .(cl))),
scales = "free")
} else {
graph <- graph + ggplot2::facet_grid(p_vs + p_i + p_j ~ G + cl,
labeller = ggplot2::label_bquote(
cols = paste("Cluster ", .(G), ", ", italic(S)[.(G)], " = ", .(cl)),
rows = paste("Dimension ", .(p_j), " vs. ", "Dimension ", .(p_i))),
scales = "free")
}
if (label_objects) {
graph <- graph +
ggrepel::geom_text_repel(ggplot2::aes(label = lbl), size = geom_args$size_lbl,
segment.size = .25, nudge_x = nudge_x, nudge_y = nudge_y)
}
graph <- graph +
ggplot2::scale_color_manual(G, values = choose_colors(n_G))
graph <- graph +
bayesplot::legend_move("none") +
bayesplot::xaxis_title(on = (p == 2)) +
bayesplot::yaxis_title(on = (p == 2)) +
force_axes_in_facets()
if (p == 2) {
graph <- graph + ggplot2::xlab(expression(italic(p)[1]))
graph <- graph + ggplot2::ylab(expression(italic(p)[2]))
}
graph
}
)
#' Extract the final cluster memberships from a \code{dmbc_config} class instance.
#'
#' @param object An object of class \code{\link{dmbc_config}}.
#' @param newdata An object of no explicit specification (currently ignored).
#' @param ... Further arguments to pass on (currently ignored).
#'
#' @author Sergio Venturini \email{sergio.venturini@unicatt.it}
#'
#' @aliases clusters,dmbc_config-method
#' @aliases dmbc_config-clusters
#'
#' @importFrom modeltools clusters
#' @exportMethod clusters
setMethod("clusters",
signature(object = "dmbc_config", newdata = "ANY"),
function(object, newdata = NULL, ...) {
object@cluster
}
)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.