Introduction

JAGS (“Just Another Gibbs Sampler”; Plummer, 2003) is a free software package for the analysis of Bayesian models. JAGS uses a suite of Markov chain Monte Carlo methods—that is, general-purpose stochastic simulation methods—to draw samples from the joint posterior distribution of the parameters of a Bayesian model. These samples can then be used to draw inferences regarding the model parameters and model fit. JAGS uses a dialect of the BUGS language (Lunn, Jackson, Best, Thomas, & Spiegelhalter, 2012; Thomas, Spiegelhalter, & Gilks, 1992) to express directed acyclic graphs, a mathematical formalism to define joint densities. What makes JAGS particularly interesting, however, is that it is designed to be extensible with user-defined functions, monitors, distributions, and samplers.Footnote 1

The goal of the present article is to provide a how-to guide to writing custom extensions for JAGS. Although JAGS is envisioned as a flexible and extensible (modular) framework and contains an infrastructure for customization, at the time of writing, no tutorials, technical manuals, or other sources on how to do so were available.

In the following section, we describe in detail how custom distributions may be implemented in JAGS. This section assumes some knowledge of the C++ programming language (and of object-oriented programming) on the part of the reader. Throughout, we will use the Bernoulli distribution as a didactic example, but JAGS can accommodate all other types of probability distributions (continuous and discrete, univariate and multivariate) with little extra effort. In the second part of the article, we present a module implementing a distribution of particular interest to the cognitive science community. The JAGS Wiener module (JWM) adds the first-passage time distribution of a drift diffusion process to JAGS. We also provide two sanity checks of the JWM: an extensive simulation study, showing good recovery, and an application to a previously analyzed data set, showing parallel results.

Steps to extending JAGS with a new module

In this section, we detail the steps required to write a new module for JAGS. JAGS modules are library files that are located in the /modules directory of JAGS. These modules can be loaded during JAGS runtime in order to extend its functionality with more functions and distributions, or even sampling algorithms or monitors.

Modules are written in C++. In order to write a custom module that provides a new distribution, we need to define two new C++ classes: one for the module itself, and one for the distribution that we will use. Throughout, we will display much of the required code, as well as templates for files used in building and compiling. The code can be copied from this text, or it can be found in our SourceForge archive (and can be used as a template for new modules). The Appendix contains a quick reference table of the required steps.

We use the Bernoulli distribution as an example to illustrate the basics of extending JAGS because the functions that define this distribution are relatively easy to write, without the need for calling advanced functions from extra libraries.Footnote 2 Using the naming schemeFootnote 3 from JAGS, we will name our module class BERNModule and the distribution class DBern.

Step 1: Defining a module class

We start by creating a specific module class that is a child of the base Module class, which JAGS provides. To make a new module class, create a new file like the one displayed in Box 1—this is our Bernoulli.cc module file. Place it in src/, a subfolder of your working directory for this project.

Box 1 The Bernoulli.cc module definition file

figure a

In the first lines, we reference two header files containing premade functions that we will use. Apart from the JAGS header file Module.h (which is available through JAGS), we will include the Bernoulli distribution class file, which we will create at a later point.

The new functions and classes that we create need to be defined within a single namespace—a named environment that is used to group programmatic entities that are frequently used together. The third line of Box 1 names the namespace.

Then, we define the module class, which we will call BERNModule. BERNModule is a subclass of Module with public inheritance, and needs to define two functions—both are required.

Constructor and destructor for the BERNModule class In Box 1, the constructor function BERNModule() is called whenever an object of our Bernoulli class is instantiated.

The constructor function BERNModule() needs to do two things. First, it needs to instantiate a generic module by calling the constructor function of the parent Module class. It calls that generic constructor with a single argument: the name of the module (as a string). Second, our new module’s constructor needs to call the insert function with, as input, a new instance of the class describing the new distribution (in our case, DBern, which is detailed below).

The destructor function ~BERNModule() needs to remove the instance. To do so, we use a for loop over all instances of the module.Footnote 4

The last line of Box 1 actually calls the new module’s constructor function and instantiates one BERNModule.

The module class BERNModule is now complete. The second step in creating the JAGS module is to define the DBern scalar distribution, which is inserted with the constructor call in Box 1.

Step 2: Defining a scalar distribution

A new scalar distribution for JAGS is implemented through a new C++ class. The new class will need, in addition to its own constructor function, the following four specific functions:

  • logDensity, a function to calculate and return the log-density;

  • randomSample, a function to draw and return random samples;

  • typicalValue, a function to return typical values; and

  • checkParameterValue, a function to check whether the given parameter values are in the allowed parameter space

We now create two new files in a subdirectory called src/distributions/. Our files will be src/distributions/DBern.h, which will contain the class prototype, and src/distributions/DBern.cc, which will contain the actual computational implementation.

DBern.h Box 2 contains our file src/distributions/DBern.h. Near the top of the code, we include the (JAGS) parent class ScalarDist in order to be able to use it as a base class and inherit from it. Since we will be adding this to the namespace of our new module, we use the same name for it here as in Box 1. What follows in Box 2 are the prototypes of (1) the constructor function (which has the same name as the class), (2) the four required functions (which will be identical for other scalar distributionsFootnote 5), and (3) the function isDiscreteValued, which will tell JAGS that the function has a discrete domain.

Box 2 The Bernoulli scalar distribution class header file DBern.h

figure b

The logDensity function takes five input arguments:

  1. 1.

    The first argument (double x) contains the data point at which the function is to be calculated.

  2. 2.

    The second argument, of the JAGS-specific type PDFType (defined in Distribution.h, which is a parent of ScalarDist.h), defines different ways of calculating the log density. During JAGS runtime, depending on how the node is used, this argument can take different values. The possible values are PDF_FULL (for when the node is used for the full evaluation of the likelihood, when the parameters and sampled values are not constant), PDF_PRIOR (for when parameters are constant, so that terms that depend on them can be omitted from the calculations), and PDF_LIKELIHOOD (used when the sampled value is constant, so that terms that depend on it can be omitted from the calculations). The PDFType can simply be ignored (as it is in our example).

  3. 3.

    The third argument (std::vector<double const *> const &parameters) is a reference to a vector that contains pointers for the parameters of the distribution. In our case, the vector consists of a pointer to pi (the probability parameter of the Bernoulli distribution; see Eq. 1).

  4. 4.

    The final two arguments (double const *lower and double const *upper) are pointers to the lower and upper boundaries of the distribution, in the case of truncated sampling (see lib/graph/ScalarStochasticNode.cc in the JAGS code to find the call to the logDensity function). In our example, these are not used, and these arguments are ignored.

The arguments used in the other functions have similar definitions. One final input argument, which appears only in the randomSample function, is the JAGS-specific RNG type argument (RNG *rng). It points to an rng object that provides the following functions to generate random numbers:

  • rng- > uniform()

  • rng- > exponential()

  • rng- > normal()

These random-number-generating functions can be used to drive samplers for any distribution.

DBern.cc Now we need to implement the four functions and the constructor function prototyped in Box 2. The distribution function of a Bernoulli distributed random variable X, used as the basis for the computations, is in Eq. 1.

$$ \begin{array}{l}P\left(X=0\right)=\left(1-\pi \right)\hfill \\ {}P\kern-0.44em \left(X=1\right)\kern-0.46em =\pi \hfill \end{array} $$
(1)

Note that, although these four functions need to be present in the code, not all are strictly speaking required to run an analysis. If drawing random samples from the distribution is not a functionality your analyses require, it is possible to let randomSample return JAGS_NAN. This will cause JAGS to throw an error if a model is defined that requires random sampling from the new distribution, but it does not affect using the new node as a likelihood. Similarly, if typicalValue is not implemented, JAGS will not run a sampling process if no starting values are supplied, but will work otherwise. The checkParameterValue function is a prerequisite for logDensity, however.

Furthermore, as the Bernoulli distribution is a discrete-valued distribution, we redefine the function isDiscreteValued and let it return true (implementation of that particular function is not necessary for continuous distributions: it will be inherited from the base class and return false by default).

The implementations of the required functions are provided, with comments, in Boxes 3, 4, and 5.

Box 3 The DBern.cc file. Note that we need to include rng/RNG.h and util/nainf.h from the JAGS library, to provide the RNG struct and the JAGS_* constants, as well as the jags_* functions. cmath.h is needed for standard math operations.

figure c

Box 4 Functions from the DBern.cc file

figure d

Box 5 Further functions from the DBern.cc file

figure e

Constructor and destructor for the DBern class Finally, we need to set up constructor and destructor functions for our ScalarDist class. The constructor function calls the constructor for ScalarDist with three input arguments (as defined in the JAGS header file distribution/ScalarDist.h). The first input argument is std::string const &name and is used to give the distribution node a name, which can later be used in the JAGS model file to define a node with this distribution (in this example, the distribution node will be called “dbern2”). The second input argument of the constructor is unsigned int npar and is used to define the number of arguments that the new distribution node can take (that is one parameter, the probability π, for this example). The third argument of the constructor is Support support, a JAGS-specific variable that defines the support of our new distribution. In our case, support is DIST_PROPORTION, indicating a distribution that spans from 0 to 1. Other valid options for this argument are DIST_POSITIVE (for a distribution that has support only on the positive real half-line), DIST_UNBOUNDED (for a distribution that spans the entire real line), and DIST_SPECIAL (for other domains).

Step 3: Building the module

To configure and build our module, we used the GNU Autotools and libtoolFootnote 6 (these are the same tools used in building the JAGS project), according to the steps below.

We will give a brief overview, together with all the information necessary to organize the build process of the module described in this article. However, these building methods are not guaranteed to stay the same over time,Footnote 7 so the build process may be slightly different on newer machines and operating systems. However, build processes are central to software development and are typically very well-documented for new architectures and operating systems.

Step 3.1: Creating a file structure

First, make sure that a project directory exists that contains within it a subdirectory called src/, and that all of the source files (which contain the C++ code for the module) are present in the src/ subdirectory. For the given example, this would be

  • the module main file src/Bernoulli.cc, and

  • in the src/distributions/ sub-subdirectory, the distribution class files DBern.h and DBern.cc

Now create one configure.ac file in the main directory of the project and one Makefile.am in every directory of the project (including all subdirectories). The configure.ac file will contain information on how to correctly configure the project. The Makefile.am files will contain the necessary arguments for building, including libraries to link to and the paths to the include directories (and more).

Since the newly created module uses functions from the JAGS library, it needs to be linked to the appropriate libraries, and JAGS’s include files need to be available (i.e., JAGS needs to be installed properly).

Step 3.2: Creating the configure.ac file

Boxes 6 and 7 together show the configure.ac file. All of the functions used in this file are documented in the documentation of Autoconf (see note 6). This file contains instructions on how to prepare the module for building (e.g., which compiler to use). Much in the file can remain unchanged.

Box 6 The first part of the configure.ac file (file continues in Box 7). Lines beginning with dnl (“delete to new line”) denote comments.

figure f

Box 7 The second part of the configure.ac file (continued from Box 6)

figure g

The AC_INIT directive takes as its first argument the package name, as the second argument the package version, as the third argument a contact e-mail address (for bug reports), and as the fourth argument the name for the .tar file. JAGS_MAJOR and JAGS_MINOR should be edited to the current JAGS version that is being used. AC_CONFIG_SRCDIR should contain a path to any file of the package. The AC_CONFIG_FILES directive at the end tells the configuration process where to create Makefiles, using the Makefile.am files as a template.Footnote 8 The libltdl directory and its Makefile.am will be created automatically. If one follows all directions given here, nothing else has to be changed.

Step 3.3: Creating several Makefile.am files

Box 8 shows the Makefile.am file in the main directory. The configuration process will use a subdirectory called m4/, which needs to be created manually. The file also needs to reference the building subdirectories: the src/ directory containing the actual source, and a libtool/ directory, which will be created automatically.

Box 8 The Makefile.am file

figure h

Box 9 shows the Makefile.am file in the src/ directory. It contains the instructions to produce the module library.

Box 9 The src/Makefile.am file

figure i

The first directive, SUBDIRS, points to the subdirectories to be used; for JAGS modules, this will typically be distributions. The next directive defines the JAGSMODULE libtool library to be created. The library file will be named <longname>.la (e.g., Bernoulli.la). The next directive defines where to find the library source code. The following directive gives the C++ compiler additional options. The -I path option passes to the compiler the path of JAGS’s include directory (a subdirectory of the system’s include path indicated by $(includedir)).

The next set of directives tell the compiler which object code (i.e., compiled code) should be linked to the libraries. For the Bernoulli module, this is the object code for the Bernoulli distribution and the JAGS library. The -ljags argument searches for the JAGS library in the system library directory (note that since Windows organizes its libraries differently, and the JAGS library will be in a slightly different location, we need to add a Windows-specific part here). Had we used additional external libraries, we could add these here as well. For example, we could link to the JRmath.h library, with -ljrmath (and -ljrmath-0 in the Windows part).

Box 10 shows the Makefile.am file of the src/distributions/ directory. It contains the instructions to build the Bernoulli distribution code, which will then be included in the package library. The sublibrary distributions/Bernoullidist.la is named on line 1 of Box 10. This name must be matched in the src/Makefile.am (line 9 of Box 9). In Box 10, the directive noinst_LTLIBRARIES will cause the sublibrary to be built, but not installed (i.e., the code will be compiled but not copied to a system location). The sublibrary will be incorporated in the module library by src/Makefile.am. As before, the Bernoullidist_la_CPPFLAGS directive takes arguments to be passed to the compiler. In this case, it passes the include directories where the header files for the code can be found.

Box 10 The src/distributions/Makefile.am file

figure j

With all of the files organized as described in this section, we can now proceed to configure and build the module using the steps below. We first provide the steps for building on a Unix-like (i.e., Mac or Linux) environment, in which building is much easier than under Windows systems. Unfortunately, building on Windows systems is a tedious and involved affair. In particular, the autoreconf -fvi command on a Windows system is not possible with standard emulators like MinGW. It is, however, possible to perform the first two steps below on a Mac/Linux system, then to create a source .tar file, copy and extract that onto a Windows machine, and use the MinGW environment with msys to configure and build the source. We will outline the full building process for Mac/Linux systems, and then describe the steps needed to compile for Windows systems.

Step 3.4a: Building the module ( Mac / Linux )

To proceed in a Mac/Linux system, follow these steps:

  • autoreconf -fvi: This command will generate a number of auxiliary files that are necessary for the configuring and building process.

  • ./configure: This command configures the source package for building on your system.

  • make: This command compiles the source code into system-specific object code.

  • make install: This command creates local copies of the object code, placing it in the correct locations in the file system (this step typically requires administrator/superuser privileges). In our case, this command will copy the module library to an appropriate location in the system where JAGS can find and load it.Unfortunately, this command does not work under Windows. It is possible, however, to create installers for Windows machines that include precompiled binaries and to copy them to the correct locations. One easy-to-use third-party tool to create such installers is NSIS (http://nsis.sourceforge.net).Footnote 9 Alternatively, this can be done manually by copying the libraries to the JAGS directory that contains the other module libraries.

Additionally, after running the first two commands (the configuration process), one can easily create source .tar files with make dist-gzip or make dist-bzip2 (or whichever format one prefers and is supported by the Makefile routines).

Step 3.4b: Building the module (Windows)

Building under Windows follows a similar process, with some added steps. First, create the Windows libraries. We start by installing MinGW (the MinGW installer includes msys) and the TDM-GCC Compiler Suite, which can be obtained from www.mingw.org and http://tdm-gcc.tdragon.net.Footnote 10 In the rest of this paragraph, it will be assumed that MinGW and TDM-GCC are installed in their default directories.

Second, delete all *.dll.a files in the TDM-GCC directory, to force the compiler to link to the static libraries (the *.dll files). This is necessary to build libraries that will work on systems that do not have TDM-GCC.

Third, change the path in the file C:\mingw\msys\1.0\etc.\fstab from C:\mingw /mingw to C:\MinGW64 /mingw.Footnote 11 This is necessary in order to use the TDM-GCC compilers instead of the standard MinGW compilers.

Now, the actual configuration and building process can commence. The module needs to locate the JAGS include files and the JAGS libraries. Because Windows has no standard path in which to look for these files, this needs to be done manually. Edit the code in Box 11 to reflect the paths to the JAGS libraries (with the -L option) and the JAGS include files (with the -I option). The code is given for a standard installation of JAGS 3.3.0.

Box 11 Building the source code under Windows. Note that if you build libraries for both 32-bit and 64-bit systems, you need to run make clean between the two building processes.

figure k

Start msys (the MinGW shell), extract the source .tar file in any directory, navigate to that directory by using the cd command, and run the appropriate commands given in Box 11. Executing these commands will create a number of files in the src/.libs/ directory. Now copy those files that are named after your module and end in .dll, .dll.a, and .la to your JAGS modules directory (where all of the other module libraries are located). You can also use these files to create an installer (e.g., with NSIS; see note 9).

After installing the new module, JAGS will recognize a new stochastic node, to be used as follows: x ~ dbern2(pi), where the string “dbern2” is defined in Box 3 (line 23) and its parameter is interpreted through line 13.

Extending the module with logical nodes

Adding a custom distribution is only one of several possible extensions that a JAGS module can provide. One other possible extension is the creation of logical nodes (i.e., deterministic functions, rather than stochastic ones; logical nodes are useful for transforming variables within a JAGS model file, and are always preceded by an assignment operator <−).

Here we demonstrate how to supplement our Module example with a logical node that calculates the Bernoulli log density at a given data point for a given parameter set.Footnote 12 To do so, we will write a scalar function. Much like the scalar distribution, the ScalarFunction class requires the implementation of two functions:

  • evaluate: the function that does the calculations and returns the result

  • checkParameterValue: a function to check whether the parameter values are in the domain of the function

In the src/ directory of our module, add a subdirectory functions/, and in it the two files LogBernFun.h and LogBernFun.cc. Boxes 12 and 13 show examples of such files.

Box 12 The src/functions/LogBernFun.h file

figure l

Box 13 The src/functions/LogBernFun.cc file

figure m

After the function class is ready, we need to add it to the Bernoulli.cc module by loading it with the constructor and deleting it with the destructor (as in Box 1). See Box 14 for the lines to add.

Box 14 The code to add to the src/Bernoulli.cc file

figure n

Finally, add an appropriate Makefile.am to the newly created src/functions/ directory. In src/functions/Makefile.am, the sublibrary is named (here, Bernoullifunc.la). Also add the name of that sublibrary to src/Makefile.am like this: bernoulli_la_LIBADD += functions/Bernoulifunc.la (in order to incorporate the sublibrary in the common Bernoulli library). Additionally, add the name of the subdirectory (/functions) to the SUBDIRS directive in the first line of src/Makefile.am, and add the appropriate directives to generate the Makefile file in the configure.ac file. The required Makefile.am will be similar to that in Box 10, with only small edits (replacing the file names and the name for the sublibrary in the first line and in the following lines, where it is part of the directive). In the configure.ac file, merely add a line to the last directive to tell Autoconf to create a Makefile in the src/functions/ directory.

Now, once the module is installed, JAGS will recognize a new logical node that can be used as follows: p <− logbern(x,pi), where the name of the node is defined in Box 13 (line 18), and the order of the parameters is defined in lines 10 and 11. This line, when used in a JAGS model file, will store in the variable p the log likelihood under a Bernoulli distribution of data point x given parameter pi. Storing the log likelihood at every iteration of the Gibbs sampler can be useful for convergence checks or to compute model fit statistics in the Bayesian framework.

Note that still more ways exist to extend JAGS. One can, for example, write custom sampling algorithms, or multivariate distributions—fully describing the installation procedure of every possible extension would be beyond the scope of a single article. However, the open-source nature of JAGS allows enterprising researchers to study the implementation of various components of the programFootnote 13 (we recommend the multivariate normal distribution, DMNorm.cc, as an example of a multivariate distribution, and the code for the slice sampler, Slicer.cc, as a sampler example). The generic method of creating modules presented here will work for more sophisticated components as well. The main difference between the workflow outlined here and that for a new component will be in the time that needs to be invested in studying the JAGS framework for the appropriate function class.

The next section focuses on the installation—from an end-user point of view—of a custom stochastic node that we believe to be of use for the cognitive science community.

The JAGS Wiener module (JWM)

Few cognitive models have had as much success as the diffusion model for two-choice response times (see Wagenmakers, 2009, for a review). Recently, research actually applying the diffusion model has greatly increased, since practical, user-friendly software has become available (Vandekerckhove & Tuerlinckx, 2007; Voss & Voss, 2007; Wagenmakers, van der Maas, & Grasman, 2007; Wiecki, Sofer, & Frank, 2011). The extension into hierarchical Bayesian modeling provides the most flexible modeling framework yet (Vandekerckhove, Tuerlinckx, & Lee, 2011). However, flexible implementations of the hierarchical diffusion model (HDM) are currently limited to WinBUGS (Lunn, Thomas, Best, & Spiegelhalter, 2000), which is not truly free software (although it is gratis Footnote 14) and only runs natively on Windows systems. The present article hence addresses a collateral need by providing an HDM implementation in an open-source, platform-independent framework.

The JWM is an extension for JAGS and is designed to integrate seamlessly with the existing JAGS platform. It extends JAGS to recognize dwiener, the first-passage time density of a drift diffusion process, as a new stochastic node with four parameters: the boundary separation α, the nondecision time τ, the initial bias β, and the drift rate δ. The psychological interpretations of the four Wiener distribution parameters are summarized in Table 1, but for a more detailed description of the assumptions of the Wiener diffusion model, we refer the reader to Vandekerckhove et al. (2011).

Table 1 The four main parameters of the Wiener diffusion model, with their substantive interpretations

Installing and using the module is straightforward, as we describe below.

Installation

In order to install and run the JWM, a recent version of JAGS is required.Footnote 15 JAGS can be freely downloaded from http://mcmc-jags.sourceforge.net/, but is also available as a package for various Linux distributions that use the RPM Package Manager, as well as Debian and Ubuntu. For the purposes of the installation procedure, we will assume that a recent version of JAGS is already available.

Because JAGS is designed with the capacity for extension in mind and is capable of dynamically loading libraries, installing a module does not compromise the regular functioning of JAGS. To activate the extensions, modules have to be loaded explicitly by JAGS. To allow the program to locate the new libraries, they need to be installed according to the instructions below.

We have split up the installation instructions between a procedure for Windows systems and one for Linux and Mac systems.

Windows systems

We provide precompiled libraries for Windows that can be readily downloaded from the SourceForge repository. The installation proceeds according to the following steps:

  1. 1.

    Obtain the Windows installer from SourceForge, at https://sourceforge.net/projects/jags-wiener/files.

  2. 2.

    Execute the installer, and make sure to save the libraries at the correct position so that JAGS will be able to find the module. The correct location will be C:\Program Files\JAGS\JAGS-3.3.0\i386\modules\ (for JAGS Version 3.3.0; the version number in the path changes for older or newer releases) on 32-bit systems, and C:\Program Files\JAGS\JAGS-3.3.0\x64\modules\ (for JAGS Version 3.3.0) on 64-bit systems (possibly replacing the JAGS root directory with its actual installation directory). The installer just needs the JAGS root directory to install the libraries at the correct position.

Linux and Mac systems, compiling from source

The general way to install our module on Linux and Mac systems is to compile the source for your system and then to install them. For this operation, some general knowledge of GNU Tools and using a console interface (i.e., the terminal) is required. The following instructions are Linux- and Mac-specific. The installation instructions will assume that JAGS is already installed, so that it is possible to link correctly to the JAGS library. (For those interested in compiling under Windows, please note that the MinGW environment with mysys needs to be installed, as well as the TDM-GCC compiler suite. Please see the previous instructions in the Step 3: Building the Module subsection and the README.win file, present in the win source directory, for further instructions on compiling under Windows.)

  1. 1.

    Get the JWM tar archive, containing the source code, from SourceForge at https://sourceforge.net/projects/jags-wiener/files, and extract it.Footnote 16

  2. 2.

    Use cd to change to the directory containing the source code.

  3. 3.

    If you are compiling the source from a cloned repository, run autoreconf -fvi before configuring (not necessary if you use the stable .tar file release).

  4. 4.

    Configure and compile the source code for your system with the following commands in a terminal window: ./configure && make. When this is done, install the libraries on your system with the following command, which usually requires root privileges: sudo make install.

The JWM should now be installed on your system.

Ubuntu Linux, with the Advanced Packaging Tool

If you use Ubuntu Linux, you can alternatively install the JWM with the Advanced Packaging Tool (APT). The authors maintain an online repository (a personal package archive, or PPA, called cidlab) from which the module can be downloaded. Adding the PPA, updating APT, and installing the module are achieved with this code:

  • sudo apt-add-repository ppa:cidlab/jwm

  • sudo apt-get update

  • sudo apt-get install jags-wiener-module

Testing the installation

The successful installation can be confirmed by loading the Wiener module in JAGS. In order to do that, bring up a terminal or command window and start JAGS. Then type the following at the JAGS interface:

  • $ load wiener

  • Loading module: wiener: ok

Alternatively, R users can use the rjags package to interface with JAGS. Loading the module in R works like this:

  • > library(rjags)

  • Loading required package: coda

  • Loading required package: lattice

  • linking to JAGS 3.3.0

  • module basemod loaded

  • module bugs loaded

  • > load.module(“wiener”)

  • module wiener loaded

We have also made available a version of the MATJAGS MATLAB-to-JAGS interface that loads extra modules, as specified. The MATJAGS interface is a single MATLAB file that can be downloaded from the SourceForge repository. Documentation is provided within that file.

Using the JAGS Wiener module

With the JWM installed, the dwiener stochastic node is now available for use in a model definition. To define that a certain node is distributed according to a Wiener distribution with the four parameters given in Table 1, the following stochastic node can be used:

  • x ~ dwiener(alpha,tau,beta,delta)

This syntax is almost identical to that used in the wiener.odc extension to WinBUGS, as presented in the supplemental material to Vandekerckhove et al. (2011).

The main differences to that implementation are that (a) the diffusion coefficient s is set to 1 in our JAGS implementation (rather than 0.1 in the WinBUGS one), and (b) the JAGS implementation takes as a third input argument the relative bias parameter β (rather than the absolute starting point ζ 0 = αβ in the WinBUGS version). Setting the diffusion coefficient s to 1 instead of 0.1 is computationally more efficient and further yields a more natural interpretation of the drift rate parameter (because it now has a range similar to that of a standard normally distributed variate). To convert to a different evidence scale (i.e., a different value for s), multiply the obtainedFootnote 17 drift rate and boundary separation parameters by s. We use the β parameter instead of the ζ 0 starting point primarily because β has an interpretation on an absolute scale (the unit scale), whereas ζ 0 can only be interpreted relative to the boundary separation α.

Some of the parameters have limited domains: α, τ > 0 and 0 < β < 1. Note that the distribution is implemented as a univariate distribution. To use the distribution, choice response time data should be coded in such a way that error response times (or whichever response type is associated with the lower boundary) are given negative values. Specifically,

$$ x=\left\{\begin{array}{ll}\mathrm{RT}\hfill & \mathrm{if}\ \mathrm{correct}\hfill \\ {}-\mathrm{RT}\hfill & \mathrm{if}\ \mathrm{error}\hfill \end{array}\right. $$

For a worked example on how the module can be used, see the examples directory in the JWM repository on SourceForge, or download a zip file of the example from SourceForge at https://sourceforge.net/projects/jags-wiener/files.

Testing of the JAGS Wiener module

In this section, we report results of some of the extensive testing of our JAGS Wiener module, which we created and built according to the instructions above. The first part describes results from a numerical simulation experiment, whereas the second part describes an application to a benchmark data set.

Parameter recovery simulations

To affirm the accuracy of our implementation of the Wiener distribution, we ran a comprehensive numerical experiment. Using the DMAT toolbox (Vandekerckhove & Tuerlinckx, 2008), we generated data from known parameter sets and then used JAGS to recover the parameters. In all cases, we generated 1,000 data sets with three conditions, which differed in their drift rates only. The three drift rates were always (−3.0, 1.0, 4.0). The boundary separation was either 0.6, 1.2, or 3.0, the bias was either 0.3, 0.5, or 0.7, and the nondecision time was always 0.4 s, yielding a total of nine distinct parameter sets. Using these parameters, we generated data sets with 200 data points per condition. We then used JAGS to estimate the parameters. Our estimate was the posterior mean obtained after running four chains for 2,000 iterations and discarding the first 1,000 as burn-in and using no thinning, for a total of 4,000 samples.

In only one case (0.01 %) was the \( \widehat{R} \) < 1.05 criterion for convergence (see Gelman, Carlin, Stern, & Rubin, 2004) violated. This case was discarded. In general, recovery was good. Figure 1 shows the distributions of each parameter’s estimates by parameter set. Though some parameters in some conditions—particularly the drift rate parameter—show some variability, no systematic biases are evident for any parameter.

Fig. 1
figure 1

Results of the parameter recovery simulations. In all panels, the nine levels refer to the nine different parameter sets. The true parameter values are given in the text. Values are recovered as posterior means, and the plots display the population of parameter sets over 1,000 simulated data sets. Circles indicate the mean recovered values, whiskers span 1.5 times the interquartile range, and the edges of the boxes are the 25th and 75th percentiles. The results indicate very good recovery, with slightly more variability in the drift rate estimates than in the other parameters, especially in the conditions with low boundary separations

Benchmark data

Data set and theoretical framework

To further illustrate the functionality and utility of our module, we applied a Bayesian hierarchical diffusion model analysis to benchmark data from Vandekerckhove, Panis, and Wagemans (2007, data used with permission). The data came from nine participants, who performed a visual detection task (i.e., they reported whether or not a change occurred in a figure, in a temporal two-alternative forced choice task). The difficulty of this task was manipulated in a 2 × 2 factorial design, resulting in four experimental conditions plus one control condition (in which no change in the figure occurred). The dependent variables X (binary: 1 if the correct response was given, 0 otherwise) and T (the response time) were recoded into a single variable Y = (2X − 1)T, preserving all the information in the bivariate data. For a more detailed description on the data and the research question, we refer the reader to Vandekerckhove et al. (2007).

Vandekerckhove et al. (2011) have previously used the same data to show the advantages of a hierarchical diffusion model in a Bayesian framework. Here the focus lay on demonstrating how we used the same theoretical framework with the presented, open-source software. Moreover, we did our analysis with two different models that differed in the exact specification of the hierarchical structure. For a more detailed description of the model framework, see Vandekerckhove et al. (2011).

Model definitions

Following Vandekerckhove et al. (2011), we assumed an unbiased diffusion process and set β = .5 accordingly. We further allowed the boundary separation parameter α to differ between persons as a random effect. In other words, every person p was allowed to have their own boundary separation parameter α (p), which was a drawn from a joint population distribution: α (p)N(μ α , σ 2 α ).

The nondecision time τ (pij) was also allowed to differ between persons p, conditions i, and trials j. It was a random variable with mean θ (p) and variance χ 2 (p), both of which were themselves seen as random variables, differing between persons, with respective population means μ χ and μ θ and variances σ 2 χ and σ 2 θ : τ (pij)N(θ (p), χ 2 (p)) and θ (p)N(μ θ , σ 2 θ ), and χ (p)N(μ χ , σ 2 χ ).

In the first model, M 1, we did not estimate μ χ and σ χ , but set them to the fixed values μ χ = .35 and σ χ = .125, casting these as fixed rather than random effects.Footnote 18

The drift rate parameter δ (pij) was also allowed to differ between trials, conditions, and persons. Moreover, it was cast as a random variable with a condition-by-person-specific mean ν (pi) and person-specific variance η 2 (p). The mean ν (pi) was itself a random parameter with mean μ ν(i) and variance σ 2 ν(i), which both differed between conditions as fixed effects. The variance η 2 (p) differed between persons and was a random variable with the population mean μ η and variance σ η . In the M 1 model, we did not estimate μ η and σ η , but set them to 3.5 and 3.5, respectively.

The second model (M 2) differed from the first (M 1) only in that μ χ , σ χ , μ η , and σ η were free parameters whose values were estimated by JAGS. Figure 2 shows a graphical model representation of the second model. A graphical representation of the first model would look identical, but for the removal of the four now-fixed nodes.

Fig. 2
figure 2

Graphical model for the M 2 model. The nodes that are indicated with stars and connected with dashed arrows—μ χ , σ χ , μ η , and σ η —are set to fixed values in the M 1 model

Results

We ran three chains with 50,000 iterations each for M 1, and three chains with 250,000 iterations each for M 2. For no parameters did the potential scale reduction factor \( \widehat{R} \) exceed 1.02 in M 1. In M 2, μ η and σ η exceeded acceptable values, with 1.8 and 2.88, respectively, whereas all other parameters and the deviance (i.e., the log posterior) had \( \widehat{R} \) values of approximately 1.

The Monte Carlo chains of M 1 showed high autocorrelation and poor mixing, prompting us to subsample the chains by a factor of 1,000 in addition to using a burn-in period of 1,000, resulting in 150 independent samples. The overall shape of the thinned chains was satisfactory, on the basis of visual examination. Figure 3 shows the deviance chain of the easy model before and after burning and thinning, with the chain after burning and thinning showing good mixing.

Fig. 3
figure 3

Examples of a good and a bad chain. The upper left panel shows the deviance chain, as it was sampled. The lower left panel shows the deviance chain after thinning and burning. The two panels on the right show the autocorrelation plots before and after thinning, respectively

Similar results held for the more complex model M 2: High autocorrelation and unconverged chains at the beginning were removed by burning the chains with 20,000 and thinning them by 2,000, resulting in 348 independent samples. Figure 4 shows a typical chain of M 2. Two parameters showed poor mixing even after burning and thinning: μ η and σ η . However, the total deviance of the model did not show convergence issues, suggesting that these two parameters converged poorly due to the modest impact that they had on the model’s fit to the data (i.e., they were poorly identified by the data).

Fig. 4
figure 4

Trace and autocorrelation plots for the first drift rate mean of the M 2 model

Finally, Table 2 shows the summary statistics of the drift rate population distributions for the M 1 and M 2 models. As the table shows, the parameters differed between conditions. The outcome of our analysis with JAGS did not differ substantially from the results reported by Vandekerckhove et al. (2011).

Table 2 Estimates (posterior means) of the five population means (μ) and standard deviations (σ) of the drift rate parameters for the two models

Summary

JAGS is an open-source software package for the analysis of graphical models that is written with extensibility in mind. Additionally, open-source software is advantageous for academic research because of its permanency and accessibility. We provide step-by-step instructions on how to implement custom distributions and logical nodes in JAGS. It is hoped that our documenting the process of extending JAGS will contribute to the formation of a collaborative community that will extend the usefulness of JAGS even more.

We have implemented the first-passage time distribution of a Wiener diffusion model as a worked example, and provide the resulting module as a freely downloadable package.