Clinician's
corner

Back to main page

Programmer's
corner
So, you use WinBUGS a lot? Want more?
Patrick Blisle
Division of Clinical Epidemiology
McGill University Health Center
Montreal, Quebec CANADA
patrick.belisle@rimuhc.ca

Last modification: 12 oct 2016















Version 1.5 (October 2016)
wb.stats
[R] Scanning node statistics from WinBUGS output
[ wb.stats is an R function developped to read node statistics from a WinBUGS text output file. ]


Menu



Top
Syntax

wb.stats(file, nodes=character(0), node.statistics=c("mean", "sd", "MC error", "2.5%", "median", "97.5%", "start", "sample"), return.as.list=length(node.statistics)>1)


Top
wb.stats arguments

Argument Value Comment
file The full path to WinBUGS text output file from which you wish to scan node statistics. Note that wb.stats expects the WinBUGS output to be in txt; it will NOT work with an odc file.
nodes The list of nodes for which node statistics are to be returned. Argument is case sensitive.
Default is to return node statistics for each node found in file.
node.statistics The list of node statistics to be returned. Argument is case insensitive.
Default is to return all node statistics.
return.as.list Whether the output should be returned as a list or a matrix. Logical.
Default is to return as a list when more than one node statistic is read.


Top
Example

The code below reproduces the input and output of an R session involving calls to wb.stats.

> source('c:/users/pbelisle/My Documents/Home/bin/R/fcts/WBStats.R') # load file where wb.stats R function is defined
> u <- wb.stats('c:/users/pbelisle/My Documents/Home/MyProject/WinBUGS/log/mywinbugsoutput.txt', return.as.list=F) # read (all) the node statistics from the cited WinBUGS text output file
> u

meansdMC error2.5%median97.5%startsample
x.p.null[1]6.5e-016.7e-021.1e-035.2e-016.6e-017.8e-015.0e+032.0e+04
x.p.null[2]5.0e-022.7e-025.7e-046.9e-034.7e-021.1e-015.0e+032.0e+04
x.mu[1]5.5e+002.8e-011.2e-025.0e+005.5e+006.1e+005.0e+032.0e+04
x.mu[2]8.2e+002.4e-013.3e-037.7e+008.2e+008.7e+005.0e+032.0e+04
x.mudiff2.7e+003.6e-011.1e-022.0e+002.7e+003.4e+005.0e+032.0e+04
x.sd.between[1]7.0e-013.8e-012.5e-025.2e-027.2e-011.4e+005.0e+032.0e+04
x.sd.between[2]1.7e+001.9e-012.1e-031.4e+001.7e+002.1e+005.0e+032.0e+04
x.sd.within[1]7.5e-013.3e-011.6e-021.2e-017.6e-011.4e+005.0e+032.0e+04
x.sd.within[2]9.3e-018.9e-021.3e-037.8e-019.3e-011.1e+005.0e+032.0e+04
prev8.2e-016.2e-021.5e-036.9e-018.3e-019.3e-015.0e+032.0e+04


# Read the same WinBUGS output file as above, but select only a subset of nodes found therein
u <- wb.stats('c:/users/pbelisle/My Documents/Home/MyProject/WinBUGS/log/mywinbugsoutput.txt', nodes=c('x.mudiff', 'prev'), return.as.list=F)
> u

meansdMC error2.5%median97.5%startsample
x.mudiff2.7e+003.6e-011.1e-022.0e+002.7e+003.4e+005.0e+032.0e+04
prev8.2e-016.2e-021.5e-036.9e-018.3e-019.3e-015.0e+032.0e+04


# Read the same WinBUGS output file as above, but select only a subset of nodes found therein and only a few node statistics AND RETURN OUTPUT AS A LIST
u <- wb.stats('c:/users/pbelisle/My Documents/Home/MyProject/WinBUGS/log/mywinbugsoutput.txt', nodes=c('prev', 'x.mu'), node.statistics=c('mean', 'median'))
> u
$prev
$prev$mean
[1] 0.82

$prev$median
[1] 0.83


$x.mu
$x.mu$mean
x.mu[1] x.mu[2]
5.5 8.2

$x.mu$median
x.mu[1] x.mu[2]
5.5 8.2

Note that the first line of code in the block above reads in the wb.stats function code, if it was previously saved to the file referenced by that command. Alternatively, you can cut and paste the function (see code in section below) directly into R. Of course, if you save your session, this step needs to be done only the first time you use wb.stats.


Top
Code
wb.stats <- function(file, nodes=character(0), node.statistics=c("mean", "sd", "MC error", "2.5%", "median", "97.5%", "start", "sample"), return.as.list=length(node.statistics)>1)
{
   wblog <- scan(file, sep='\n', what="", quiet=T)
   tab.found <- grep('\t', wblog)
   wblog <- wblog[tab.found]

   # Remove node stats header line(s) -- if found twice, ignore lines below second occurence
   median.found <- grep('\\smedian', wblog, perl=T)
   if (length(median.found) > 1) wblog <- wblog[seq(median.found[2]-1)]
   wblog <- wblog[-median.found[1]]

   # Exhaustive list of node stats returned by WinBUGS
   wb.node.statistics <- c("mean", "sd", "MC error", "2.5%", "median", "97.5%", "start", "sample")

   # remove lines with comments (possible if output was modified)
   found.comment <- grep('#', wblog)
   if (length(found.comment)) wblog <- wblog[-found.comment]

   wblog <- gsub('\\s+', ' ', wblog, perl=T)
   wblog <- sub('\\s', '', wblog, perl=T) # remove leading space
   wblog <- wblog[nchar(wblog)>0] # drop empty lines
     # drop lines with /odds ratio/
     slashes <- grep("/", wblog)
     if (length(slashes)) wblog <- wblog[-slashes]
   wblog <- matrix(unlist(strsplit(wblog, ' ')), byrow=T, nrow=length(wblog))
   nodes.names <- wblog[,1]
   wblog <- matrix(suppressWarnings(as.numeric(wblog[,-1])), nrow=nrow(wblog))
   colnames(wblog) <- wb.node.statistics
   rownames(wblog) <- nodes.names

   # Keep only required node statistics
   w <- match(tolower(node.statistics), tolower(wb.node.statistics))
   wblog <- wblog[,w,drop=F]

   if (length(nodes))
   {
     nodes.names0 <- sub('\\[.*', '', nodes.names, perl=T)
     w <- !is.na(match(nodes.names0, nodes)) | !is.na(match(nodes.names, nodes))
     wblog <- wblog[w,,drop=F]
   }


   if (return.as.list)
   {
     stats <- sort(colnames(wblog))
     n.stats <- length(stats)

     dim <- rownames(wblog)
     dim <- sub("\\[.*", "", dim, perl=T)
     dimensions <- unique(sort(dim))
     n.dim <- length(dimensions)

     out <- list()

     for (i in seq(n.dim))
     {
       d <- dimensions[i]
       w <- which(dim==d)

       my.out <- list()
       my.i <- 0
       for (s in stats)
       {
         my.i <- my.i + 1
         my.out[[my.i]] <- wblog[w, s]
       }

       if (n.stats == 1)
       {
         my.out <- unlist(my.out)
       }
       else
       {
         names(my.out) <- stats
       }

       out[[i]] <- my.out
     }

     names(out) <- dimensions
   }
   else
   {
     out <- wblog
   }

   out
} # end of wb.stats


Top
Download

wb.stats is a free R function. Download version 1.5 now.