xfim / ggmcmc

Graphical tools for analyzing Markov Chain Monte Carlo simulations from Bayesian inference
111 stars 31 forks source link

Can something like the " gelman.plot" from the "coda" package be added? Simple function below #69

Closed ltrainstg closed 4 years ago

ltrainstg commented 4 years ago

I found this package recently and I like how it uses ggplot instead of base R for graphs. One of the only graphs that seemed to be missing that was in "coda" that I used frequently was 'gelman.plot'. The current 'ggs_Rhat' function only shows the end result. It is useful to see how this changes over time to see if the parameter has truly converged and not simply lucked into an early convergence.

I attempted to create a function that mimics the gelman.plot from coda as much as possible below. Could something similar be added to this package?

library(coda)
library(ggmcmc)
library(dplyr)
data(s.binary)
gelman.plot(s.binary)
gg_gelmanPlot <- function(x){

  x <- as.mcmc.list(x)
  y <- coda:::gelman.preplot(x, bin.width = 10, max.bins = 50, confidence = 0.95, 
                             transform = FALSE, autoburnin = TRUE)

  df <- NULL
  for(i in 1:2){
    df.piece <- cbind(
      y$last.iter,
      y$shrink[, , 1]) %>%
      as.data.frame()
    baseParams <- colnames(y$shrink[,,i])
    names( df.piece) <- c("iterations", baseParams)
    df.piece <- 
      gather(df.piece, condition, measurement,
             baseParams[1]:rev(baseParams)[1],
             factor_key=TRUE)
    df.piece$i <- i
    df <- rbind(df, df.piece)
  }

  df <- df %>% mutate(Variable = ifelse(i==1, "Median", "97.5%"))
  ggplot(dat = df,
         aes(x=iterations, y=measurement, col= Variable)) +
    geom_line(aes(linetype = Variable)) + 
    geom_hline(yintercept = 1.05, linetype = "dashed") +
    facet_grid(vars(condition), scales = "free") +
    xlab("last iteration in chain") + 
    ylab("shrink factor")

}
gg_gelmanPlot(s.binary)
xfim commented 4 years ago

Dear @ltrainstg ,

Thank you very much for your suggestion and the code provided.

I adapted it and produced a new function ggs_grb() (the name includes the whole group of contributors to the plot idea).

Let me know what you think, see commit f6dddd455116d07e909dfe1676c3f526073889d0.

It is already packed in version 1.5.0.

I would love to hear suggestions specifically for: