Title: | 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>. |
Authors: | 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 <[email protected]> |
License: | GPL (>= 3) |
Version: | 2.1.0 |
Built: | 2024-11-14 04:30:10 UTC |
Source: | https://github.com/ampl-psych/emc2 |
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.
chain_n(emc)
chain_n(emc)
emc |
A list, the output of |
A matrix
chain_n(samples_LNR)
chain_n(samples_LNR)
Runs a series of convergence checks, prints statistics to the console, and makes traceplots of the worst converged parameter per selection.
## S3 method for class 'emc' check( emc, selection = c("mu", "sigma2", "alpha"), digits = 3, plot_worst = TRUE, ... ) check(emc, ...)
## S3 method for class 'emc' check( emc, selection = c("mu", "sigma2", "alpha"), digits = 3, plot_worst = TRUE, ... ) check(emc, ...)
emc |
An emc object |
selection |
A Character vector. Indicates which parameter types to check (e.g., |
digits |
Integer. How many digits to round the ESS and Rhat to in the plots |
plot_worst |
Boolean. If |
... |
Optional arguments that can be passed to |
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
.
a list with the statistics for the worst converged parameter per selection
check(samples_LNR)
check(samples_LNR)
Returns the BPIC/DIC or marginal deviance (-2*marginal likelihood) for a list of samples objects.
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, ... )
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, ... )
sList |
List of samples objects |
stage |
A string. Specifies which stage the samples are to be taken from |
filter |
An integer or vector. If it's an integer, iterations up until the value set by |
use_best_fit |
Boolean, defaults to |
BayesFactor |
Boolean, defaults to |
cores_for_props |
Integer, how many cores to use for the Bayes factor calculation, here 4 is the default for the 4 different proposal densities to evaluate, only 1, 2 and 4 are sensible. |
cores_per_prop |
Integer, how many cores to use for the Bayes factor calculation if you have more than 4 cores available. Cores used will be cores_for_props * cores_per_prop. Best to prioritize cores_for_props being 4 or 2 |
print_summary |
Boolean (default |
digits |
Integer, significant digits in printed table for information criteria |
digits_p |
Integer, significant digits in printed table for model weights |
... |
Additional, optional arguments |
Matrix of effective number of parameters, mean deviance, deviance of
mean, DIC, BPIC, Marginal Deviance (if BayesFactor=TRUE
) and associated weights.
## Not run: # 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) ## End(Not run)
## Not run: # 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) ## End(Not run)
Returns the BPIC/DIC based model weights for each participant in a list of samples objects
compare_subject( sList, stage = "sample", filter = 0, use_best_fit = TRUE, print_summary = TRUE, digits = 3 )
compare_subject( sList, stage = "sample", filter = 0, use_best_fit = TRUE, print_summary = TRUE, digits = 3 )
sList |
List of samples objects |
stage |
A string. Specifies which stage the samples are to be taken from |
filter |
An integer or vector. If it's an integer, iterations up until the value set by |
use_best_fit |
Boolean, defaults to |
print_summary |
Boolean (defaults to |
digits |
Integer, significant digits in printed table |
List of matrices for each subject of effective number of parameters, mean deviance, deviance of mean, DIC, BPIC and associated weights.
## Not run: # 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, cores_for_chains = 1) null_model <- make_emc(forstmann, design_null, cores_for_chains = 1) null_model <- fit(null_model) sList <- list(full_model, null_model) compare_subject(sList) # prints a set of weights for each model for the different participants # And returns the DIC and BPIC for each participant for each model. ## End(Not run)
## Not run: # 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, cores_for_chains = 1) null_model <- make_emc(forstmann, design_null, cores_for_chains = 1) null_model <- fit(null_model) sList <- list(full_model, null_model) compare_subject(sList) # prints a set of weights for each model for the different participants # And returns the DIC and BPIC for each participant for each model. ## End(Not run)
Similar to contr.helmert
, but then scaled to estimate differences between conditions. Use in design()
.
contr.anova(n)
contr.anova(n)
n |
An integer. The number of items for which to create the contrast |
A contrast matrix.
{ 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))) }
{ 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))) }
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).
contr.bayes(n)
contr.bayes(n)
n |
An integer. The number of items for which to create the contrast |
A contrast matrix.
{ 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))) }
{ 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))) }
Each level will be estimated as a reduction from the previous level
contr.decreasing(n)
contr.decreasing(n)
n |
an integer. The number of items for which to create the contrast. |
a contrast matrix.
{ 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))) }
{ 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))) }
Each level will be estimated additively from the previous level
contr.increasing(n)
contr.increasing(n)
n |
an integer. The number of items for which to create the contrast. |
a contrast matrix.
{ 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))) }
{ 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))) }
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.
## 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, ...)
## 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, ...)
x |
An emc object |
x_name |
A character string. Name of the parameter to be tested for |
x_fun |
Function applied to the MCMC chains to create variable to be tested. |
x_fun_name |
Name to give to quantity calculated by |
selection |
A character string designating parameter type (e.g. |
y |
A second emc object |
y_name |
A character string. Name of the parameter to be tested for |
y_fun |
Function applied to the MCMC chains to create variable to be tested. |
y_fun_name |
Name to give to quantity calculated by |
x_subject |
Integer or name selecting a subject |
y_subject |
Integer or name selecting a subject |
mu |
Numeric. |
alternative |
|
probs |
Vector defining quantiles to return. |
digits |
Integer, significant digits for estimates in printed results |
p_digits |
Integer, significant digits for probability in printed results |
print_table |
Boolean (defaults to |
... |
Additional optional arguments that can be passed to |
Invisible results table with no rounding.
## Not run: # Run a credible interval test (Bayesian ''t-test'') # Here the full model is an emc object with 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))) full_model <- make_emc(forstmann, design_full) full_model <- fit(full_model) credible(full_model, x_name = "v") # 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"]) ## End(Not run)
## Not run: # Run a credible interval test (Bayesian ''t-test'') # Here the full model is an emc object with 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))) full_model <- make_emc(forstmann, design_full) full_model <- fit(full_model) credible(full_model, x_name = "v") # 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"]) ## End(Not run)
Model file to estimate the Diffusion Decision Model (DDM) in EMC2.
DDM()
DDM()
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
.
Parameter | Transform | Natural scale | Default | Mapping | Interpretation |
v | - | [-Inf, Inf] | 1 | Mean evidence-accumulation rate (drift rate) | |
a | log | [0, Inf] | log(1) | Boundary separation | |
t0 | log | [0, Inf] | log(0) | Non-decision time | |
s | log | [0, Inf] | log(1) | Within-trial standard deviation of drift rate | |
Z | probit | [0, 1] | qnorm(0.5) | z = Z x a | Relative start point (bias) |
SZ | probit | [0, 1] | qnorm(0) | sz = 2 x SZ x min(a x Z, a x (1-Z)) | Relative between-trial variation in start point |
sv | log | [0, Inf] | log(0) | Between-trial standard deviation of drift rate | |
st0 | log | [0, Inf] | log(0) | Between-trial variation (range) in non-decision time | |
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, sv
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.
A model list with all the necessary functions for EMC2 to sample
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).
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).
This function combines information regarding the data, type of model, and the model specification.
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, ... )
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, ... )
formula |
A list. Contains the design formulae in the
format |
factors |
A named list containing all the factor variables that span
the design cells and that should be taken into account by the model.
The name Example: |
Rlevels |
A character vector. Contains the response factor levels.
Example: |
model |
A function, specifies the model type.
Choose from the drift diffusion model ( |
data |
A data frame. |
contrasts |
Optional. A named list specifying a design matrix.
Example for supplying a customized design matrix:
|
matchfun |
A function. Only needed for race models. Specifies whether a
response was correct or not. Example: |
constants |
A named vector that sets constants. Any parameter in
|
covariates |
Names of numeric covariates. |
functions |
List of functions to create new factors based on those in
the factors argument. These new factors can then be used in |
report_p_vector |
Boolean. If TRUE (default), it returns the vector of parameters to be estimated. |
custom_p_vector |
A character vector. If specified, a custom likelihood function can be supplied. |
... |
Additional, optional arguments |
A design list.
# 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)))
# 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)))
Returns the effective sample size (ESS) of the selected parameter type.
Full range of possible samples manipulations described in get_pars
.
## S3 method for class 'emc' ess_summary( emc, selection = "mu", stat = "min", stat_only = FALSE, digits = 1, ... ) ess_summary(emc, ...)
## S3 method for class 'emc' ess_summary( emc, selection = "mu", stat = "min", stat_only = FALSE, digits = 1, ... ) ess_summary(emc, ...)
emc |
An emc object |
selection |
A Character vector. Indicates which parameter types to check (e.g., |
stat |
A string. Should correspond to a function that can be applied to a vector, which will be performed on the vector/rows or columns of the matrix of the parameters |
stat_only |
Boolean. If |
digits |
Integer. How many digits to round the output to |
... |
Optional additional arguments that can be passed to |
A matrix or vector of ESS values for the selected parameter type.
ess_summary(samples_LNR, selection = "alpha")
ess_summary(samples_LNR, selection = "alpha")
General purpose function to estimate models specified in EMC2.
## S3 method for class 'emc' fit( emc, stage = NULL, iter = 1000, stop_criteria = NULL, report_time = TRUE, p_accept = 0.8, step_size = 100, verbose = TRUE, verboseProgress = FALSE, fileName = NULL, particles = NULL, particle_factor = 50, cores_per_chain = 1, cores_for_chains = length(emc), max_tries = 20, n_blocks = 1, ... ) fit(emc, ...)
## S3 method for class 'emc' fit( emc, stage = NULL, iter = 1000, stop_criteria = NULL, report_time = TRUE, p_accept = 0.8, step_size = 100, verbose = TRUE, verboseProgress = FALSE, fileName = NULL, particles = NULL, particle_factor = 50, cores_per_chain = 1, cores_for_chains = length(emc), max_tries = 20, n_blocks = 1, ... ) fit(emc, ...)
emc |
An emc object created with |
stage |
A string. Indicates which stage to start the run from, either |
iter |
An integer. Indicates how many iterations to run in the sampling stage. |
stop_criteria |
A list. Defines the stopping criteria and for which types of parameters these should hold. See the details and examples section. |
report_time |
Boolean. If |
p_accept |
A double. The target acceptance probability of the MCMC process. This fine-tunes the width of the search space to obtain the desired acceptance probability. Defaults to .8 |
step_size |
An integer. After each step, the stopping requirements as specified
by |
verbose |
Logical. Whether to print messages between each step with the current status regarding the |
verboseProgress |
Logical. Whether to print a progress bar within each step or not.
Will print one progress bar for each chain and only if |
fileName |
A string. If specified, will auto-save emc object at this location on every iteration. |
particles |
An integer. How many particles to use, default is |
particle_factor |
An integer. |
cores_per_chain |
An integer. How many cores to use per chain. Parallelizes across
participant calculations. Only available on Linux or Mac OS. For Windows, only
parallelization across chains ( |
cores_for_chains |
An integer. How many cores to use across chains.
Defaults to the number of chains. The total number of cores used is equal to |
max_tries |
An integer. How many times should it try to meet the finish
conditions as specified by |
n_blocks |
An integer. Number of blocks. Will block the parameter chains such that they are updated in blocks. This can be helpful in extremely tough models with a large number of parameters. |
... |
Additional optional arguments |
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
An emc object
## Not run: # 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 make the emc object, we've omitted a prior here for brevity so default priors will be used. emc_forstmann <- make_emc(forstmann, design) # With the emc object we can start sampling by simply calling fit emc_forstmann <- fit(emc_forstmann, fileName = "intermediate_save_location.RData") # For particularly hard models it pays off to increase the ``particle_factor`` # and, although to a lesser extent, lower ``p_accept``. emc_forstmann <- fit(emc_forstmann, particle_factor = 100, p_accept = .6) # Example of how to use the stop_criteria: emc_forstmann <- fit(emc_forstmann, stop_criteria = list(mean_gd = 1.1, max_gd = 1.5, selection = c('alpha', 'sigma2'), omit_mpsrf = TRUE, min_es = 1000)) # In this case the stop_criteria are set for the sample stage, which will be # run until the mean_gd < 1.1, the max_gd < 1.5 (omitting the multivariate psrf) # and the effective sample size > 1000, # for both the individual-subject parameters ("alpha") # and the group-level variance parameters. # For the unspecified stages in the ``stop_criteria`` the default values # are assumed which are found in Stevenson et al. 2024 <doi.org/10.31234/osf.io/2e4dq> # Alternatively, you can also specify the stop_criteria for specific stages by creating a # nested list emc_forstmann <- fit(emc_forstmann, stop_criteria = list("burn" = list(mean_gd = 1.1, max_gd = 1.5, selection = c('alpha')), "adapt" = list(min_unique = 100))) ## End(Not run)
## Not run: # 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 make the emc object, we've omitted a prior here for brevity so default priors will be used. emc_forstmann <- make_emc(forstmann, design) # With the emc object we can start sampling by simply calling fit emc_forstmann <- fit(emc_forstmann, fileName = "intermediate_save_location.RData") # For particularly hard models it pays off to increase the ``particle_factor`` # and, although to a lesser extent, lower ``p_accept``. emc_forstmann <- fit(emc_forstmann, particle_factor = 100, p_accept = .6) # Example of how to use the stop_criteria: emc_forstmann <- fit(emc_forstmann, stop_criteria = list(mean_gd = 1.1, max_gd = 1.5, selection = c('alpha', 'sigma2'), omit_mpsrf = TRUE, min_es = 1000)) # In this case the stop_criteria are set for the sample stage, which will be # run until the mean_gd < 1.1, the max_gd < 1.5 (omitting the multivariate psrf) # and the effective sample size > 1000, # for both the individual-subject parameters ("alpha") # and the group-level variance parameters. # For the unspecified stages in the ``stop_criteria`` the default values # are assumed which are found in Stevenson et al. 2024 <doi.org/10.31234/osf.io/2e4dq> # Alternatively, you can also specify the stop_criteria for specific stages by creating a # nested list emc_forstmann <- fit(emc_forstmann, stop_criteria = list("burn" = list(mean_gd = 1.1, max_gd = 1.5, selection = c('alpha')), "adapt" = list(min_unique = 100))) ## End(Not run)
A dataset containing the speed or accuracy manipulation for a Random Dot Motion experiment.
forstmann
forstmann
A data frame with 15818 rows and 5 variables:
Factor with 3 levels for Speed, Accuracy and Neutral
Factor with 2 levels for Left and Right responses
Factor with 2 levels for Left and Right trials
reaction time for each trial as a double
integer ID for each subject
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
https://www.pnas.org/doi/10.1073/pnas.0805903105
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.
## S3 method for class 'emc' gd_summary( emc, selection = "mu", omit_mpsrf = TRUE, stat = "max", stat_only = FALSE, digits = 3, ... ) gd_summary(emc, ...)
## S3 method for class 'emc' gd_summary( emc, selection = "mu", omit_mpsrf = TRUE, stat = "max", stat_only = FALSE, digits = 3, ... ) gd_summary(emc, ...)
emc |
An emc object |
selection |
A Character vector. Indicates which parameter types to check (e.g., |
omit_mpsrf |
Boolean. If |
stat |
A string. Should correspond to a function that can be applied to a vector, which will be performed on the vector/rows or columns of the matrix of the parameters |
stat_only |
Boolean. If |
digits |
Integer. How many digits to round the output to |
... |
Optional additional arguments that can be passed to |
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
.
A matrix or vector of R-hat values for the selected parameter type.
gd_summary(samples_LNR, selection = "correlation", stat = "mean", flatten = TRUE)
gd_summary(samples_LNR, selection = "correlation", stat = "mean", flatten = TRUE)
returns the Bayes Factor for two models
get_BayesFactor(MLL1, MLL2)
get_BayesFactor(MLL1, MLL2)
MLL1 |
Numeric. Marginal likelihood of model 1. Obtained with |
MLL2 |
Numeric. Marginal likelihood of model 2. Obtained with |
The BayesFactor for model 1 over model 2
## Not run: # First get the marginal likelihood for two_models # Here the full model is an emc object with the hypothesized effect # The null model is an emc object without the hypothesized effect MLL_full <- run_bridge_sampling(full_model, cores_per_prop = 3) MLL_null <- run_bridge_sampling(null_model, cores_per_prop = 3) # Now we can calculate their Bayes factor get_BayesFactor(MLL_full, MLL_null) ## End(Not run)
## Not run: # First get the marginal likelihood for two_models # Here the full model is an emc object with the hypothesized effect # The null model is an emc object without the hypothesized effect MLL_full <- run_bridge_sampling(full_model, cores_per_prop = 3) MLL_null <- run_bridge_sampling(null_model, cores_per_prop = 3) # Now we can calculate their Bayes factor get_BayesFactor(MLL_full, MLL_null) ## End(Not run)
Extracts data from an emc object
## S3 method for class 'emc' get_data(emc) get_data(emc)
## S3 method for class 'emc' get_data(emc) get_data(emc)
emc |
an emc object |
emc adds columns and rows to a dataframe in order to facilitate efficient likelihood calculations. This function will return the data as provided originally.
A dataframe of the original data
get_data(samples_LNR)
get_data(samples_LNR)
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.
get_pars( emc, selection = "mu", stage = "sample", 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 )
get_pars( emc, selection = "mu", stage = "sample", 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 )
emc |
an emc object. |
selection |
A Character string. Indicates which parameter type to select (e.g., |
stage |
A character string. Indicates from which sampling stage(s) to take the samples from (i.e. |
thin |
An integer. By how much to thin the chains |
filter |
Integer or numeric vector. If an integer is supplied, iterations up until that integer are removed. If a vector is supplied, the iterations within the range are kept. |
map |
Boolean. If |
add_recalculated |
Boolean. If |
length.out |
Integer. Alternatively to thinning, you can also select a desired length of the MCMC chains, which will be thinned appropriately. |
by_subject |
Boolean. If |
return_mcmc |
Boolean. If |
merge_chains |
Boolean. If |
subject |
Integer (vector) or character (vector). If an integer will select the 'x'th subject(s), if a character it should match subject names in the data which will be selected. |
flatten |
Boolean. If |
remove_dup |
Boolean. If |
remove_constants |
Boolean. If |
use_par |
Character (vector). If specified, only these parameters are returned. Should match the parameter names
(i.e. these are collapsed when |
type |
Character indicating the group-level model selected. Only necessary if sampler isn't specified. |
true_pars |
Set of |
chain |
Integer. Which of the chain(s) to return |
covariates |
Only needed with |
An mcmc.list object of the selected parameter types with the specified manipulations
# 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)
# 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)
Approximates the Bayes factor for parameter effects using the savage-dickey ratio.
## 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_plot_args = list(), ... ) hypothesis(emc, ...)
## 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_plot_args = list(), ... ) hypothesis(emc, ...)
emc |
An emc object |
parameter |
A string. A parameter which you want to compare to H0. Will not be used if a FUN is specified. |
H0 |
An integer. The H0 value which you want to compare to |
fun |
A function. Specifies an operation to be performed on the sampled or mapped parameters. |
selection |
A Character string. Indicates which parameter type to use (e.g., |
do_plot |
Boolean. If |
use_prior_lim |
Boolean. If |
N |
Integer. How many prior samples to draw |
prior_plot_args |
A list. Optional additional arguments to be passed to plot.default for the plotting of the prior density (see |
... |
Optional arguments that can be passed to |
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).
The Bayes factor for the hypothesis against H0.
# 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)
# 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)
Adds a set of start points to each chain. These start points are sampled from a user-defined multivariate normal across subjects.
init_chains( emc, start_mu = NULL, start_var = NULL, particles = 1000, cores_per_chain = 1, cores_for_chains = length(emc) )
init_chains( emc, start_mu = NULL, start_var = NULL, particles = 1000, cores_per_chain = 1, cores_for_chains = length(emc) )
emc |
An emc object made by |
start_mu |
A vector. Mean of multivariate normal used in proposal distribution |
start_var |
A matrix. Variance covariance matrix of multivariate normal used in proposal distribution. Smaller values will lead to less deviation around the mean. |
particles |
An integer. Number of starting values |
cores_per_chain |
An integer. How many cores to use per chain. Parallelizes across participant calculations. |
cores_for_chains |
An integer. How many cores to use to parallelize across chains. Default is the number of chains. |
An emc object
## Not run: # 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, Z~1, sv~1, SZ~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),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5)) # 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 = p_vector, start_var = var, cores_per_chain = 4) # Afterwards we can just use fit DDMaE <- fit(DDMaE, cores_per_chain = 4) ## End(Not run)
## Not run: # 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, Z~1, sv~1, SZ~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),Z=qnorm(.5),sv=log(.5),SZ=qnorm(.5)) # 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 = p_vector, start_var = var, cores_per_chain = 4) # Afterwards we can just use fit DDMaE <- fit(DDMaE, cores_per_chain = 4) ## End(Not run)
Model file to estimate the Linear Ballistic Accumulator (LBA) in EMC2.
LBA()
LBA()
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
.
Parameter | Transform | Natural scale | Default | Mapping | Interpretation |
v | - | [-Inf, Inf] | 1 | Mean evidence-accumulation rate | |
A | log | [0, Inf] | log(0) | Between-trial variation (range) in start point | |
B | log | [0, Inf] | log(1) | b = B+A | Distance from A to b (response threshold) |
t0 | log | [0, Inf] | log(0) | Non-decision time | |
sv | log | [0, Inf] | log(1) | Between-trial variation in evidence-accumulation rate | |
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
A model list with all the necessary functions for EMC2 to sample
# 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).
# 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).
Model file to estimate the Log-Normal Race Model (LNR) in EMC2.
LNR()
LNR()
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
.
Parameter | Transform | Natural scale | Default | Mapping | Interpretation |
m | - | [-Inf, Inf] | 1 | Scale parameter | |
s | log | [0, Inf] | log(1) | Shape parameter | |
t0 | log | [0, Inf] | log(0) | Non-decision time | |
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
A model list with all the necessary functions for EMC2 to sample
# 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).
# 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).
Simulates data based on a model design and a parameter vector (p_vector
) by one of two methods:
Creating a fully crossed and balanced design specified by the design,
with number of trials per cell specified by the n_trials
argument
Using the design of a data frame supplied, which allows creation of unbalanced and other irregular designs, and replacing previous data with simulated data
make_data( parameters, design = NULL, n_trials = NULL, data = NULL, expand = 1, mapped_p = FALSE, hyper = FALSE, ... )
make_data( parameters, design = NULL, n_trials = NULL, data = NULL, expand = 1, mapped_p = FALSE, hyper = FALSE, ... )
parameters |
parameter vector used to simulate data.
Can also be a matrix with one row per subject (with corresponding row names)
or an emc object with sampled parameters
(in which case posterior medians of |
design |
Design list created by |
n_trials |
Integer. If |
data |
Data frame. If supplied, the factors are taken from the data. Determines the number of trials per level of the design factors and can thus allow for unbalanced designs |
expand |
Integer. Replicates the |
mapped_p |
If |
hyper |
If |
... |
Additional optional arguments |
To create data for multiple subjects see ?make_random_effects()
.
A data frame with simulated data
# 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)
# 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)
Creates an emc object by combining the data, prior,
and model specification into a emc
object that is needed in fit()
.
make_emc( data, design, model = NULL, type = "standard", n_chains = 3, compress = TRUE, rt_resolution = 0.02, prior_list = NULL, grouped_pars = NULL, par_groups = NULL, ... )
make_emc( data, design, model = NULL, type = "standard", n_chains = 3, compress = TRUE, rt_resolution = 0.02, prior_list = NULL, grouped_pars = NULL, par_groups = NULL, ... )
data |
A data frame, or a list of data frames. Needs to have the variable |
design |
A list with a pre-specified design, the output of |
model |
A model list. If none is supplied, the model specified in |
type |
A string indicating whether to run a |
n_chains |
An integer. Specifies the number of mcmc chains to be run (has to be more than 1 to compute |
compress |
A Boolean, if |
rt_resolution |
A double. Used for compression, response times will be binned based on this resolution. |
prior_list |
A named list containing the prior. Default prior created if |
grouped_pars |
An integer vector. Parameters on this location of the vector of parameters are treated as constant across subjects |
par_groups |
A vector. Only to be specified with type |
... |
Additional, optional arguments. |
An uninitialized emc object
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)
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)
Simulates subject-level parameters in the format required by make_data()
.
make_random_effects( design, group_means, n_subj = NULL, variance_proportion = 0.2, covariances = NULL )
make_random_effects( design, group_means, n_subj = NULL, variance_proportion = 0.2, covariances = NULL )
design |
A design list. The design as specified by |
group_means |
A numeric vector. The group level means for each parameter, in the same order as |
n_subj |
An integer. The number of subjects to generate parameters for. If |
variance_proportion |
A double. Optional. If |
covariances |
A covariance matrix. Optional. Specify the intended covariance matrix. |
A matrix of subject-level parameters.
# 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 = 5) # We can also define a covariance matrix to simulate from subj_pars <- make_random_effects(design_DDMaE, group_means, n_subj = 5, 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)
# 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 = 5) # We can also define a covariance matrix to simulate from subj_pars <- make_random_effects(design_DDMaE, group_means, n_subj = 5, 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)
Maps a parameter vector that corresponds to sampled parameters
of the cognitive model back to the experimental design. The parameter vector
can be created using sampled_p_vector()
. The returned matrix shows whether/how parameters
differ across the experimental factors.
mapped_par( p_vector, design, model = NULL, digits = 3, remove_subjects = TRUE, covariates = NULL, ... )
mapped_par( p_vector, design, model = NULL, digits = 3, remove_subjects = TRUE, covariates = NULL, ... )
p_vector |
A parameter vector. Must be in the form of |
design |
A design list. Created by |
model |
Optional model type (if not already specified in |
digits |
Integer. Will round the output parameter values to this many decimals |
remove_subjects |
Boolean. Whether to include subjects as a factor in the design |
covariates |
Covariates specified in the design can be included here. |
... |
optional arguments |
Matrix with a column for each factor in the design and for each model parameter type (p_type
).
# 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 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_par(p_vector,design_DDMaE)
# 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 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_par(p_vector,design_DDMaE)
Merges samples from all chains as one unlisted object.
merge_chains(emc)
merge_chains(emc)
emc |
An emc object, commonly the output of |
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.
An unlisted emc object with all chains merged
Plots within-chain parameter correlations (upper triangle) and corresponding scatterplots (lower triangle) to visualize parameter sloppiness.
pairs_posterior( emc, selection = "alpha", scale_subjects = TRUE, do_plot = TRUE, N = 500, ... )
pairs_posterior( emc, selection = "alpha", scale_subjects = TRUE, do_plot = TRUE, N = 500, ... )
emc |
An emc object |
selection |
A Character string. Indicates which parameter type to
plot ( |
scale_subjects |
Boolean. To standardize each participant with |
do_plot |
Boolean. Whether to plot the pairs plot, if |
N |
Integer for maximum number of iterations used (defaults to 500). If number of samples in stage or selection exceeds N, a random subset will be taken of size N |
... |
Optional arguments that can be passed to |
If selection = alpha
the parameter chains are concatenated across participants,
(after standardizing if scale_subjects = TRUE
) and then correlated.
Invisibly returns a matrix with the correlations between the parameters.
# 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")
# 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")
Returns a parameter type from an emc object as a data frame.
## S3 method for class 'emc' parameters(emc, selection = "mu", N = NULL, resample = FALSE, ...) parameters(emc, ...)
## S3 method for class 'emc' parameters(emc, selection = "mu", N = NULL, resample = FALSE, ...) parameters(emc, ...)
emc |
An emc object |
selection |
String designating parameter type (e.g. mu, sigma2, correlation, alpha) |
N |
Integer. How many samples to take from the posterior. If |
resample |
Boolean. If |
... |
Optional arguments that can be passed to |
A data frame with one row for each sample (with a subjects column if selection = "alpha")
Plots panels that contain a set of densities for each response option in the data. These densities are defective; their areas are relative to the respective response proportion. Across all responses, the area sums to 1.
plot_defective_density( data, subject = NULL, factors = NULL, layout = NA, correct_fun = NULL, rt_pos = "top", accuracy = "topright", ... )
plot_defective_density( data, subject = NULL, factors = NULL, layout = NA, correct_fun = NULL, rt_pos = "top", accuracy = "topright", ... )
data |
A data frame. The experimental data in EMC2 format with at least |
subject |
An integer or character string selecting a subject from the data.
If specified, only that subject is plotted. Per default (i.e., |
factors |
A character vector of the factor names in the design to aggregate across
Defaults to all (i.e., |
layout |
A vector indicating which layout to use as in par(mfrow = layout). If NA, will automatically generate an appropriate layout. |
correct_fun |
If specified, the accuracy for each subject is calculated, using the supplied function and an accuracy vector for each subject is returned invisibly. |
rt_pos |
legend function position character string for the mean response time (defaults to |
accuracy |
legend function position character string for accuracy (defaults to |
... |
Optional arguments that can be passed to |
If correct_fun
is specified, a subject accuracy vector is returned invisibly
# First for each subject and the factor combination in the design: plot_defective_density(forstmann) # Now collapsing across subjects: plot_defective_density(forstmann, factors = c("S", "E")) # If the data is response coded, it generally makes sense to include the "S" factor # because EMC2 will plot the "R" factor automatically. This way, choice accuracy can # be examined # Each subject's accuracy can be returned using a custom function: print(plot_defective_density(forstmann, correct_fun = function(d) d$R == d$S))
# First for each subject and the factor combination in the design: plot_defective_density(forstmann) # Now collapsing across subjects: plot_defective_density(forstmann, factors = c("S", "E")) # If the data is response coded, it generally makes sense to include the "S" factor # because EMC2 will plot the "R" factor automatically. This way, choice accuracy can # be examined # Each subject's accuracy can be returned using a custom function: print(plot_defective_density(forstmann, correct_fun = function(d) d$R == d$S))
Plot (defective) cumulative density functions of the observed data and data from the posterior predictive distribution: the probability of a response, p(R) as a function of response time for the experimental data and posterior predictive simulations.
plot_fit( data, pp, subject = NULL, factors = NULL, functions = NULL, stat = NULL, stat_name = "", adjust = 1, quants = c(0.025, 0.5, 0.975), do_plot = TRUE, xlim = NULL, ylim = NULL, layout = NULL, mfcol = FALSE, probs = c(1:99)/100, data_lwd = 2, fit_lwd = 1, q_points = c(0.1, 0.3, 0.5, 0.7, 0.9), qp_cex = 1, pqp_cex = 0.5, lpos = "right", main = "" )
plot_fit( data, pp, subject = NULL, factors = NULL, functions = NULL, stat = NULL, stat_name = "", adjust = 1, quants = c(0.025, 0.5, 0.975), do_plot = TRUE, xlim = NULL, ylim = NULL, layout = NULL, mfcol = FALSE, probs = c(1:99)/100, data_lwd = 2, fit_lwd = 1, q_points = c(0.1, 0.3, 0.5, 0.7, 0.9), qp_cex = 1, pqp_cex = 0.5, lpos = "right", main = "" )
data |
A data frame. The experimental data in EMC2 format with at least |
pp |
A data frame. Posterior predictives created by |
subject |
Integer or string selecting a subject from the data. If specified only that subject
is plotted. |
factors |
Character vector of factors in data to display separately. If NULL (default) use names of all columns in data except "trials","R", and "rt". Omitted factors are aggregated over. If NA treats entire data set as a single cell. Must be NA or NULL when using stat argument. |
functions |
A named list of functions that create new factors which can then be used by the factors and stat arguments. |
stat |
A function that takes the data/the posterior predictives and returns a single value. For the posterior predictives it will use a single value per replicate, which are then plotted as a density. |
stat_name |
A string naming what the stat argument calculates, used in labeling the x-axis of the plot. |
adjust |
Numeric. Density function bandwidth adjust parameter. See “?density' |
quants |
A vector. Quantiles of the posterior predictives to return when stat argument is supplied. |
do_plot |
Boolean. Set to |
xlim |
A numeric vector. x-axis plot limit. |
ylim |
A numeric vector. y-axis plot limit. |
layout |
A vector specifying the layout as in |
mfcol |
Boolean. If |
probs |
Vector of probabilities at which to calculate cumulative density function |
data_lwd |
Integer. Line width for data |
fit_lwd |
Integer. Line width for posterior predictives |
q_points |
Vector. Quantile points to plot |
qp_cex |
Numeric. Cex for data quantile points |
pqp_cex |
Numeric. Cex for predicted quantile points |
lpos |
Character. Legend position, see |
main |
Character. Pasted before the plot title, especially useful when specifying a stat argument. |
The data is plotted in black. Large grey points show the average quantiles across the posterior predictives. The small grey points represent the predicted quantile of an individual replicate, providing a representation of uncertainty in the model predictions.
If the stat argument is supplied (which calculates a statistic based on the data), the posterior predictives are plotted as a density over the different replicates. A vertical line is plotted at the value of that statistic for the experimental data.
If more than one subject is included, the data and fits are aggregated across subjects by default.
Also see ?plot_defective_density()
for more details.
If stat argument is provided, a vector of observed values and predicted quantiles is returned
# First generate posterior predictives based on an emc object run with run_emc pp <- predict(samples_LNR, n_cores = 1, n_post = 10) # Then visualize the model fit plot_fit(forstmann, pp, factors = c("S", "E"), layout = c(2,3)) # Specific statistics on the posterior predictives can also be specified # This function calculates the difference in rt between two S levels. # It takes the data (or the posterior predictives) as an argument drt <- function(data) diff(tapply(data$rt,data[,c("S")],mean)) plot_fit(forstmann, pp, stat=drt,stat_name="Rt difference", main=("Left vs Right"))
# First generate posterior predictives based on an emc object run with run_emc pp <- predict(samples_LNR, n_cores = 1, n_post = 10) # Then visualize the model fit plot_fit(forstmann, pp, factors = c("S", "E"), layout = c(2,3)) # Specific statistics on the posterior predictives can also be specified # This function calculates the difference in rt between two S levels. # It takes the data (or the posterior predictives) as an argument drt <- function(data) diff(tapply(data$rt,data[,c("S")],mean)) plot_fit(forstmann, pp, stat=drt,stat_name="Rt difference", main=("Left vs Right"))
Plots the posterior and prior density for selected parameters of a model.
Full range of samples manipulations described in get_pars
.
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_plot_args = list(), true_plot_args = list(), ... )
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_plot_args = list(), true_plot_args = list(), ... )
emc |
An emc object |
layout |
A vector indicating which layout to use as in par(mfrow = layout). If NA, will automatically generate an appropriate layout. |
selection |
A Character string. Indicates which parameter type to use (e.g., |
show_chains |
Boolean (defaults to |
plot_prior |
Boolean. If |
N |
Integer. How many prior samples to draw |
use_prior_lim |
Boolean. If |
lpos |
Character. Where to plot the contraction statistic. |
true_pars |
A vector or emc object. Can be used to visualize recovery. If a vector will plot a vertical line for each parameter at the appropriate place. If an emc object will plot the densities of the object as well, assumed to be the data-generating posteriors. |
all_subjects |
Boolean. Will plot the densities of all (selected) subjects overlaid with the group-level distribution |
prior_plot_args |
A list. Optional additional arguments to be passed to plot.default for the plotting of the prior density (see |
true_plot_args |
A list. Optional additional arguments to be passed to plot.default for the plotting of the true parameters (see |
... |
Optional arguments that can be passed to |
An invisible return of the contraction statistics for the selected parameter type
# 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_plot_args = list(col = 'blue'), adjust = 2)
# 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_plot_args = list(col = 'blue'), adjust = 2)
Title
plot_prior( prior, design, selection = "mu", do_plot = TRUE, covariates = NULL, layout = NA, N = 50000, ... )
plot_prior( prior, design, selection = "mu", do_plot = TRUE, covariates = NULL, layout = NA, N = 50000, ... )
prior |
A prior list created with |
design |
A design list created with |
selection |
A Character string. Indicates which parameter type to use (e.g., |
do_plot |
Boolean. If |
covariates |
dataframe/functions as specified by the design |
layout |
A vector indicating which layout to use as in par(mfrow = layout). If NA, will automatically generate an appropriate layout. |
N |
Integer. How many prior samples to draw |
... |
Optional arguments that can be passed to get_pars, histogram, plot.default (see par()), or arguments required for the types of models e.g. n_factors for type = "factor" |
An mcmc.list object with prior samples of the selected type
# 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(prior_DDMaE, design_DDMaE, selection = "mu", N = 1e3) plot_prior(prior_DDMaE, design_DDMaE, selection = "mu", mapped = FALSE, N=1e3) # We can also plot the implied prior on the participant level effects. plot_prior(prior_DDMaE, design_DDMaE, selection = "alpha", col = "green", N = 1e3)
# 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(prior_DDMaE, design_DDMaE, selection = "mu", N = 1e3) plot_prior(prior_DDMaE, design_DDMaE, selection = "mu", mapped = FALSE, N=1e3) # We can also plot the implied prior on the participant level effects. plot_prior(prior_DDMaE, design_DDMaE, selection = "alpha", col = "green", N = 1e3)
An adjusted version of the corrplot
package function corrplot()
tailored
to EMC2
and the plotting of estimated correlations.
plot_relations( emc = NULL, stage = "sample", plot_cred = TRUE, plot_means = TRUE, only_cred = FALSE, nice_names = NULL, ... )
plot_relations( emc = NULL, stage = "sample", plot_cred = TRUE, plot_means = TRUE, only_cred = FALSE, nice_names = NULL, ... )
emc |
An EMC2 object, commonly the output of |
stage |
Character. The stage from which to take the samples, defaults to
the sampling stage |
plot_cred |
Boolean. Whether to plot the 95 percent credible intervals or not |
plot_means |
Boolean. Whether to plot the means or not |
only_cred |
Boolean. Whether to only plot credible values |
nice_names |
Character string. Alternative names to give the parameters |
... |
Optional additional arguments |
No return value, creates a plot of group-level relations
# 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)
# 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)
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.
plot_sbc_ecdf(ranks, layout = NA)
plot_sbc_ecdf(ranks, layout = NA)
ranks |
A list of named dataframes of the rank statistic |
layout |
Optional. A numeric vector specifying the layout using |
No returns
Note that this plot is dependent on the number of bins, and a more general
visualization is to use plot_sbc_ecdf
plot_sbc_hist(ranks, bins = 10, layout = NA)
plot_sbc_hist(ranks, bins = 10, layout = NA)
ranks |
A list of named dataframes of the rank statistic |
bins |
An integer specifying the number of bins to use when plotting the histogram |
layout |
Optional. A numeric vector specifying the layout using |
No returns
Makes trace plots for model parameters.
## S3 method for class 'emc' plot( x, stage = "sample", selection = c("mu", "sigma2", "alpha"), layout = NA, ... )
## S3 method for class 'emc' plot( x, stage = "sample", selection = c("mu", "sigma2", "alpha"), layout = NA, ... )
x |
An object of class |
stage |
A character string indicating the sampling stage to be summarized.
Can be |
selection |
A character vector indicating the parameter group(s).
Defaults to |
layout |
A vector indicating which layout to use as in par(mfrow = layout). If NA, will automatically generate an appropriate layout. |
... |
Optional arguments that can be passed to |
A trace/acf plot of the selected MCMC chains
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(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)
Returns the quantiles of the selected parameter type.
Full range of possible samples manipulations described in get_pars
.
## S3 method for class 'emc' posterior_summary( emc, selection = "mu", probs = c(0.025, 0.5, 0.975), digits = 3, ... ) posterior_summary(emc, ...)
## S3 method for class 'emc' posterior_summary( emc, selection = "mu", probs = c(0.025, 0.5, 0.975), digits = 3, ... ) posterior_summary(emc, ...)
emc |
An emc object |
selection |
A Character vector. Indicates which parameter types to check (e.g., |
probs |
A vector. Indicates which quantiles to return from the posterior. |
digits |
Integer. How many digits to round the output to |
... |
Optional additional arguments that can be passed to |
A list of posterior quantiles for each parameter group in the selected parameter type.
posterior_summary(samples_LNR)
posterior_summary(samples_LNR)
Simulate n_post
data sets using the posterior parameter estimates
## S3 method for class 'emc' predict( object, hyper = FALSE, n_post = 100, n_cores = 1, stat = c("random", "mean", "median")[1], ... )
## S3 method for class 'emc' predict( object, hyper = FALSE, n_post = 100, n_cores = 1, stat = c("random", "mean", "median")[1], ... )
object |
An emc object from which posterior predictives should be generated |
hyper |
Boolean. Defaults to |
n_post |
Integer. Number of generated datasets |
n_cores |
Integer. Number of cores across which there should be parallellized |
stat |
Character. Can be |
... |
Optional additional arguments passed to |
A list of simulated data sets of length n_post
# based on an emc object ran by fit() we can generate posterior predictives predict(samples_LNR, n_cores = 1, n_post = 10)
# 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. These values are entered manually by default but can be
recycled from another prior (given in the update
argument).
prior( design, type = "standard", update = NULL, do_ask = NULL, fill_default = TRUE, ... )
prior( design, type = "standard", update = NULL, do_ask = NULL, fill_default = TRUE, ... )
design |
Design list for which a prior is constructed, typically the output of |
type |
Character. What type of group-level model you plan on using i.e. |
update |
Prior list from which to copy values |
do_ask |
Character. For which parameter types or hyperparameters to ask for prior specification,
i.e. |
fill_default |
Boolean, If |
... |
Either values to prefill, i.e. |
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'.
A prior list object
# 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)
# 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)
Creates likelihood profile plots from a design and the experimental data by varying one model parameter while holding all others constant.
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_plot_args = list(), ... )
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_plot_args = list(), ... )
data |
A dataframe. Experimental data used, needed for the design mapping |
design |
A design list. Created using |
p_vector |
Named vector of parameter values (typically created with |
range |
Numeric. The max and min will be p_vector + range/2 and p_vector - range/2, unless specified in p_min or p_max. |
layout |
A vector indicating which layout to use as in par(mfrow = layout). If NA, will automatically generate an appropriate layout. |
p_min |
Named vector. If specified will instead use these values for minimum range of the selected parameters. |
p_max |
Named vector. If specified will instead use these values for maximum range of the selected parameters. |
use_par |
Character vector. If specified will only plot the profiles for the specified parameters. |
n_point |
Integer. Number of evenly spaced points at which to calculate likelihood |
n_cores |
Number of likelihood points evenly spaced between the minimum and maximum likelihood range. |
round |
Integer. To how many digits will the output be rounded. |
true_plot_args |
A list. Optional additional arguments that can be passed to plot.default for the plotting of the true vertical line. |
... |
Optional additional arguments that can be passed to plot.default. |
Vector with highest likelihood point, input and mismatch between true and highest point
# 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)
# 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)
Model file to estimate the Racing Diffusion Model (RDM), also known as the Racing Wald Model.
RDM()
RDM()
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
.
Parameter | Transform | Natural scale | Default | Mapping | Interpretation |
v | log | [0, Inf] | log(1) | Evidence-accumulation rate (drift rate) | |
A | log | [0, Inf] | log(0) | Between-trial variation (range) in start point | |
B | log | [0, Inf] | log(1) | b = B + A | Distance from A to b (response threshold) |
t0 | log | [0, Inf] | log(0) | Non-decision time | |
sv | log | [0, Inf] | log(1) | Within-trial standard deviation of drift rate | |
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
A list defining the cognitive model
# 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).
# 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).
Plots recovery of data generating parameters/samples.
Full range of samples manipulations described in get_pars
## 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, ...)
## 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, ...)
emc |
An emc object |
true_pars |
A vector of data-generating parameters or an emc object with data-generating samples |
selection |
A Character vector. Indicates which parameter types to plot (e.g., |
layout |
A vector indicating which layout to use as in par(mfrow = layout). If NA, will automatically generate an appropriate layout. |
do_CI |
Boolean. If |
correlation |
Character. Which correlation to include in the plot. Options are either |
stat |
Character. Which statistic to include in the plot. Options are either |
digits |
Integer. How many digits to round the statistic and correlation in the plot to |
CI |
Numeric. The size of the credible intervals. Default is .95 (95%). |
ci_plot_args |
A list. Optional additional arguments to be passed to plot.default for the plotting of the credible intervals (see |
... |
Optional arguments that can be passed to |
Invisible list with RMSE, coverage, and Pearson and Spearman correlations.
# 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"))
# 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"))
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.
run_bridge_sampling( emc, stage = "sample", filter = NULL, repetitions = 1, cores_for_props = 4, cores_per_prop = 1, both_splits = TRUE, ... )
run_bridge_sampling( emc, stage = "sample", filter = NULL, repetitions = 1, cores_for_props = 4, cores_per_prop = 1, both_splits = TRUE, ... )
emc |
An emc object with a set of converged samples |
stage |
A character indicating which stage to use, defaults to |
filter |
An integer or vector. If integer, it will exclude up until that integer. If vector it will include everything in that range. |
repetitions |
An integer. How many times to repeat the bridge sampling scheme. Can help get an estimate of stability of the estimate. |
cores_for_props |
Integer. Warp-III evaluates the posterior over 4 different proposal densities. If you have the CPU, 4 cores will do this in parallel, 2 is also already helpful. |
cores_per_prop |
Integer. Per density we can also parallelize across subjects. Eventual cores will be |
both_splits |
Boolean. Bridge sampling uses a proposal density and a target density. We can estimate the stability of our samples and therefore MLL estimate, by running 2 bridge sampling iterations
The first one uses the first half of the samples as the proposal and the second half as the target, the second run uses the opposite. If this is is set to |
... |
Additional, optional more in-depth hyperparameters |
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
A vector of length repetitions which contains the marginal log likelihood estimates per repetition
## Not run: # 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(list(samples_LNR), cores_per_prop = 2) # This will run on 2*4 cores (since 4 is the default for ``cores_for_props``) ## End(Not run)
## Not run: # 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(list(samples_LNR), cores_per_prop = 2) # This will run on 2*4 cores (since 4 is the default for ``cores_for_props``) ## End(Not run)
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
).
run_emc( emc, stage, stop_criteria, p_accept = 0.8, step_size = 100, verbose = FALSE, verboseProgress = FALSE, fileName = NULL, particles = NULL, particle_factor = 50, cores_per_chain = 1, cores_for_chains = length(emc), max_tries = 20, n_blocks = 1 )
run_emc( emc, stage, stop_criteria, p_accept = 0.8, step_size = 100, verbose = FALSE, verboseProgress = FALSE, fileName = NULL, particles = NULL, particle_factor = 50, cores_per_chain = 1, cores_for_chains = length(emc), max_tries = 20, n_blocks = 1 )
emc |
An emc object |
stage |
A string. Indicates which stage is to be run, either |
stop_criteria |
A list. Defines the stopping criteria and for which types of parameters these should hold. See |
p_accept |
A double. The target acceptance probability of the MCMC process. This fine-tunes the width of the search space to obtain the desired acceptance probability. Defaults to .8 |
step_size |
An integer. After each step, the stopping requirements as
specified by |
verbose |
Logical. Whether to print messages between each step with the current status regarding the stop_criteria. |
verboseProgress |
Logical. Whether to print a progress bar within each step or not. Will print one progress bar for each chain and only if cores_for_chains = 1. |
fileName |
A string. If specified will autosave emc at this location on every iteration. |
particles |
An integer. How many particles to use, default is |
particle_factor |
An integer. |
cores_per_chain |
An integer. How many cores to use per chain.
Parallelizes across participant calculations. Only available on Linux or Mac OS.
For Windows, only parallelization across chains ( |
cores_for_chains |
An integer. How many cores to use across chains.
Defaults to the number of chains. the total number of cores used is equal to |
max_tries |
An integer. How many times should it try to meet the finish conditions as specified by stop_criteria? Defaults to 20. max_tries is ignored if the required number of iterations has not been reached yet. |
n_blocks |
An integer. Number of blocks. Will block the parameter chains such that they are updated in blocks. This can be helpful in extremely tough models with a large number of parameters. |
An emc object
## Not run: # 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 make the emc, we've omitted a prior here for brevity so default priors will be used. emc <- make_emc(forstmann, design) # Now for example we can specify that we only want to run the "preburn" phase # for MCMC 200 iterations emc <- run_emc(emc, stage = "preburn", stop_criteria = list(iter = 200)) ## End(Not run)
## Not run: # 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 make the emc, we've omitted a prior here for brevity so default priors will be used. emc <- make_emc(forstmann, design) # Now for example we can specify that we only want to run the "preburn" phase # for MCMC 200 iterations emc <- run_emc(emc, stage = "preburn", stop_criteria = list(iter = 200)) ## End(Not run)
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.
run_sbc( design_in, prior_in, replicates = 250, trials = 100, n_subjects = 30, plot_data = FALSE, verbose = TRUE, fileName = NULL, ... )
run_sbc( design_in, prior_in, replicates = 250, trials = 100, n_subjects = 30, plot_data = FALSE, verbose = TRUE, fileName = NULL, ... )
design_in |
An emc design list. The design of the model to be used in SBC |
prior_in |
An emc prior list. The prior for the design to be used in SBC |
replicates |
Integer. The number of samples to draw from the prior |
trials |
Integer. The number of trials of the simulated data (per subject) |
n_subjects |
Integer. Only used for hierarchical models. The number of subjects to be used in data generation of each replicate |
plot_data |
Boolean. Whether to plot the data simulated (aggregated across subjects) |
verbose |
Verbose. Whether to print progress related messages |
fileName |
Character. Highly recommended, saves temporary results to the fileName |
... |
A list of optional additional arguments that can be passed to |
The ranks and prior samples. For hierarchical models also the prior-generated subject-level parameters.
Makes a vector with zeroes, with names and length corresponding to the model parameters of the design.
sampled_p_vector( design, model = NULL, doMap = TRUE, add_da = FALSE, all_cells_dm = FALSE )
sampled_p_vector( design, model = NULL, doMap = TRUE, add_da = FALSE, all_cells_dm = FALSE )
design |
a list of the design made with |
model |
a model list. Defaults to the model specified in the design list. |
doMap |
logical. If |
add_da |
Boolean. Whether to include the relevant data columns in the map attribute |
all_cells_dm |
Boolean. Whether to include all levels of a factor in the mapping attribute, even when one is dropped in the design |
Named vector.
# 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_p_vector(design_DDMaE)
# 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_p_vector(design_DDMaE)
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.
samples_LNR
samples_LNR
An emc object. An emc object is a list with a specific structure and elements, as outlined below.
A list of dataframes, one for each subject included
A character vector containing the model parameter names
The number of parameters in the model
The number of unique subject ID's in the data
A vector containing the unique subject ID's
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
)
The log likelihood function used by pmwg for model estimation
A list with defined structure containing the samples, see the Samples Element section for more detail
Which parameters are grouped across subjects, in this case none
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.
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.
samples used for estimating the model parameters (group level), an array of size (n_pars x n_samples)
samples used for estimating the parameter covariance matrix, an array of size (n_pars x n_pars x n_samples)
samples used for estimating the subject random effects, an array of size (n_pars x n_subjects x n_samples)
A vector containing what PMwG stage each sample was drawn in
The winning particles log-likelihood for each subject and sample
Mixing weights used during the Gibbs step when creating a new sample for the covariance matrix
The inverse of the last samples covariance matrix
The index of the last sample drawn
The scaling parameter of the proposal distributions for each subject array of size (n_subjects x n_samples)
From which propoosal distribution the accepted samples in the MCMC chain came, an array of size (n_subjects x n_samples)
https://www.pnas.org/doi/10.1073/pnas.0805903105
Shorten an emc object
## S3 method for class 'emc' subset( x, stage = "sample", filter = NULL, thin = 1, keep_stages = FALSE, length.out = NULL, ... )
## S3 method for class 'emc' subset( x, stage = "sample", filter = NULL, thin = 1, keep_stages = FALSE, length.out = NULL, ... )
x |
an emc object |
stage |
A character string. Indicates from which sampling stage(s) to take the samples from (i.e. |
filter |
Integer or numeric vector. If an integer is supplied, iterations up until that integer are removed. If a vector is supplied, the iterations within the range are kept. |
thin |
An integer. By how much to thin the chains |
keep_stages |
Boolean. If |
length.out |
Integer. Alternatively to thinning, you can also select a desired length of the MCMC chains, which will be thinned appropriately. |
... |
additional optional arguments |
A shortened emc object
subset(samples_LNR, length.out = 10)
subset(samples_LNR, length.out = 10)
Computes quantiles, Rhat
and ESS
for selected model parameters.
## S3 method for class 'emc' summary( object, selection = c("mu", "sigma2", "alpha"), probs = c(0.025, 0.5, 0.975), digits = 3, ... )
## S3 method for class 'emc' summary( object, selection = c("mu", "sigma2", "alpha"), probs = c(0.025, 0.5, 0.975), digits = 3, ... )
object |
An object of class |
selection |
A character string indicating the parameter type
Defaults to |
probs |
The quantiles to be computed. Defaults to the the 2.5%, 50% and 97.5% quantiles. |
digits |
An integer specifying rounding of output. |
... |
Optional arguments that can be passed to |
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.
A list of summary output.