Nothing
#' Summarize and extract posterior chains from MCMC output while preserving parameter structure
#'
#' Extract summary information and posterior chains from MCMC output (specific function specified) for specific parameters of interest while preserving original parameter structure (i.e., scalar, vector, matrix, array). Function outputs a \code{list} with calculated values or posterior chains for each specified parameter.
#'
#' @param object Object containing MCMC output. See DETAILS below.
#'
#' @param params Character string (or vector of character strings) denoting parameters to be returned in output.
#'
#' Default \code{'all'} returns all parameters in output.
#'
#' @param excl Character string (or vector of character strings) denoting parameters to exclude. Used in conjunction with \code{params} argument to select parameters of interest.
#'
#' @param ISB Ignore Square Brackets (ISB). Logical specifying whether square brackets should be ignored in the \code{params} and \code{excl} arguments. If \code{TRUE}, square brackets are ignored. If \code{FALSE}, square brackets are not ignored. This allows partial names to be used when specifying parameters of interest. Use \code{exact} argument to specify whether input from \code{params} and \code{excl} arguments should be matched exactly.
#'
#' @param exact Logical specifying whether input from \code{params} and \code{excl} arguments should be matched exactly (after ignoring square brackets if \code{ISB = FALSE}). If \code{TRUE}, input from \code{params} and \code{excl} are matched exactly (after taking \code{ISB} argument into account). If \code{FALSE}, input from \code{params} and \code{excl} are matched using regular expression format (after taking \code{ISB} argument into account).
#'
#' @param func Function to be performed on MCMC output. When output of specified function is greater than length 1, an extra dimension is added. For instance, output of length 3 for a parameter with dimensions 2x2 results in a 2x2x3 output. Functions that produce output with dimensionality greater than 1 are not permitted. \code{func} is ignored when \code{type = 'chains'}.
#'
#' @param type Character string specifying whether to return summary information (calculated based on \code{func} argument) or posterior chains. Valid options are \code{'summary'} and \code{'chains'}. When \code{type = 'chains'}, the \code{'func'} argument is ignored. When \code{type = 'chains'}, posterior chains are concatenated and stored in the last dimension in the array for each element (parameter) of the list.
#'
#' @section Details:
#' \code{object} argument can be a \code{stanfit} object (\code{rstan} package), a \code{CmdStanMCMC} object (\code{cmdstanr} package), a \code{stanreg} object (\code{rstanarm} package), a \code{brmsfit} object (\code{brms} package), an \code{mcmc.list} object (\code{coda} and \code{rjags} packages), \code{mcmc} object (\code{coda} and \code{nimble} packages), \code{list} object (\code{nimble} package), an \code{R2jags} model object (\code{R2jags} package), a \code{jagsUI} model object (\code{jagsUI} package), or a matrix containing MCMC chains (each column representing MCMC output for a single parameter, rows representing iterations in the chain). The function automatically detects the object type and proceeds accordingly.
#'
#' @examples
#' #Load data
#' data(MCMC_data)
#'
#' MCMCpstr(MCMC_data, func = function(x) quantile(x, probs = c(0.01, 0.99)))
#'
#' @export
MCMCpstr <- function(object,
params = 'all',
excl = NULL,
ISB = TRUE,
exact = TRUE,
func = mean,
type = 'summary')
{
#SORTING BLOCK
if (methods::is(object, 'matrix'))
{
object2 <- MCMCchains(object, params, excl, ISB, exact = exact, mcmc.list = FALSE)
} else {
object2 <- MCMCchains(object, params, excl, ISB, exact = exact, mcmc.list = TRUE)
}
if (coda::is.mcmc.list(object2) == TRUE)
{
temp_in <- object2
cti <- colnames(temp_in[[1]])
if (methods::is(object, 'brmsfit') & length(grep('Intercept]', cti)) > 1)
{
#remove Intercept and close bracket
i_idx <- grep('Intercept]', cti)
cti[i_idx] <- gsub(',Intercept', '', cti[i_idx])
colnames(temp_in[[1]]) <- cti
}
if (ISB == TRUE)
{
names <- vapply(strsplit(cti,
split = "[", fixed = TRUE), `[`, 1, FUN.VALUE=character(1))
} else {
names <- cti
}
np <- NCOL(object2[[1]])
if (np > 1)
{
ch_bind <- do.call('rbind', object2)
} else {
ch_bind <- as.matrix(object2)
}
#how many elements will be in the list
un <- unique(names)
onames <- colnames(temp_in[[1]])
}
if (methods::is(object2, 'matrix'))
{
temp_in <- object2
if (ISB == TRUE)
{
names <- vapply(strsplit(colnames(temp_in),
split = "[", fixed = TRUE), `[`, 1, FUN.VALUE=character(1))
} else {
names <- colnames(temp_in)
}
np <- NCOL(object2)
ch_bind <- object2
#how many elements will be in the list
un <- unique(names)
onames <- colnames(temp_in)
}
#create empty list
out_list <- vector('list', length(un))
#iterate through each element
for (i in 1:length(un))
{
#i <- 1
ind <- which(un[i] == names)
#determine how many ',' and therefore how many dimensions for parameter
#if only one ind, there should only be one dim
if (length(ind) == 1)
{
dims <- 1
} else {
dims <- length(strsplit(onames[ind[1]], split = ',', fixed = TRUE)[[1]])
}
#scalar or vector
if (dims == 1)
{
if (type == 'summary')
{
func_out <- func(ch_bind[,ind[1]])
if (length(func_out) > 1)
{
if (!is.null(dim(func_out)))
{
stop("Output from 'func' argument must be of dimension 1.")
}
temp_obj <- matrix(NA, nrow = length(ind), ncol = length(func_out))
dimnames(temp_obj)[[2]] <- names(func_out)
dimnames(temp_obj)[[1]] <- onames[ind]
for (j in 1:length(ind))
{
temp_obj[j,] <- func(ch_bind[,ind[j]])
}
} else {
temp_obj <- rep(NA, length(ind))
for (j in 1:length(ind))
{
temp_obj[j] <- func(ch_bind[,ind[j]])
}
}
}
if (type == 'chains')
{
temp_obj <- matrix(NA, nrow = length(ind), ncol = NROW(ch_bind))
if (NROW(temp_obj) > 1)
{
dimnames(temp_obj)[[1]] <- onames[ind]
} else {
rownames(temp_obj) <- onames[ind]
}
for (j in 1:length(ind))
{
temp_obj[j,] <- ch_bind[,ind[j]]
}
}
if (type != 'summary' & type != 'chains')
{
stop("Invalid input for argument 'type'. Valid options are 'summary' and 'chains'.")
}
#fill list
out_list[[i]] <- temp_obj
}
#2 dimensions
if (dims == 2)
{
pnames <- vapply(strsplit(onames[ind],
'[', fixed = TRUE), `[`, 2, FUN.VALUE=character(1))
pnames2 <- vapply(strsplit(pnames,
']', fixed = TRUE), `[`, 1, FUN.VALUE=character(1))
#determine how large matrix should be
RI <- c()
CI <- c()
for (j in 1:length(ind))
{
#j <- 1
tt <- strsplit(pnames2[j], ',')
ri <- as.numeric(tt[[1]][1])
ci <- as.numeric(tt[[1]][2])
RI <- c(RI, ri)
CI <- c(CI, ci)
}
if (type == 'summary')
{
func_out <- func(ch_bind[,ind[1]])
if (length(func_out) > 1)
{
if (!is.null(dim(func_out)))
{
stop("Output from 'func' argument must be of dimension 1.")
}
#create blank array
temp_obj <- array(NA, dim = c(max(RI), max(CI), length(func_out)))
dimnames(temp_obj)[[3]] <- names(func_out)
for (j in 1:length(ind))
{
#j <- 1
temp_obj[RI[j], CI[j], ] <- func(ch_bind[,ind[j]])
}
} else {
#create blank matrix
temp_obj <- matrix(NA, nrow = max(RI), ncol = max(CI))
for (j in 1:length(ind))
{
temp_obj[RI[j], CI[j]] <- func(ch_bind[,ind[j]])
}
}
}
if (type == 'chains')
{
#create blank array
temp_obj <- array(NA, dim = c(max(RI), max(CI), NROW(ch_bind)))
for (j in 1:length(ind))
{
temp_obj[RI[j], CI[j], ] <- ch_bind[,ind[j]]
}
}
if (type != 'summary' & type != 'chains')
{
stop("Invalid input for argument 'type'. Valid options are 'summary' and 'chains'.")
}
#fill list
out_list[[i]] <- temp_obj
}
#3 dimensions
if (dims == 3)
{
pnames <- vapply(strsplit(onames[ind],
'[', fixed = TRUE), `[`, 2, FUN.VALUE=character(1))
pnames2 <- vapply(strsplit(pnames,
']', fixed = TRUE), `[`, 1, FUN.VALUE=character(1))
#determine how large array should be
RI <- c()
CI <- c()
TI <- c()
for (j in 1:length(ind))
{
#j <- 1
tt <- strsplit(pnames2[j], ',')
ri <- as.numeric(tt[[1]][1])
ci <- as.numeric(tt[[1]][2])
ti <- as.numeric(tt[[1]][3])
RI <- c(RI, ri)
CI <- c(CI, ci)
TI <- c(TI, ti)
}
if (type == 'summary')
{
func_out <- func(ch_bind[,ind[1]])
if (length(func_out) > 1)
{
if (!is.null(dim(func_out)))
{
stop("Output from 'func' argument must be of dimension 1.")
}
#create blank array
temp_obj <- array(NA, dim = c(max(RI), max(CI), max(TI), length(func_out)))
dimnames(temp_obj)[[4]] <- names(func_out)
for (j in 1:length(ind))
{
#j <- 1
temp_obj[RI[j], CI[j], TI[j], ] <- func(ch_bind[,ind[j]])
}
} else {
#create blank array
temp_obj <- array(NA, dim = c(max(RI), max(CI), max(TI)))
for (j in 1:length(ind))
{
temp_obj[RI[j], CI[j], TI[j]] <- func(ch_bind[,ind[j]])
}
}
}
if (type == 'chains')
{
#create blank array
temp_obj <- array(NA, dim = c(max(RI), max(CI), max(TI), NROW(ch_bind)))
for (j in 1:length(ind))
{
temp_obj[RI[j], CI[j], TI[j], ] <- ch_bind[,ind[j]]
}
}
if (type != 'summary' & type != 'chains')
{
stop("Invalid input for argument 'type'. Valid options are 'summary' and 'chains'.")
}
#fill list
out_list[[i]] <- temp_obj
}
#4 dimensions
if (dims == 4)
{
pnames <- vapply(strsplit(onames[ind],
'[', fixed = TRUE), `[`, 2, FUN.VALUE=character(1))
pnames2 <- vapply(strsplit(pnames,
']', fixed = TRUE), `[`, 1, FUN.VALUE=character(1))
#determine how large matrix should be
RI <- c()
CI <- c()
TI <- c()
FI <- c()
for (j in 1:length(ind))
{
#j <- 1
tt <- strsplit(pnames2[j], ',')
ri <- as.numeric(tt[[1]][1])
ci <- as.numeric(tt[[1]][2])
ti <- as.numeric(tt[[1]][3])
fi <- as.numeric(tt[[1]][4])
RI <- c(RI, ri)
CI <- c(CI, ci)
TI <- c(TI, ti)
FI <- c(FI, fi)
}
if (type == 'summary')
{
func_out <- func(ch_bind[,ind[1]])
if (length(func_out) > 1)
{
if (!is.null(dim(func_out)))
{
stop("Output from 'func' argument must be of dimension 1.")
}
#create blank array
temp_obj <- array(NA, dim = c(max(RI), max(CI), max(TI), max(FI), length(func_out)))
dimnames(temp_obj)[[5]] <- names(func_out)
for (j in 1:length(ind))
{
#j <- 1
temp_obj[RI[j], CI[j], TI[j], FI[j], ] <- func(ch_bind[,ind[j]])
}
} else {
#create blank array
temp_obj <- array(NA, dim = c(max(RI), max(CI), max(TI), max(FI)))
for (j in 1:length(ind))
{
temp_obj[RI[j], CI[j], TI[j], FI[j]] <- func(ch_bind[,ind[j]])
}
}
}
if (type == 'chains')
{
#create blank array
temp_obj <- array(NA, dim = c(max(RI), max(CI), max(TI), max(FI), NROW(ch_bind)))
for (j in 1:length(ind))
{
temp_obj[RI[j], CI[j], TI[j], FI[j], ] <- ch_bind[,ind[j]]
}
}
if (type != 'summary' & type != 'chains')
{
stop("Invalid input for argument 'type'. Valid options are 'summary' and 'chains'.")
}
#fill list
out_list[[i]] <- temp_obj
}
if (dims > 4)
{
stop('This function does not currently support parameters with > 4 dimensions. If you have a need for this functionality, please create an "issue" at https://github.com/caseyyoungflesh/MCMCvis.')
}
}
#name elements in list
names(out_list) <- un
return(out_list)
}
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.