minimal_example

This vignette fits the longitudinal drift-diffusion mixed model to the bundled data, checks the fit, and points to a simulated example in which the true parameters are known. The chains below are kept very short so that the vignette builds quickly; for an analysis use something like Niter = 5000, burnin = 2000, thin = 5.

Data

One row per trial: participant, training block, input stimulus s, decision d, response time in seconds and a censoring indicator. For numerical stability, most response times should lie in \([0, 10]\).

data('data')
head(data)
#> # A tibble: 6 × 6
#>   subject block     s     d r_time  cens
#>     <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>
#> 1       1     1     2     1  1.63      0
#> 2       1     1     2     2  0.802     0
#> 3       1     1     1     1  1.10      0
#> 4       1     1     2     1  2.36      0
#> 5       1     1     1     2  1.73      0
#> 6       1     1     2     1  1.00      0
plot_accuracy(data)

plot_RT(data)

Fit

set.seed(123)
fit = LDDMM(data = data,
            hypers = list(s_sigma_mu = 0.1, s_sigma_b = 0.1),
            Niter = 200, burnin = 100, thin = 1,
            save_z = TRUE, verbose = FALSE)

Population curves

Posterior means and 90% credible bands of the drifts \(\mu_{d,s}(t)\) and boundaries \(b_{d,s}(t)\); panels are input stimuli, colours are decisions.

plot_post_pars(data, fit, par = 'drift')

plot_post_pars(data, fit, par = 'boundary')

Posterior predictive checks

For every saved draw, the fit simulates one response of every participant to every stimulus in every block. Compare them with the data: the proportion of each decision, and the 10%, 50% and 90% quantiles of the response times.

plot_predictive(data, fit, type = 'response')

plot_predictive(data, fit, type = 'rt')

Convergence

Traces should fluctuate around a stable level (with a chain this short they will not yet), and the acceptance rates should be roughly between 0.2 and 0.6.

traces = data.frame(draw = seq_along(fit$sigma2_1_mu),
                    sigma2_1_mu = fit$sigma2_1_mu,
                    sigma2_mu_ua = fit$sigma2_mu_ua,
                    drift_block1 = fit$post_mean_mu[1, 1, 1, ],
                    boundary_block1 = fit$post_mean_b[1, 1, 1, ])
ggplot(tidyr::pivot_longer(traces, -draw), aes(draw, value)) +
  geom_line() +
  facet_wrap(~ name, scales = 'free_y')

round(sapply(fit$acceptance, mean), 2)
#>        delta     beta1_mu         U_mu          U_b  sigma2_1_mu sigma2_mu_ua 
#>         0.74         0.15         0.62         0.47         0.78         0.53 
#> sigma2_mu_us  sigma2_b_ua     scale_mu      scale_b    ridge_pop   ridge_subj 
#>         0.53         0.45         0.37         0.18         0.00         0.18 
#>      beta1_b   sigma2_1_b  sigma2_b_us      alpha_C      alpha_I 
#>         0.15         0.88         0.17         0.80         0.80

Co-clustering probabilities of curves that are only beginning to separate carry more Monte Carlo error than the curves themselves; when they matter, run two or more chains with different seeds and compare them. The same holds for the level of the boundaries and the offsets, which can trade off against each other.

Results

head(extract_post_mean(data, fit, par = 'drift'))
#> # A tibble: 6 × 4
#>   block subject drift_1_1 drift_2_2
#>   <int>   <dbl>     <dbl>     <dbl>
#> 1     1       1     0.568     0.568
#> 2     2       1     0.569     0.576
#> 3     3       1     0.593     0.599
#> 4     4       1     0.673     0.692
#> 5     5       1     0.765     0.789
#> 6     1       2     1.12      1.12
compute_WAIC(fit)
#> [1] 31733.49
# posterior probability that the correct-response curves of the two stimuli
# share a cluster, for each spline coefficient
round(apply(fit$Z[1, 1, , ] == fit$Z[2, 2, , ], 1, mean), 2)
#> [1] 1.00 1.00 0.00 0.99 0.00 0.95

The boundary model is chosen with boundaries ("flexible", "constant", "fixed" or "fixed-constant") and the random effects with random_effects ("correct-incorrect", as in the article, or "shared"); models can be compared with compute_WAIC().

A simulated example

The script toy_simulation.R, installed with the package, simulates data from the model’s likelihood with known parameters and compares the posterior with the truth. Three tones are trained: tone 2 is learned, tone 3 is learned faster, and tone 1 is “anti-learned”, so that it is answered with category 2 more and more often. The script recovers the population and subject curves, the offsets, the variance parameters and the local clustering, and runs the same posterior predictive checks as above.

file.edit(system.file("examples", "toy_simulation.R", package = "lddmm"))

toy_example.R repeats the analysis of this vignette with a longer chain and adds subject-level curves, co-clustering probabilities by block and a model comparison.