In jeff324/powder: Power posterior sampling via differential evolution

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This vignette describes the hierarchical examples in the paper "Thermodynamic integration and Steppingstone sampling methods for computing Bayes factors: A tutorial" by Annis, Evans, Miller, and Palmeri (submitted). The paper comes bundled with the powder package.

vignette("powder-paper",package = 'powder')

Additionally, the examples from Evans and Annis (submitted), "Thermodynamic integration via differential evolution: A method for approximating marginal likelihoods", can be run by supplying the method='parallel' argument to powder. The two methods from the paper, TIDE and TIDEz, can then be run by supplying de_params = list(randomize_phi=FALSE) and de_params = list(randomize_phi=TRUE), respectively to powder.

vignette("powder-paper-tide",package = 'powder')

In the paper, we used two different simulated datasets, each with 10 simulated subjects in two conditions with 300 trials per condition. The "null" dataset contained no differences in parameters over conditions. The "drift" dataset contained changes in drift rates across conditions. In addition to the two simulated datasets we also used an empirical dataset from Rae et al. (2014).

We fit 4 different models to each dataset. We used two models that are commonly of theoretical interest: one that contained different drift rates for each condition, and one that contained different thresholds for each condition. In addition, we included a simple variant – a model with no parameters varying over condition – and a complex variant – with drift rate, threshold, and non-decision time all varying over condition – in order to see how different levels of dimensionality would affect the estimates.

To load the "null" dataset run the following. If you want to run the drift dataset, set dataset = 'drift'. For the Rae et al. dataset, set dataset = 'rae'. For the individual (single-subject) dataset, set dataset = 'individual.

library(powder)

dataset = 'null'

if(dataset=='null'){
     data('null')
     dat = null
}
if(dataset=='drift'){
     data('drift')
     dat = drift
}
if(dataset=='rae'){
     data('rae')
     dat = rae
}

To define the simple hierarchical LBA model in the paper do the following:

simple.model = LBA$new(contaminant=list(pct=2,upper.bound=5))

Next, define the priors in the paper. Note that LBA priors are all truncated normals with mean mu, standard deviation sigma, lower bound 0, and upper bound Inf.

simple.model$prior$A$mu = c(mu=1,sigma=1) 
simple.model$prior$A$sigma = c(mu=1,sigma=1)

#note that b is the relative threshold (denoted as k in the paper)
simple.model$prior$b$mu = c(mu=.4,sigma=.4) 
simple.model$prior$b$sigma = c(mu=.4,sigma=.4)

simple.model$prior$t0$mu = c(mu=.3,sigma=.3)
simple.model$prior$t0$sigma = c(mu=.3,sigma=.3)

simple.model$prior$vc$mu = c(mu=3,sigma=3)
simple.model$prior$vc$sigma = c(mu=3,sigma=3)

simple.model$prior$ve$mu = c(mu=1,sigma=1)
simple.model$prior$ve$sigma = c(mu=1,sigma=1)

simple.model$prior$sve$mu = c(mu=1,sigma=1)
simple.model$prior$sve$sigma = c(mu=1,sigma=1)

We can create the other models quite easily because they share the same priors. Note the priors are the same across conditions. To create the complex model:

complex.model = LBA$new(vc = TRUE, b = TRUE, t0 = TRUE, prior = simple.model$prior,contaminant=list(pct=2,upper.bound=5))

Drift model

drift.model = LBA$new(vc = TRUE, prior = simple.model$prior,contaminant=list(pct=2,upper.bound=5))

Threshold model

threshold.model = LBA$new(b = TRUE, prior = simple.model$prior,contaminant=list(pct=2,upper.bound=5))

We can then run each model on the chosen data and compute marginal likelihood. This will take a while at these settings.

n.samples = 300
burnin = 1000
meltin = 200
num.temps = 35

simple.out = powder(data=dat, model=simple.model, num.temps = num.temps, 
                    n.samples = n.samples, burnin = burnin, meltin = meltin)
simple.ml = summary(simple.out)

complex.out = powder(data=dat, model=complex.model, num.temps = num.temps, 
                    n.samples = n.samples, burnin = burnin, meltin = meltin)
complex.ml = summary(complex.out)

drift.out = powder(data=dat, model=drift.model, num.temps = num.temps, 
                    n.samples = n.samples, burnin = burnin, meltin = meltin)
drift.ml = summary(drift.out)

threshold.out = powder(data=dat, model=threshold.model, num.temps = num.temps, 
                    n.samples = n.samples, burnin = burnin, meltin = meltin)
threshold.ml = summary(threshold.out)