Title:	Bayesian Hierarchical Analysis of Cognitive Models of Choice
Version:	3.2.0
Description:	Fit Bayesian (hierarchical) cognitive models using a linear modeling language interface using particle Metropolis Markov chain Monte Carlo sampling with Gibbs steps. The diffusion decision model (DDM), linear ballistic accumulator model (LBA), racing diffusion model (RDM), and the lognormal race model (LNR) are supported. Additionally, users can specify their own likelihood function and/or choose for non-hierarchical estimation, as well as for a diagonal, blocked or full multivariate normal group-level distribution to test individual differences. Prior specification is facilitated through methods that visualize the (implied) prior. A wide range of plotting functions assist in assessing model convergence and posterior inference. Models can be easily evaluated using functions that plot posterior predictions or using relative model comparison metrics such as information criteria or Bayes factors. References: Stevenson et al. (2024) <doi:10.31234/osf.io/2e4dq>.
License:	GPL (≥ 3)
URL:	https://ampl-psych.github.io/EMC2/, https://github.com/ampl-psych/EMC2
BugReports:	https://github.com/ampl-psych/EMC2/issues
Encoding:	UTF-8
RoxygenNote:	7.3.2
VignetteBuilder:	knitr
Suggests:	testthat (≥ 3.0.0), vdiffr, knitr, rmarkdown
Config/testthat/edition:	3
Imports:	abind, coda, graphics, grDevices, magic, MASS, matrixcalc, methods, msm, mvtnorm, parallel, stats, Matrix, Rcpp, Brobdingnag, corrplot, colorspace, psych, utils, lpSolve, WienR
LinkingTo:	Rcpp, RcppArmadillo
Depends:	R (≥ 3.5.0)
LazyData:	true
Config/testthat/parallel:	true
NeedsCompilation:	yes
Packaged:	2025-06-27 11:48:17 UTC; nstevenson
Author:	Niek Stevenson [aut, cre], Michelle Donzallaz [aut], Andrew Heathcote [aut], Steven Miletić [ctb], Raphael Hartmann [ctb], Karl C. Klauer [ctb], Steven G. Johnson [ctb], Jean M. Linhart [ctb], Brian Gough [ctb], Gerard Jungman [ctb], Rudolf Schuerer [ctb], Przemyslaw Sliwa [ctb], Jason H. Stover [ctb]
Maintainer:	Niek Stevenson <niek.stevenson@gmail.com>
Repository:	CRAN
Date/Publication:	2025-06-27 12:10:05 UTC

EMC2: Bayesian Hierarchical Analysis of Cognitive Models of Choice

Description

Fit Bayesian (hierarchical) cognitive models using a linear modeling language interface using particle Metropolis Markov chain Monte Carlo sampling with Gibbs steps. The diffusion decision model (DDM), linear ballistic accumulator model (LBA), racing diffusion model (RDM), and the lognormal race model (LNR) are supported. Additionally, users can specify their own likelihood function and/or choose for non-hierarchical estimation, as well as for a diagonal, blocked or full multivariate normal group-level distribution to test individual differences. Prior specification is facilitated through methods that visualize the (implied) prior. A wide range of plotting functions assist in assessing model convergence and posterior inference. Models can be easily evaluated using functions that plot posterior predictions or using relative model comparison metrics such as information criteria or Bayes factors. References: Stevenson et al. (2024) doi:10.31234/osf.io/2e4dq.

Author(s)

Maintainer: Niek Stevenson niek.stevenson@gmail.com (ORCID)

Authors:

• Michelle Donzallaz

• Andrew Heathcote

Other contributors:

• Steven Miletić [contributor]

• Raphael Hartmann [contributor]

• Karl C. Klauer [contributor]

• Steven G. Johnson [contributor]

• Jean M. Linhart [contributor]

• Brian Gough [contributor]

• Gerard Jungman [contributor]

• Rudolf Schuerer [contributor]

• Przemyslaw Sliwa [contributor]

• Jason H. Stover [contributor]

See Also

Useful links:

• https://ampl-psych.github.io/EMC2/

• https://github.com/ampl-psych/EMC2

• Report bugs at https://github.com/ampl-psych/EMC2/issues

The Diffusion Decision Model

Description

Model file to estimate the Diffusion Decision Model (DDM) in EMC2.

Usage

DDM()

Details

Model files are almost exclusively used in design().

Default values are used for all parameters that are not explicitly listed in the formula argument of design().They can also be accessed with DDM()$p_types.

a, t0, sv, st0, s are sampled on the log scale because these parameters are strictly positive, Z, SZ and DP are sampled on the probit scale because they should be strictly between 0 and 1.

Z is estimated as the ratio of bias to one boundary where 0.5 means no bias. DP comprises the difference in non-decision time for each response option.

Conventionally, s is fixed to 1 to satisfy scaling constraints.

See Ratcliff, R., & McKoon, G. (2008). The diffusion decision model: theory and data for two-choice decision tasks. Neural computation, 20(4), 873-922. doi:10.1162/neco.2008.12-06-420.

Value

A model list with all the necessary functions for EMC2 to sample

Examples

design_DDMaE <- design(data = forstmann,model=DDM,
                           formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                           constants=c(s=log(1)))
# For all parameters that are not defined in the formula, default values are assumed
# (see Table above).

The GNG (go/nogo) Diffusion Decision Model

Description

In the GNG paradigm one of the two possible choices results in a response being withheld (a non-response), which is indicated in the data by an NA for the rt, with the corresponding level of the R (response) factor still being specified. For example, suppose the go response is coded as "yes" and nogo is coded as "no", then for a non-response (R,rt) = ("no",NA) and for a response e.g., (R,rt) = ("yes",1.36). The GNG paradigm must also have a response

Usage

DDMGNG()

Details

The model used is described in the following paper, with the addition of modeling the TIMEOUT (which is considered but not used in this paper).

Gomez, P., Ratcliff, R., & Perea, M. (2007). A Model of the Go/No-Go Task. Journal of Experimental Psychology: General, 136(3), 389–413. https://doi.org/10.1037/0096-3445.136.3.389

The likelihood of non-responses requires and evaluation of the DDM cdf, specifically 1 - p(hitting the yes boundary before TIMEOUT).

To use these models three functions must be supplied in the design's function argument with the names TIMEOUT, Rnogo and Rgo. For example, assuming a 2.5 second timeout, and R factor with levels c("no","yes") and "no" mapping to a non-response.

TIMEOUT=function(d)rep(2.5,nrow(d)) Rnogo=function(d)factor(rep("no",nrow(d)),levels=c("no","yes")) Rgo=function(d)factor(rep("yes",nrow(d)),levels=c("no","yes")))

See the help for DDM for further details. At present this model is not fully implemented in C, so is a little slower to use than the DDM, but not greatly.

Value

A model list with all the necessary functions to sample

window (i.e., a length of time, TIMEOUT period, after which withholding is

assumed).

Examples

dGNG <- design(Rlevels = c("left","right"),
               factors=list(subjects=1,S=c("left","right")),
               functions=list(
               TIMEOUT=function(d)rep(2.5,nrow(d)),
               # no go response level
               Rnogo=function(d)factor(rep("left",nrow(d)),levels=c("left","right")),
               # go response level
               Rgo=function(d)factor(rep("right",nrow(d)),levels=c("left","right"))),
               formula=list(v~S,a~1, Z~1, t0~1),
               model=DDMGNG)

p_vector <- sampled_pars(dGNG)

The Linear Ballistic Accumulator model

Description

Model file to estimate the Linear Ballistic Accumulator (LBA) in EMC2.

Usage

LBA()

Details

Model files are almost exclusively used in design().

Default values are used for all parameters that are not explicitly listed in the formula argument of design().They can also be accessed with LBA()$p_types.

All parameters are estimated on the log scale, except for the drift rate which is estimated on the real line.

Conventionally, sv is fixed to 1 to satisfy scaling constraints.

The b = B + A parameterization ensures that the response threshold is always higher than the between trial variation in start point of the drift rate.

Because the LBA is a race model, it has one accumulator per response option. EMC2 automatically constructs a factor representing the accumulators lR (i.e., the latent response) with level names taken from the R column in the data.

The lR factor is mainly used to allow for response bias, analogous to Z in the DDM. For example, in the LBA, response thresholds are determined by the B parameters, so B~lR allows for different thresholds for the accumulator corresponding to left and right stimuli (e.g., a bias to respond left occurs if the left threshold is less than the right threshold). For race models, the design() argument matchfun can be provided, a function that takes the lR factor (defined in the augmented data (d) in the following function) and returns a logical defining the correct response. In the example below, the match is simply such that the S factor equals the latent response factor: matchfun=function(d)d$S==d$lR. Then matchfun is used to automatically create a latent match (lM) factor with levels FALSE (i.e., the stimulus does not match the accumulator) and TRUE (i.e., the stimulus does match the accumulator). This is added internally and can also be used in model formula, typically for parameters related to the rate of accumulation.

Brown, S. D., & Heathcote, A. (2008). The simplest complete model of choice response time: Linear ballistic accumulation. Cognitive Psychology, 57(3), 153-178. https://doi.org/10.1016/j.cogpsych.2007.12.002

Value

A model list with all the necessary functions for EMC2 to sample

Examples

# When working with lM it is useful to design  an "average and difference"
# contrast matrix, which for binary responses has a simple canonical from:
ADmat <- matrix(c(-1/2,1/2),ncol=1,dimnames=list(NULL,"d"))
# We also define a match function for lM
matchfun=function(d)d$S==d$lR
# We now construct our design, with v ~ lM and the contrast for lM the ADmat.
design_LBABE <- design(data = forstmann,model=LBA,matchfun=matchfun,
                       formula=list(v~lM,sv~lM,B~E+lR,A~1,t0~1),
                       contrasts=list(v=list(lM=ADmat)),constants=c(sv=log(1)))
# For all parameters that are not defined in the formula, default values are assumed
# (see Table above).

The Log-Normal Race Model

Description

Model file to estimate the Log-Normal Race Model (LNR) in EMC2.

Usage

LNR()

Details

Model files are almost exclusively used in design().

Default values are used for all parameters that are not explicitly listed in the formula argument of design().They can also be accessed with LNR()$p_types.

Because the LNR is a race model, it has one accumulator per response option. EMC2 automatically constructs a factor representing the accumulators lR (i.e., the latent response) with level names taken from the R column in the data.

In design(), matchfun can be used to automatically create a latent match (lM) factor with levels FALSE (i.e., the stimulus does not match the accumulator) and TRUE (i.e., the stimulus does match the accumulator). This is added internally and can also be used in the model formula, typically for parameters related to the rate of accumulation (see the example below).

Rouder, J. N., Province, J. M., Morey, R. D., Gomez, P., & Heathcote, A. (2015). The lognormal race: A cognitive-process model of choice and latency with desirable psychometric properties. Psychometrika, 80, 491-513. https://doi.org/10.1007/s11336-013-9396-3

Value

A model list with all the necessary functions for EMC2 to sample

Examples

# When working with lM it is useful to design  an "average and difference"
# contrast matrix, which for binary responses has a simple canonical from:
ADmat <- matrix(c(-1/2,1/2),ncol=1,dimnames=list(NULL,"d"))
# We also define a match function for lM
matchfun=function(d)d$S==d$lR
# We now construct our design, with v ~ lM and the contrast for lM the ADmat.
design_LNRmE <- design(data = forstmann,model=LNR,matchfun=matchfun,
                       formula=list(m~lM + E,s~1,t0~1),
                       contrasts=list(m=list(lM=ADmat)))
# For all parameters that are not defined in the formula, default values are assumed
# (see Table above).

GLM model for fMRI data

Description

Creates a model specification for fMRI data using a normal distribution. This model assumes that the observed BOLD signal follows a normal distribution with a mean determined by the design matrix and betas, and a standard deviation parameter for noise.

Usage

MRI()

Details

The model uses a normal distribution to model fMRI BOLD signals. Beta parameters represent the effect sizes for different conditions, and the sd parameter represents the standard deviation of the noise.

The log-likelihood function centers the predicted values by subtracting the mean, which helps with model identifiability.

Value

A list containing model specification

Examples

# Create a normal MRI model specification
model_spec <- MRI()

# Access model parameters
model_spec$p_types

Create an AR(1) GLM model for fMRI data

Description

This function creates a model specification for MRI data with an AR(1) error structure. The model includes beta parameters for the design matrix, a rho parameter for the autocorrelation, and a standard deviation parameter for the noise.

Usage

MRI_AR1()

Details

The AR(1) model accounts for temporal autocorrelation in the data, where each timepoint is correlated with the previous timepoint according to the rho parameter.

Value

A list containing the model specifications

Examples

# Create an AR(1) GLM model for fMRI data
model_spec <- MRI_AR1()

# Access model parameters
model_spec$p_types

The Racing Diffusion Model

Description

Model file to estimate the Racing Diffusion Model (RDM), also known as the Racing Wald Model.

Usage

RDM()

Details

Model files are almost exclusively used in design().

Default values are used for all parameters that are not explicitly listed in the formula argument of design().They can also be accessed with RDM()$p_types.

All parameters are estimated on the log scale.

The parameterization b = B + A ensures that the response threshold is always higher than the between trial variation in start point.

Conventionally, s is fixed to 1 to satisfy scaling constraints.

Because the RDM is a race model, it has one accumulator per response option. EMC2 automatically constructs a factor representing the accumulators lR (i.e., the latent response) with level names taken from the R column in the data.

The lR factor is mainly used to allow for response bias, analogous to Z in the DDM. For example, in the RDM, response thresholds are determined by the B parameters, so B~lR allows for different thresholds for the accumulator corresponding to "left" and "right" stimuli, for example, (e.g., a bias to respond left occurs if the left threshold is less than the right threshold).

For race models in general, the argument matchfun can be provided in design(). One needs to supply a function that takes the lR factor (defined in the augmented data (d) in the following function) and returns a logical defining the correct response. In the example below, this is simply whether the S factor equals the latent response factor: matchfun=function(d)d$S==d$lR. Using matchfun a latent match factor (lM) with levels FALSE (i.e., the stimulus does not match the accumulator) and TRUE (i.e., the stimulus does match the accumulator). This is added internally and can also be used in model formula, typically for parameters related to the rate of accumulation.

Tillman, G., Van Zandt, T., & Logan, G. D. (2020). Sequential sampling models without random between-trial variability: The racing diffusion model of speeded decision making. Psychonomic Bulletin & Review, 27(5), 911-936. https://doi.org/10.3758/s13423-020-01719-6

Value

A list defining the cognitive model

Examples

# When working with lM it is useful to design  an "average and difference"
# contrast matrix, which for binary responses has a simple canonical from:
ADmat <- matrix(c(-1/2,1/2),ncol=1,dimnames=list(NULL,"d"))
# We also define a match function for lM
matchfun=function(d)d$S==d$lR
# We now construct our design, with v ~ lM and the contrast for lM the ADmat.
design_RDMBE <- design(data = forstmann,model=RDM,matchfun=matchfun,
                       formula=list(v~lM,s~lM,B~E+lR,A~1,t0~1),
                       contrasts=list(v=list(lM=ADmat)),constants=c(s=log(1)))
# For all parameters that are not defined in the formula, default values are assumed
# (see Table above).

Gaussian Signal Detection Theory Model for Binary Responses

Description

Discrete binary choice based on continuous Gaussian latent, with no rt (rt must be set to NA in data).

Usage

SDT()

Details

Model parameters are: mean (unbounded) sd (log scale) and threshold (unbounded).

For identifiability in one condition two parameters must be fixed (conventionally mean=0 and sd = 1). When used with data that records only accuracy (so reponse bias cannot be evaluated) a single threshold must be assumed and fixed (e.g., threshold = 0).

At present this model is not fully implemented in C, but as its likelihood requires only pnorm evaluation it is quite fast.

Value

A model list with all the necessary functions to sample

Examples

dprobit <- design(Rlevels = c("left","right"),
           factors=list(subjects=1,S=c("left","right")),
           formula=list(mean ~ 0+S, sd ~ 1,threshold ~ 1),
           matchfun=function(d)d$S==d$lR,
           constants=c(sd=log(1),threshold=0),
           model=SDT)

p_vector <- sampled_pars(dprobit)

Reorder MCMC Samples of Factor Loadings

Description

This function reorders MCMC samples of factor loadings to address the label switching problem in Bayesian factor analysis. It implements a parallelized version of the code and algorithm proposed by Papastamoulis and Ntzoufras (2022)

Usage

align_loadings(
  emc = NULL,
  lambda = NULL,
  n_cores = 1,
  verbose = TRUE,
  rotate_fun = NULL
)

Arguments

Value

A list containing:

References

Papastamoulis, P., & Ntzoufras, I. (2022). On the identifiability of Bayesian factor analytic models. Statistical Computing, 32(2), 1-29. doi: 10.1007/s11222-022-10084-4

Examples

# This function works natively with emc objects, but also factor arrays:
# Simulate a small example with 5 variables, 2 factors, and 10 MCMC iterations
set.seed(123)
p <- 5  # Number of variables
q <- 2  # Number of factors
n <- 10 # Number of MCMC iterations

# Create random factor loadings with label switching
lambda <- array(0, dim = c(p, q, n))
for (i in 1:n) {
  # Generate base loadings
  base_loadings <- matrix(rnorm(p*q, 0, 0.5), p, q)
  base_loadings[1:3, 1] <- abs(base_loadings[1:3, 1]) + 0.5  # Strong loadings on factor 1
  base_loadings[4:5, 2] <- abs(base_loadings[4:5, 2]) + 0.5  # Strong loadings on factor 2

  # Randomly switch labels and signs
  if (runif(1) > 0.5) {
    # Switch factor order
    base_loadings <- base_loadings[, c(2, 1)]
  }
  if (runif(1) > 0.5) {
    # Switch sign of factor 1
    base_loadings[, 1] <- -base_loadings[, 1]
  }
  if (runif(1) > 0.5) {
    # Switch sign of factor 2
    base_loadings[, 2] <- -base_loadings[, 2]
  }

  lambda[,,i] <- base_loadings
}

# Align the loadings
result <- align_loadings(lambda = lambda, verbose = TRUE, n_cores = 1)

# Examine the aligned loadings
print(result)

Automatically Thin an emc Object

Description

Uses the effective sample size of selection to determine how much to optimally thin an emc object

Usage

## S3 method for class 'emc'
auto_thin(emc, stage = "sample", selection = c("alpha", "mu"), ...)

auto_thin(emc, stage = "sample", selection = c("alpha", "mu"), ...)

Arguments

MCMC Chain Iterations

Description

Returns a matrix with the number of samples per chain for each stage that is present in the emc object (i.e., preburn, burn, adapt, sample). The number of rows of the matrix reflects the number of chains and the number of columns the number of sampling stages.

Usage

chain_n(emc)

Arguments

Value

A matrix

Examples

chain_n(samples_LNR)

Convergence Checks for an emc Object

Description

Runs a series of convergence checks, prints statistics to the console, and makes traceplots of the worst converged parameter per selection.

Usage

## S3 method for class 'emc'
check(
  emc,
  selection = c("mu", "sigma2", "alpha"),
  digits = 3,
  plot_worst = TRUE,
  ...
)

check(emc, ...)

Arguments

Details

Note that the Rhat is calculated by doubling the number of chains by first splitting chains into first and second half, so it also a test of stationarity.

Efficiency of sampling is indicated by the effective sample size (ESS) (from the coda R package). Full range of possible samples manipulations described in get_pars.

Value

a list with the statistics for the worst converged parameter per selection

Examples

check(samples_LNR)

Information Criteria and Marginal Likelihoods

Description

Returns the BPIC/DIC or marginal deviance (-2*marginal likelihood) for a list of samples objects.

Usage

compare(
  sList,
  stage = "sample",
  filter = NULL,
  use_best_fit = TRUE,
  BayesFactor = TRUE,
  cores_for_props = 4,
  cores_per_prop = 1,
  print_summary = TRUE,
  digits = 0,
  digits_p = 3,
  ...
)

Arguments

Value

Matrix of effective number of parameters, mean deviance, deviance of mean, DIC, BPIC, Marginal Deviance (if BayesFactor=TRUE) and associated weights.

Examples

compare(list(samples_LNR), cores_for_props = 1)
# Typically we would define a list of two (or more) different models:
# # Here the full model is an emc object with the hypothesized effect
# # The null model is an emc object without the hypothesized effect
# design_full <- design(data = forstmann,model=DDM,
#                            formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
#                            constants=c(s=log(1)))
# # Now without a ~ E
# design_null <- design(data = forstmann,model=DDM,
#                            formula =list(v~0+S,a~1, t0~1, s~1, Z~1, sv~1, SZ~1),
#                            constants=c(s=log(1)))
#
# full_model <- make_emc(forstmann, design_full)
# full_model <- fit(full_model)
#
# null_model <- make_emc(forstmann, design_null)
# null_model <- fit(null_model)
# sList <- list(full_model, null_model)
# # By default emc uses 4 cores to parallelize marginal likelihood estimation across proposals
# # So cores_per_prop = 3 results in 12 cores used.
# compare(sList, cores_per_prop = 3)

Information Criteria For Each Participant

Description

Returns the BPIC/DIC based model weights for each participant in a list of samples objects

Usage

compare_subject(
  sList,
  stage = "sample",
  filter = 0,
  use_best_fit = TRUE,
  print_summary = TRUE,
  digits = 3
)

Arguments

Value

List of matrices for each subject of effective number of parameters, mean deviance, deviance of mean, DIC, BPIC and associated weights.

Examples

# For a broader illustration see `compare`.
# Here we just take two times the same model, but normally one would compare
# different models
compare_subject(list(m0 = samples_LNR, m1 = samples_LNR))

Anova Style Contrast Matrix

Description

Similar to contr.helmert, but then scaled to estimate differences between conditions. Use in design().

Usage

contr.anova(n)

Arguments

Value

A contrast matrix.

Examples

{
design_DDMaE <- design(data = forstmann,model=DDM, contrasts = list(E = contr.anova),
formula =list(v~S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
constants=c(s=log(1)))
}

Contrast Enforcing Equal Prior Variance on each Level

Description

Typical contrasts impose different levels of marginal prior variance for the different levels. This contrast can be used to ensure that each level has equal marginal priors (Rouder, Morey, Speckman, & Province; 2012).

Usage

contr.bayes(n)

Arguments

Value

A contrast matrix.

Examples

{
design_DDMaE <- design(data = forstmann,model=DDM, contrasts = list(E = contr.bayes),
formula =list(v~S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
constants=c(s=log(1)))
}

Contrast Enforcing Decreasing Estimates

Description

Each level will be estimated as a reduction from the previous level

Usage

contr.decreasing(n)

Arguments

Value

a contrast matrix.

Examples

{
design_DDMaE <- design(data = forstmann,model=DDM, contrasts = list(E = contr.decreasing),
formula =list(v~S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
constants=c(s=log(1)))
}

Contrast Enforcing Increasing Estimates

Description

Each level will be estimated additively from the previous level

Usage

contr.increasing(n)

Arguments

Value

a contrast matrix.

Examples

{
design_DDMaE <- design(data = forstmann,model=DDM, contrasts = list(E = contr.increasing),
formula =list(v~S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
constants=c(s=log(1)))
}

Convolve Events with HRF to Construct Design Matrices

Description

This function convolves events with the HRF to construct design matrices for fMRI analysis.

Usage

convolve_design_matrix(
  timeseries,
  events,
  factors = NULL,
  contrasts = NULL,
  covariates = NULL,
  add_constant = TRUE,
  hrf_model = "glover",
  cell_coding = NULL,
  scale = TRUE,
  high_pass = TRUE,
  high_pass_model = "cosine",
  cut_off = 1e-12
)

Arguments

Value

A list containing the design matrices

Examples

# Generate a simple example timeseries
ts <- data.frame(
  subjects = rep(1, 100),
  run = rep(1, 100),
  time = seq(0, 99),
  ROI1 = rnorm(100)
)

# Generate example events
events <- data.frame(
  subjects = rep(1, 4),
  run = rep(1, 4),
  onset = c(10, 30, 50, 70),
  duration = rep(0.5, 4),
  event_type = c("hard", "easy", "hard", "easy"),
  modulation = c(1, 1, 1, 1)
)

# Build design matrices
design_matrices <-  convolve_design_matrix(
  timeseries = ts,
  events = events,
  factors = list(difficulty = c("hard", "easy")),
  contrasts = list(difficulty = matrix(c(-1, 1)))
)

Posterior Credible Interval Tests

Description

Modeled after t.test, returns the credible interval of the parameter or test and what proportion of the posterior distribution (or the difference in posterior distributions in case of a two sample test) overlaps with mu. For a one sample test provide x and for two sample also provide y. Note that for comparisons within one model, we recommend using hypothesis() if the priors were well chosen.

Usage

## S3 method for class 'emc'
credible(
  x,
  x_name = NULL,
  x_fun = NULL,
  x_fun_name = "fun",
  selection = "mu",
  y = NULL,
  y_name = NULL,
  y_fun = NULL,
  y_fun_name = "fun",
  x_subject = NULL,
  y_subject = NULL,
  mu = 0,
  alternative = c("less", "greater")[1],
  probs = c(0.025, 0.5, 0.975),
  digits = 2,
  p_digits = 3,
  print_table = TRUE,
  ...
)

credible(x, ...)

Arguments

Value

Invisible results table with no rounding.

Examples

{
# Run a credible interval test (Bayesian ''t-test'')
credible(samples_LNR, x_name = "m")
# We can also compare between two sets of emc objects

# # Now without a ~ E
# design_null <- design(data = forstmann,model=DDM,
#                            formula =list(v~0+S,a~1, t0~1, s~1, Z~1, sv~1, SZ~1),
#                            constants=c(s=log(1)))
#
# null_model <- make_emc(forstmann, design_null)
# null_model <- fit(null_model)
# credible(x = null_model, x_name = "a", y = full_model, y_name = "a")
#
# # Or provide custom functions:
# credible(x = full_model, x_fun = function(d) d["a_Eaccuracy"] - d["a_Eneutral"])
}

Posterior Quantiles

Description

Returns the quantiles of the selected parameter type. Full range of possible samples manipulations described in get_pars.

Usage

## S3 method for class 'emc.prior'
credint(
  x,
  selection = "mu",
  probs = c(0.025, 0.5, 0.975),
  digits = 3,
  N = 1000,
  covariates = NULL,
  ...
)

## S3 method for class 'emc'
credint(x, selection = "mu", probs = c(0.025, 0.5, 0.975), digits = 3, ...)

credint(x, ...)

Arguments

Value

A list of posterior quantiles for each parameter group in the selected parameter type.

Examples

credint(samples_LNR)

Cut Factors Based on Credible Loadings

Description

This function removes factors that do not have more than one credible loading based on the specified confidence interval.

Usage

cut_factors(emc, CI = 95)

Arguments

Value

An 'emc' object with factors that don't meet the credibility criterion removed

Specify a Design and Model

Description

This function combines information regarding the data, type of model, and the model specification.

Usage

design(
  formula = NULL,
  factors = NULL,
  Rlevels = NULL,
  model,
  data = NULL,
  contrasts = NULL,
  matchfun = NULL,
  constants = NULL,
  covariates = NULL,
  functions = NULL,
  report_p_vector = TRUE,
  custom_p_vector = NULL,
  transform = NULL,
  bound = NULL,
  ...
)

Arguments

Value

A design list.

Examples

# load example dataset
dat <- forstmann

# create a function that takes the latent response (lR) factor (d) and returns a logical
# defining the correct response for each stimulus. Here the match is simply
# such that the S factor equals the latent response factor
matchfun <- function(d)d$S==d$lR

# When working with lM and lR, it can be useful to design  an
# "average and difference" contrast matrix. For binary responses, it has a
# simple canonical form
ADmat <- matrix(c(-1/2,1/2),ncol=1,dimnames=list(NULL,"diff"))

# Create a design for a linear ballistic accumulator model (LBA) that allows
# thresholds to be a function of E and lR. The final result is a 9 parameter model.
design_LBABE <- design(data = dat,model=LBA,matchfun=matchfun,
                            formula=list(v~lM,sv~lM,B~E+lR,A~1,t0~1),
                            contrasts=list(v=list(lM=ADmat)),
                            constants=c(sv=log(1)))

Create fMRI Design for EMC2 Sampling

Description

This function takes the output from convolve_design_matrix and transforms it into a design suitable for sampling with EMC2. It properly configures parameter types, bounds, and transformations for the specified model.

Usage

design_fmri(design_matrix, model = MRI_AR1, ...)

Arguments

Value

An object of class 'emc.design' suitable for EMC2 sampling

Examples

# Generate a simple example timeseries
ts <- data.frame(
 subjects = rep(1, 100),
 run = rep(1, 100),
 time = cumsum(rep(1.38, 100)),
 ROI1 = rnorm(100)
)

# Generate example events
events <- data.frame(
 subjects = rep(1, 4),
 run = rep(1, 4),
 onset = c(10, 30, 50, 70),
 duration = rep(0.5, 4),
 event_type = c("A", "B", "A", "B"),
 modulation = c(1, 1, 1, 1)
)

# Create convolved design matrix
design_matrix <- convolve_design_matrix(
 timeseries = ts,
 events = events,
 factors = list(condition = c("A", "B")),
 hrf_model = "glover"
)

# Create fMRI design for EMC2
fmri_design <- design_fmri(design_matrix, model = MRI_AR1)

Effective Sample Size

Description

Returns the effective sample size (ESS) of the selected parameter type. Full range of possible samples manipulations described in get_pars.

Usage

## S3 method for class 'emc'
ess_summary(
  emc,
  selection = "mu",
  stat = "min",
  stat_only = FALSE,
  digits = 1,
  ...
)

ess_summary(emc, ...)

Arguments

Value

A matrix or vector of ESS values for the selected parameter type.

Examples

ess_summary(samples_LNR, selection = "alpha")

Factor diagram plot #Makes a factor diagram plot. Heavily based on the fa.diagram function of the psych package.

Description

Factor diagram plot #Makes a factor diagram plot. Heavily based on the fa.diagram function of the psych package.

Usage

factor_diagram(
  emc = NULL,
  stage = "sample",
  loadings = NULL,
  standardize = TRUE,
  simple = FALSE,
  only_cred = TRUE,
  cut = 0,
  nice_names = NULL,
  factor_names = NULL,
  sort = TRUE,
  adj = 1,
  main = NULL,
  cex = NULL
)

Arguments

Model Estimation in EMC2

Description

General purpose function to estimate models specified in EMC2.

Usage

## S3 method for class 'emc'
fit(
  emc,
  stage = NULL,
  iter = 1000,
  stop_criteria = NULL,
  search_width = 1,
  step_size = 100,
  verbose = TRUE,
  fileName = NULL,
  particle_factor = 50,
  cores_per_chain = 1,
  cores_for_chains = length(emc),
  max_tries = 20,
  thin = FALSE,
  ...
)

fit(emc, ...)

Arguments

Details

stop_criteria is either a list of lists with names of the stages, or a single list in which case its assumed to be for the sample stage (see examples). The potential stop criteria to be set are:

selection (character vector): For which parameters the stop_criteria should hold

mean_gd (numeric): The mean Gelman-Rubin diagnostic across all parameters in the selection

max_gd (numeric): The max Gelman-Rubin diagnostic across all parameters in the selection

min_unique (integer): The minimum number of unique samples in the MCMC chains across all parameters in the selection

min_es (integer): The minimum number of effective samples across all parameters in the selection

omit_mpsrf (Boolean): Whether to include the multivariate point-scale reduction factor in the Gelman-Rubin diagnostic. Default is FALSE.

iter (integer): The number of MCMC samples to collect.

The estimation is performed using particle-metropolis within-Gibbs sampling. For sampling details see:

Gunawan, D., Hawkins, G. E., Tran, M.-N., Kohn, R., & Brown, S. (2020). New estimation approaches for the hierarchical linear ballistic accumulator model. Journal of Mathematical Psychology ,96, 102368. doi.org/10.1016/j.jmp.2020.102368

Stevenson, N., Donzallaz, M. C., Innes, R. J., Forstmann, B., Matzke, D., & Heathcote, A. (2024). EMC2: An R Package for cognitive models of choice. doi.org/10.31234/osf.io/2e4dq

Value

An emc object

Examples

# Define a design first
ADmat <- matrix(c(-1/2,1/2),ncol=1,dimnames=list(NULL,"d"))
# We also define a match function for lM
matchfun=function(d)d$S==d$lR

# Drop most subjects for speed
dat <- forstmann[forstmann$subjects %in% unique(forstmann$subjects)[1:2],]
dat$subjects <- droplevels(dat$subjects)

design_LNR <- design(data = dat,model=LNR,matchfun=matchfun,
                     formula=list(m~lM,s~1,t0~1),
                     contrasts=list(m=list(lM=ADmat)))
# Before fit can be called, we first need to make an emc object
LNR_s <- make_emc(dat, design_LNR, rt_resolution = 0.05, n_chains = 2)
# Run fit, here illustrating how to use stop_criteria (also for speed purposes)
LNR_s <- fit(LNR_s, cores_for_chains = 1, stop_criteria = list(
  preburn = list(iter = 10), burn = list(mean_gd = 2.5), adapt = list(min_unique = 20),
  sample = list(iter = 25, max_gd = 2)), verbose = FALSE, particle_factor = 30, step_size = 25)

Forstmann et al.'s Data

Description

A dataset containing the speed or accuracy manipulation for a Random Dot Motion experiment.

Usage

forstmann

Format

A data frame with 15818 rows and 5 variables:

E

Factor with 3 levels for Speed, Accuracy and Neutral

R

Factor with 2 levels for Left and Right responses

S

Factor with 2 levels for Left and Right trials

rt

reaction time for each trial as a double

subjects

integer ID for each subject

Details

Details on the dataset can be found in the following paper:

Striatum and pre-SMA facilitate decision-making under time pressure

Birte U. Forstmann, Gilles Dutilh, Scott Brown, Jane Neumann, D. Yves von Cramon, K. Richard Ridderinkhof, Eric-Jan Wagenmakers.

Proceedings of the National Academy of Sciences Nov 2008, 105 (45) 17538-17542; DOI: 10.1073/pnas.0805903105

Source

https://www.pnas.org/doi/10.1073/pnas.0805903105

Gelman-Rubin Statistic

Description

Returns the Gelman-Rubin diagnostics (otherwise known as the R-hat) of the selected parameter type; i.e. the ratio of between to within MCMC chain variance.

Usage

## S3 method for class 'emc'
gd_summary(
  emc,
  selection = "mu",
  omit_mpsrf = TRUE,
  stat = "max",
  stat_only = FALSE,
  digits = 3,
  ...
)

gd_summary(emc, ...)

Arguments

Details

See: Gelman, A and Rubin, DB (1992) Inference from iterative simulation using multiple sequences, Statistical Science, 7, 457-511.

Full range of possible samples manipulations described in get_pars.

Value

A matrix or vector of R-hat values for the selected parameter type.

Examples

gd_summary(samples_LNR, selection = "correlation", stat = "mean", flatten = TRUE)

Bayes Factors

Description

returns the Bayes Factor for two models

Usage

get_BayesFactor(MLL1, MLL2)

Arguments

Value

The BayesFactor for model 1 over model 2

Examples

# Normally one would compare two different models
# Here we use two times the same model:
M1 <- M0 <- run_bridge_sampling(samples_LNR, both_splits = FALSE, cores_for_props = 1)
get_BayesFactor(M1, M0)

Get Data

Description

Extracts data from an emc object

Usage

## S3 method for class 'emc'
get_data(emc)

get_data(emc)

Arguments

Details

emc adds columns and rows to a dataframe in order to facilitate efficient likelihood calculations. This function will return the data as provided originally.

Value

A dataframe of the original data

Examples

get_data(samples_LNR)

Get Design

Description

Extracts design from an emc object

Usage

## S3 method for class 'emc.prior'
get_design(x)

## S3 method for class 'emc'
get_design(x)

get_design(x)

Arguments

Value

A design with class emc.design

Examples

get_design(samples_LNR)

Filter/Manipulate Parameters from emc Object

Description

Underlying function used in most plotting and object handling functions in EMC2. Can for example be used to filter/thin a parameter type (i.e, group-level means mu) and convert to an mcmc.list.

Usage

get_pars(
  emc,
  selection = "mu",
  stage = get_last_stage(emc),
  thin = 1,
  filter = 0,
  map = FALSE,
  add_recalculated = FALSE,
  length.out = NULL,
  by_subject = FALSE,
  return_mcmc = TRUE,
  merge_chains = FALSE,
  subject = NULL,
  flatten = FALSE,
  remove_dup = FALSE,
  remove_constants = TRUE,
  use_par = NULL,
  type = NULL,
  true_pars = NULL,
  chain = NULL,
  covariates = NULL
)

Arguments

Value

An mcmc.list object of the selected parameter types with the specified manipulations

Examples

# E.g. get the group-level mean parameters mapped back to the design
get_pars(samples_LNR, stage = "sample", map = TRUE, selection = "mu")

# Or return the flattened correlation, with 10 iterations per chain
get_pars(samples_LNR, stage = "sample", selection = "correlation", flatten = TRUE, length.out = 10)

Get Prior

Description

Extracts prior from an emc object

Usage

## S3 method for class 'emc'
get_prior(emc)

get_prior(emc)

Arguments

Value

A prior with class emc.prior

Examples

get_prior(samples_LNR)

Get parameter types from trend object

Description

Get parameter types from trend object

Usage

get_trend_pnames(trend)

Arguments

Value

A character vector of parameter names used in the trend

Examples

trend <- make_trend(par_names = "v", cov_names = "trial", kernels = "exp_incr")
get_trend_pnames(trend)

Create Group-Level Design Matrices

Description

Creates design matrices for group-level parameters based on subject-level design and formulas. This function is used for hierarchical modeling to specify how subject-level parameters vary across groups or conditions.

Usage

group_design(formula, data, subject_design, contrasts = NULL)

Arguments

Details

Here it is important to consider the interpretation of the group-level mean. This allows one to add covariates/group-level factors to the model. However, mu, the group-level mean, is still included for all parameters. Mu represents the intercept in the design matrix, this intercept is always added to the group-level model. Therefore, to keep the interpretation of mu as the group-level mean, it is important to ensure that the design matrix has a mean of zero. If not, this function will throw a warning. For some unbalanced designs, this is unavoidable and the warning can be ignored.

Value

A list of design matrices, one for each parameter specified in the formula. The intercept is automatically included as the group-level mean and is omitted from the design matrices.

Examples

# Create subject-level design 
subj_design <- design(data = forstmann, model = DDM,
                      formula = list(v ~ S, a ~ E, t0 ~ 1),
                      contrasts = list(S = contr.helmert))
# Add some age covariate and roughly demeans
# Demeaning is important to ensure that the interpretation of the group-level intercept
# is the mean of the group (i.e., 'mu' still represents the group-level mean)
forstmann$age <- as.numeric(forstmann$subjects) -mean(as.numeric(forstmann$subjects))
# Create fake group column
forstmann$group <- ifelse(forstmann$subjects %in%
              unique(forstmann$subjects)[seq(1, 19, 2)], "A", "B")

# Create group-level design matrices
group_des <- group_design(
  formula = list(v_S1 ~ age + group, a ~ age),
  data = forstmann,
  subject_design = subj_design,
  contrasts = list(group = contr.bayes)
)
# Then you can make the emc object with
emc <- make_emc(forstmann, subj_design, compress = FALSE, group_design = group_des)

Apply High-Pass Filtering to fMRI Data

Description

This function applies high-pass filtering to fMRI data to remove low-frequency noise and drift. It supports two filtering methods: cosine basis functions and polynomial regressors.

Usage

high_pass_filter(X, high_pass_model = "cosine", frame_times = NULL, ...)

Arguments

Value

A data frame or matrix with the same structure as X, but with high-frequency components removed from the data columns.

Examples

# Create a simple example data frame with drift
set.seed(123)
n_frames <- 100
time <- seq(0, 99)

# Create a signal with low-frequency drift
drift <- 0.1 * time
signal <- sin(2 * pi * 0.1 * time) + drift
noise <- rnorm(n_frames, 0, 0.5)
data <- signal + noise

# Create a data frame
df <- data.frame(
  time = time,
  signal = data
)

# Apply high-pass filtering using cosine basis functions
filtered_df <- high_pass_filter(df, high_pass_model = "cosine")

Within-Model Hypothesis Testing

Description

Approximates the Bayes factor for parameter effects using the savage-dickey ratio.

Usage

## S3 method for class 'emc'
hypothesis(
  emc,
  parameter = NULL,
  H0 = 0,
  fun = NULL,
  selection = "mu",
  do_plot = TRUE,
  use_prior_lim = TRUE,
  N = 10000,
  prior_args = list(),
  ...
)

hypothesis(emc, ...)

Arguments

Details

Note this is different to the computation of the marginal deviance in compare since it only considers the group level effect and not the whole model (i.e. subject-level parameters). For details see: Wagenmakers, Lodewyckx, Kuriyal, & Grasman (2010).

Value

The Bayes factor for the hypothesis against H0.

Examples

# Here the emc object has an effect parameter (e.g. m),
# that maps onto a certain hypothesis.
# The hypothesis here is that m is different from zero.
# We can test whether there's a group-level effect on m:
hypothesis(samples_LNR, parameter = "m")
# Alternatively we can also test whether two parameters differ from each other
mdiff <- function(p)diff(p[c("m","m_lMd")])
hypothesis(samples_LNR,fun=mdiff)

Initialize Chains

Description

Adds a set of start points to each chain. These start points are sampled from a user-defined multivariate normal across subjects.

Usage

init_chains(
  emc,
  start_mu = NULL,
  start_var = NULL,
  particles = 1000,
  cores_per_chain = 1,
  cores_for_chains = length(emc)
)

Arguments

Value

An emc object

Examples

# Make a design and an emc object
design_DDMaE <- design(data = forstmann,model=DDM,
                           formula =list(v~0+S,a~E, t0~1, s~1),
                           constants=c(s=log(1)))

DDMaE <- make_emc(forstmann, design_DDMaE)
# set up our mean starting points (same used across subjects).
mu <- c(v_Sleft=-2,v_Sright=2,a=log(1),a_Eneutral=log(1.5),a_Eaccuracy=log(2),
       t0=log(.2))
# Small variances to simulate start points from a tight range
var <- diag(0.05, length(mu))
# Initialize chains, 4 cores per chain, and parallelizing across our 3 chains as well
# so 4*3 cores used.
DDMaE <- init_chains(DDMaE, start_mu = mu, start_var = var,
                     cores_per_chain = 1, cores_for_chains = 1, particles = 10)
# Afterwards we can just use fit
# DDMaE <- fit(DDMaE, cores_per_chain = 4)

Simulate Data

Description

Simulates data based on a model design and a parameter vector (p_vector) by one of two methods:

1. Creating a fully crossed and balanced design specified by the design, with number of trials per cell specified by the n_trials argument

2. Using the design of a data frame supplied, which allows creation of unbalanced and other irregular designs, and replacing previous data with simulated data

Usage

make_data(
  parameters,
  design = NULL,
  n_trials = NULL,
  data = NULL,
  expand = 1,
  staircase = NULL,
  functions = NULL,
  ...
)

Arguments

Details

To create data for multiple subjects see ?make_random_effects().

Value

A data frame with simulated data

Examples

# First create a design
design_DDMaE <- design(factors = list(S = c("left", "right"),
                                           E = c("SPD", "ACC"),
                                           subjects = 1:30),
                            Rlevels = c("left", "right"), model = DDM,
                            formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                            constants=c(s=log(1)))
# Then create a p_vector:
parameters <- c(v_Sleft=-2,v_Sright=2,a=log(1),a_EACC=log(2), t0=log(.2),
              Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5))

# Now we can simulate data
data <- make_data(parameters, design_DDMaE, n_trials = 30)

# We can also simulate data based on a specific dataset
design_DDMaE <- design(data = forstmann,model=DDM,
                            formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                            constants=c(s=log(1)))
parameters <- c(v_Sleft=-2,v_Sright=2,a=log(1),a_Eneutral=log(1.5),a_Eaccuracy=log(2),
              t0=log(.2),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5))

data <- make_data(parameters, design_DDMaE, data = forstmann)

Make an emc Object

Description

Creates an emc object by combining the data, prior, and model specification into a emc object that is needed in fit().

Usage

make_emc(
  data,
  design,
  model = NULL,
  type = "standard",
  n_chains = 3,
  compress = TRUE,
  rt_resolution = 0.02,
  prior_list = NULL,
  group_design = NULL,
  par_groups = NULL,
  ...
)

Arguments

Value

An uninitialized emc object

Examples

dat <- forstmann

# function that takes the lR factor (named diff in the following function) and
# returns a logical defining the correct response for each stimulus. In this
# case the match is simply such that the S factor equals the latent response factor.
matchfun <- function(d)d$S==d$lR

# design an "average and difference" contrast matrix
ADmat <- matrix(c(-1/2,1/2),ncol=1,dimnames=list(NULL,"diff"))

# specify design
design_LBABE <- design(data = dat,model=LBA,matchfun=matchfun,
formula=list(v~lM,sv~lM,B~E+lR,A~1,t0~1),
contrasts=list(v=list(lM=ADmat)),constants=c(sv=log(1)))

# specify priors
pmean <- c(v=1,v_lMdiff=1,sv_lMTRUE=log(.5), B=log(.5),B_Eneutral=log(1.5),
           B_Eaccuracy=log(2),B_lRright=0, A=log(0.25),t0=log(.2))
psd <- c(v=1,v_lMdiff=0.5,sv_lMTRUE=.5,
         B=0.3,B_Eneutral=0.3,B_Eaccuracy=0.3,B_lRright=0.3,A=0.4,t0=.5)
prior_LBABE <- prior(design_LBABE, type = 'standard',pmean=pmean,psd=psd)

# create emc object
LBABE <- make_emc(dat,design_LBABE,type="standard",  prior=prior_LBABE)

Generate Subject-Level Parameters

Description

Simulates subject-level parameters in the format required by make_data().

Usage

make_random_effects(
  design,
  group_means,
  n_subj = NULL,
  variance_proportion = 0.2,
  covariances = NULL
)

Arguments

Value

A matrix of subject-level parameters.

Examples

# First create a design
design_DDMaE <- design(data = forstmann,model=DDM,
                            formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                            constants=c(s=log(1)))
# Then create a group-level means vector:
group_means =c(v_Sleft=-2,v_Sright=2,a=log(1),a_Eneutral=log(1.5),a_Eaccuracy=log(2),
               t0=log(.2),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5))
# Now we can create subject-level parameters
subj_pars <- make_random_effects(design_DDMaE, group_means, n_subj = 19)

# We can also define a covariance matrix to simulate from
subj_pars <- make_random_effects(design_DDMaE, group_means, n_subj = 19,
             covariances = diag(.1, length(group_means)))

# The subject level parameters can be used to generate data
make_data(subj_pars, design_DDMaE, n_trials = 10)

Create a trend specification for model parameters

Description

Create a trend specification for model parameters

Usage

make_trend(
  par_names,
  cov_names,
  kernels,
  bases = NULL,
  shared = NULL,
  trend_pnames = NULL,
  premap = TRUE,
  pretransform = FALSE
)

Arguments

Value

A list containing the trend specifications for each parameter

Examples

# Put trend on B and v parameters
trend <- make_trend(
  par_names = c("B", "v"),
  cov_names = "strial",
  kernels = c("exp_incr", "poly3"),
  premap = TRUE,
  shared = list(shrd = list("B.B0", "v.d1"))
)
get_trend_pnames(trend)

Parameter Mapping Back to the Design Factors

Description

Maps parameters of the cognitive model back to the experimental design. If p_vector is left unspecified will print a textual description of the mapping. Otherwise the p_vector can be created using sampled_pars(). The returned matrix shows whether/how parameters differ across the experimental factors.

Usage

mapped_pars(
  x,
  p_vector = NULL,
  model = NULL,
  digits = 3,
  remove_subjects = TRUE,
  covariates = NULL,
  ...
)

## S3 method for class 'emc.design'
mapped_pars(
  x,
  p_vector = NULL,
  model = NULL,
  digits = 3,
  remove_subjects = TRUE,
  covariates = NULL,
  ...
)

## S3 method for class 'emc.prior'
mapped_pars(
  x,
  p_vector = NULL,
  model = NULL,
  digits = 3,
  remove_subjects = TRUE,
  covariates = NULL,
  ...
)

## S3 method for class 'emc'
mapped_pars(
  x,
  p_vector = NULL,
  model = NULL,
  digits = 3,
  remove_subjects = TRUE,
  covariates = NULL,
  ...
)

Arguments

Value

Matrix with a column for each factor in the design and for each model parameter type (p_type).

Examples

# First define a design:
design_DDMaE <- design(data = forstmann,model=DDM,
                           formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                           constants=c(s=log(1)))
mapped_pars(design_DDMaE)
# Then create a p_vector:
p_vector=c(v_Sleft=-2,v_Sright=2,a=log(1),a_Eneutral=log(1.5),a_Eaccuracy=log(2),
          t0=log(.2),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5))
# This will map the parameters of the p_vector back to the design
mapped_pars(design_DDMaE, p_vector)

Merge Samples

Description

Merges samples from all chains as one unlisted object.

Usage

merge_chains(emc)

Arguments

Details

Note that all sampling stages are included in the merged output, including iterations from the preburn, burn, and adapt stages. merge_chains(emc)$samples$stage shows the corresponding sampling stages.

Value

An unlisted emc object with all chains merged

Model Averaging

Description

Computes model weights and a Bayes factor by comparing two groups of models based on their Information Criterion (IC) values. The function works with either numeric vectors or data frames containing multiple IC measures (e.g., MD, BPIC, DIC).

Usage

model_averaging(IC_for, IC_against)

Arguments

Details

When provided with numeric vectors, it computes the weights for the two groups by first converting the IC values into relative weights and then normalizing them. When provided with a data frame, it assumes that the data frame is the output of a call to compare and applies averaging to each IC metric

Value

A data.frame with the following columns:

wFor

The aggregated weight of the models in favor.

wAgainst

The aggregated weight of the models against.

Factor

The Bayes factor (ratio of wFor to wAgainst).

If IC_for is a data frame, a matrix with rows corresponding to each IC measure is returned.

Examples

# First set up some example models (normally these would be alternative models)
samples_LNR2 <- subset(samples_LNR, length.out = 45)
samples_LNR3 <- subset(samples_LNR, length.out = 40)
samples_LNR4 <- subset(samples_LNR, length.out = 35)

# Run compare on them, BayesFactor = F is set for speed.
ICs <- compare(list(S1 = samples_LNR, S2 = samples_LNR2,
                    S3 = samples_LNR3, S4 = samples_LNR4), BayesFactor = FALSE)

# Model averaging can either be done with a vector of ICs:
model_averaging(ICs$BPIC[1:2], ICs$BPIC[2:4])

# Or the output of compare:
model_averaging(ICs[1:2,], ICs[3:4,])

Plot Within-Chain Correlations

Description

Plots within-chain parameter correlations (upper triangle) and corresponding scatterplots (lower triangle) to visualize parameter sloppiness.

Usage

pairs_posterior(
  emc,
  selection = "alpha",
  scale_subjects = TRUE,
  do_plot = TRUE,
  N = 500,
  ...
)

Arguments

Details

If selection = alpha the parameter chains are concatenated across participants, (after standardizing if scale_subjects = TRUE) and then correlated.

Value

Invisibly returns a matrix with the correlations between the parameters.

Examples

# Plot the sloppiness for the individual-level subjects
pairs_posterior(samples_LNR, selection = "alpha")

# We can also choose group-level parameters and subsets of the parameter space
pairs_posterior(samples_LNR, use_par = c("m", "t0"), selection = "sigma2")

Return Data Frame of Parameters

Description

Return Data Frame of Parameters

Usage

## S3 method for class 'emc.prior'
parameters(x, selection = "mu", N = 1000, covariates = NULL, ...)

## S3 method for class 'emc'
parameters(x, selection = "mu", N = NULL, resample = FALSE, ...)

parameters(x, ...)

Arguments

Value

A data frame with one row for each sample (with a subjects column if selection = "alpha" and using draws from the posterior)

Examples

# For prior inference:
# First set up a prior
design_DDMaE <- design(data = forstmann,model=DDM,
                       formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                       constants=c(s=log(1)))
# Then set up a prior using make_prior
p_vector=c(v_Sleft=-2,v_Sright=2,a=log(1),a_Eneutral=log(1.5),a_Eaccuracy=log(2),
           t0=log(.2),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5))
psd <- c(v_Sleft=1,v_Sright=1,a=.3,a_Eneutral=.3,a_Eaccuracy=.3,
         t0=.4,Z=1,sv=.4,SZ=1)
# Here we left the variance prior at default
prior_DDMaE <- prior(design_DDMaE,mu_mean=p_vector,mu_sd=psd)
# Get our prior samples
parameters(prior_DDMaE, N = 100)
# For posterior inference:
# Get 100 samples of the group-level mean (the default)
parameters(samples_LNR, N = 100)
# or from the individual-level parameters and mapped
parameters(samples_LNR, selection = "alpha", map = TRUE)

Plot Function for emc Objects

Description

Makes trace plots for model parameters.

Usage

## S3 method for class 'emc'
plot(
  x,
  stage = "sample",
  selection = c("mu", "sigma2", "alpha"),
  layout = NA,
  ...
)

Arguments

Details

If an emc object that has not been run with fit yet is supplied prior plots will be returned.

Value

A trace/acf plot of the selected MCMC chains

Examples

plot(samples_LNR)
# Or trace autocorrelation for the second subject:
plot(samples_LNR, subject = 2, selection = "alpha")

# Can also plot the trace of for example the group-level correlation:
plot(samples_LNR, selection = "correlation", col = c("green", "purple", "orange"), lwd = 2)

Plot method for emc.design objects

Description

Makes design illustration by plotting simulated data based on the design

Usage

## S3 method for class 'emc.design'
plot(
  x,
  p_vector,
  data = NULL,
  factors = NULL,
  plot_factor = NULL,
  n_data_sim = 10,
  functions = NULL,
  ...
)

Arguments

Value

No return value, called for side effect of plotting

Plot a prior

Description

Takes a prior object and plots the selected implied prior

Usage

## S3 method for class 'emc.prior'
plot(
  x,
  selection = "mu",
  do_plot = TRUE,
  covariates = NULL,
  layout = NA,
  N = 50000,
  ...
)

Arguments

Value

An invisible mcmc.list object with prior samples of the selected type

Examples

# First define a design for the model
design_DDMaE <- design(data = forstmann,model=DDM,
                           formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                           constants=c(s=log(1)))
# Then set up a prior using make_prior
p_vector=c(v_Sleft=-2,v_Sright=2,a=log(1),a_Eneutral=log(1.5),a_Eaccuracy=log(2),
          t0=log(.2),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5))
psd <- c(v_Sleft=1,v_Sright=1,a=.3,a_Eneutral=.3,a_Eaccuracy=.3,
          t0=.4,Z=1,sv=.4,SZ=1)
# Here we left the variance prior at default
prior_DDMaE <- prior(design_DDMaE,mu_mean=p_vector,mu_sd=psd)
# Now we can plot all sorts of (implied) priors
plot(prior_DDMaE, selection = "mu", N = 1e3)
plot(prior_DDMaE, selection = "mu", mapped = FALSE, N=1e3)
# We can also plot the implied prior on the participant level effects.
plot(prior_DDMaE, selection = "alpha", col = "green", N = 1e3)

Plot Defective Cumulative Distribution Functions

Description

Plots panels of cumulative distribution functions (CDFs) for each level of the specified defective factor in the data. The CDFs are defective; each factor level's CDF scales only up to that level's proportion. Summed across levels, the maximum is 1. Optionally, posterior and/or prior predictive CDFs can be overlaid.

Usage

plot_cdf(
  input,
  post_predict = NULL,
  prior_predict = NULL,
  subject = NULL,
  quants = c(0.025, 0.975),
  functions = NULL,
  factors = NULL,
  defective_factor = "R",
  n_cores = 1,
  n_post = 50,
  layout = NA,
  to_plot = c("data", "posterior", "prior")[1:2],
  use_lim = c("data", "posterior", "prior")[1:2],
  legendpos = c("top", "topright"),
  posterior_args = list(),
  prior_args = list(),
  ...
)

Arguments

Value

Returns NULL invisibly.

Examples

# Plot defective CDF for data only
# plot_cdf(forstmann, to_plot = "data")
#
# Plot with posterior predictions
# plot_cdf(samples_LNR, to_plot = c("data","posterior"), n_post=10)
#
# Or a list of multiple emc objects ...

Plot Defective Densities

Description

Plots panels that contain a set of densities for each level of the specified defective factor in the data. These densities are defective; their areas are relative to the respective proportions of the defective factor levels. Across all levels, the area sums to 1. Optionally, posterior/prior predictive densities can be overlaid.

Usage

plot_density(
  input,
  post_predict = NULL,
  prior_predict = NULL,
  subject = NULL,
  quants = c(0.025, 0.975),
  functions = NULL,
  factors = NULL,
  defective_factor = "R",
  n_cores = 1,
  n_post = 50,
  layout = NA,
  to_plot = c("data", "posterior", "prior")[1:2],
  use_lim = c("data", "posterior", "prior")[1:2],
  legendpos = c("topright", "top"),
  posterior_args = list(),
  prior_args = list(),
  ...
)

Arguments

Examples

# Plot defective densities for each subject and the factor combination in the design:
plot_density(forstmann)
# or for one subject:
plot_density(forstmann, subject = 1)
# Now collapsing across subjects and using a different defective factor:
plot_density(forstmann, factors = "S", defective_factor = "E")
# Or plot posterior predictives
plot_density(samples_LNR, n_post = 10)

Plot Design

Description

Makes design illustration by plotting simulated data based on the design

Usage

## S3 method for class 'emc.design'
plot_design(
  x,
  data = NULL,
  factors = NULL,
  plot_factor = NULL,
  n_data_sim = 10,
  p_vector = NULL,
  functions = NULL,
  ...
)

## S3 method for class 'emc.prior'
plot_design(
  x,
  data = NULL,
  factors = NULL,
  plot_factor = NULL,
  n_data_sim = 10,
  p_vector = NULL,
  functions = NULL,
  ...
)

plot_design(
  x,
  data = NULL,
  factors = NULL,
  plot_factor = NULL,
  n_data_sim = 10,
  p_vector = NULL,
  functions = NULL,
  ...
)

## S3 method for class 'emc'
plot_design(
  x,
  data = NULL,
  factors = NULL,
  plot_factor = NULL,
  n_data_sim = 10,
  p_vector = NULL,
  functions = NULL,
  ...
)

Arguments

Value

No return value. Just plots the design

Plot fMRI Design Matrix

Description

This function creates a visualization of an fMRI design matrix, showing the temporal evolution of regressors over time. It can handle various input formats and provides options to customize the visualization.

Usage

plot_design_fmri(
  design_matrix,
  TRs = 100,
  events = NULL,
  remove_nuisance = TRUE,
  subject = 1,
  legend_pos = "bottomleft",
  ...
)

Arguments

Value

A plot showing the design matrix regressors over time.

Examples

# Example time series
ts <- data.frame(
  subjects = rep(1, 100),
  run      = rep(1, 100),
  time     = seq(0, 99),
  ROI      = rnorm(100)
)
# Create a simple events data frame
events <- data.frame(
  subjects = rep(1, 10),
  run = rep(1, 10),
  onset = seq(0, 90, by = 10),
  condition = rep(c("A", "B"), 5),
  rt = runif(10, 0.5, 1.5),
  accuracy = sample(0:1, 10, replace = TRUE)
)
# Reshape with custom duration for each event_type
reshaped <- reshape_events(events,
                           event_types = c("condition", "accuracy", "rt"),
                           duration = list(condition = 0.5,
                                          accuracy = 0.2,
                                          rt = function(x) x$rt))
design_matrices <-  convolve_design_matrix(
                    timeseries = ts,
                    events = reshaped,
                    covariates = c('accuracy', 'rt'),
                    factors = list(cond = c("condition_A", "condition_B")),
                    contrasts = list(cond = matrix(c(-1, 1))))

# Plot the design matrix
plot_design_fmri(design_matrices)

Plot fMRI peri-stimulus time courses

Description

This function plots average BOLD response around specified events for a single ROI by using FIR based event estimation, all event_types in events are taken into account in the FIR. Posterior predictives can be overlaid via the post_predict argument.

Usage

plot_fmri(
  timeseries,
  post_predict = NULL,
  events,
  event_type,
  high_pass = TRUE,
  high_pass_model = "cosine",
  posterior_args = list(),
  legend_pos = "topleft",
  layout = NA,
  n_cores = 1,
  ...
)

Arguments

Value

NULL. Produces plots as a side-effect.

Examples

ts <- data.frame(
  subjects = rep(1, 100),
  run      = rep(1, 100),
  time     = seq(0, 99),
  ROI      = rnorm(100)
)
events <- data.frame(
  subjects   = rep(1, 5),
  run        = rep(1, 5),
  onset      = c(10, 30, 50, 70, 90),
  event_type = rep("A", 5),
  modulation = rep(1, 5),
  duration   = rep(0.5, 5)
)
plot_fmri(ts, events = events, event_type = "A")

Plots Density for Parameters

Description

Plots the posterior and prior density for selected parameters of a model. Full range of samples manipulations described in get_pars.

Usage

plot_pars(
  emc,
  layout = NA,
  selection = "mu",
  show_chains = FALSE,
  plot_prior = TRUE,
  N = 10000,
  use_prior_lim = !all_subjects,
  lpos = "topright",
  true_pars = NULL,
  all_subjects = FALSE,
  prior_args = list(),
  true_args = list(),
  ...
)

Arguments

Value

An invisible return of the contraction statistics for the selected parameter type

Examples

# Full range of possibilities described in get_pars
plot_pars(samples_LNR)
# Or plot all subjects
plot_pars(samples_LNR, all_subjects = TRUE, col = 'purple')
# Or plot recovery
true_emc <- samples_LNR # This would normally be the data-generating samples
plot_pars(samples_LNR, true_pars = true_emc, true_args = list(col = 'blue'), adjust = 2)

Plot Group-Level Relations

Description

An adjusted version of the corrplot package function corrplot() tailored to EMC2 and the plotting of estimated correlations.

Usage

plot_relations(
  emc = NULL,
  stage = "sample",
  plot_cred = FALSE,
  plot_means = TRUE,
  only_cred = TRUE,
  nice_names = NULL,
  selection = "correlation",
  use_par = NULL,
  ...
)

Arguments

Value

No return value, creates a plot of group-level relations

Examples

# For a given set of hierarchical model samples we can make a
# correlation matrix plot.
plot_relations(samples_LNR, only_cred = TRUE, plot_cred = TRUE)
# We can also only plot the correlations where the credible interval does not include zero
plot_relations(samples_LNR, plot_means = TRUE, only_cred = TRUE)

Plot the ECDF Difference in SBC Ranks

Description

Plots the difference in observed cumulative rank statistics and the expected cumulative distribution of a uniform distribution. The blue shaded areas indicate the 95% credible interval.

Usage

plot_sbc_ecdf(ranks, layout = NA)

Arguments

Value

No returns

Plot the Histogram of the Observed Rank Statistics of SBC

Description

Note that this plot is dependent on the number of bins, and a more general visualization is to use plot_sbc_ecdf

Usage

plot_sbc_hist(ranks, bins = 10, layout = NA)

Arguments

Value

No returns

Plot Statistics on Data

Description

Plots panels that contain a set of densities for each level of the specified factor The densities represent the predicted data across the posterior, the vertical lines represent the real data.

Usage

plot_stat(
  input,
  post_predict = NULL,
  prior_predict = NULL,
  stat_fun,
  stat_name = NULL,
  subject = NULL,
  factors = NULL,
  n_cores = 1,
  n_post = 50,
  quants = c(0.025, 0.5, 0.975),
  functions = NULL,
  layout = NA,
  to_plot = c("data", "posterior", "prior")[1:2],
  use_lim = c("data", "posterior", "prior")[1:2],
  legendpos = c("topleft", "top"),
  posterior_args = list(),
  prior_args = list(),
  ...
)

Arguments

Value

an invisible data frame with the stat applied to the real data, posterior predictives and/or prior predictives

Examples

# For example plot the observed and predicted response accuracy
# Can also apply more sophisticated statistics
drt <- function(data) diff(tapply(data$rt,data[,c("E")],mean))
plot_stat(samples_LNR, stat_fun = drt, n_post = 10, stat_name = "RT diff Speed - A/N")

Generate Posterior/Prior Predictives

Description

Simulate n_post data sets using the posterior/prior parameter estimates

Usage

## S3 method for class 'emc.prior'
predict(object, data = NULL, n_post = 50, n_cores = 1, n_trials = NULL, ...)

## S3 method for class 'emc'
predict(
  object,
  hyper = FALSE,
  n_post = 50,
  n_cores = 1,
  stat = c("random", "mean", "median")[1],
  ...
)

Arguments

Value

A list of simulated data sets of length n_post

Examples

# based on an emc object ran by fit() we can generate posterior predictives
predict(samples_LNR, n_cores = 1, n_post = 10)

Specify Priors for the Chosen Model

Description

These values are entered manually by default but can be recycled from another prior (given in the update argument).

Usage

prior(
  design,
  type = NULL,
  group_design = NULL,
  update = NULL,
  do_ask = NULL,
  fill_default = TRUE,
  ...
)

Arguments

Details

Where a value is not supplied, the user is prompted to enter numeric values (or functions that evaluate to numbers).

To get the prior help use prior_help(type). With type e.g. 'diagonal'.

Value

A prior list object

Examples

# First define a design for the model
design_DDMaE <- design(data = forstmann,model=DDM,
                           formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                           constants=c(s=log(1)))
# Then set up a prior using prior
p_vector=c(v_Sleft=-2,v_Sright=2,a=log(1),a_Eneutral=log(1.5),a_Eaccuracy=log(2),
                     t0=log(.2),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5))
psd <- c(v_Sleft=1,v_Sright=1,a=.3,a_Eneutral=.3,a_Eaccuracy=.3,
                     t0=.4,Z=1,sv=.4,SZ=1)
# Here we left the variance prior at default
prior_DDMaE <- prior(design_DDMaE,mu_mean=p_vector,mu_sd=psd)
# Also add a group-level variance prior:
pscale <- c(v_Sleft=.6,v_Sright=.6,a=.3,a_Eneutral=.3,a_Eaccuracy=.3,
                             t0=.2,Z=.5,sv=.4,SZ=.3)
df <- .4
prior_DDMaE <- prior(design_DDMaE,mu_mean=p_vector,mu_sd=psd, A = pscale, df = df)
# If we specify a new design
design_DDMat0E <- design(data = forstmann,model=DDM,
                           formula =list(v~0+S,a~E, t0~E, s~1, Z~1, sv~1, SZ~1),
                           constants=c(s=log(1)))
# We can easily update the prior
prior_DDMat0E <- prior(design_DDMat0E, update = prior_DDMaE)

Prior Specification Information

Description

Prints information associated with the prior for certain 'type'

Usage

prior_help(type)

Arguments

Value

Invisible return with a list of all the information that is also printed

Examples

prior_help('diagonal')

Likelihood Profile Plots

Description

Creates likelihood profile plots from a design and the experimental data by varying one model parameter while holding all others constant.

Usage

profile_plot(
  data,
  design,
  p_vector,
  range = 0.5,
  layout = NA,
  p_min = NULL,
  p_max = NULL,
  use_par = NULL,
  n_point = 100,
  n_cores = 1,
  round = 3,
  true_args = list(),
  ...
)

Arguments

Value

Vector with highest likelihood point, input and mismatch between true and highest point

Examples

# First create a design
design_DDMaE <- design(data = forstmann,model=DDM,
                      formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                      constants=c(s=log(1)))
# Then create a p_vector:
p_vector=c(v_Sleft=-2,v_Sright=2,a=log(.95),a_Eneutral=log(1.5),a_Eaccuracy=log(2),
          t0=log(.25),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5))
# Make a profile plot for some parameters. Specifying a custom range for t0.
profile_plot(p_vector = p_vector, p_min = c(t0 = -1.35),
             p_max = c(t0 = -1.45), use_par = c("a", "t0", "SZ"),
             data = forstmann, design = design_DDMaE, n_point = 10)

Recovery Plots

Description

Plots recovery of data generating parameters/samples. Full range of samples manipulations described in get_pars

Usage

## S3 method for class 'emc'
recovery(
  emc,
  true_pars,
  selection = "mu",
  layout = NA,
  do_CI = TRUE,
  correlation = "pearson",
  stat = "rmse",
  digits = 3,
  CI = 0.95,
  ci_plot_args = list(),
  ...
)

recovery(emc, ...)

Arguments

Value

Invisible list with RMSE, coverage, and Pearson and Spearman correlations.

Examples

# Make up some values that resemble posterior samples
# Normally this would be true values that were used to simulate the data
# Make up some values that resemble posterior samples
# Normally this would be true values that were used to simulate the data
pmat <- matrix(rnorm(12, mean = c(-1, -.6, -.4, -1.5), sd = .01), ncol = 4, byrow = TRUE)
# Conventionally this would be created before one makes data with true values
recovery(samples_LNR, pmat, correlation = "pearson", stat = "rmse", selection = "alpha")
# Similarly we can plot recovery of other parameters with a set of true samples
true_samples <- samples_LNR # Normally this would be data-generating samples
recovery(samples_LNR, true_samples, correlation = "pearson", stat = "rmse",
         selection = "correlation", cex = 1.5,
         ci_plot_args = list(lty = 3, length = .2, lwd = 2, col = "brown"))

Reshape events data for fMRI analysis

Description

This function reshapes event data into a format suitable for fMRI analysis by converting specified event_types into separate event types with appropriate modulation values.

Usage

reshape_events(events, event_types, duration = 0.001, modulation = NULL)

Arguments

Value

A data frame with columns 'subjects', 'onset', 'run', 'modulation', 'duration', and 'event_type'

Examples

# Create a simple events data frame
events <- data.frame(
  subjects = rep(1, 10),
  run = rep(1, 10),
  onset = seq(0, 90, by = 10),
  condition = rep(c("A", "B"), 5),
  rt = runif(10, 0.5, 1.5),
  accuracy = sample(0:1, 10, replace = TRUE)
)

# Reshape with default duration
reshaped1 <- reshape_events(events, event_types = c("condition", "accuracy"))

# Reshape with custom duration for each event_type
reshaped2 <- reshape_events(events,
                           event_types = c("condition", "accuracy", "rt"),
                           duration = list(condition = 0.5,
                                          accuracy = 0.2,
                                          rt = function(x) x$rt))

Estimating Marginal Likelihoods Using WARP-III Bridge Sampling

Description

Uses bridge sampling that matches a proposal distribution to the first three moments of the posterior distribution to get an accurate estimate of the marginal likelihood. The marginal likelihood can be used for computing Bayes factors and posterior model probabilities.

Usage

run_bridge_sampling(
  emc,
  stage = "sample",
  filter = NULL,
  repetitions = 1,
  cores_for_props = 4,
  cores_per_prop = 1,
  both_splits = TRUE,
  ...
)

Arguments

Details

If not enough posterior samples were collected using fit(), bridge sampling can be unstable. It is recommended to run run_bridge_sampling() several times with the repetitions argument and to examine how stable the results are.

It can be difficult to converge bridge sampling for exceptionally large models, because of a large number of subjects (> 100) and/or cognitive model parameters.

For a practical introduction:

Gronau, Q. F., Heathcote, A., & Matzke, D. (2020). Computing Bayes factors for evidence-accumulation models using Warp-III bridge sampling. Behavior research methods, 52(2), 918-937. doi.org/10.3758/s13428-019-01290-6

For mathematical background:

Meng, X.-L., & Wong, W. H. (1996). Simulating ratios of normalizing constants via a simple identity: A theoretical exploration. Statistica Sinica, 6, 831-860. http://www3.stat.sinica.edu.tw/statistica/j6n4/j6n43/j6n43.htm

Meng, X.-L., & Schilling, S. (2002). Warp bridge sampling. Journal of Computational and Graphical Statistics, 11(3), 552-586. doi.org/10.1198/106186002457

Value

A vector of length repetitions which contains the marginal log likelihood estimates per repetition

Examples

# After `fit` has converged on a specific model
# We can take those samples and calculate the marginal log-likelihood for them
MLL <- run_bridge_sampling(samples_LNR, cores_for_props = 1, both_splits = FALSE)
# This will run on 2*4 cores (since 4 is the default for ``cores_for_props``)

Fine-Tuned Model Estimation

Description

Although typically users will rely on fit, this function can be used for more fine-tuned specification of estimation needs. The function will throw an error if a stage is skipped, the stages have to be run in order ("preburn", "burn", "adapt", "sample"). More details can be found in the fit help files (?fit).

Usage

run_emc(
  emc,
  stage,
  stop_criteria,
  search_width = 1,
  step_size = 100,
  verbose = FALSE,
  verboseProgress = FALSE,
  fileName = NULL,
  particle_factor = 50,
  cores_per_chain = 1,
  cores_for_chains = length(emc),
  max_tries = 20,
  n_blocks = 1,
  thin = FALSE,
  trim = TRUE
)

Arguments

Value

An emc object

Examples

# First define a design
design_in <- design(data = forstmann,model=DDM,
                           formula =list(v~0+S,a~E, t0~1, s~1, Z~1),
                           constants=c(s=log(1)))
# Then make the emc, we've omitted a prior here for brevity so default priors will be used.
emc <- make_emc(forstmann, design_in)

# Now for example we can specify that we only want to run the "preburn" phase
# for MCMC 10 iterations
emc <- run_emc(emc, stage = "preburn", stop_criteria = list(iter = 10), cores_for_chains = 1)

Simulation-Based Calibration

Description

Runs SBC for an EMC2 model and associated design. Returns normalized rank (between 0 and 1) and prior samples. For hierarchical models the group-level mean and the (implied) group-level (co-)variance are returned. For non-hierarchical models only the subject-level parameters rank is returned.

Usage

run_sbc(
  design_in,
  prior_in,
  replicates = 250,
  trials = 100,
  n_subjects = 30,
  plot_data = FALSE,
  verbose = TRUE,
  fileName = NULL,
  ...
)

Arguments

Value

The ranks and prior samples. For hierarchical models also the prior-generated subject-level parameters.

Get Model Parameters from a Design

Description

Makes a vector with zeroes, with names and length corresponding to the model parameters of the design.

Usage

sampled_pars(
  x,
  group_design = NULL,
  doMap = FALSE,
  add_da = FALSE,
  all_cells_dm = FALSE,
  data = NULL
)

## S3 method for class 'emc.design'
sampled_pars(
  x,
  group_design = NULL,
  doMap = FALSE,
  add_da = FALSE,
  all_cells_dm = FALSE,
  data = NULL
)

## S3 method for class 'emc.group_design'
sampled_pars(
  x,
  group_design = NULL,
  doMap = FALSE,
  add_da = FALSE,
  all_cells_dm = FALSE,
  data = NULL
)

## S3 method for class 'emc.prior'
sampled_pars(
  x,
  group_design = NULL,
  doMap = FALSE,
  add_da = FALSE,
  all_cells_dm = FALSE,
  data = NULL
)

## S3 method for class 'emc'
sampled_pars(
  x,
  group_design = NULL,
  doMap = FALSE,
  add_da = FALSE,
  all_cells_dm = FALSE,
  data = NULL
)

Arguments

Value

Named vector.

Examples

# First define a design
design_DDMaE <- design(data = forstmann,model=DDM,
                           formula =list(v~0+S,a~E, t0~1, s~1, Z~1, sv~1, SZ~1),
                           constants=c(s=log(1)))
# Then for this design get which cognitive model parameters are sampled:
sampled_pars(design_DDMaE)

LNR Model of Forstmann Data (First 3 Subjects)

Description

An emc object with a limited number of samples and subjects of the Forstmann dataset. The object is a nested list of lenght three, each list containing the MCMC samples of the respective chain. The MCMC samples are stored in the samples element.

Usage

samples_LNR

Format

An emc object. An emc object is a list with a specific structure and elements, as outlined below.

data

A list of dataframes, one for each subject included

par_names

A character vector containing the model parameter names

n_pars

The number of parameters in the model

n_subjects

The number of unique subject ID's in the data

model

A list containing the model functions

nuisance

A logical vector indicating which parameters are nuisance parameters

subjects

A vector containing the unique subject ID's

type

The type of model e.g., "standard" or "diagonal"

prior

A list that holds the prior for theta_mu (the model parameters). Contains the mean (theta_mu_mean), covariance matrix (theta_mu_var), degrees of freedom (v), and scale (A) and inverse covariance matrix (theta_mu_invar)

samples

A list with defined structure containing the samples, see the Samples Element section for more detail

sampler_nuis

A sampler list for nuisance parameters (in this case there are none), similarly structured to the overall samples list of one of the MCMC chains.

Samples Element

The samples element of a emc object contains the different types of samples estimated by EMC2. These include the three main types of samples theta_mu, theta_var and alpha as well as a number of other items which are detailed here.

theta_mu

samples used for estimating the model parameters (group level), an array of size (n_pars x n_samples)

theta_var

samples used for estimating the parameter covariance matrix, an array of size (n_pars x n_pars x n_samples)

alpha

samples used for estimating the subject random effects, an array of size (n_pars x n_subjects x n_samples)

stage

A vector containing what PMwG stage each sample was drawn in

subj_ll

The winning particles log-likelihood for each subject and sample

a_half

Mixing weights used during the Gibbs step when creating a new sample for the covariance matrix

last_theta_var_inv

The inverse of the last samples covariance matrix

idx

The index of the last sample drawn

Source

https://www.pnas.org/doi/10.1073/pnas.0805903105

Split fMRI Timeseries Data by ROI Columns

Description

This function splits a timeseries data frame containing multiple ROI columns into a list of data frames, where each data frame contains the common columns (subjects, run, time) and one ROI column.

Usage

split_timeseries(timeseries, columns = NULL)

Arguments

Value

A named list of data frames, where each data frame contains the common columns (subjects, run, time) and one ROI column. The names of the list elements correspond to the ROI column names.

Examples

# Create a simple example timeseries with multiple ROIs
set.seed(123)
n_frames <- 100

# Create a data frame with multiple ROIs
timeseries <- data.frame(
  subjects = rep(1, n_frames),
  run = rep(1, n_frames),
  time = seq(0, n_frames-1),
  ROI1 = rnorm(n_frames),
  ROI2 = rnorm(n_frames),
  ROI3 = rnorm(n_frames)
)

# Split the timeseries by all ROI columns
split_data <- split_timeseries(timeseries)

Shorten an emc Object

Description

Shorten an emc Object

Usage

## S3 method for class 'emc'
subset(
  x,
  stage = "sample",
  filter = NULL,
  thin = 1,
  keep_stages = FALSE,
  length.out = NULL,
  ...
)

Arguments

Value

A shortened emc object

Examples

subset(samples_LNR, length.out = 10)

Summary Statistics for emc Objects

Description

Computes quantiles, Rhat and ESS for selected model parameters.

Usage

## S3 method for class 'emc'
summary(
  object,
  selection = c("mu", "sigma2", "alpha"),
  probs = c(0.025, 0.5, 0.975),
  digits = 3,
  ...
)

Arguments

Details

Note that if selection = alpha and by_subject = TRUE (default) is used, summary statistics are computed at the individual level. to the console but summary statistics for all subjects are returned by the function.

If an emc object that has not been run with fit yet is supplied, summary of the design will be returned.

Value

A list of summary output.

Summary method for emc.design objects

Description

Prints a summary of the design object, including sampled parameters and design matrices. For continuous covariates just prints one row, instead of all covariates.

Usage

## S3 method for class 'emc.design'
summary(object, ...)

Arguments

Value

Invisibly returns the design matrices

Summary method for emc.prior objects

Description

Prints a summary of the prior specification, including descriptions of the prior types and their associated hyperparameters.

Usage

## S3 method for class 'emc.prior'
summary(object, ...)

Arguments

Value

Invisibly returns NULL. Called for its side effect of printing the summary.

See Also

prior for creating prior objects

Examples

# Take a prior object
prior <- get_prior(samples_LNR)
summary(prior)

Get help information for trend kernels and bases

Description

Get help information for trend kernels and bases

Usage

trend_help(kernel = NULL, base = NULL, ...)

Arguments

Value

Formatted trend information

Examples

# Get information about exponential increasing kernel
trend_help(kernel = "exp_incr")

# Get information about linear base
trend_help(base = "lin")

# Return available kernel and base types
trend_help()

Update EMC Objects to the Current Version

Description

This function updates EMC objects created with older versions of the package to be compatible with the current version.

Usage

update2version(emc)

Arguments

Value

An updated EMC object compatible with the current version

Examples

# Update the model to current version
updated_model <- update2version(samples_LNR)