R/plot_functions.R

In abnormally-distributed/abdisttools: Random Functions for Personal Use

#' Plot Posterior Densities in a Faceted Density Plot
#'
#'
#' @param D The model, data frame, or MCMC object.
#' @param cols Number of columns in facet plot
#' @param interval should credible intervals be plotted
#' @param cred.level the confidence level. defaults to  .95
#' @param greek should names be parsed for greek symbols
#' @param plot should the results be plotted? Defaults to FALSE.
#' @param extract.post Should the posterior be plotted? Defaults to TRUE.
#' @param extract.prior Should the prior be plotted? Defaults to FALSE.
#' @param ... other
#'
#' @export
#' @examples
#' post_dens_facet()

post_dens_facet <-
  function(D, cols = 2, interval = TRUE, greek = TRUE, cred.level = 0.95,
           extract.post = TRUE, extract.prior = FALSE, ...) {
    if (extract.post == TRUE) {
      D <- abdisttools::extract_post(D, ...)
    }
    if (extract.prior == TRUE) {
      D <- abdisttools::extract_prior(D, ...)
    }
    D <- as.data.frame(D)
    D <- D %>% gather(key = Parameter, value = value, factor_key = TRUE)
    f <- ggplot(D, aes(x = value))
    f <- f + stat_density(
      alpha = 0.15, size = 0.49, adjust = 4, bw = "SJ", n = 2048, trim = FALSE)
    
    if (!greek) {
      f <- f + facet_wrap(~Parameter, ncol = cols, scales = "free") +
        aes(fill = as.factor(Parameter), color = as.factor(Parameter)) +
        guides(fill = FALSE, color = FALSE) + theme(strip.text.x = element_text(margin = margin(
          0.088,
          0, 0.088, 0, "cm"
        ))) + xlab(label = NULL) + ylab(label = NULL)
    }
    else {
      f <- f + facet_wrap(~Parameter,
                          ncol = cols, scales = "free",
                          labeller = label_parsed
      ) + aes(
        fill = as.factor(Parameter),
        color = as.factor(Parameter)
      ) + guides(
        fill = FALSE,
        color = FALSE
      ) + theme(strip.text.x = element_text(margin = margin(
        0.088,
        0, 0.088, 0, "cm"
      ))) + xlab(label = NULL) + ylab(label = NULL)
    }
    if (interval == TRUE) {
      f <- f + tidybayes::stat_pointintervalh(aes(y = 0),
                                              point_interval = mean_qi,
                                              .width = cred.level
      ) + scale_fill_manual(values =c("#FA004F", "#00FF88", "#794CFF", "#FF9100", "#47BCFF",
                                      "#CCFF00", "#FF0000", "#0066FF", "#7B00FF", "#64E02A", 
                                      "#FFEE00", "#00FFE1", "#FF3300", "#3300FF", "#0047E0", 
                                      "#FFBB00", "#52BFAF", "#E100FF", "#70BAFF", "#00C227"), 
                            aesthetics = c("fill", "color"))
    }
    return(f)
  }


#' Generic plotting of intervals from a tidy summary
#'
#' @param terms the column containing the terms
#' @param estimates the column containing the estimates
#' @param H0 the null hypothesis. defaults to 0.
#' @param lower the column containing the lower limit of the intervals
#' @param upper the column containing the upper limit of the intervals
#' @param point.shape default to 18 (diamond)
#' @param h0.line.color standard error bar color
#' @param point.color point color
#' @param bar.color bar color
#' @param point.size point size
#' @param bar.size bar size
#' @export
#' @examples
#' interval_plot()
#'
interval_plot = function (terms, estimates, lower, upper, H0 = 0, point.color = "firebrick2", 
                          point.size = 4.75, point.shape = 18, bar.color = "dodgerblue2", 
                          bar.size = 0.85, h0.line.color = "white") 
{
  sum = cbind.data.frame(terms = factor(terms, levels = terms), 
                         estimates = estimates, lower = lower, upper = upper)
  sum %>% 
    
    ggplot(aes(y = forcats::fct_rev(terms), x = estimates, xmin = lower, xmax = upper)) + 
    labs(y = "term", x = "\u03B8") + 
    geom_vline(xintercept = H0, linetype = "dotted", size = 1, color = h0.line.color) +
    ggplot2::geom_errorbarh(color = bar.color, height = 0.18, size = bar.size, alpha = 1) + 
    geom_point(size = point.size, color = point.color, shape = point.shape) 
}


#' Plot Credible Intervals of a Bayesian Model 
#'
#' @param out the model fit
#' @param cred.level Confidence level. Defaults to .95.
#' @param H0 the null hypothesis. defaults to 0.
#' @param SEbars Should standard error bars (68.2 pct credible intervals) be superimposed on the wider interval? Defaults to TRUE.
#' @param ROPE If you would like, ROPE limits added to the plot.
#' @param point.shape default to 18 (diamond)
#' @param rope.color rope color
#' @param SEbar.color standard error bar color
#' @param point.color point color
#' @param bar.color bar color
#' @param point.size point size
#' @param bar.size bar size
#' @param droppars c("lp__", "hs_c2", "g", "g_inv", "log_lik", "inv_g", "eta", "K", "lambda", "inv_lambda")
#' @param ... other arguments to pass to post_summary, which is called internally
#' @export
#' @examples
#' bayes_intervals()
#'
bayes_intervals <- function(out, cred.level = 0.95, ROPE = NULL, rope.color = "white", SEbars = TRUE,
           SEbar.color = "dodgerblue4", point.color = "firebrick2", point.size = 3.75, point.shape = 18,
           bar.color = "dodgerblue2", bar.size = 0.85, H0 = 0, droppars = c(
             "lp__", "hs_c2", "g", "g_inv", "log_lik", "inv_g", "eta", "K",
             "lambda", "inv_lambda"), ...) {
    summary <- post_summary(out,
                            droppars = droppars, cred.level = cred.level,
                            ...
    )
    SEbarlower <- post_summary(out,
                               droppars = droppars, cred.level = 0.682,
                               ...
    )$cred.low
    SEbarupper <- post_summary(out,
                               droppars = droppars, cred.level = 0.682,
                               ...
    )$cred.high
    sum <- cbind.data.frame(
      terms = factor(summary$term, levels = summary$term),
      estimates = summary$estimate, lower = summary$cred.low,
      upper = summary$cred.high, SElower = SEbarlower, SEupper = SEbarupper
    )
    
    if (is.null(ROPE) == FALSE) {
      plot <- sum %>%
        ggplot(aes(
          y = forcats::fct_rev(terms), x = estimates, xmin = lower, xmax = upper)) + 
        labs(y = "term", x = "\u03B8") +
        geom_vline(xintercept = ROPE[1], linetype = "dotted", size = 0.7, color = rope.color) + 
        geom_vline(xintercept = ROPE[2], linetype = "dotted", size = 0.7, color = rope.color) 
    } else if (is.null(ROPE) == TRUE){
      plot <- sum %>%
        ggplot(aes(
          y = forcats::fct_rev(terms), x = estimates,
          xmin = lower, xmax = upper)) + labs(y = "term", x = "\u03B8") +
        geom_vline(xintercept = H0, linetype = "dotted", size = 0.75, color = rope.color)
    }
    
    if (SEbars == TRUE) {
      plot <- plot + ggplot2::geom_errorbarh(
        color = bar.color,
        height = 0.18, size = bar.size, alpha = 1,
      ) + ggstance::geom_linerangeh(
        xmin = sum$SElower,
        xmax = sum$SEupper, color = SEbar.color, size = bar.size * 3,
        alpha = 1
      ) + geom_point(size = point.size, color = point.color , shape = point.shape)
    }
    else {
      plot <- plot + ggplot2::geom_errorbarh(
        color = bar.color,
        height = 0.18, size = bar.size, alpha = 1
      ) + geom_point(
        size = point.size,
        color = point.color,
        shape = point.shape
      )
    }
    plot(plot)
}

#' Plot posterior distributions from model fit
#'
#' @param mcmc.list. an MCMC list
#' @param family. group of parameters
#' with the same name but different numerical value
#' between square brackets (as beta[1], beta[2])
#' @param description. Character vector giving a short descriptive text that identifies the model.
#' @param par_labels. data frame with two colums. One named "Parameter"
#' with the same names of the parameters of the model.
#' @param sort. Logical. When TRUE (the default), parameters are sorted first by
#' family name and then by numerical value.
#' @param inc_warmup. Logical. When dealing with stanfit objects from rstan,
#' logical value whether the warmup samples are included. Defaults to FALSE.
#' @param stan_include_auxiliar.	Logical value to include "lp__" parameter in rstan, and "lp__", "treedepth__" and
#' "stepsize__" in stan running without rstan. Defaults to FALSE.
#' @param droppars variables to exclude. Defaults to droppars = c("lp__","hs_c2","nu","sigma").
#' @export
#' @examples
#' tidyGg()
#'
#'
ridge.plots = function(mcmc,family. = NA, description. = NA, burnin. = FALSE, par_labels. = NA,
                       sort. = TRUE, inc_warmup. = FALSE, stan_include_auxiliar. = FALSE,
                       droppars = c("lp__","hs_c2","nu","sigma"), alpha=.65, hline.col="#95bdd6", hline.alpha=.90)
{
  
  
  mcmc.list = as.mcmc(as.matrix(mcmc))
  
  t_col <- function(color, percent = 50, name = NULL) {
    #	  color = color name
    #	percent = % transparency
    #	   name = an optional name for the color
    ## Get RGB values for named color
    rgb.val <- col2rgb(color)
    ## Make new color using input color as base and alpha set by transparency
    t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3],
                 max = 255,
                 alpha = (100-percent)*255/100,
                 names = name)
    ## Save the color
    invisible(t.col)
    
  }
  color.transparent = t_col(hline.col, perc = hline.alpha, name = "transparent.white")
  GG = mcmc.list
  fit.gg = tidyGg(GG)
  ggplot(fit.gg, aes(x = `value`, y = forcats::fct_rev(`Parameter`), fill=0.5 - abs(0.5-..ecdf..))) +
    stat_density_ridges(geom = "density_ridges_gradient", calc_ecdf = TRUE,color=color.transparent) +
    scico::scale_fill_scico(direction = 1,palette="oslo",alpha= alpha) +
    labs(y="term", x="\u03B8") +
    guides(fill=FALSE)
  
}


#' Check the sensitivity of credible intervals for the same data fit to different models to different priors
#'
#' @param list a list of tidy summaries where the summaries are all of the same time with the same contents
#' @param H0 the null hypothesis, defaults to 0
#' @param adj the amount of space between intervals, deaults to .75
#' @param plot set to TRUE to plot
#' @export
#' @examples
#' compare_intervals()
#'
compare_intervals = function(list,H0=0,model.names="default",adj=.75){
  
  if (model.names=="default") {
    Model = as.character(c(sapply(1:length(list), function(n) rep(paste0("model", n), nrow(list[[1]])))))
  } else {
    Model = as.character(c(sapply(1:length(model.names), function(n) rep(model.names[n], nrow(list[[1]])))))
  }
  
  posts = as.data.frame(dplyr::bind_rows(list))
  low = which(colnames(posts)=="hdi.low"|colnames(posts)=="LPL.low"|colnames(posts)=="hdi.low"|colnames(posts)=="hdp.low"|colnames(posts)=="cred.low")
  high = which(colnames(posts)=="hdi.high"|colnames(posts)=="LPL.high"|colnames(posts)=="hdi.high"|colnames(posts)=="hdp.high"|colnames(posts)=="cred.high")
  colnames(posts)[low]  = "lower"
  colnames(posts)[high] = "upper"
  
  posts = cbind.data.frame(Model=Model,posts)
  posts$Model = as.character(posts$Model)
  
  posts %>%
    ggplot(aes(y = forcats::fct_rev(term), x = estimate, xmin = lower, xmax = upper, color = Model)) +
    labs(y="term", x="\u03B8") +
    geom_point(size=1.6,position=ggstance::position_dodgev(height = adj+.001, preserve="total")) +
    geom_vline(xintercept=H0, linetype="dotted",size=1) +
    geom_errorbarh(height = .5, size=.75, position=ggstance::position_dodgev(height = adj, preserve="total")) 
}

#' Plot bayes factors from the output of adv_summary
#'
#' @param summary the summary
#' @param vjust the amount of vertical adjustment
#' @export
#' @examples
#' bf_plot()
#'
bf_plot = function(summary, vjust=0.6, sig.bf = 4){
  
  summary =cbind.data.frame(summary, color=ifelse(10*log10(summary$BF10) < 0, "NULL", "ALT"))
  
  ggplot(summary, aes(y=forcats::fct_rev(term), x=10*log10(BF10))) +
    ggplot2::scale_color_manual(values = c("steelblue4", "coral2")) +
    guides(color=guide_legend(title="Evidence (in dB) in favor of")) +
    geom_point(size=3.5,aes(color=color)) +
    geom_vline(xintercept = 10*log10(sig.bf), linetype="dashed") +
    geom_vline(xintercept = 10*log10(1/sig.bf), linetype="dashed") +
    geom_vline(xintercept = 0) +
    geom_segment(aes(y=term,
                     yend=term,
                     x=0,
                     xend=10*log10(BF10),color=color),size=1.5) +
    theme(axis.text.x = element_text(vjust=vjust)) +
    scale_x_continuous(name="10*log10(BF10)") +
    scale_y_discrete(name="term")
}

#' Check the sensitivity of bayes factors for the same data fit to different models to different priors
#'
#' @param list a list of advanced summaries where the summaries are all of different model fit to the same data
#' @param H0 the null hypothesis, defaults to 0
#' @param adj the amount of space between intervals, deaults to .75
#' @param plot set to TRUE to plot
#' @export 
#' @examples
#' bf_robust()
#'
bf_robust =  function(list,H0=0,model.names="default",adj=.75){
  
  if (model.names=="default") {
    Model = as.character(c(sapply(1:length(list), function(n) rep(paste0("model", n), nrow(list[[1]])))))
  } else {
    Model = as.character(c(sapply(1:length(model.names), function(n) rep(model.names[n], nrow(list[[1]])))))
  }
  
  posts = as.data.frame(dplyr::bind_rows(list))
  low = which(colnames(posts)=="hdi.low"|colnames(posts)=="LPL.low"|colnames(posts)=="hdi.low"|colnames(posts)=="hdp.low"|colnames(posts)=="cred.low")
  high = which(colnames(posts)=="hdi.high"|colnames(posts)=="LPL.high"|colnames(posts)=="hdi.high"|colnames(posts)=="hdp.high"|colnames(posts)=="cred.high")
  colnames(posts)[low]  = "lower"
  colnames(posts)[high] = "upper"
  
  posts = cbind.data.frame(Model=Model,posts)
  posts$Model = as.character(posts$Model)
  posts$term =  forcats::fct_rev(posts$term)
  posts = cbind.data.frame(posts, Direction=ifelse(log(posts$BF10) < 0, "NULL", "ALT"))
  
  
  posts %>%
    ggplot(aes(y = forcats::fct_rev(term), x = log(BF10),fill=Model,color=Direction)) +
    geom_segment(aes(y=term,
                     yend=term,
                     x=0,
                     xend=log(BF10)),size=2,alpha=1) +
    geom_point(aes(fill=Model, color=Direction),shape=21,size=5,alpha=.9)+
    ggplot2::scale_color_manual(values = c("steelblue4", "coral2")) +
    guides(color=guide_legend(title="Evidence in favor of")) +
    scale_shape_discrete(solid=F) +
    #guides(fill = guide_legend(override.aes=list(shape=21))) +
    geom_vline(xintercept = 2.302585, linetype="dashed") +
    geom_vline(xintercept = -2.302585, linetype="dashed") +
    geom_vline(xintercept = 0) +
    ggplot2::labs(y="term") +
    scale_fill_aaas()
  
}