The EZ-Diffusion Model (EZDM)

1 Introduction to the model

The diffusion model (Ratcliff 1978) is a widely used model in speeded decision-making. It describes how participants accumulate noisy evidence over time until reaching one of two decision boundaries. Usually, the diffusion model is fit to trial-level reaction time distributions, which can be computationally intensive, especially for hierarchical Bayesian estimation.

The EZ-Diffusion Model (Wagenmakers et al. 2007) provides a simplified alternative building on closed-form equations that relate the parameters of the diffusion model to easily computed summary statistics: mean reaction time (MRT), reaction time variance (VRT), and accuracy (proportion correct, Pc). This makes it possible to calculate diffusion model parameters from aggregated data without fitting to individual trials.

1.1 The Diffusion Model

In the standard diffusion model, evidence accumulates from a starting point toward one of two boundaries. The key parameters are:

Drift rate (\(v\)): The average rate of evidence accumulation, reflecting task difficulty or ability
Boundary separation (\(a\)): The distance between the two decision boundaries, reflecting the speed-accuracy tradeoff
Non-decision time (\(T_{er}\)): Time for processes outside the decision (encoding, motor response)
Starting point (\(z\)): A response bias indicating where evidence accumulation begins relative to the boundaries

Illustration of the diffusion model. Evidence accumulates from the starting point until reaching the upper (correct) or lower (error) boundary.

Figure 1.1: Illustration of the diffusion model. Evidence accumulates from the starting point until reaching the upper (correct) or lower (error) boundary.

1.2 The EZ Equations

Wagenmakers et al. (2007) provided closed-form solutions relating the diffusion model parameters to observable summary statistics. For the 3-parameter version (assuming symmetric starting point, \(z = a/2\)):

Given the proportion correct (\(P_c\)), mean RT of correct responses (\(MRT\)), and variance of correct RTs (\(VRT\)), the parameters can be computed as:

\[ v = \text{sign}(P_c - 0.5) \cdot s \cdot \left[ L(P_c) \cdot \frac{P_c^2 L(P_c) - P_c \cdot L(P_c) + P_c - 0.5}{VRT} \right]^{1/4} \]

\[ a = s^2 \cdot \frac{L(P_c)}{v} \]

\[ T_{er} = MRT - \frac{a}{2v} \cdot \frac{1 - e^{-\frac{va}{s^2}}}{1 + e^{-\frac{va}{s^2}}} \]

where \(L(P_c) = \ln\left(\frac{P_c}{1 - P_c}\right)\) is the logit function and \(s\) is the diffusion constant (typically fixed to 1 for identification).

1.3 The 4-Parameter Extension

The original EZ model assumes a symmetric starting point. Srivastava et al. (2016) derived the moments of first-passage times for the Wiener diffusion process with arbitrary starting points, enabling a 4-parameter version that can estimate starting point bias (\(z_r\), the relative starting point). This 4-parameter version, however, requires separate summary statistics for responses to each boundary:

Mean and variance of RTs for upper boundary responses (typically correct)
Mean and variance of RTs for lower boundary responses (typically errors)

With sufficient responses at both the upper and lower bound, this allows estimation of response bias, which is important when participants favor one response over another.

1.4 Bayesian Estimation

So far, studies using the EZ-Diffusion model have adopted step-wise analyses approaches. First, diffusion model parameters are calculated for summary statistics for each individual in each design cell. Second, the calculated parameters are submitted to statistical tests evaluating differences between experimental conditions. This approach has two downsides: 1) uncertainty in parameter estimation is not sufficiently propagated between analysis steps, and 2) all parameters of the diffusion model were allowed to vary over all design cells, bearing the risk of over fitting. To address these issues, Chávez De la Peña and Vandekerckhove (2025) introduced a Bayesian hierarchical approach to the EZ-diffusion model. This approach combines the derivation of summary statistics from diffusion model parameters with a probabilistic model of the sampling distributions of the summary statistics. Specifically, this approach:

Treats the summary statistics as observed data generated from a known sampling distribution
Uses likelihood-based inference to estimate parameters with proper uncertainty quantification
Enables hierarchical modeling with random effects across participants

Chávez De la Peña and Vandekerckhove (2025) shared JAGS code alongside their model. The downside of this implementation is that the model code has to be adapted to different experimental designs which likely leads to errors and increases the bar for using the implementation. The bmm package implements this Bayesian approach, allowing researchers to:

Estimate uncertainty in parameter estimates
Include random effects for participants
Compare conditions using Bayesian hypothesis testing
Leverage the full power of the brms modeling framework

1.5 The sampling distribution of the summary statistics

Because the summary statistics are the data, the likelihood is a statement about how they vary from one sample of trials to the next. The number of upper-boundary responses is binomial. The mean and the variance of the response times are the sample mean and sample variance of \(n\) decision times, and their joint sampling distribution is what the model needs.

It is tempting to treat them the way one would for normally distributed data: \(\overline{\mathrm{RT}}\) normal around the predicted mean with variance \(\mathrm{VRT}/n\), \(s^2\) scaled chi-square, and the two independent. All three parts of that hold only if the response times themselves are normal, and diffusion decision times are conspicuously right-skewed. With a symmetric starting point their kurtosis is 8.8 at zero drift and falls towards 3 as drift grows. This inflates the sampling variance of \(s^2\) by a factor of about 3.8 with 100 trials per cell and 3.5 to 3.6 with 10 (less at accuracies above .95; with an asymmetric starting point it is more at the boundary nearer the starting point, up to 5.7 for a relative starting point between .3 and .7, and less at the other), and the skew makes the mean and the variance correlate at about 0.7. A likelihood that ignores both treats each cell as more informative than it is, and the resulting posteriors are too narrow — the credible intervals are not wrong about where the parameter is, but about how sure one should be.

bmm therefore matches the first four cumulants of the first-passage time distribution. Writing \(\kappa_3\) and \(\kappa_4\) for its third and fourth cumulants and \(W = \kappa_4 / n + 2\,\mathrm{VRT}^2 / (n - 1)\) for the exact variance of a sample variance, the observed variance follows \(\mathrm{Gamma}(\mathrm{VRT}^2 / W,\ \mathrm{VRT} / W)\) and the observed mean is normal conditional on it, centred at \(\mathrm{ndt} + \mathrm{MDT} + (\kappa_3 / n)\,W^{-1}(s^2 - \mathrm{VRT})\). Both reduce to the familiar normal and chi-square terms when \(\kappa_3\) and \(\kappa_4\) vanish, so nothing is lost when the decision times happen to be symmetric.

One consequence is worth keeping in mind when designing a study: with a symmetric starting point the decision time does not depend on which boundary was reached, so all \(n\) trials in a cell inform one set of moments. With a free starting point they do differ, which is why the 4-parameter version needs the summaries split by boundary — and why a boundary reached fewer than twice in a cell has no variance to contribute. Such a cell still contributes through the response counts, and so does a cell whose summaries at one boundary are missing, as ezdm_summary_stats() returns them for a boundary reached fewer than min_trials times. bmm() keeps these cells and warns how many there are. Dropping them instead would remove exactly the cells with the most extreme accuracy, and with them information about drift and starting point.

2 Parametrization in the `bmm` package

The bmm package implements two versions of the EZ-diffusion model:

2.1 3-Parameter Version (`version = "3par"`)

Estimates three parameters with a fixed symmetric starting point (\(z_r = 0.5\)):

Parameter	Description	Link Function
`drift`	Drift rate (\(v\)) - rate of evidence accumulation	identity
`bound`	Boundary separation (\(a\)) - response caution	log
`ndt`	Non-decision time (\(T_{er}\)) - encoding + motor time	log

The diffusion constant s is fixed to 1 for identification.

2.2 4-Parameter Version (`version = "4par"`)

Adds estimation of starting point bias:

Parameter	Description	Link Function
`drift`	Drift rate (\(v\))	identity
`bound`	Boundary separation (\(a\))	log
`ndt`	Non-decision time (\(T_{er}\))	log
`zr`	Relative starting point (0-1 scale)	logit

Note that zr = 0.5 indicates no bias (symmetric starting point), zr > 0.5 indicates bias toward the upper boundary, and zr < 0.5 indicates bias toward the lower boundary.

3 Data requirements

Unlike other models in bmm that work with trial-level data, the EZ-diffusion model requires aggregated summary statistics:

Variable	Description
`mean_rt`	Mean reaction time (in seconds)
`var_rt`	Variance of reaction times
`n_upper`	Number of responses hitting the upper boundary
`n_trials`	Total number of trials

For the 3-parameter version, mean_rt and var_rt should be computed from all correct responses pooled together.

For the 4-parameter version, you need separate statistics for upper and lower boundary responses: - mean_rt_upper, var_rt_upper for upper boundary (correct) responses - mean_rt_lower, var_rt_lower for lower boundary (error) responses

4 Preparing data with `ezdm_summary_stats()`

The bmm package provides the ezdm_summary_stats() function to compute the required aggregated statistics from trial-level RT data in the correct format. RT contaminants (fast guesses, lapses, attentional failures) can severely distort mean and variance estimates, leading to biased parameter estimates (Wagenmakers et al. 2008). The function offers three methods to address this:

Method	Description	Use when
`"simple"`	Standard mean/variance	Data is pre-cleaned
`"robust"`	Median + IQR/MAD-based variance	Moderate contamination
`"mixture"`	EM mixture modeling (default)	General use; provides contamination estimates

For detailed coverage of RT contamination handling, including mixture model theory, distribution options, contaminant bounds, and diagnostics, see vignette("bmm_rt_contamination").

4.1 Example with Simulated Data

To illustrate the impact of contamination, we generate trial-level data with known parameters and introduce 5% contaminant trials:

library(bmm)
library(rtdists)
library(dplyr)

set.seed(42)

# Generate clean trial-level data with known parameters
n_trials <- 200
true_drift <- 2.5
true_bound <- 1.2
true_ndt <- 0.25

# Simulate diffusion process
trials_clean <- rdiffusion(
  n = n_trials,
  a = true_bound,
  v = true_drift,
  t0 = true_ndt,
  z = 0.5 * true_bound
)

# Add contaminants (5% fast guesses and lapses)
n_contaminants <- round(n_trials * 0.05)
contaminant_idx <- sample(n_trials, n_contaminants)

trials_contaminated <- trials_clean
trials_contaminated$rt[contaminant_idx] <- runif(n_contaminants, 0.1, 3)
trials_contaminated$response[contaminant_idx] <- sample(c("upper", "lower"),
                                                         n_contaminants,
                                                         replace = TRUE)

trials_contaminated <- trials_contaminated |>
  mutate(
    correct = as.numeric(response == "upper"),
    subject = 1,
    condition = "standard"
  )

We can now compare the three methods. The function takes rt and response as vectors. For grouped operations, use dplyr::group_by() with reframe():

summary_simple <- ezdm_summary_stats(
  trials_contaminated$rt, trials_contaminated$correct,
  method = "simple", version = "3par"
)

summary_robust <- ezdm_summary_stats(
  trials_contaminated$rt, trials_contaminated$correct,
  method = "robust", version = "3par"
)

summary_mixture <- ezdm_summary_stats(
  trials_contaminated$rt, trials_contaminated$correct,
  method = "mixture", version = "3par"
)

comparison <- data.frame(
  Method = c("Simple", "Robust", "Mixture", "TRUE"),
  Mean_RT = c(
    summary_simple$mean_rt,
    summary_robust$mean_rt,
    summary_mixture$mean_rt,
    mean(trials_clean$rt[trials_clean$response == "upper"])
  ),
  Var_RT = c(
    summary_simple$var_rt,
    summary_robust$var_rt,
    summary_mixture$var_rt,
    var(trials_clean$rt[trials_clean$response == "upper"])
  ),
  Contaminant = c(NA, NA, summary_mixture$contaminant_prop, 0.05)
)

print(comparison)
#>    Method   Mean_RT     Var_RT Contaminant
#> 1  Simple 0.5409176 0.10605857          NA
#> 2  Robust 0.4568634 0.02440439          NA
#> 3 Mixture 0.4964177 0.03642890  0.03532347
#> 4    TRUE 0.4755037 0.02050638  0.05000000

The mixture method identifies contaminants and provides more accurate moment estimates.

4.2 Accuracy and Trial Counts

The mixture method also corrects n_upper and n_trials, so that the counts describe the responses the moments were computed from. The simple and robust methods estimate no contaminant proportion and return the raw counts:

summary_raw <- ezdm_summary_stats(
  trials_contaminated$rt, trials_contaminated$correct,
  method = "simple", version = "3par"
)

cat("Raw counts:", summary_raw$n_upper, "of", summary_raw$n_trials,
    "- accuracy", round(summary_raw$n_upper / summary_raw$n_trials, 3), "\n")
#> Raw counts: 186 of 200 - accuracy 0.93
cat("Contaminant-free:", summary_mixture$n_upper, "of", summary_mixture$n_trials,
    "- accuracy", round(summary_mixture$n_upper / summary_mixture$n_trials, 3), "\n")
#> Contaminant-free: 182 of 193 - accuracy 0.943
cat("True accuracy (no contaminants):",
    round(mean(trials_clean$response == "upper"), 3), "\n")
#> True accuracy (no contaminants): 0.945

Splitting the estimated contaminants between the two boundaries needs an assumption about how accurate a contaminant response is. The guess_rate argument carries it and defaults to 0.5, the chance level of a two-choice task. The 4-parameter version estimates contamination per boundary and ignores it.

4.3 Working with the 4-Parameter Version

The 4-parameter version estimates starting point bias (zr) and requires separate statistics for upper and lower boundary responses:

# Generate data with response bias (zr = 0.6, favoring upper boundary)
trials_biased <- rdiffusion(
  n = 200,
  a = 1.2,
  v = 2.0,
  t0 = 0.25,
  z = 0.7 * 1.2
)

trials_biased <- trials_biased |>
  mutate(
    correct = as.numeric(response == "upper"),
    subject = 1
  )

# Compute summary statistics for 4-parameter model
summary_4par <- ezdm_summary_stats(
  trials_biased$rt, trials_biased$correct, version = "4par"
)

summary_4par
#>    mean_rt_upper mean_rt_lower var_rt_upper var_rt_lower n_upper n_trials
#> mu     0.3950143            NA   0.01629836           NA     193      198
#>    contaminant_prop_upper contaminant_prop_lower
#> mu             0.01143807                     NA

Notice the output now includes: - mean_rt_upper, var_rt_upper: Statistics for upper boundary responses - mean_rt_lower, var_rt_lower: Statistics for lower boundary responses

This format is required for fitting the 4-parameter EZ-DM model.

4.4 Complete Workflow: From Trial-Level to Model Fitting

Here’s the complete pipeline from trial-level data to parameter estimation:

# Step 1: Generate or load trial-level data
trial_data <- rdiffusion(
  n = 100 * 20,  # 20 subjects, 100 trials each
  a = 1.5, v = 2.5, t0 = 0.3, z = 0.75
) |>
  mutate(
    subject = rep(1:20, each = 100),
    correct = as.numeric(response == "upper")
  )

# Step 2: Compute summary statistics with robust estimation
summary_data <- trial_data |>
  group_by(subject) |>
  reframe(ezdm_summary_stats(
    rt, correct,
    method = "mixture", distribution = "exgaussian", version = "3par"
  ))

# Step 3: Inspect contaminant proportions (quality control)
summary(summary_data$contaminant_prop)

# Step 4: Fit model with bmm
model <- ezdm("mean_rt", "var_rt", "n_upper", "n_trials", version = "3par")
formula <- bmf(drift ~ 1, bound ~ 1, ndt ~ 1)

fit <- bmm(
  formula = formula,
  data = summary_data,
  model = model,
  backend = "cmdstanr"
)

# Step 5: Examine results
summary(fit)

This workflow ensures robust parameter estimation by properly handling RT contaminants throughout the analysis pipeline.

5 Fitting the model with `bmm`

We start by loading the bmm package:

library(bmm)

5.1 Generating simulated data

For illustration, we will generate simulated data with known parameters using the rezdm() function. As for all other models, bmm provides a density and random generation function for all implemented models. In the case of the ezdm there is no single dependent variable for the density function, which makes efficient illustration challenging. However, the rezdm function allows for evaluating parameter recovery simulations and predicitions based on the ezdm.

set.seed(123)

# Define true parameters for two conditions
true_params <- data.frame(
  condition = c("easy", "hard"),
  drift = c(3.0, 1.5),
  bound = c(1.5, 1.5),
  ndt = c(0.3, 0.3)
)

# Number of subjects and trials per condition
n_subjects <- 20
n_trials <- 100

# Generate data
sim_data <- data.frame()
for (subj in 1:n_subjects) {
  for (cond in 1:2) {
    # Add some random variation for each subject
    drift_subj <- true_params$drift[cond] + rnorm(1, 0, 0.3)
    bound_subj <- true_params$bound[cond] + rnorm(1, 0, 0.1)
    ndt_subj <- true_params$ndt[cond] + rnorm(1, 0, 0.02)

    # Generate summary statistics
    stats <- rezdm(
      n = 1,
      n_trials = n_trials,
      drift = max(drift_subj, 0.1),
      bound = max(bound_subj, 0.5),
      ndt = max(ndt_subj, 0.1),
      version = "3par"
    )

    stats$subject <- subj
    stats$condition <- true_params$condition[cond]
    sim_data <- rbind(sim_data, stats)
  }
}

# Convert condition to factor
sim_data$condition <- factor(sim_data$condition, levels = c("easy", "hard"))

head(sim_data)
#>     mean_rt     var_rt n_upper n_trials subject condition
#> 1 0.5849755 0.02704076      96      100       1      easy
#> 2 0.7349579 0.05481245      89      100       1      hard
#> 3 0.5064378 0.01015706     100      100       2      easy
#> 4 0.6327151 0.07406956      91      100       2      hard
#> 5 0.5354869 0.01897598      98      100       3      easy
#> 6 0.6624182 0.07254991      85      100       3      hard

5.2 Specifying the formula

Next, we specify the model formula using bmmformula() (or bmf()). For this example, we want to estimate separate drift rates for each condition while keeping boundary and non-decision time constant:

model_formula <- bmf(
  drift ~ 0 + condition,
  bound ~ 1,
  ndt ~ 1
)

5.3 Specifying the model

Next, we specify the bmmodelobject for the EZDM model by calling ezdm and passing the respective variable names from our dataset:

model <- ezdm(
  mean_rt = "mean_rt",
  var_rt = "var_rt",
  n_upper = "n_upper",
  n_trials = "n_trials",
  version = "3par"
)

As for all other bmmodel this provides the necessary information for bmm to properly set up the model for estimation.

5.4 Fitting the model

Now we can fit the model using bmm():

fit <- bmm(
  formula = model_formula,
  data = sim_data,
  model = model,
  cores = 4,
  backend = "cmdstanr",
  file = "assets/bmmfit_ezdm_vignette"
)

5.5 Inspecting results

First, check the model summary:

summary(fit)
Loading required namespace: rstan

  Model: ezdm(mean_rt = "mean_rt",
              var_rt = "var_rt",
              n_upper = "n_upper",
              n_trials = "n_trials",
              version = "3par") 
  Links: drift = identity; bound = log; ndt = log; s = log 
Formula: drift ~ 0 + condition
         bound ~ 1
         ndt ~ 1
         s = 0 
   Data: sim_data (Number of observations: 40)
  Draws: 4 chains, each with iter = 1000; warmup = 500; thin = 1;
         total post-warmup draws = 2000

Regression Coefficients:
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
bound_Intercept         0.43      0.01     0.41     0.46 1.00     1231     1134
ndt_Intercept          -1.29      0.02    -1.33    -1.25 1.00     1230     1195
drift_conditioneasy     2.74      0.06     2.63     2.85 1.00     1358     1196
drift_conditionhard     1.37      0.04     1.29     1.44 1.00     1424     1365

Constant Parameters:
                Value
s_Intercept      0.00

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

To provide stable sampling all parameters use log link functions, so we need to exponentiate to get estimates on the natural scale:

# Extract fixed effects
fixef_est <- brms::fixef(fit)

# Transform to natural scale
cat("Drift rate (easy):", exp(fixef_est["drift_conditioneasy", "Estimate"]), "\n")
#> Drift rate (easy): 15.42381
cat("Drift rate (hard):", exp(fixef_est["drift_conditionhard", "Estimate"]), "\n")
#> Drift rate (hard): 3.92211
cat("Boundary:", exp(fixef_est["bound_Intercept", "Estimate"]), "\n")
#> Boundary: 1.543554
cat("Non-decision time:", exp(fixef_est["ndt_Intercept", "Estimate"]), "\n")
#> Non-decision time: 0.2750842

As you can see, these match the generating values well:

print(true_params)
#>   condition drift bound ndt
#> 1      easy   3.0   1.5 0.3
#> 2      hard   1.5   1.5 0.3

5.6 Visualizing posterior distributions

library(tidybayes)
library(dplyr)
library(tidyr)
library(ggplot2)

# Extract posterior draws
draws <- tidybayes::tidy_draws(fit) |>
  select(starts_with("b_")) |>
  mutate(across(everything(), exp))

# Rename for plotting
colnames(draws) <- gsub("b_", "", colnames(draws))

# True values for comparison
true_vals <- data.frame(
  par = c("drift_conditioneasy", "drift_conditionhard", "bound_Intercept", "ndt_Intercept"),
  value = c(3.0, 1.5, 1.5, 0.3)
)

# Plot
draws |>
  pivot_longer(everything(), names_to = "par", values_to = "value") |>
  filter(par != "Intercept") |>
  ggplot(aes(x = value, y = par)) +
  stat_halfeye(normalize = "groups", fill = "steelblue", alpha = 0.7) +
  geom_point(data = true_vals, aes(x = value, y = par),
             color = "red", shape = 18, size = 4) +
  labs(x = "Parameter value", y = "",
       title = "Posterior distributions of EZDM parameters",
       subtitle = "Red diamonds indicate true generating values") +
  coord_cartesian(xlim = c(0,5)) +
  theme_minimal()

6 Posterior predictive checks

The EZ-diffusion likelihood constrains three statistics per cell — the mean RT, the RT variance, and the number of upper-boundary (correct) responses — but by default pp_check() only shows the mean RT. The resp_var argument selects the other observables of the likelihood; pp_check_vars() lists the available checks:

pp_check_vars(fit)
#>   resp_var                         label default_type          slot default
#> 1  mean_rt            Mean response time dens_overlay             Y    TRUE
#> 2   var_rt        Response time variance dens_overlay        vreal1   FALSE
#> 3  mean_pc Proportion of upper responses  bars_binned vint1, trials   FALSE

The mean_pc check is the proportion of correct responses (n_upper / n_trials), which stays comparable across cells even when the number of trials differs. With resp_var = "all" all checks are drawn from one shared joint simulation. Each panel compares how a statistic is distributed across all cells: the RT panels overlay the density of the observed statistic (dark line) on the densities of the simulated datasets, and the accuracy panel is a histogram of the observed proportions (bars) with the predicted number of cells per bin (points with 90% intervals):

pp_check(fit, resp_var = "all", ndraws = 50)

The observed mean RTs are spread more evenly than the simulated ones, which cluster around one value per condition. The data were simulated with subject-level variation in all three parameters, which this model omits; the model with random effects below includes it. To compare each cell with its own predictive interval instead, select a single check with type = "intervals", e.g. pp_check(fit, resp_var = "mean_rt", type = "intervals").

7 Comparing conditions

As bmm feeds seamlessly into brms, we can use Bayesian hypothesis testing via Savage-Dickey density ratios to compare drift rates between conditions:

# Test if drift rate is higher in easy vs hard condition
brms::hypothesis(fit, "exp(drift_conditioneasy) > exp(drift_conditionhard)")
#> Hypothesis Tests for class b:
#>                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (exp(drift_condit... > 0    11.52      0.88    10.14    13.05        Inf
#>   Post.Prob Star
#> 1         1    *
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

8 Adding random effects

For hierarchical models with random effects across subjects:

# Formula with random effects
model_formula_re <- bmf(
  drift ~ 0 + condition + (0 + condition | subject),
  bound ~ 1 + (1 | subject),
  ndt ~ 1 + (1 | subject)
)

fit_re <- bmm(
  formula = model_formula_re,
  data = sim_data,
  model = model,
  cores = 4,
  backend = "cmdstanr"
)

References

Chávez De la Peña, Adriana F., and Joachim Vandekerckhove. 2025. “An EZ Bayesian Hierarchical Drift Diffusion Model for Response Time and Accuracy.” Psychonomic Bulletin & Review, ahead of print. https://doi.org/10.3758/s13423-025-02729-y.

Ratcliff, Roger. 1978. “A Theory of Memory Retrieval.” Psychological Review 85 (2): 59–108. https://doi.org/10.1037/0033-295X.85.2.59.

Srivastava, Vaibhav, Philip Holmes, and Patrick Simen. 2016. “A Martingale Analysis of First Passage Times of Time-Dependent Wiener Diffusion Models.” Journal of Mathematical Psychology 77: 94–110. https://doi.org/10.1016/j.jmp.2016.10.001.

Wagenmakers, Eric-Jan, Han L. J. van der Maas, Conor V. Dolan, and Raoul P. P. P. Grasman. 2008. “EZ Does It! Extensions of the EZ-diffusion Model.” Psychonomic Bulletin & Review 15 (6): 1229–35. https://doi.org/10/bwq9dw.

Wagenmakers, Eric-Jan, Han L. J. Van Der Maas, and Raoul P. P. P. Grasman. 2007. “An EZ-diffusion Model for Response Time and Accuracy.” Psychonomic Bulletin & Review 14 (1): 3–22. https://doi.org/10.3758/BF03194023.

Gidon Frischkorn

Ven Popov

1 Introduction to the model

1.1 The Diffusion Model

1.2 The EZ Equations

1.3 The 4-Parameter Extension

1.4 Bayesian Estimation

1.5 The sampling distribution of the summary statistics

2 Parametrization in the `bmm` package

2.1 3-Parameter Version (`version = "3par"`)

2.2 4-Parameter Version (`version = "4par"`)

3 Data requirements

4 Preparing data with `ezdm_summary_stats()`

4.1 Example with Simulated Data

4.2 Accuracy and Trial Counts

4.3 Working with the 4-Parameter Version

4.4 Complete Workflow: From Trial-Level to Model Fitting

5 Fitting the model with `bmm`

5.1 Generating simulated data

5.2 Specifying the formula

5.3 Specifying the model

5.4 Fitting the model

5.5 Inspecting results

5.6 Visualizing posterior distributions

6 Posterior predictive checks

7 Comparing conditions

8 Adding random effects

References

The EZ-Diffusion Model (EZDM)

Gidon Frischkorn

Ven Popov

1 Introduction to the model

1.1 The Diffusion Model

1.2 The EZ Equations

1.3 The 4-Parameter Extension

1.4 Bayesian Estimation

1.5 The sampling distribution of the summary statistics

2 Parametrization in the bmm package

2.1 3-Parameter Version (version = "3par")

2.2 4-Parameter Version (version = "4par")

3 Data requirements

4 Preparing data with ezdm_summary_stats()

4.1 Example with Simulated Data

4.2 Accuracy and Trial Counts

4.3 Working with the 4-Parameter Version

4.4 Complete Workflow: From Trial-Level to Model Fitting

5 Fitting the model with bmm

5.1 Generating simulated data

5.2 Specifying the formula

5.3 Specifying the model

5.4 Fitting the model

5.5 Inspecting results

5.6 Visualizing posterior distributions

6 Posterior predictive checks

7 Comparing conditions

8 Adding random effects

References

2 Parametrization in the `bmm` package

2.1 3-Parameter Version (`version = "3par"`)

2.2 4-Parameter Version (`version = "4par"`)

4 Preparing data with `ezdm_summary_stats()`

5 Fitting the model with `bmm`