SHARE
TWEET

Untitled

a guest Oct 17th, 2019 88 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #BAMM run evaluation function
  2.  
  3. # libraries required: coda, data.table, BAMMtools
  4.  
  5. # if mcmcFile or chainFile are NULL, then it will look for the most likely files
  6.  
  7. # returns effective size and Geweke z-score for the log likelihood and the number of shifts
  8. # as well as frequency of hot/cold chain swapping, if multi-chain analysis.
  9.  
  10. # ideal:
  11. # effective size > 200
  12. # Geweke z-score between -2 and 2
  13. # chain swap frequency between 20 and 60%.
  14.  
  15.  
  16. bammCheck <- function(expectedNumberOfShifts = 1, burnin=0.15, mcmcFile = NULL, chainFile = NULL, plotToPDF = FALSE, title = NULL) {
  17.    
  18.     if (is.null(mcmcFile)) {
  19.         files <- list.files(pattern='mcmc_out')
  20.         if (length(files) == 1) {
  21.             mcmcFile <- files
  22.         } else {
  23.             stop('Please specify an mcmc_out file. Several were detected.')
  24.         }
  25.     }
  26.     if (is.null(chainFile)) {
  27.         files <- list.files(pattern='chain_swap')
  28.         if (length(files) == 1) {
  29.             chainFile <- files
  30.         } else if (length(files) == 0) {
  31.             chainFile <- 'chain_swap'
  32.             cat('\n\tNo chain swap file detected. Assuming single chain analysis.\n')
  33.         } else {
  34.             stop('Please specify an mcmc_out file. Several were detected.')
  35.         }
  36.     }
  37.    
  38.     if (class(mcmcFile) == 'character') {
  39.         mcmc <- data.table::fread(mcmcFile, data.table = FALSE)
  40.     }
  41.     mcmc2 <- mcmc[floor(burnin * nrow(mcmc)):nrow(mcmc), ]
  42.    
  43.     # if BAMM was running with a single chain, there will be no chain swap file.
  44.     if (chainFile %in% list.files()) {
  45.         chainswap <- data.table::fread(chainFile, data.table = FALSE)
  46.         chainswap <- chainswap[(burnin * nrow(chainswap)):nrow(chainswap),]
  47.         chainswap <- chainswap[which(chainswap$rank_1 == 1), ]
  48.         success <- length(which(chainswap$swapAccepted == 1)) / nrow(chainswap)
  49.         success <- round(success, digits = 2)
  50.     }
  51.    
  52.     #get prior distribution of shifts
  53.     obsK <- seq(from=0, to = max(mcmc2[, "N_shifts"]), by = 1)
  54.     prior <- sapply(obsK, prob.k, poissonRatePrior=1 / expectedNumberOfShifts)
  55.     prior <- data.frame(N_shifts = obsK, prob = prior)
  56.    
  57.     #get posterior distribution of shifts
  58.     posterior <- sapply(obsK, function(x) length(which(mcmc2[,'N_shifts'] == x))) / nrow(mcmc2)
  59.     names(posterior) <- obsK
  60.     posterior <- data.frame(N_shifts = names(posterior), prob=posterior)
  61.    
  62.     if (!plotToPDF) {
  63.         dev.new(width = 10, height = 7)
  64.     }
  65.     par(mfrow = c(2, 2), mar = c(4, 4, 3, 3))
  66.     plot(mcmc[, 1], mcmc[, 4], type = 'l', lwd = 0.5, xlab = 'generations', ylab = 'loglik')
  67.     abline(v = max(mcmc[,1]) * burnin, lty = 2)
  68.    
  69.     if (!is.null(title)) {
  70.         title(main = title)
  71.     }
  72.  
  73.     #barplot(table(mcmc2$N_shifts), xlab='n shifts')
  74.     barplot(prior[,2], names.arg = prior[,1], ylim = c(0, max(c(prior[,2], posterior[,2]))), border = 'black', col = 'light blue', xlab = 'n shifts')
  75.     barplot(posterior[,2], add = TRUE, border = 'black', col = BAMMtools::transparentColor('red', 0.4), axes=FALSE)
  76.     legend('topright', legend = c('prior','posterior'), fill = c('light blue', BAMMtools::transparentColor('red', 0.4)), bty = 'n', cex=1.5)
  77.  
  78.     plot(mcmc[,1], mcmc[,2], type='p', pch=20, cex=0.5, xlab='generations', ylab='n Events')
  79.     abline(v=max(mcmc[,1]) * burnin, lty=2)
  80.    
  81.     plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n', xpd=NA)
  82.     text(x=0.5, y=0.9, paste("effective size for log lik: ", round(coda::effectiveSize(mcmc2[,4]), 2), sep=''))
  83.     text(x=0.5, y=0.7, paste("effective size for n Events: ", round(coda::effectiveSize(mcmc2[,2]), 2), sep=''))
  84.     text(x=0.5, y=0.5, paste("Geweke z-score for log lik: ", round(coda::geweke.diag(coda::as.mcmc(mcmc2[,4]))$z, 2), sep=''))
  85.     text(x=0.5, y=0.3, paste("Geweke z-score for n Events: ", round(coda::geweke.diag(coda::as.mcmc(mcmc2[,2]))$z, 2), sep=''))
  86.  
  87.     if (chainFile %in% list.files()) {
  88.         text(x=0.5, y=0.1, paste("successful swap frequency with cold chain: ", success, sep=''))
  89.     }  
  90. }
  91.  
  92. prob.k <- function(k, poissonRatePrior) {
  93.     Denom <- (poissonRatePrior + 1) ^ (k + 1)
  94.     Prob <- poissonRatePrior / Denom
  95.     return(Prob)
  96. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top