| Type: | Package |
| Title: | Cognitive Models |
| Version: | 0.2.6.2 |
| Date: | 2025-04-13 |
| Maintainer: | Yi-Shin Lin <yishinlin001@gmail.com> |
| Description: | Hierarchical Bayesian models. The package provides tools to fit two response time models, using the population-based Markov Chain Monte Carlo. |
| License: | GPL-2 |
| URL: | https://github.com/yxlin/ggdmc |
| BugReports: | https://github.com/yxlin/ggdmc/issues |
| Imports: | Rcpp (≥ 1.0.7), coda, stats, utils, ggplot2, matrixStats, data.table (≥ 1.10.4) |
| Depends: | R (≥ 3.3.0) |
| LinkingTo: | Rcpp (≥ 1.0.7), RcppArmadillo (≥ 0.7.100.3.0) |
| Suggests: | testthat |
| RoxygenNote: | 7.3.2 |
| Encoding: | UTF-8 |
| NeedsCompilation: | yes |
| Packaged: | 2025-04-13 11:29:48 UTC; yslin |
| Repository: | CRAN |
| Author: | Yi-Shin Lin [aut, cre], Andrew Heathcote [aut] |
| Date/Publication: | 2025-04-13 11:50:02 UTC |
Bind data and models
Description
Binding a data set with a model object. The function also checks whether they are compatible and adds attributes on a data model instance.
Usage
BuildDMI(x, model)
Arguments
x |
data as in data frame |
model |
a model object |
Value
a data model instance
Create a model object
Description
A model object consists of arraies with model attributes.
Usage
BuildModel(
p.map,
responses,
factors = list(A = "1"),
match.map = NULL,
constants = numeric(0),
type = "norm",
posdrift = TRUE,
verbose = TRUE
)
## S3 method for class 'model'
print(x, p.vector = NULL, ...)
## S3 method for class 'dmi'
print(x, ...)
Arguments
p.map |
parameter map. This option maps a particular factorial design to model parameters |
responses |
specifying the response names and levels |
factors |
specifying a list of factors and their levels |
match.map |
match map. This option matches stimuli and responses |
constants |
specifying the parameters with fixed values |
type |
specifying model type, either "rd" or "norm". |
posdrift |
a Boolean, switching between enforcing strict postive drift rates by using truncated normal distribution. This option is only useful in "norm" model type. |
verbose |
Print p.vector, constants and model type |
x |
a model object |
p.vector |
parameter vector |
... |
other arguments |
Examples
model <- BuildModel(
p.map = list(a = "1", v = "1", z = "1", d = "1", t0 = "1",
sv = "1", sz = "1", st0 = "1"),
constants = c(st0 = 0, d = 0, sz = 0, sv = 0),
match.map = list(M = list(s1 = "r1", s2 = "r2")),
factors = list(S = c("s1", "s2")),
responses = c("r1", "r2"),
type = "rd")
Specifying Parameter Prior Distributions
Description
BuildPrior sets up parameter prior distributions for each model
parameter. p1 and p2 refer to the first and second parameters
a prior distribution.
Usage
BuildPrior(
p1,
p2,
lower = rep(NA, length(p1)),
upper = rep(NA, length(p1)),
dists = rep("tnorm", length(p1)),
untrans = rep("identity", length(p1)),
types = c("tnorm", "beta", "gamma", "lnorm", "unif", "constant", "tnorm2", NA)
)
Arguments
p1 |
the first parameter of a distribution |
p2 |
the second parameter of a distribution |
lower |
lower support (boundary) |
upper |
upper support (boundary) |
dists |
a vector of character string specifying a distribution. |
untrans |
whether to do log transformation. Default is not |
types |
available distribution types |
Details
Four distribution types are implemented:
Normal and truncated normal, where: p1 = mean, p2 = sd. It specifies a normal distribution when bounds are set -Inf and Inf,
Beta, where: p1 = shape1 and p2 = shape2 (see pbeta). Note the uniform distribution is a special case of the beta with p1 and p2 = 1),
Gamma, where p1 = shape and p2 = scale (see pgamma). Note p2 is scale, not rate,
Lognormal, where p1 = meanlog and p2 = sdlog (see plnorm).
Value
a list of list
Prepare posterior samples for plotting functions version 1
Description
Convert MCMC chains to a data frame for plotting functions
Usage
ConvertChains(x, start = 1, end = NA, pll = TRUE)
Arguments
x |
posterior samples |
start |
which iteration to start |
end |
end at which iteration |
pll |
a Boolean switch to make posterior log likelihood |
Deviance information criteria
Description
Calculate DIC and BPIC.
Usage
DIC(object, ...)
BPIC(object, ...)
Arguments
object |
posterior samples |
... |
other plotting arguments passing through dot dot dot. |
Get a n-cell matrix
Description
Constructs a matrix, showing how many responses to in each
cell. The function checks whether the format of n and ns
conform.
Usage
GetNsim(model, n, ns)
Arguments
model |
a model object. |
n |
number of trials. |
ns |
number of subjects. |
Details
n can be:
an integer for a balanced design,
a matrix for an unbalanced design, where rows are subjects and columns are cells. If the matrix is a row vector, all subjects have the same
nin each cell. If it is a column vector, all cells have the samen. Otherwise each entry specifies thenfor a particular subject x cell combination. See below for concrete examples.
Examples
model <- BuildModel(
p.map = list(A = "1", B = "R", t0 = "1", mean_v = "M", sd_v = "M",
st0 = "1"),
match.map = list(M = list(s1 = 1, s2 = 2)),
constants = c(sd_v.false = 1, st0 = 0),
factors = list(S = c("s1","s2")),
responses = c("r1", "r2"),
type = "norm")
#######################30
## Example 1
#######################30
GetNsim(model, ns = 2, n = 1)
# [,1] [,2]
# [1,] 1 1
# [2,] 1 1
#######################30
## Example 2
#######################30
n <- matrix(c(1:2), ncol = 1)
# [,1]
# [1,] 1 ## subject 1 has 1 response for each cell
# [2,] 2 ## subject 2 has 2 responses for each cell
GetNsim(model, ns = 2, n = n)
# [,1] [,2]
# [1,] 1 1
# [2,] 2 2
#######################30
## Example 3
#######################30
n <- matrix(c(1:2), nrow = 1)
# [,1] [,2]
# [1,] 1 2
GetNsim(model, ns = 2, n = n)
# [,1] [,2]
# [1,] 1 2 ## subject 1 has 1 response for cell 1 and 2 responses for cell 2
# [2,] 1 2 ## subject 2 has 1 response for cell 1 and 2 responses for cell 2
#######################30
## Example 4
#######################30
n <- matrix(c(1:4), nrow=2)
# [,1] [,2]
# [1,] 1 3
# [2,] 2 4
ggdmc::GetNsim(model, ns = 2, n = n)
# [,1] [,2]
# [1,] 1 3 ## subject 1 has 1 response for cell 1 and 3 responses for cell 2
# [2,] 2 4 ## subject 2 has 2 responses for cell 1 and 4 responses for cell 2
Extract parameter names from a model object
Description
Extract parameter names from a model object
Usage
GetPNames(x)
Arguments
x |
a model object |
Constructs a ns x npar matrix,
Description
The matrix is used to simulate data. Each row represents one set of parameters for a participant.
Usage
GetParameterMatrix(x, nsub, prior = NA, ps = NA, seed = NULL)
Arguments
x |
a model object |
nsub |
number of subjects. |
prior |
a prior object |
ps |
a vector or a matirx. |
seed |
an integer specifying a random seed. |
Details
One must enter either a vector or a matrix as true parameters
to the argument, ps, when presuming to simulate data based on a
fixed-effect model. When the assumption is to simulate data based on a
random-effect model, one must enter a prior object to the argument,
prior to first randomly generate a true parameter matrix.
Value
a ns x npar matrix
Examples
model <- BuildModel(
p.map = list(a ="1", v = "1",z = "1", d = "1", sz = "1", sv = "1",
t0 = "1", st0 = "1"),
match.map = list(M = list(s1 = "r1", s2 ="r2")),
factors = list(S = c("s1", "s2")),
constants = c(st0 = 0, d = 0),
responses = c("r1", "r2"),
type = "rd")
p.prior <- BuildPrior(
dists = c("tnorm", "tnorm", "beta", "beta", "tnorm", "beta"),
p1 = c(a = 1, v = 0, z = 1, sz = 1, sv = 1, t0 = 1),
p2 = c(a = 1, v = 2, z = 1, sz = 1, sv = 1, t0 = 1),
lower = c(0, -5, NA, NA, 0, NA),
upper = c(2, 5, NA, NA, 2, NA))
## Example 1: Randomly generate 2 sets of true parameters from
## parameter priors (p.prior)
GetParameterMatrix(model, 2, p.prior)
## a v z sz sv t0
## [1,] 1.963067 1.472940 0.9509158 0.5145047 1.344705 0.0850591
## [2,] 1.512276 -1.995631 0.6981290 0.2626882 1.867853 0.1552828
## Example 2: Use a user-selected true parameters
true.vector <- c(a=1, v=1, z=0.5, sz=0.2, sv=1, t0=.15)
GetParameterMatrix(model, 2, NA, true.vector)
## a v z sz sv t0
## [1,] 1 1 0.5 0.2 1 0.15
## [2,] 1 1 0.5 0.2 1 0.15
GetParameterMatrix(model, 2, ps = true.vector)
## Example 3: When a user enter arbritary sequence of parameters.
## Note sv is before sz. It should be sz before sv
## See correct sequence, by entering "attr(model, 'p.vector')"
## GetParameterMatrix will rearrange the sequence.
true.vector <- c(a=1, v=1, z=0.5, sv=1, sz = .2, t0=.15)
GetParameterMatrix(model, 2, NA, true.vector)
## a v z sz sv t0
## [1,] 1 1 0.5 0.2 1 0.15
## [2,] 1 1 0.5 0.2 1 0.15
Which chains get stuck
Description
Calculate each chain separately for the mean (across many MCMC iterations)
of posterior log-likelihood. If the difference of the means and
the median (across chains) of the mean of posterior is greater than the
cut, chains are considered stuck. The default value for cut
is 10. unstick manually removes stuck chains from posterior samples.
Usage
PickStuck(
x,
hyper = FALSE,
cut = 10,
start = 1,
end = NA,
verbose = FALSE,
digits = 2
)
Arguments
x |
posterior samples |
hyper |
whether x are hierarhcial samples |
cut |
a criterion deciding if a chain is stuck. |
start |
start to evaluate from which iteration. |
end |
end at which iteration for evaeuation. |
verbose |
a boolean switch to print more information |
digits |
print how many digits. Default is 2 |
Value
PickStuck gives an index vector; unstick gives a DMC
sample.
Examples
model <- BuildModel(
p.map = list(A = "1", B = "1", t0 = "1", mean_v = "M", sd_v = "1", st0 = "1"),
match.map = list(M = list(s1 = 1, s2 = 2)),
factors = list(S = c("s1", "s2")),
constants = c(st0 = 0, sd_v = 1),
responses = c("r1", "r2"),
type = "norm")
p.vector <- c(A = .75, B = .25, t0 = .2, mean_v.true = 2.5, mean_v.false = 1.5)
p.prior <- BuildPrior(
dists = c("tnorm", "tnorm", "beta", "tnorm", "tnorm"),
p1 = c(A = .3, B = .3, t0 = 1, mean_v.true = 1, mean_v.false = 0),
p2 = c(1, 1, 1, 3, 3),
lower = c(0, 0, 0, NA, NA),
upper = c(NA,NA, 1, NA, NA))
## Not run:
dat <- simulate(model, 30, ps = p.vector)
dmi <- BuildDMI(dat, model)
sam <- run(StartNewsamples(dmi, p.prior))
bad <- PickStuck(sam)
## End(Not run)
Start new model fits
Description
Fit a hierarchical or a fixed-effect model, using Bayeisan optimisation. We use a specific type of pMCMC algorithm, the DE-MCMC. This particular sampling method includes crossover and two different migration operators. The migration operators are similar to random-walk algorithm. They wouold be less efficient to find the target parameter space, if been used alone.
Usage
StartNewsamples(
data,
prior = NULL,
nmc = 200,
thin = 1,
nchain = NULL,
report = 100,
rp = 0.001,
gammamult = 2.38,
pm_Hu = 0.05,
pm_BT = 0.05,
block = TRUE,
ncore = 1,
add = FALSE,
is_old = FALSE
)
run(
samples,
nmc = 500,
thin = 1,
report = 100,
rp = 0.001,
gammamult = 2.38,
pm_Hu = 0,
pm_BT = 0,
block = TRUE,
ncore = 1,
add = FALSE,
is_old = TRUE
)
Arguments
data |
data model instance(s) |
prior |
prior objects. For hierarchical model, this must be a list with three sets of prior distributions. Each is respectively named, "pprior", "location", and "scale". |
nmc |
number of Monte Carlo samples |
thin |
thinning length |
nchain |
number of chains |
report |
progress report interval |
rp |
tuning parameter 1 |
gammamult |
tuning parameter 2. This is the step size. |
pm_Hu |
probability of migration type 0 (Hu & Tsui, 2010) |
pm_BT |
probability of migration type 1 (Turner et al., 2013) |
block |
Only for hierarchical modeling. A Boolean switch for update one parameter at a time |
ncore |
Only for non-hierarchical, fixed-effect models with many subjects. |
add |
Boolean whether to add new samples |
is_old |
start sampling from a DMI or fit samples |
samples |
posterior samples. |
Table response and parameter
Description
TableParameters arranges the values in a parameter
vector and creates a response x parameter matrix. The matrix is used
by the likelihood function, assigning a trial to a cell for calculating
probability densities.
Usage
TableParameters(p.vector, cell, model, n1order)
Arguments
p.vector |
a parameter vector |
cell |
a string or an integer indicating a design cell, e.g.,
|
model |
a model object |
n1order |
a Boolean switch, indicating using node 1 ordering. This is only for LBA-like models and its n1PDF likelihood function. |
Value
each row corresponding to the model parameter for a response.
When n1.order is FALSE, TableParameters returns a martix without
rearranging into node 1 order. For example, this is used in
the simulate function. By default n1.order is TRUE.
Examples
m1 <- BuildModel(
p.map = list(a = "1", v = "F", z = "1", d = "1", sz = "1", sv = "F",
t0 = "1", st0 = "1"),
match.map = list(M = list(s1 = "r1", s2 = "r2")),
factors = list(S = c("s1", "s2"), F = c("f1","f2")),
constants = c(st0 = 0, d = 0),
responses = c("r1","r2"),
type = "rd")
m2 <- BuildModel(
p.map = list(A = "1", B = "1", mean_v = "M", sd_v = "1",
t0 = "1", st0 = "1"),
constants = c(st0 = 0, sd_v = 1),
match.map = list(M = list(s1 = 1, s2 = 2)),
factors = list(S = c("s1", "s2")),
responses = c("r1", "r2"),
type = "norm")
pvec1 <- c(a = 1.15, v.f1 = -0.10, v.f2 = 3, z = 0.74, sz = 1.23,
sv.f1 = 0.11, sv.f2 = 0.21, t0 = 0.87)
pvec2 <- c(A = .75, B = .25, mean_v.true = 2.5, mean_v.false = 1.5,
t0 = .2)
print(m1, pvec1)
print(m2, pvec2)
accMat1 <- TableParameters(pvec1, "s1.f1.r1", m1, FALSE)
accMat2 <- TableParameters(pvec2, "s1.r1", m2, FALSE)
## a v t0 z d sz sv st0
## 1.15 -0.1 0.87 0.26 0 1.23 0.11 0
## 1.15 -0.1 0.87 0.26 0 1.23 0.11 0
## A b t0 mean_v sd_v st0
## 0.75 1 0.2 2.5 1 0
## 0.75 1 0.2 1.5 1 0
Does a model object specify a correct p.vector
Description
Check a parameter vector
Usage
check_pvec(p.vector, model)
Arguments
p.vector |
parameter vector |
model |
a model object |
A modified dbeta function
Description
A modified dbeta function
Usage
dbeta_lu(x, p1, p2, lower, upper, lg = FALSE)
Arguments
x |
quantile |
p1 |
shape1 parameter |
p2 |
shape2 parameter |
lower |
lower bound |
upper |
upper bound |
lg |
logical; if TRUE, return log density. |
A modified dcauchy functions
Description
A modified dcauchy functions
Usage
dcauchy_l(x, p1, p2, lg = FALSE)
Arguments
x |
quantile |
p1 |
location parameter |
p2 |
scale parameter |
lg |
log density? |
A pseudo constant function to get constant densities
Description
Used with constant prior
Usage
dconstant(x, p1, p2, lower, upper, lg = FALSE)
Arguments
x |
quantile |
p1 |
constant value |
p2 |
unused argument |
lower |
dummy varlable |
upper |
dummy varlable |
lg |
log density? |
Calculate the statistics of model complexity
Description
Calculate deviance for a model object for which a log-likelihood value can be obtained, according to the formula -2*log-likelihood.
Usage
## S3 method for class 'model'
deviance(object, ...)
Arguments
object |
posterior samples |
... |
other plotting arguments passing through dot dot dot. |
A modified dgamma function
Description
A modified dgamma function
Usage
dgamma_l(x, p1, p2, lower, upper, lg = FALSE)
Arguments
x |
quantile |
p1 |
shape parameter |
p2 |
scale parameter |
lower |
lower bound |
upper |
upper bound |
lg |
log density? |
A modified dlnorm functions
Description
A modified dlnorm functions
Usage
dlnorm_l(x, p1, p2, lower, upper, lg = FALSE)
Arguments
x |
quantile |
p1 |
meanlog parameter |
p2 |
sdlog parameter |
lower |
lower bound |
upper |
upper bound |
lg |
log density? |
Truncated Normal Distribution
Description
Random number generation, probability density and cumulative density functions for truncated normal distribution.
Usage
dtnorm(x, p1, p2, lower, upper, lg = FALSE)
rtnorm(n, p1, p2, lower, upper)
ptnorm(q, p1, p2, lower, upper, lt = TRUE, lg = FALSE)
Arguments
x, q |
vector of quantiles; |
p1 |
mean (must be scalar). |
p2 |
standard deviation (must be scalar). |
lower |
lower truncation value (must be scalar). |
upper |
upper truncation value (must be scalar). |
lg |
log probability. If TRUE (default is FALSE) probabilities p are
given as |
n |
number of observations. n must be a scalar. |
lt |
lower tail. If TRUE (default) probabilities are |
Value
a column vector.
Examples
## rtn example
dat1 <- rtnorm(1e5, 0, 1, 0, Inf)
hist(dat1, breaks = "fd", freq = FALSE, xlab = "",
main = "Truncated normal distributions")
## dtn example
x <- seq(-5, 5, length.out = 1e3)
dat1 <- dtnorm(x, 0, 1, -2, 2, 0)
plot(x, dat1, type = "l", lwd = 2, xlab = "", ylab= "Density",
main = "Truncated normal distributions")
## ptn example
x <- seq(-10, 10, length.out = 1e2)
mean <- 0
sd <- 1
lower <- 0
upper <- 5
dat1 <- ptnorm(x, 0, 1, 0, 5, lg = TRUE)
Calculate effective sample sizes
Description
effectiveSize calls effectiveSize in coda package to
calculate sample sizes.
Usage
effectiveSize_hyper(x, start, end, digits, verbose)
effectiveSize_many(x, start, end, verbose)
effectiveSize_one(x, start, end, digits, verbose)
effectiveSize(
x,
hyper = FALSE,
start = 1,
end = NA,
digits = 0,
verbose = FALSE
)
Arguments
x |
posterior samples |
start |
starting iteration |
end |
ending iteraton |
digits |
printing how many digits |
verbose |
printing more information |
hyper |
a Boolean switch to extract hyper attribute |
Examples
#################################40
## effectiveSize example
#################################40
## Not run:
es1 <- effectiveSize_one(hsam[[1]], 1, 100, 2, TRUE)
es2 <- effectiveSize_one(hsam[[1]], 1, 100, 2, FALSE)
es3 <- effectiveSize_many(hsam, 1, 100, TRUE)
es4 <- effectiveSize_many(hsam, 1, 100, FALSE)
es5 <- effectiveSize_hyper(hsam, 1, 100, 2, TRUE)
es6 <- effectiveSize(hsam, TRUE, 1, 100, 2, TRUE)
es7 <- effectiveSize(hsam, TRUE, 1, 100, 2, FALSE)
es8 <- effectiveSize(hsam, FALSE, 1, 100, 2, TRUE)
es9 <- effectiveSize(hsam, FALSE, 1, 100, 2, FALSE)
es10 <- effectiveSize(hsam[[1]], FALSE, 1, 100, 2, TRUE)
## End(Not run)
Potential scale reduction factor
Description
gelman function calls the function, gelman.diag in the
coda package to calculates PSRF.
Usage
gelman(
x,
hyper = FALSE,
start = 1,
end = NA,
confidence = 0.95,
transform = TRUE,
autoburnin = FALSE,
multivariate = TRUE,
split = TRUE,
subchain = FALSE,
nsubchain = 3,
digits = 2,
verbose = FALSE,
...
)
hgelman(
x,
start = 1,
end = NA,
confidence = 0.95,
transform = TRUE,
autoburnin = FALSE,
split = TRUE,
subchain = FALSE,
nsubchain = 3,
digits = 2,
verbose = FALSE,
...
)
Arguments
x |
posterior samples |
hyper |
a Boolean switch, indicating posterior samples are from hierarchical modeling |
start |
start iteration |
end |
end iteration |
confidence |
confident inteval |
transform |
turn on transform |
autoburnin |
turn on auto burnin |
multivariate |
multivariate Boolean switch |
split |
split whether split mcmc chains; When split is TRUE, the function doubles the number of chains by spliting into 1st and 2nd halves. |
subchain |
whether only calculate a subset of chains |
nsubchain |
indicate how many chains in a subset |
digits |
print out how many digits |
verbose |
print more information |
... |
arguments passing to |
Examples
## Not run:
rhat1 <- hgelman(hsam); rhat1
rhat2 <- hgelman(hsam, end = 51); rhat2
rhat3 <- hgelman(hsam, confidence = .90); rhat3
rhat4 <- hgelman(hsam, transform = FALSE); rhat4
rhat5 <- hgelman(hsam, autoburnin = TRUE); rhat5
rhat6 <- hgelman(hsam, split = FALSE); rhat6
rhat7 <- hgelman(hsam, subchain = TRUE); rhat7
rhat8 <- hgelman(hsam, subchain = TRUE, nsubchain = 4);
rhat9 <- hgelman(hsam, subchain = TRUE, nsubchain = 4,
digits = 1, verbose = TRUE);
hat1 <- gelman(hsam[[1]], multivariate = FALSE); hat1
hat2 <- gelman(hsam[[1]], hyper = TRUE, verbose = TRUE); hat2
hat3 <- gelman(hsam, hyper = TRUE, verbose = TRUE); hat3
hat4 <- gelman(hsam, multivariate = TRUE, verbose = FALSE);
hat5 <- gelman(hsam, multivariate = FALSE, verbose = FALSE);
hat6 <- gelman(hsam, multivariate = FALSE, verbose = TRUE);
hat7 <- gelman(hsam, multivariate = T, verbose = TRUE);
## End(Not run)
Retrieve information of operating system
Description
A wrapper function to extract system information from Sys.info
and .Platform
Usage
get_os()
Examples
get_os()
## sysname
## "linux"
Bayeisan computation of response time models
Description
ggdmc uses the population-based Markov chain Monte Carlo to conduct Bayesian computation on cognitive models.
Author(s)
Yi-Shin Lin <yishinlin001@gmail.com>
Andrew Heathcote <andrew.heathcote@utas.edu.au>
References
Heathcote, A., Lin, Y.-S., Reynolds, A., Strickland, L., Gretton, M. &
Matzke, D., (2018). Dynamic model of choice.
Behavior Research Methods.
https://doi.org/10.3758/s13428-018-1067-y.
Turner, B. M., & Sederberg P. B. (2012). Approximate Bayesian computation
with differential evolution, Journal of Mathematical Psychology, 56,
375–385.
Ter Braak (2006). A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces. Statistics and Computing, 16, 239-249.
Model checking functions
Description
The function tests whether we have drawn enough samples.
Usage
iseffective(x, minN, nfun, verbose = FALSE)
Arguments
x |
posterior samples |
minN |
specify the size of minimal effective samples |
nfun |
specify to use the |
verbose |
print more information |
Model checking functions
Description
The function tests whether Markov chains converge prematurelly:
Usage
isflat(
x,
p1 = 1/3,
p2 = 1/3,
cut_location = 0.25,
cut_scale = Inf,
verbose = FALSE
)
Arguments
x |
posterior samples |
p1 |
the range of the head of MCMC chains |
p2 |
the range of the tail of the MCMC chains |
cut_location |
how far away a location chains been considered as stuck |
cut_scale |
how far away a scale chains been considered as stuck |
verbose |
print more information |
Model checking functions
Description
The function tests whether Markov chains are mixed well.
Usage
ismixed(x, cut = 1.01, split = TRUE, verbose = FALSE)
Arguments
x |
posterior samples |
cut |
psrf criterion for well mixed |
split |
whether to split MCMC chains. This is an argument passing to
|
verbose |
print more information |
See Also
Model checking functions
Description
The function tests whether Markov chains encounter a parameter
region that is difficult to search. CheckConverged is
a wrapper function running the four checking functions,
isstuck, isflat, ismixed and iseffective.
Usage
isstuck(x, hyper = FALSE, cut = 10, start = 1, end = NA, verbose = FALSE)
CheckConverged(x)
Arguments
x |
posterior samples |
hyper |
a Boolean switch, extracting hyper attribute. |
cut |
the criteria for suggesting abnormal chains found |
start |
start iteration |
end |
end iteration |
verbose |
print more information |
Calculate likelihoods
Description
These function calculate likelihoods. likelihood_rd implements
the equations in Voss, Rothermund, and Voss (2004). These equations
calculate diffusion decision model (Ratcliff & Mckoon, 2008). Specifically,
this function implements Voss, Rothermund, and Voss's (2004) equations A1
to A4 (page 1217) in C++.
Usage
likelihood(p_vector_r, data, min_lik = 1e-10)
Arguments
p_vector_r |
a parameter vector |
data |
data model instance |
min_lik |
minimal likelihood. |
Value
a vector
References
Voss, A., Rothermund, K., & Voss, J. (2004). Interpreting the
parameters of the diffusion model: An empirical validation.
Memory & Cognition, 32(7), 1206-1220.
Ratcliff, R. (1978). A theory of memory retrival. Psychological
Review, 85, 238-255.
Examples
model <- BuildModel(
p.map = list(A = "1", B = "1", t0 = "1", mean_v = "M", sd_v = "1",
st0 = "1"),
match.map = list(M = list(s1 = 1, s2 = 2)),
factors = list(S = c("s1", "s2")),
constants = c(st0 = 0, sd_v = 1),
responses = c("r1", "r2"),
type = "norm")
p.vector <- c(A = .25, B = .35, t0 = .2, mean_v.true = 1,
mean_v.false = .25)
dat <- simulate(model, 1e3, ps = p.vector)
dmi <- BuildDMI(dat, model)
den <- likelihood(p.vector, dmi)
model <- BuildModel(
p.map = list(a = "1", v = "1", z = "1", d = "1", t0 = "1", sv = "1",
sz = "1", st0 = "1"),
constants = c(st0 = 0, d = 0),
match.map = list(M = list(s1 = "r1", s2 = "r2")),
factors = list(S = c("s1", "s2")),
responses = c("r1", "r2"),
type = "rd")
p.vector <- c(a = 1, v = 1, z = 0.5, sz = 0.25, sv = 0.2, t0 = .15)
dat <- simulate(model, 1e2, ps = p.vector)
dmi <- BuildDMI(dat, model)
den <- likelihood (p.vector, dmi)
Create a MCMC list
Description
Create a MCMC list
Usage
mcmc_list.model(x, start = 1, end = NA, pll = TRUE)
Arguments
x |
posterior samples |
start |
start from which iteration |
end |
end at which iteration |
pll |
a Boolean switch for calculating posterior log-likelihood |
Plot prior distributions
Description
plot_prior plots one member in a prior object. plot.prior
plots all members in a prior object.
Usage
plot_prior(
i,
prior,
xlim = NA,
natural = TRUE,
npoint = 100,
trans = NA,
save = FALSE,
...
)
## S3 method for class 'prior'
plot(x, save = FALSE, ps = NULL, ...)
Arguments
i |
an integer or a character string indicating which parameter to plot |
prior |
a prior object |
xlim |
set the range of on x axis. This is usually the range for each parameter. |
natural |
default TRUE. |
npoint |
default to plot 100 |
trans |
default NA. trans can be a scalar or vector. |
save |
whether to save the data out |
... |
other plotting arguments passing through dot dot dot. |
x |
a prior object |
ps |
true parameter vectors or matrix in the case of many observation units |
Examples
p.prior <- BuildPrior(
dists = rep("tnorm", 7),
p1 = c(a = 2, v.f1 = 4, v.f2 = 3, z = 0.5, sv = 1,
sz = 0.3, t0 = 0.3),
p2 = c(a = 0.5, v.f1 = .5, v.f2 = .5, z = 0.1, sv = .3,
sz = 0.1, t0 = 0.05),
lower = c(0,-5, -5, 0, 0, 0, 0),
upper = c(5, 7, 7, 1, 2, 1, 1))
plot_prior("a", p.prior)
plot_prior(2, p.prior)
plot(p.prior)
Print Prior Distribution
Description
a convenient function to rearrange p.prior or an element in a
pp.prior as a data frame for inspection.
Usage
## S3 method for class 'prior'
print(x, ...)
Arguments
x |
a list of prior distributions list, usually created by
|
... |
other arguments |
Value
a data frame listing prior distributions and their settings
Examples
pop.mean <- c(a=1, v.f1=1, v.f2=.2, z=.5, sz=.3, sv.f1=.25, sv.f2=.23,
t0=.3)
pop.scale <- c(a=.2, v.f1=.2, v.f2=.2, z=.1, sz=.05, sv.f1=.05, sv.f2=.05,
t0=.05)
p.prior <- BuildPrior(
dists = rep("tnorm", 8),
p1 = pop.mean,
p2 = pop.scale,
lower = c(0, -5, -5, 0, 0, 0, 0, 0),
upper = c(2, 5, 5, 1, 2, 2, 1, 1))
print(p.prior)
Generate random numbers
Description
A wrapper function for generating random numbers of either
the model type, rd, or norm.
Usage
random(type, pmat, n, seed = NULL)
Arguments
type |
a character string of the model type |
pmat |
a matrix of response x parameter |
n |
number of observations |
seed |
an integer specifying a random seed |
Generate Random Responses of the LBA Distribution
Description
rlba_norm used the LBA process to generate response times and
responses.
Usage
rlba_norm(n, A, b, mean_v, sd_v, t0, st0, posdrift)
Arguments
n |
is the numbers of observation. |
A |
start point upper bound, a vector of a scalar. |
b |
decision threshold, a vector or a scalar. |
mean_v |
mean drift rate vector |
sd_v |
standard deviation of drift rate vector |
t0 |
non-decision time, a vector. |
st0 |
non-decision time variation, a vector. |
posdrift |
if exclude negative drift rates |
Value
a n x 2 matrix of response time (first column) and responses (second column).
Parameter Prior Distributions
Description
Probability density functions and random generation for parameter prior distributions.
Usage
rprior(prior, n = 1)
Arguments
prior |
a list of list usually created by BuildPrior to store the information about parameter prior distributions. |
n |
number of observations/random draws |
Examples
p.prior <- BuildPrior(
dists = c("tnorm", "tnorm", "beta", "tnorm", "beta", "beta"),
p1 = c(a = 1, v = 0, z = 1, sz = 1, sv = 1, t0 = 1),
p2 = c(a = 1, v = 2, z = 1, sz = 1, sv = 1, t0 = 1),
lower = c(0,-5, NA, NA, 0, NA),
upper = c(2, 5, NA, NA, 2, NA))
rprior(p.prior, 9)
## a v z sz sv t0
## [1,] 0.97413686 0.78446178 0.9975199 -0.5264946 0.5364492 0.55415052
## [2,] 0.72870190 0.97151662 0.8516604 1.6008591 0.3399731 0.96520848
## [3,] 1.63153685 1.96586939 0.9260939 0.7041254 0.4138329 0.78367440
## [4,] 1.55866180 1.43657110 0.6152371 0.1290078 0.2957604 0.23027759
## [5,] 1.32520281 -0.07328408 0.2051155 2.4040387 0.9663111 0.06127237
## [6,] 0.49628528 -0.19374770 0.5142829 2.1452972 0.4335482 0.38410626
## [7,] 0.03655549 0.77223432 0.1739831 1.4431507 0.6257398 0.63228368
## [8,] 0.71197612 -1.15798082 0.8265523 0.3813370 0.4465184 0.23955415
## [9,] 0.38049166 3.32132034 0.9888108 0.9684292 0.8437480 0.13502154
pvec <- c(a=1, v=1, z=0.5, sz=0.25, sv=0.2,t0=.15)
p.prior <- BuildPrior(
dists = rep("tnorm", 6),
p1 = c(a=2, v=2.5, z=0.5, sz=0.3, sv=1, t0=0.3),
p2 = c(a=0.5, v=.5, z=0.1, sz=0.1, sv=.3, t0=0.05) * 5,
lower = c(0,-5, 0, 0, 0, 0),
upper = c(5, 7, 2, 2, 2, 2))
Simulate response time data
Description
Simulate response time data either for one subject or multiple subjects.
The simulation is based on a model object. For one subject, one must supply
a true parameter vector to the ps argument.
Usage
## S3 method for class 'model'
simulate(object, nsim = NA, seed = NULL, nsub = NA, prior = NA, ps = NA, ...)
Arguments
object |
a model object. |
nsim |
number of trials / responses. |
seed |
a user specified random seed. |
nsub |
number of subjects |
prior |
a prior object |
ps |
a true parameter vector or matrix. |
... |
additional optional arguments. |
Details
For multiple subjects, one can enter a matrix (or a row vector) as true
parameters. Each row is to generate data separately for a subject. This is
the fixed-effect model. To generate data based on a random-effect
model, one must supply a prior object. In this case, ps argument
is unused. Note in some cases, a random-effect model may fail to draw data
from the model, because true parameters are randomly drawn from
a prior object. This would happen sometimes in diffusion model, because
certain parameter combinations are considered invalid.
ps can be a row vector, in which case each subject has identical
parameters. It can also be a matrix with one row per subject, in which
case it must have ns rows. The true values will be saved as
parameters attribute in the output object.
Value
a data frame
Summarise posterior samples
Description
This calls seven different variants of summary function to summarise posterior samples
Usage
## S3 method for class 'model'
summary(
object,
hyper = FALSE,
start = 1,
end = NA,
hmeans = FALSE,
hci = FALSE,
prob = c(0.025, 0.25, 0.5, 0.75, 0.975),
recovery = FALSE,
ps = NA,
type = 1,
verbose = FALSE,
digits = 2,
...
)
Arguments
object |
posterior samples |
hyper |
whether to summarise hyper parameters |
start |
start from which iteration. |
end |
end at which iteration. For example, set
|
hmeans |
a boolean switch indicating to calculate mean of hyper parameters |
hci |
boolean switch; whether to calculate credible intervals of hyper parameters |
prob |
a numeric vector, indicating the quantiles to calculate |
recovery |
a boolean switch indicating if samples are from a recovery study |
ps |
true parameter values. This is only for recovery studies |
type |
calculate type 1 or 2 hyper parameters |
verbose |
print more information |
digits |
printing digits |
... |
other arguments |
Examples
## Not run:
est1 <- summary(hsam[[1]], FALSE)
est2 <- summary(hsam[[1]], FALSE, 1, 100)
est3 <- summary(hsam)
est4 <- summary(hsam, verbose = TRUE)
est5 <- summary(hsam, verbose = FALSE)
hest1 <- summary(hsam, TRUE)
## End(Not run)
Summary statistic for posterior samples
Description
Calculate summary statistics for posterior samples
Usage
summary_mcmc_list(object, prob = c(0.025, 0.25, 0.5, 0.75, 0.975), ...)
Arguments
object |
posterior samples |
prob |
summary quantile summary |
... |
other arguments passing in |
Convert theta to a mcmc List
Description
Extracts the parameter array (ie theta) from posterior samples of a partiipant and convert it to a coda mcmc.list.
Usage
theta2mcmclist(
x,
start = 1,
end = NA,
split = FALSE,
subchain = FALSE,
nsubchain = 3,
thin = NA
)
phi2mcmclist(
x,
start = 1,
end = NA,
split = FALSE,
subchain = FALSE,
nsubchain = 3
)
Arguments
x |
posterior samples |
start |
start iteration |
end |
end iteraton |
split |
whether to divide one MCMC sequence into two sequences. |
subchain |
boolean swith convert only a subset of chains |
nsubchain |
indicate the number of chains in the subset |
thin |
thinning lenght of the posterior samples |
Details
phi2mcmclist extracts the phi parameter array, which stores
the location and scale parameters at the hyper level.
Examples
## Not run:
model <- BuildModel(
p.map = list(a = "RACE", v = c("S", "RACE"), z = "RACE", d = "1",
sz = "1", sv = "1", t0 = c("S", "RACE"), st0 = "1"),
match.map = list(M = list(gun = "shoot", non = "not")),
factors = list(S = c("gun", "non"), RACE = c("black", "white")),
constants = c(st0 = 0, d = 0, sz = 0, sv = 0),
responses = c("shoot", "not"),
type = "rd")
pnames <- GetPNames(model)
npar <- length(pnames)
pop.mean <- c(1, 1, 2.5, 2.5, 2.5, 2.5, .50, .50, .4, .4, .4, .4)
pop.scale <- c(.15, .15, 1, 1, 1, 1, .05, .05, .05, .05, .05, .05)
names(pop.mean) <- pnames
names(pop.scale) <- pnames
pop.prior <- BuildPrior(
dists = rep("tnorm", npar),
p1 = pop.mean,
p2 = pop.scale,
lower = c(rep(0, 2), rep(-5, 4), rep(0, 6)),
upper = c(rep(5, 2), rep(7, 4), rep(2, 6)))
p.prior <- BuildPrior(
dists = rep("tnorm", npar),
p1 = pop.mean,
p2 = pop.scale*10,
lower = c(rep(0, 2), rep(-5, 4), rep(0, 6)),
upper = c(rep(10, 2), rep(NA, 4), rep(5, 6)))
mu.prior <- BuildPrior(
dists = rep("tnorm", npar),
p1 = pop.mean,
p2 = pop.scale*10,
lower = c(rep(0, 2), rep(-5, 4), rep(0, 6)),
upper = c(rep(10, 2), rep(NA, 4), rep(5, 6)))
sigma.prior <- BuildPrior(
dists = rep("beta", npar),
p1 = rep(1, npar),
p2 = rep(1, npar),
upper = rep(2, npar))
names(sigma.prior) <- GetPNames(model)
priors <- list(pprior=p.prior, location=mu.prior, scale=sigma.prior)
dat <- simulate(model, nsim = 10, nsub = 10, prior = pop.prior)
dmi <- BuildDMI(dat, model)
ps <- attr(dat, "parameters")
fit0 <- StartNewsamples(dmi, priors)
fit <- run(fit0)
tmp1 <- theta2mcmclist(fit[[1]])
tmp2 <- theta2mcmclist(fit[[2]], start = 10, end = 90)
tmp3 <- theta2mcmclist(fit[[3]], split = TRUE)
tmp4 <- theta2mcmclist(fit[[4]], subchain = TRUE)
tmp5 <- theta2mcmclist(fit[[5]], subchain = TRUE, nsubchain = 4)
tmp6 <- theta2mcmclist(fit[[6]], thin = 2)
## End(Not run)
Unstick posterios samples (One subject)
Description
Unstick posterios samples (One subject)
Usage
unstick_one(x, bad)
Arguments
x |
posterior samples |
bad |
a numeric vector, indicating which chains to remove |