#' Posteriors Sampling Diagnostic #' #' Extract diagnostic metrics (Effective Sample Size (ESS), Rhat and Monte #' Carlo Standard Error MCSE). #' #' @param posteriors A stanreg or brms model. #' @param diagnostic Diagnostic metrics to compute. Character (vector) or list #' with one or more of these options: "ESS", "Rhat", "MCSE" or "all". #' #' @details #' **Effective Sample (ESS)** should be as large as possible, although for #' most applications, an effective sample size greater than 1000 is sufficient #' for stable estimates (Bürkner, 2017). The ESS corresponds to the number of #' independent samples with the same estimation power as the N autocorrelated #' samples. It is is a measure of \dQuote{how much independent information #' there is in autocorrelated chains} (\cite{Kruschke 2015, p182-3}). #' \cr \cr #' **Rhat** should be the closest to 1. It should not be larger than 1.1 #' (\cite{Gelman and Rubin, 1992}) or 1.01 (\cite{Vehtari et al., 2019}). The #' split Rhat statistic quantifies the consistency of an ensemble of Markov #' chains. #' \cr \cr #' **Monte Carlo Standard Error (MCSE)** is another measure of accuracy of the #' chains. It is defined as standard deviation of the chains divided by their #' effective sample size (the formula for mcse() is from Kruschke 2015, p. #' 187). The MCSE \dQuote{provides a quantitative suggestion of how big the #' estimation noise is}. #' #' #' @examples #' \dontrun{ #' # rstanarm models #' # ----------------------------------------------- #' if (require("rstanarm", quietly = TRUE)) { #' model <- stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0) #' diagnostic_posterior(model) #' } #' #' # brms models #' # ----------------------------------------------- #' if (require("brms", quietly = TRUE)) { #' model <- brms::brm(mpg ~ wt + cyl, data = mtcars) #' diagnostic_posterior(model) #' } #' } #' @references #' \itemize{ #' \item Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical science, 7(4), 457-472. #' \item Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., \& Bürkner, P. C. (2019). Rank-normalization, folding, and localization: An improved Rhat for assessing convergence of MCMC. arXiv preprint arXiv:1903.08008. #' \item Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press. #' } #' @export diagnostic_posterior <- function(posteriors, diagnostic = c("ESS", "Rhat"), ...) { UseMethod("diagnostic_posterior") } #' @export diagnostic_posterior.numeric <- function(posteriors, diagnostic = c("ESS", "Rhat"), ...) { stop("diagnostic_posterior only works with rstanarm or brms models.") } #' @export diagnostic_posterior.data.frame <- diagnostic_posterior.numeric #' @export diagnostic_posterior.BFBayesFactor <- diagnostic_posterior.numeric #' @inheritParams insight::get_parameters.BFBayesFactor #' @inheritParams insight::get_parameters #' @rdname diagnostic_posterior #' @export diagnostic_posterior.stanreg <- function(posteriors, diagnostic = "all", effects = c("fixed", "random", "all"), component = c("location", "all", "conditional", "smooth_terms", "sigma", "distributional", "auxiliary"), parameters = NULL, ...) { # Find parameters effects <- match.arg(effects) component <- match.arg(component) params <- insight::find_parameters( posteriors, effects = effects, component = component, parameters = parameters, flatten = TRUE ) # If no diagnostic if (is.null(diagnostic)) { return(data.frame("Parameter" = params)) } diagnostic <- match.arg(diagnostic, c("ESS", "Rhat", "MCSE", "all"), several.ok = TRUE) if ("all" %in% diagnostic) { diagnostic <- c("ESS", "Rhat", "MCSE", "khat") } else { diagnostic <- c(diagnostic) if ("Rhat" %in% diagnostic) diagnostic <- c(diagnostic, "khat") } # Get indices and rename diagnostic_df <- as.data.frame(posteriors$stan_summary) diagnostic_df$Parameter <- row.names(diagnostic_df) if ("n_eff" %in% names(diagnostic_df)) { diagnostic_df$ESS <- diagnostic_df$n_eff } # special handling for MCSE, due to some parameters (like lp__) missing in rows MCSE <- mcse(posteriors, effects = "all") diagnostic_df <- merge(diagnostic_df, MCSE, by = "Parameter", all = FALSE) # Select columns available_columns <- intersect(colnames(diagnostic_df), c("Parameter", diagnostic)) diagnostic_df <- diagnostic_df[available_columns] names(diagnostic_df)[available_columns == "khat"] <- "Khat" row.names(diagnostic_df) <- NULL # Remove columns with all Nans diagnostic_df <- diagnostic_df[!sapply(diagnostic_df, function(x) all(is.na(x)))] # Select rows diagnostic_df[diagnostic_df$Parameter %in% params, ] } #' @inheritParams insight::get_parameters #' @rdname diagnostic_posterior #' @export diagnostic_posterior.stanmvreg <- function(posteriors, diagnostic = "all", effects = c("fixed", "random", "all"), parameters = NULL, ...) { # Find parameters effects <- match.arg(effects) all_params <- insight::find_parameters( posteriors, effects = effects, parameters = parameters, flatten = FALSE ) params <- unlist(lapply(names(all_params), function(i) { all_params[[i]]$sigma <- NULL unlist(all_params[[i]]) })) # If no diagnostic if (is.null(diagnostic)) { return(data.frame("Parameter" = params)) } diagnostic <- match.arg(diagnostic, c("ESS", "Rhat", "MCSE", "all"), several.ok = TRUE) if ("all" %in% diagnostic) { diagnostic <- c("ESS", "Rhat", "MCSE", "khat") } else { diagnostic <- c(diagnostic) if ("Rhat" %in% diagnostic) diagnostic <- c(diagnostic, "khat") } # Get indices and rename diagnostic_df <- as.data.frame(posteriors$stan_summary) diagnostic_df$Parameter <- row.names(diagnostic_df) if ("n_eff" %in% names(diagnostic_df)) { diagnostic_df$ESS <- diagnostic_df$n_eff } # special handling for MCSE, due to some parameters (like lp__) missing in rows MCSE <- mcse(posteriors, effects = effects) diagnostic_df <- merge(diagnostic_df, MCSE, by = "Parameter", all = FALSE) # Select columns available_columns <- intersect(colnames(diagnostic_df), c("Parameter", diagnostic)) diagnostic_df <- diagnostic_df[available_columns] names(diagnostic_df)[available_columns == "khat"] <- "Khat" row.names(diagnostic_df) <- NULL # Remove columns with all Nans diagnostic_df <- diagnostic_df[!sapply(diagnostic_df, function(x) all(is.na(x)))] diagnostic_df$Response <- gsub("(b\\[)*(.*)\\|(.*)", "\\2", diagnostic_df$Parameter) for (i in unique(diagnostic_df$Response)) { diagnostic_df$Parameter <- gsub(sprintf("%s|", i), "", diagnostic_df$Parameter, fixed = TRUE) } # Select rows diagnostic_df[diagnostic_df$Parameter %in% params, ] } #' @inheritParams insight::get_parameters #' @rdname diagnostic_posterior #' @export diagnostic_posterior.brmsfit <- function(posteriors, diagnostic = "all", effects = c("fixed", "random", "all"), component = c("conditional", "zi", "zero_inflated", "all"), parameters = NULL, ...) { # Find parameters effects <- match.arg(effects) component <- match.arg(component) params <- insight::find_parameters(posteriors, effects = effects, component = component, parameters = parameters, flatten = TRUE ) # If no diagnostic if (is.null(diagnostic)) { return(data.frame("Parameter" = params)) } # Get diagnostic diagnostic <- match.arg(diagnostic, c("ESS", "Rhat", "MCSE", "all"), several.ok = TRUE) if ("all" %in% diagnostic) { diagnostic <- c("ESS", "Rhat", "MCSE", "khat") # Add MCSE } else { if ("Rhat" %in% diagnostic) diagnostic <- c(diagnostic, "khat") } insight::check_if_installed("rstan") # Get indices and rename diagnostic_df <- as.data.frame(rstan::summary(posteriors$fit)$summary) diagnostic_df$Parameter <- row.names(diagnostic_df) diagnostic_df$ESS <- diagnostic_df$n_eff # special handling for MCSE, due to some parameters (like lp__) missing in rows MCSE <- mcse(posteriors, effects = "all", component = "all") diagnostic_df <- merge(diagnostic_df, MCSE, by = "Parameter", all = FALSE) # Select columns available_columns <- intersect(colnames(diagnostic_df), c("Parameter", diagnostic)) diagnostic_df <- diagnostic_df[available_columns] names(diagnostic_df)[available_columns == "khat"] <- "Khat" row.names(diagnostic_df) <- NULL # Remove columns with all Nans diagnostic_df <- diagnostic_df[!sapply(diagnostic_df, function(x) all(is.na(x)))] # Select rows diagnostic_df[diagnostic_df$Parameter %in% params, ] } #' @inheritParams insight::get_parameters #' @export diagnostic_posterior.stanfit <- function(posteriors, diagnostic = "all", effects = c("fixed", "random", "all"), parameters = NULL, ...) { # Find parameters effects <- match.arg(effects) params <- insight::find_parameters(posteriors, effects = effects, parameters = parameters, flatten = TRUE) # If no diagnostic if (is.null(diagnostic)) { return(data.frame("Parameter" = params)) } # Get diagnostic diagnostic <- match.arg(diagnostic, c("ESS", "Rhat", "MCSE", "all"), several.ok = TRUE) if ("all" %in% diagnostic) { diagnostic <- c("ESS", "Rhat", "MCSE") } insight::check_if_installed("rstan") all_params <- insight::find_parameters(posteriors, effects = effects, flatten = TRUE ) diagnostic_df <- data.frame( Parameter = all_params, stringsAsFactors = FALSE ) if ("ESS" %in% diagnostic) { diagnostic_df$ESS <- effective_sample(posteriors, effects = effects)$ESS } if ("MCSE" %in% diagnostic) { diagnostic_df$MCSE <- mcse(posteriors, effects = effects)$MCSE } if ("Rhat" %in% diagnostic) { s <- as.data.frame(rstan::summary(posteriors)$summary) diagnostic_df$Rhat <- s[rownames(s) %in% all_params, ]$Rhat } # Remove columns with all Nans diagnostic_df <- diagnostic_df[!sapply(diagnostic_df, function(x) all(is.na(x)))] # Select rows diagnostic_df[diagnostic_df$Parameter %in% params, ] } #' @export diagnostic_posterior.blavaan <- function(posteriors, diagnostic = "all", ...) { # Find parameters params <- suppressWarnings(insight::find_parameters(posteriors, flatten = TRUE)) out <- data.frame("Parameter" = params) # If no diagnostic if (is.null(diagnostic)) { return(out) } diagnostic <- match.arg(diagnostic, c("ESS", "Rhat", "MCSE", "all"), several.ok = TRUE) if ("all" %in% diagnostic) { diagnostic <- c("ESS", "Rhat", "MCSE") } else { diagnostic <- c(diagnostic) if ("Rhat" %in% diagnostic) diagnostic <- c(diagnostic, "khat") } # Get indices if ("Rhat" %in% diagnostic) { insight::check_if_installed("blavaan") Rhat <- blavaan::blavInspect(posteriors, what = "psrf") Rhat <- data.frame( Parameter = colnames(insight::get_parameters(posteriors)), Rhat = Rhat ) out <- merge(out, Rhat, by = "Parameter", all = TRUE) } if ("ESS" %in% diagnostic) { ESS <- effective_sample(posteriors) out <- merge(out, ESS, by = "Parameter", all = TRUE) } if ("MCSE" %in% diagnostic) { MCSE <- mcse(posteriors) out <- merge(out, MCSE, by = "Parameter", all = TRUE) } unique(out) }