Functions to assess convergence

[lmsimp]

We are missing some helper functions used to help assess convergence which are in pRoloc e.g. mcmc_get_outliers and friends. If you're happy for me to add these I will add to my to-do list and then open a PR.

[lmsimp]

@ococrook I have taken some functionality from pRoloc and also created some new plots which I think help with assessing convergence of the bandle chains.

Before running bandleProcess or summaries and straight after the bandle run, the bandleParams object is created and I think it would be nice to enable users to visualise the convergence of the chains at this initial stage to allow users to make decisions about throwing away bad chains.

Three functions I'd like to add:

1. bandle_get_outliers which is essentially the same function as mcmc_get_outliers from pRoloc. It takes a bandleParams object and extracts the outliers using coda::mcmc something like [1]
2. plotOutliers which would again take a bandleParams object and generate trace and density plots for all chains and produce something like this:
   
   using code from coda e.g. see [2]
3. A helper function [3] for computing Gelman diagnostics. This would work on bandleParams before we have processed to allow us to check the chains and remove any at this stage that are bad. The idea is to run a function say calculateGelman and output all pairwise Gelman stats for each chain, output would be something like this:

Condition1
            comb_12  comb_13   comb_14  comb_15  comb_23   comb_24  comb_25  comb_34  comb_35  comb_45
Point_Est 0.9993745 1.000839 0.9980787 1.001300 1.003825 0.9987998 1.006697 1.001453 1.000449 1.002392
Upper_CI  1.0003297 1.005529 0.9980988 1.001796 1.005075 0.9994918 1.009668 1.005490 1.009232 1.003116

$Condition2
            comb_12  comb_13   comb_14  comb_15  comb_23   comb_24  comb_25   comb_34  comb_35  comb_45
Point_Est 0.9982833 1.001201 0.9992053 1.008020 1.001212 0.9998492 1.010451 0.9988859 1.003818 1.002878
Upper_CI  0.9993789 1.002797 0.9992706 1.018785 1.001277 1.0014875 1.028742 1.0010528 1.022219 1.011496

Look forward to hearing your thoughts. Will open a PR and update the vignettes once we've decided the best way forward. I realise some of this may already be in bandle but perhaps not explicitly accessible?

Example code
[1]

## ================================================
## Taken and modified from mcmc_get_outliers
## compute the number of outliers for each iteration for each chain
## ================================================
bandle_get_outliers <- function(params) {
  ## check for validity
  stopifnot(inherits(params, "bandleParams"))
  ## get outliers
  out_cond1 <- lapply(params@chains@chains, function(mc) 
    coda::mcmc(colSums(1 - mc@outlier$cond1)))
  out_cond2 <- lapply(params@chains@chains, function(mc) 
    coda::mcmc(colSums(1 - mc@outlier$cond2)))
  out <- list(cond1 = out_cond1, 
              cond2 = out_cond2)
  return(out)
}

[2]

## ================================================
## Taken and modified from pRoloc TAGM MCMC code
## ================================================

plotOutliers <- function(params, 
                         auto.layout = TRUE) {
  

  ## check for validity
  stopifnot(inherits(params, "bandleParams")) 
  
  ## get outliers
  out <- bandle_get_outliers(params)
  
  ## set plotting parameters
  oldpar <- NULL
  on.exit(par(oldpar))
  if (auto.layout) {
    mfrow <- c(length(out$cond1), 4)
    oldpar <- par(mfrow = mfrow)
  }
  
  ## set names for plots
  chain_id <- seq_along(out$cond1)
  t_names1 <- paste0("Trace - condition ", 1, ", chain ", chain_id)
  d_names1 <- paste0("Density - condition ", 1, ", chain ", chain_id)
  t_names2 <- paste0("Trace - condition ", 2, ", chain ", chain_id)
  d_names2 <- paste0("Density - condition ", 2, ", chain ", chain_id)
  
  ## colours
  cols <- brewer.pal(n = 9, "Set1")[seq_along(chain_id)]
  
  ## use coda for plots
  for (i in seq_along(chain_id)) {
    traceplot(out$cond1[[i]], main = t_names1[i], col = cols[i])
    densplot(out$cond1[[i]], main = d_names1[i], col = cols[i])
    traceplot(out$cond2[[i]], main = t_names2[i], col = cols[i])
    densplot(out$cond2[[i]], main = d_names2[i], col = cols[i])
  }
}

[3]

## ================================================
## Calculate the Gelman diagnostics for all pairwise 
## combinations of chains in the bandle run
## ================================================

calculateGelman <- function(params, all = FALSE) {
  ## get number of conditions (default is 2 but add here in case of pkg developemnt to allow more conditions)
  n_conds <- length(params@chains@chains[[1]]@outlier)
  ## initiate list
  gelmanStats <- vector("list", n_conds)
  ## Look over both conditions (n_conds = 2)
  for (j in seq(n_conds)) {
    ## get MCMC outputs
    outliers_cond <- lapply(params@chains@chains, function(mc)
      coda::mcmc(colSums(1 - mc@outlier[[j]])))
    ## get combinations
    p <- combinations(length(params@chains), 2)
    m <- vector("list", nrow(p))
    for (i in seq(nrow(p))) {
      m[[i]] <- as.vector(p[i, ])
    }
    ## calculate gelmans
    gelmanStats[[j]] <- sapply(m, function(z) 
      gelman.diag(outliers_cond[z], autoburnin = FALSE)$psrf)
    rownames(gelmanStats[[j]]) <- c("Point_Est", "Upper_CI") 
    colnames(gelmanStats[[j]]) <- sapply(m, function(x) 
      paste0("comb_", paste0(x, collapse = "")))
  }
  names(gelmanStats) <- paste0("Condition", 1:2)
  return(gelmanStats)
}

CC @Charl-Hutchings

[ococrook]

Looks good, can you change the sapply to vapply?

[lmsimp]

Will open PR 👍