In this tutorial, I illustrated fitting multiple participants, assuming the mechanism of data generation is fixed-effect models. That is, each participant is accounted for by independent mechanisms. I also assume the LBA model is the true RT model.

I made up a two-factor factorial design. The first two-level factor is the stimulus (S). Suppose the stimuli have two types: one is low quality face photos, so people find it hard to recognize and the other is normal quality face photo. The second factor is the frequency (F), supposing one type is the celebrity photos, so people perhaps see more often, and the other is the photos of randomly selected strangers.

In the model set-up, I presume a rate model, which has its drift rates
affected by the two factors, S and F. The latent factor, *M*, is just a
LBA way to model independent accumulators. Another factor, *R*, not
explicitly in the factorial design, is an indicator factor, indicating
the response type affecting the threshold parameter (i.e., accumulators
traveling distance).

```
require(ggdmc)
model <- BuildModel(
p.map = list(A = "1", B = "R", t0 = "1",
mean_v = c("S", "F", "M"), sd_v = "M", st0 = "1"),
match.map = list(M = list(s1 = 1, s2 = 2)),
factors = list(S = c("s1", "s2"), F = c("f1", "f2")),
constants = c(sd_v.false = 1, st0 = 0),
responses = c("r1", "r2"),
type = "norm")
## [1] "A" "B.r1" "B.r2"
## [4] "t0" "mean_v.s1.f1.true" "mean_v.s2.f1.true"
## [7] "mean_v.s1.f2.true" "mean_v.s2.f2.true" "mean_v.s1.f1.false"
## [10] "mean_v.s2.f1.false" "mean_v.s1.f2.false" "mean_v.s2.f2.false"
## [13] "sd_v.true"
npar <- length(GetPNames(model))
```

To simulate many participants, I set up a population distribution, which is not in line with the assumption of fixed-effects model. That is, this way to generate data is to presume that a random-effects model at work. For the purpose of illustration, I forgo this issue for now.

```
pop.mean <- c(A = .4, B.r1 = .85, B.r2 = .8, t0 = .1,
mean_v.s1.f1.true = 2.5, mean_v.s2.f1.true = 3.5,
mean_v.s1.f2.true = 4.5, mean_v.s2.f2.true = 5.5,
mean_v.s1.f1.false = 1.00, mean_v.s2.f1.false = 1.10,
mean_v.s1.f2.false = 1.05, mean_v.s2.f2.false = 1.20,
sd_v.true = .25)
pop.scale <- c(A = .1, B.r1 = .1, B.r2 = .1, t0 = .05,
mean_v.s1.f1.true = .2, mean_v.s2.f1.true = .2,
mean_v.s1.f2.true = .2, mean_v.s2.f2.true = .2,
mean_v.s1.f1.false = .2, mean_v.s2.f1.false = .2,
mean_v.s1.f2.false = .2, mean_v.s2.f2.false = .2,
sd_v.true = .1)
pop.prior <- BuildPrior(
dists = rep("tnorm", npar),
p1 = pop.mean,
p2 = pop.scale,
lower = c(rep(0, 4), rep(NA, 8), 0),
upper = c(rep(NA, npar)))
```

We may want to check how the prior distributions look like.
*ggdmc* has a *plot* function to do just that. Note you need to load *ggdmc*
package (i.e., *require(ggdmc)*) to make *plot* function changes its default
behaviour.

```
plot(pop.prior)
```

```
## Simulate some data
dat <- simulate(model, nsim = 30, nsub = 8, prior = pop.prior)
dmi <- BuildDMI(dat, model)
dplyr::tbl_df(dat)
## # A tibble: 960 x 5
## s S F R RT
## <fct> <fct> <fct> <fct> <dbl>
## 1 1 s1 f1 r1 0.438
## 2 1 s1 f1 r1 0.517
## 3 1 s1 f1 r1 0.407
## 4 1 s1 f1 r1 0.454
## 5 1 s1 f1 r1 0.449
## 6 1 s1 f1 r1 0.463
## 7 1 s1 f1 r1 0.552
## 8 1 s1 f1 r1 0.411
## 9 1 s1 f1 r1 0.387
## 10 1 s1 f1 r1 0.486
## # ... with 950 more rows
```

The true averaged parameter vectors, which were randomly chosen based
on *pop.prior*, can be retrieved by looking up the **parameters** attribute,
attached onto the *dat* object. However, note the real true values are *pop.mean* and
*pop.scale*, because the data are generated based on random-effects model.

```
require(matrixStats)
ps <- attr(dat, "parameters")
mu <- round(colMeans2(ps), 2)
sigma <- round(colSds(ps), 2)
truevalues <- rbind(mu, sigma)
colnames(truevalues) <- GetPNames(model)
## Set up prior distributions
p.prior <- BuildPrior(
dists = rep("tnorm", npar),
p1 = pop.mean,
p2 = pop.scale*20,
lower = c(rep(0, 4), rep(NA, 8), 0),
upper = c(rep(NA, npar)))
```

Now we am ready to fit the eight participants. Ideally, this can be
done simultaneously, if a eight-core machine is available.
Here I used a four-core machine so launched 2 cores only (*ncore = 2*).

```
## Sampling -------------
fit0 <- StartNewsamples(dmi, p.prior, ncore = 2)
fit <- run(fit0, 5e2, ncore = 2)
```

Use plot to check whether posterior log-likelihood converged.

plot(fit)

*gelman* function prints the potential scale reduction
factor (psrf). A psrf value less than 1.1 suggests chains are well-mixed.

```
res <- gelman(fit, verbose = TRUE)
# Diagnosing theta for many participants separately
# 15 20 13 18 6 17 12 19 5 11 8 3 16 14 1 2 9
# 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.02 1.02
# 4 7 10
# 1.02 1.03 1.03
# Mean
# [1] 1.01
```

By setting the option, *pll* (posterior log-likelihood),
to FALSE and the option, *den* (density plot), to TRUE , we
can check the trace plots for each model parameters. Because there are
several participants, the size of the figure is considerably large.
It may be a better to plot separately in a pdf file and check it later.

```
pdf("figs/subjects-density.pdf")
lapply(fit, ggdmc:::plot.model, pll = FALSE, den = TRUE )
dev.off()
```

One specific feature in *ggdmc* is that it uses pMCMC, so occasionally,
we want to check a subset of chains. The function will randomly pick three chains
to plot

plot(sam, subchain = TRUE)

You can indicate how many subset of chains to plot, too.

plot(sam, subchain = TRUE, nsubchain = 4))

You can also indicate which chains, instead of randomly selecting a subset of chains.

plot(sam, subchain = TRUE, nsubchain = 4, chains = c(1:4))

These are a lot of checks! Finally and fortunately, because this is a parameter recovery study,
I can look up the true parameters to see if I do estimate them well. *summary* function will
do the trick by entering TRUE for the *recovery* option and entering the true parameter matrix,
which I had stored it to a *ps* object before, to *ps* option.

est <- summary(sam, recovery = TRUE, ps = ps)

```
# Summary each participant separately
# A B.r1 B.r2 t0 mean_v.s1.f1.true mean_v.s2.f1.true mean_v.s1.f2.true
# Mean 0.50 1.02 0.96 0.10 3.05 4.34 5.50
# True 0.41 0.82 0.77 0.10 2.49 3.56 4.50
# Diff -0.09 -0.20 -0.19 0.00 -0.56 -0.78 -1.00
# Sd 0.15 0.20 0.17 0.05 0.49 0.44 0.76
# True 0.11 0.10 0.09 0.05 0.21 0.19 0.22
# Diff -0.04 -0.10 -0.08 -0.01 -0.28 -0.25 -0.54
# mean_v.s2.f2.true mean_v.s1.f1.false mean_v.s2.f1.false mean_v.s1.f2.false
# Mean 6.69 1.59 1.58 0.70
# True 5.49 1.03 1.15 1.08
# Diff -1.20 -0.56 -0.43 0.38
# Sd 0.82 0.34 1.21 1.58
# True 0.22 0.15 0.23 0.15
# Diff -0.59 -0.19 -0.98 -1.43
# mean_v.s2.f2.false sd_v.true
# Mean -0.38 0.30
# True 1.11 0.24
# Diff 1.48 -0.06
# Sd 1.23 0.11
# True 0.16 0.10
# Diff -1.07 -0.01
```