## Loading required package: coda
## Loading required package: MASS
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2016 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
library(dostats) <- function(x, ...){
  print(ascii(df, include.rownames = FALSE), type = 'rest')


The harvestr package is a new approach to simulation studies that facilitates parallel execution. It builds off the structures available in the plyr, foreach and rsprng packages. What harvestr brings to the picture is abstractions of the process of performing simulation.


The theme of harvestr is that of gardening, stemming from the idea that the pseudo-random numbers generated (RNG) in replicable simulation come from initial states called seeds. Figure 1 shows the basic process for harvestr.

(workflow.png)[The basic harvestr process]

The ideas are simple.

The effect is the results can be taken in any order and independently, and the final results are the same as if each analysis was taken from start to end with setting a single seed for each stream.

Example 1 - Basic Example

Some learn best by example. Here I will show a simple example for the basic process. Here we will perform simple linear regression for 100 data sets. First off we gather the seeds. This step is separate to facilitate storing the seeds to be distributed along with research if necessary.

seeds <- gather(100, seed=12345)

Second, we generate the data.

datasets <- farm(seeds, {
  x <- rnorm(100)
  y <- rnorm(100, mean=x)

Then we analyze the data.

analyses <-  harvest(datasets, lm)

So what do we have in analyses? We have whatever lm returned. In this case we have a list of lm objects containg the results of a linear regression. Ussually we will want to do more to summarize the results.

coefs <- t(sapply(analyses, coef))
adply(coefs,2, dostats, mean, sd)
##            X1        mean        sd
## 1 (Intercept) 0.009484538 0.1030893
## 2           x 1.006747591 0.1007218

Example 2 - Stochastic Analysis

That is very nice, but rather simple as far ananalyses go. What might be more interesting is to perform an analysis with a stochastic component such as Markov Chain Monte Carlo.

posteriors <- harvest(datasets, MCMCregress, formula=y~x)
dataframes <- harvest(posteriors,
X.samples  <- harvest(dataframes, `[[`, "x")
densities  <- harvest(X.samples, density)
l_ply(densities, lines)