The Rwiener Package: an R Package Providing Distribution Functions for the Wiener Diffusion Model

We present the RWiener package that provides R functions for the Wiener diffusion model. The core of the package are the four distribution functions dwiener, pwiener, qwiener and rwiener, which use up-to-date methods, implemented in C, and provide fast and accurate computation of the density, distribution, and quantile function, as well as a random number generator for the Wiener diffusion model. We used the typical Wiener diffusion model with four parameters: boundary separation, non-decision time, initial bias and drift rate parameter. Beyond the distribution functions, we provide extended likelihood-based functions that can be used for parameter estimation and model selection. The package can be obtained via CRAN.


Introduction
The diffusion model is one of the most successful models of choice reaction time in cognitive psychology (see Wagenmakers, 2009, for an overview) and, while software packages that make diffusion model analyses possible have been made available (Vandekerckhove andTuerlinckx, 2007, 2008;Vandekerckhove et al., 2011;Voss and Voss, 2007;Wagenmakers et al., 2007;Wiecki et al., 2013), a diffusion model R package was so far missing.Due to the nature of the model, no closed analytical forms for the computation of the distribution functions exist, which has generated a need for more involved computational approximations of the various functions (Blurton et al., 2012;Navarro and Fuss, 2009;Tuerlinckx et al., 2001).
Here, we present an R package, RWiener, that provides the four functions R uses to represent a distribution, implemented for the Wiener diffusion model: the density function d (Navarro and Fuss, 2009) to compute the probability density function (PDF) at a given quantile for a given parameter set; the probability function p (Blurton et al., 2012) to compute the cumulative distribution (CDF) at a given quantile for a given parameter set; the inverse CDF or quantile function1 q to compute the quantile at a given p-value (CDF-value) for a given parameter set; and a random number generator (RNG) r (Tuerlinckx et al., 2001) to generate random samples for a given parameter set.
In addition, we extended the package to include several likelihood-based functions that can be used in combination with R's built-in estimation routines, like the nlm or optim function, to estimate the model parameters for an observed data set, assuming a Wiener diffusion process as the underlying process that created the data.
First, we will present the standard Wiener diffusion model with four parameters: the boundary separation α, the non-decision time τ, the bias parameter β and a drift rate parameter δ.Next, we will show the standard functionality of the RWiener package and how it can be used for parameter estimation.To conclude, we add a small discussion about the utility and extensibility of the here presented methods.

The Wiener diffusion model
For two-choice reaction time (2CRT) data coming from experiments in which subjects are asked to give one of two responses on a decision task, the Wiener diffusion model and its extensions have provided useful modeling approaches.In its simplest complete form, the Wiener diffusion model includes four parameters.The decision process is thought of as a continuous random walk, or diffusion, process that starts somewhere between two boundaries and ends as soon as it hits or crosses one or the other.The time of the first passage, and which boundary is hit first (and therefore which decision will be made) is probabilistic with a probability determined by the parameter set.
The standard Wiener diffusion model incorporates the following four parameters, a summary of which is provided in Table 1: (1) the boundary separation α which defines the distance between the two boundaries, with the first boundary being at 0 and the second being at α, (2) the non-decision time τ that defines the time that passes without any diffusion process going on (e.g., this incorporates the time needed to encode a stimulus and to execute a motor response), (3) the bias parameter β defining   The process terminates as soon as the accrued evidence exceeds α or deceeds 0. The decision process starts at time τ from the stimulus presentation and terminates at the reaction time.
the relative starting point of the diffusion process between the two boundaries (the absolute starting point ζ 0 can be obtained by ζ 0 = βα) and ( 4) the drift rate parameter δ which captures the tendency of the diffusion process to drift towards the upper boundary.Given these parameters, the latent information accumulation process is modeled with the stochastic process d dt X t ∼ Normal δ, s 2 , with initial condition X 0 = αβ.The decision time is modeled as the first time at which X t ≤ 0 or X t ≥ α, and the response time is the sum of the decision time and the non-decision time τ. Figure 1 shows an illustration of the Wiener diffusion model as described.
The parameters have the following restrictions: 0 < β < 1, α > 0, τ > 0. Often, the variance parameter s is fixed at 0.1, whereas we fixed it at 1, yielding a more convenient interpretation of the drift rate parameter and greater computational efficiency.Conversion to a scale with the variance parameter fixed at any s can be done by multiplying the α and δ parameters by s.
The upper boundary may be defined as corresponding to correct trials, with the lower boundary corresponding to error trials.However, it often also makes sense to map qualitatively different responses to the two boundaries (e.g., in a same-different discrimination task, the upper bound could represent the 'different' responses).For generality, we will always speak simply of upper and lower boundaries.An overview of the computational strategy for the PDF and CDF functions is provided in the Appendix.

Usage of the package
The package can be obtained via CRAN for R version 2.15.0 or higher, and includes documentation that provides useful information and examples.As a guide through this demonstration of the RWiener package functionality, we use an example data set created by the random function of the package (using a RNG seed of 0 for exact reproducibility): The R Journal Vol.6/1, June ISSN 2073-4859 set.seed(0)dat <-rwiener(n=100, alpha=2, tau=.3,beta=.5, delta=.5) The arguments for this function call are explained in the next paragraph.This command will generate a random data set, consisting of 100 observations, each generated from a random Wiener diffusion model process with the given four parameters.We used a drift rate parameter δ of 0.5, meaning that the drift diffusion process slightly tends towards the upper boundary.There was no initial bias towards either of the two boundaries, indicated by β = 0.5.The boundaries are separated by a value of α = 2 and we used a non-decision time of τ = 0.3, so there will never be a RT smaller than 0.3.
The created dat object is a dataframe with 100 observations and two variables: dat$q contains the RTs, whereas dat$resp contains a factor to indicate the response made (upper vs. lower).Note that the upper bound is assigned to the first level and the lower bound to the second one.
The four main functions of the package are described in the following subsection.All four functions are implemented in C, using the R library for C functions.These C functions are called inside end-user R functions, so that the end-user never calls the C functions directly, only indirectly through the provided R functions.

Standard distribution functions
To get the density of a specific quantile, one can use the density function dwiener: dwiener(dat$q[1], alpha=2, tau=.3,beta=.5, delta=.5, resp=dat$resp[1], give_log=FALSE) The function takes seven arguments, the first argument is the RT at which the density shall be evaluated, the next four arguments are the model parameters α, τ, β, and δ.The next argument, resp takes any of three valid values: upper, lower, or both, meaning that the function shall return the density at the given quantile for the upper boundary, the lower boundary or the sum of both densities together, respectively.Its default value is set to upper.Lastly, the function can be given an additional argument give_log that takes TRUE or FALSE as values and is set to FALSE as a default.When this argument is TRUE, the function will return the logarithm of the density instead of the density itself.
This function can also be used to draw simple plots.For example, the following code results in the plot in the left panel of Figure 2: curve(dwiener(x, 2, .3,.5, .5, rep("upper", length(x))), xlim=c(0,3), main="Density of upper responses", ylab="density", xlab="quantile") To calculate the cumulative distribution, we can use the CDF function pwiener as follows: alpha=2, tau=.3, beta=.5, delta=.5, resp=dat$resp[1]) The arguments are the same as for the previously described function, without the give_log argument.Further, this function can be used in the same manner to draw plots, like the one shown in the right panel of Figure 2.
We can also find the appropriate quantile for a given probability via the inverse CDF function, implemented as qwiener: # lookup of the .2quantile for the CDF of the lower boundary qwiener (p=.2, alpha=2, tau=.3, beta=.5, delta=.5, resp="lower") The R Journal Vol.6/1, June ISSN 2073-4859 Again, the last five arguments of this function have the same meaning as for the two previously described functions.The first argument, however, is different: instead of a quantile, the first argument needs to be a probability, indicating the CDF-value to be looked up.
To round out the four standard distribution functions, the RWiener package further provides a function named rwiener that generates random values.Again, the last four input variables correspond to the model parameters and are the same as in the other distribution functions.The first argument n takes a number, indicating the desired number of randomly generated variates.

Plot function
To easily visualize data from the Wiener diffusion model, we can use the plot function wiener_plot as follows:

wiener_plot(dat)
This draws the two densities of the observed (or, in this case, generated) lower and upper responses.Figure 3 shows an example of this plot.

Parameter estimation with R's general-purpose optimization algorithms
In addition to the four distribution functions dwiener, pwiener, qwiener, rwiener and the plot function wiener_plot, the RWiener package provides four more functions designed for use in parameter estimation and model selection: (1) the wiener_likelihood function, which calculates the logarithm of the likelihood for given parameter values and data; (2) the wiener_deviance function, which calculates the model deviance −2 ; (3) the wiener_aic function, which calculates the model AIC (Akaike's Information Criterion) −2 + 8; and (4) the wiener_bic function, which calculates the model BIC (Bayesian Information Criterion) −2 + 4 ln(N).
x <-c(2, .3,.5, .5)wiener_likelihood(x=x, dat=dat) wiener_deviance(x=x, dat=dat) wiener_aic(x=x, dat=dat) wiener_bic (x=x, dat=dat) The R Journal Vol.6/1, June ISSN 2073-4859 The first argument of these functions, x, is a vector containing the four parameter values of the model in the order: α, τ, β, δ.The second argument of these functions, dat, is a dataframe that contains the data points and the values of the two dependent variables: the RT and the boundary that was hit.This dataframe has to have the same shape as the dataframe dat as generated previously with the rwiener function.The first column of the dataframe has to contain the RTs and the second column of the dataframe has to contain the response type (upper vs. lower).

Estimating parameters for multiple conditions
The deviance function can be combined to compute a single loss value for multiple conditions.To do so, first create a new inline function: # create a second data set and a list containing both data sets dat2 <-rwiener(n=100, alpha=2, tau=.3,beta=.5, delta=1) datlist <-list(dat, dat2) # use nlm to estimate parameters nlm1 <-nlm(p=c(1, .1,.1, 1, 1), f=many_drifts, dat=datlist) If the input x is a parameter vector containing the first three parameters (α, τ, β) and then C drift rates, and dat is a list with C data sets, then the result will contain point estimates of three joint parameters in addition to C condition-specific drift rate estimates.
An alternative function can be defined that does not allow the drift rates to differ between conditions: Finally, the two competing models can be compared using the AIC criterion for model selection.3AIC1 <-wiener_aic(x=nlm1$estimate, dat=datlist, loss=many_drifts) AIC2 <-wiener_aic(x=nlm2$estimate, dat=datlist, loss=one_drift) For more details on the interpretation of the AIC criterion, see Wagenmakers and Farrell (2004).

Discussion
We presented RWiener, an R package that provides efficient algorithms to compute Wiener diffusion model functions and use these for further inference.This is the first package in R that provides full Wiener diffusion model functionality for R. We demonstrated how to use these functions and how to extend their usage to combine them with other R functions in order to do parameter estimation.
The presented functions can help researchers who work with 2CRT data to analyze their data in the Wiener diffusion model framework.

Figure 1 :
Figure1: A graphical illustration of the Wiener diffusion model for two-choice reaction times.An evidence counter starts at value αβ and evolves with random increments.The mean increment is δ.The process terminates as soon as the accrued evidence exceeds α or deceeds 0. The decision process starts at time τ from the stimulus presentation and terminates at the reaction time.

Figure 2 :
Figure 2: Example of plots made with function dwiener and pwiener.

Table 1 :
Vandekerckhove (2009)ers of the Wiener diffusion model, with their substantive interpretations.Reprinted with permission fromVandekerckhove (2009).Copyright 2009, Joachim Vandekerckhove and Department of Psychology and Educational Sciences, University of Leuven, Belgium.