Fit Bayesian MPT Models with brms & Stan

Fit a Bayesian multinomial processing tree (MPT) model using brms::brm() (so ultimately Stan) to trial-level data, potentially including a multilevel/hierarchical model structure. To fit an MPT, first create a model object using make_mpt() and then pass the model object, the data, as well as a formula providing a symbolic description of the regression (i.e., fixed-effect) and hierarchical (i.e., random-effect) structure. The returned object can be used for model post-processing.

Usage

mpt(
  formula,
  data,
  tree,
  model,
  default_prior_intercept = "normal(0, 1)",
  default_prior_coef = "normal(0, 0.5)",
  default_priors = TRUE,
  data_format = "long",
  log_p = FALSE,
  link = "probit",
  ...
)

Arguments

formula

Either a formula that is applied to each MPT model parameter or a mpt_formula object (created by mpt_formula()) specifying separate formulas for each parameter. If a formula, the LHS needs to specify the response variable.

data

data.frame containing the variables in formula. Data needs to be on an observation-level (i.e., each row is one response/observation) and cannot be aggregated in any way. TODO: change this

tree

one-sided formula or character specifying the variable in data indicating the tree (or item type) of a given observation. The values of the tree variable need to match the names of the trees in model. Can be omitted for models with only one tree.

model

mpt_model object as created by make_mpt(), ignored if formula is an object of mpt_formula.

default_prior_intercept

character string describing the prior applied to the fixed-effect intercepts for each MPT model parameter on the unconstrained scale (if default_priors = TRUE). The default, "normal(0, 1)" implies a flat prior on the MPT parameter scale.

default_prior_coef

character string describing the prior applied to the non-intercept fixed-effect parameters for each MPT model parameter on the unconstrained scale (if default_priors = TRUE).

default_priors

logical value indicating whether (the default, TRUE) or not (FALSE) the priors specified via the default_prior_intercept and default_prior_coef argument should be applied.

data_format

character string indicating whether the formula is to be generated for fitting data in long format / non-aggregated data (long, the default), where a single variable contains trial-level responses, or for data in wide format / aggregated data (wide / aggregated), where a separate column for each response category contains the respective . frequency. Only used if a formula object (not an mpt_formula) is given as input and ignored otherwise.

log_p

logical value indicating whether the likelihood should be evaluated with probabilities (the default, FALSE) or log probabilities (TRUE). Setting log_p to TRUE can help in case of convergence issues but might be slower.

link

character specifying the link function for transforming from unconstrained space to MPT model parameter (i.e., 0 to 1) space. Default is "probit".

...

Further arguments passed to brms::brm() such as prior, chains, iter, warmup, and cores.

Value

A fitted model object returned from brms::brm() of class brmsfit with additional class mpt_fit and some extra slots (i.e., call, mpt_formula, and orig_data). This object should seemlessly interact with the brms ecosystem (e.g., can be used for obtaining posterior predictive values and information criteria, see Examples).

Examples

if (FALSE) { # \dontrun{
### Step 1: Specify model using make_mpt()
## model in easy format (with model specified as text)
## unsure-extended 2-high threshold model for recognition memory
u2htm <- "
# Old Items
Do + (1 - Do) * (1 - g1) * g2
(1 - Do) * g1
(1 - Do) * (1 - g1) * (1 - g2)

# New Items
(1 - Dn) * (1 - g1) * g2
(1 - Dn) * g1
Dn + (1 - Dn) * (1 - g1) * (1 - g2)
"

# need to specify tree names and category names for easy format
u2htsm_model <- make_mpt(text = u2htm, trees = c("old", "new"),
                         categories = rep(c("old", "unsure", "new"), 2))
u2htsm_model


### Step 2: Fit model with same formula for all model parameters

## we fit to data from Singmann, Kellen, & Klauer (2013, CogSci-Proc):
# Investigating the Other-Race Effect of Germans towards Turks and Arabs using
# Multinomial Processing Tree Models
# http://singmann.org/download/publications/SKK-CogSci2013.pdf
str(skk13)

## here we use simplified syntax without any random effects for fitting speed.
fit_fast <- mpt(resp ~ race, data = skk13, model = u2htsm_model,
                tree = "type",
                cores = min(c(4, parallel::detectCores()))) ## uses multicore
## a more appropriate formula would be: resp ~ race + (race|s|id) + (1|p|stim)
## (i.e., crossed-random effects and full correlations among random terms)

### Step 3: Inspect results
fit_fast # in output, no parameter name refers to first MPT parameter (here: Dn)

## Check model fit
ppp_test(fit_fast)

## we can use mpt_emmeans() for marginal effects of all MPT model parameters
# requires package emmeans to be installed
mpt_emmeans(fit_fast, "race")

## Alternatively, we can invoke emmeans() individually for each MPT model parameter
# Default is first MPT model parameter (here: Dn)
emmeans::emmeans(fit_fast, "race", type = "response")
# using the dpar argument, we can get output for any other model parameter
emmeans::emmeans(fit_fast, "race", type = "response", dpar = "Do")
emmeans::emmeans(fit_fast, "race", type = "response", dpar = "g1x")
emmeans::emmeans(fit_fast, "race", type = "response", dpar = "g2x")

## We can also use all brms post-processing
# information criteria (e.g., loo)
(loo_model <- loo(fit_fast))
# posterior predictive checks
pp_check(fit_fast, type = "bars_grouped", group = "mpt_tree", ndraws = 100)
# get posterior mean predictions
pepred <- posterior_epred(fit_fast)
str(pepred)  ## dimensions are: samples, observations, probability (in tree)

### Alternative Step 2: First specify formula, then fit model
## Step 2a: Specify formula for each parameter separately using mpt_formula()
## (formula could be different for each model parameter)
u2htm_formula <- mpt_formula(
  Do ~ race + (race|s|id) + (1|p|stim),
  Dn ~ race + (race|s|id) + (1|p|stim),
  g1x ~ race + (race|s|id) + (1|p|stim),
  g2x ~ race + (race|s|id) + (1|p|stim),
  response = ~ resp,
  model = u2htsm_model
)
u2htm_formula

## Step 2b: Fit model using formula (takes rather long)
fit_slow <- mpt(u2htm_formula, data = skk13,
                tree = "type",
                cores = min(c(4, parallel::detectCores()))) ## uses multicore
fit_slow
} # }