The data

First we will load the pmwg package and explore the included dataset

library(pmwg)

head(forstmann)

##      subject condition stim resp     rt
## 1115       1         1    2    2 0.4319
## 1116       1         3    2    2 0.5015
## 1117       1         3    1    1 0.3104
## 1118       1         1    2    2 0.4809
## 1119       1         1    1    1 0.3415
## 1120       1         2    1    1 0.3465

The forstmann dataset is from (Forstmann et al. 2008) and is the cleaned data from an experiment that displayed a random dot motion task with one of three speed accuracy manipulations (respond accurately instructions, neutral instructions and respond quickly instructions).

We can see that there are 5 columns:

• subject - which gives an ID for each of the 19 participants
• condition - the speed accuracy instructions for the trial
• stim - whether dots were moving left or right
• resp - whether the participant responded left or right
• rt - the response time for the trial

We will use the rtdists package to look at threshold differences between the three conditions to see whether there is evidence for different thresholds (the evidence required to trigger a response) in each of the three conditions. For a full explanation of the model please see the full documentation in Chapter 3.

The log-likelihood function

Now that we have the data we will need to specify the log-likelihood function.

lba_loglike <- function(x, data) {
  x <- exp(x)
  if (any(data$rt < x["t0"])) {
    return(-1e10)
  }
  # This is faster than "paste".
  bs <- x["A"] + x[c("b1", "b2", "b3")][data$condition]

  out <- rtdists::dLBA(
    rt = data$rt, # nolint
    response = data$stim,
    A = x["A"],
    b = bs,
    t0 = x["t0"],
    mean_v = x[c("v1", "v2")],
    sd_v = c(1, 1),
    distribution = "norm",
    silent = TRUE
  )
  bad <- (out < 1e-10) | (!is.finite(out))
  out[bad] <- 1e-10
  out <- sum(log(out))
  out
}

This function contains the primary implementation of our model for the data. The pmwgs object that we run later in this example will pass two things:

• A set of parameter values x that we will be assessing.
• A subset of our dataset (data), the data for one subject.

Our job then in the function is to take the parameter values and calculate the log-likelihood of the data given the parameter values.

The overall process is as follows:

• Take the exponent of the parameter values (that are stored in log-form
• Test for improbable values of t0 (non-decision time) and return extremely low log-likelihood if they are unlikely.
• Create a vector of b (threshold) parameter values for the call to rtdists
• Use the rtdists dLBA function to calculate the likelihood of the data given our parameter estimates.
• Clean extremely small of otherwise bad values from output vector.
• Return the sum of the log of the likelihood of each row of the data.

Other necessary components

There are a couple of other necessary components to running the sampler, a list of parameters and a prior for the group level parameter values.

In this case the parameters consist of threshold parameters for each of the three conditions (b1, b2 and b3), a upper limit on the range for accumulator start points (A), the drift rates for evidence accumulation for the two stimulus types (v1 and v2) and non-decision time (t0).

The prior for this model is uninformed and is just 0 for the mean of the group parameters and standard deviation of 1 for each parameter.

pars <- c("b1", "b2", "b3", "A", "v1", "v2", "t0")

priors <- list(
  theta_mu_mean = rep(0, length(pars)),
  theta_mu_var = diag(rep(1, length(pars)))
)

Creating and running the sampler

Once we have these components we are ready to create and run the sampler.

# Create the Particle Metropolis within Gibbs sampler object
sampler <- pmwgs(
  data = forstmann,
  pars = pars,
  ll_func = lba_loglike,
  prior = priors
)

Creation of the sampler object is done through a call to the pmwgs function, which creates the object and sets numerous precomputed values based on the number of parameters and other aspects of the included function arguments. Next steps are to initialise the sampler with start values for the random effects and then run the three stages of sampling.

The stages can be long running, but once complete you will have

# Initialise (generate first random effects for sampler)
sampler <- init(sampler)  # Can also pass start points for sampler

# Run each stage of the sampler, can adjust number of particles on each
sampler <- run_stage(sampler, stage = "burn")
sampler <- run_stage(sampler, stage = "adapt")
sampler <- run_stage(sampler, stage = "sample")

The run_stage command can also be passed other arguments such as iter for number of iterations, particles for number of particles among others and more. For a full list see the description in the PMwG Tutorial Book.