Models of Cognition for Subjective Scales and Decision Making Tasks in R
This R package is dedicated to facilitate the application of computational models of cognition in R under a Bayesian framework. These are useful in the field of cognitive science and computational neuropsycholology.
If you have suggestions for improvement, please get in touch!
mu and sigma index the Gaussian component alone and tau the
exponential tail - unlike brms’s native exgaussian(), whose mu
indexes the mean of the entire distribution)
Measures from cognitive tasks, such as decision-making paradigms involving fast responses or ratings, often produce noisy, specific, and complex patterns of results. Broadly speaking, there are three ways of analysing such data.
Install the released version from CRAN:
install.packages("cogmod")
Or the development version from GitHub:
if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes")
# remotes::install_github("DominiqueMakowski/cogmod")
install.packages("cogmod", repos = c('https://dominiquemakowski.r-universe.dev', 'https://cloud.r-project.org'))
Fitting the models requires a working Stan installation through
brms. The examples use the
cmdstanr backend, which is not on CRAN: follow its installation
guide
to set up CmdStan.
Using a cogmod model requires two arguments beyond a standard brm()
call: family and stanvars. Two further helpers, cogmod_priors()
and cogmod_inits(), are strictly speaking optional but should be
treated as part of the call: they provide adapted chain-initialization
values and weakly informative priors on the sensitive parameters to
limit convergence issues and other sampling pathologies.
library(cogmod)
library(brms)
# Specify the formula using brms' bf()
f <- bf(
RT ~ Condition,
sigma ~ Condition,
ndt ~ Condition,
sigmabias = 0, # No start-point variability
family = cogmod_lognormal()
)
# Fit the model
m <- brm(
f,
data = df,
prior = cogmod_priors(f, df),
init = cogmod_inits(f, df),
stanvars = cogmod_stanvars(f),
backend = "cmdstanr"
)
We can then analyze its results, and check its predictions like with any other models. See the Subjective Ratings, RT-only Models, and Decision Making Models vignettes for more detailed examples.
log(RT - ndt) ~ Logistic. The motivation is
the hazard function: Weibull and Gamma hazards are monotone, whereas
the log-logistic rises then falls, a shape none of the current
families can produce.poutlier currently mixes in a half-Normal at zero with a fixed
scale: a contaminant with no mechanism, there so that one stray fast
response cannot drag ndt down. LATER handles the same responses with
a second accumulator - an “early” or “maverick” unit with a mean
rate near zero and a large SD - that races the main one, so an
early response is a decision that happened to win rather than a trial
the model disowns. That is the better account, and it is what the
second line on a reciprobit plot actually is. It is also a large
undertaking: a race needs the winner’s density against the loser’s
survival for every family in .OUTLIER_FAMILIES, where the present
mixture needs one scalar, and the early unit is deliberately
half-defective (with a mean rate of zero, half its trials never arrive
at all) - which our normalised densities cannot represent unless the
race is there to absorb the missing mass. Put on hold until we find a
design that does not cost every family its closed form.cogmod_loggamma() ties skew and tail weight together through one
shape, SHASH decouples them.cogmod_priors() and cogmod_inits() for
subjective scale models.cogmod_lognormal()
and cogmod_lnr() with sigmabias > 0, since the shifted LogNormal
and the LNR are that LBA with the start point removed, and adding
the start-point range back was one dpar on each rather than new
families. Two lessons from it carry over. The scaling constraint
moves with the drift distribution - fixing the drift scale
(sigmazero = 1, the LBA convention) pins the Normal and the Gamma,
but a LogNormal rate’s sigma is untouched by rescaling, so the
threshold has to be pinned instead, which is why the LNR’s offset
is fixed at 1 - and Stan’s upper-tail normal log-CDF is not accurate
enough for the fast tail (use erfc).sigmabias = 0 yields Luce / multinomial-logit choice
probabilities, bridging RT models and random-utility choice models,
but its heavy drift tail produces too many very fast responses and
its mean is undefined for shape ≤ 1. The four predicted CDFs are
nearly indistinguishable and differ only in the tails, so the drift
distribution is a modelling choice to justify theoretically and to
compare by LOO or Bayes factor, not by raw likelihood. Validate
against rtdists::dLBA(), check the sigmabias = 0 limits reduce
to the named distributions (inverse-Gamma and Weibull), and run a
cross-fitting matrix at realistic per-participant trial counts (the
published comparison used 20,000).z and the drift v -
x(t) = a(t) z + g(t) v + c(t) - crosses the threshold exactly when
the plain LBA does with an effective time g(t), an effective
start-point range a(t) * sigmabias and an effective threshold
b - c(t), so its CDF is one wrapper around the LBA core (density by
the chain rule) and it inherits every drift distribution above for
free. Three choices of (a, g, c) are worth having, in this order.
Time warp (x = z + v g(t),
e.g. g(t) = (exp(gamma t) - 1) / gamma, series-expanded near
gamma = 0): a time-dependent gain on evidence, the ballistic
counterpart of urgency signals (Hawkins et al.,
2015); it preserves
the finishing-time order on every trial, so choice probabilities are
exactly the LBA’s and only the RT shape changes - a built-in unit
test. A collapsing threshold is the additive cousin
(c(t) = b - b(t)) and does move accuracy. Ballistic
Ornstein-Uhlenbeck (dx/dt = v + gamma x, so a = exp(gamma t)):
the Leaky Competing Accumulator (Usher & McClelland,
2001) without noise or
inhibition; gamma < 0 is leak, gamma > 0 self-excitation (Bogacz
et al., 2006). The start
point is amplified or forgotten over time, an accumulator finishes
only if z + v / gamma > 0, so bias and drift interact and accuracy
changes; and with piecewise-constant input the integral stays
analytic, giving primacy (gamma > 0) or recency (gamma < 0) - the
LCA’s signature phenomenon with a closed-form likelihood.
Accumulator-specific gamma is cheap; ballistic lateral inhibition
keeps affinity but makes the winner’s density a bivariate integral
over the start-point box; a floor at zero breaks affinity and is out
of scope. Literal gravity (dx/dt = v + gamma / (b - x)^2) is not
affine in v, needs a one-dimensional quadrature and loses the
composability, so only if a story specifically needs
distance-dependent attraction. Expect gamma to trade off against
threshold and drift scale, since curvature acts on the tail exactly
where the drift distributions differ; recover before believing
anything.T* on how long the
decision should take, and trajectories heading to finish early are
slowed while late ones are sped up, so the attractor is the point
(T*, b) on the threshold rather than the threshold line -
bidirectional bending that gravity cannot produce. In the
linear-Gaussian case predictive coding is precision weighting, so the
effective rate is the sensory drift shrunk toward the rate that
reaches threshold at T*: v_eff = (1 - w) v + w v*, with w the
relative prior precision. Two readings of v* map onto the two
mechanisms above: a prior over finishing time, v* = (b - z) / T*,
is a pure time warp with g(t) = (1 - w) t / (1 - w t / T*) and
leaves accuracy untouched; a prior over rate, v* = b / T*, is an
LBA with drift scaled by (1 - w) and a linearly collapsing threshold
b (1 - w t / T*), and moves accuracy. Both are affine, so both are
wrappers on the same core; curvature proper appears when w(t)
evolves (rising as T* approaches, or falling as sensory precision
accumulates), still affine, still closed form. The distinctive
prediction is that urgency and collapsing bounds compress only the
right tail, whereas a time prior pulls the fastest responses toward
T* too - so manipulate temporal expectation (pacing, deadline,
foreperiod) and test whether the 10th percentile moves, not only the
90th. Known problems: nothing finishes after T* / w because the
prior eventually dominates (bound w(t) below 1, or give T*
trial-to-trial variability at the cost of one quadrature, as scalar
timing would suggest anyway); and T* is a new time scale competing
with boundary, sigmabias and drift scale, identified only through
RT shape in the first version. Precedents: Frazier & Yu
(2008)
and Drugowitsch et
al. (2012) let a
prior over time act on the threshold; Simen et
al. (2011) model
interval timing as an accumulator whose drift is set to b / T*. This
is a precision-weighted blend of Simen’s timing ramp with the LBA’s
evidence ramp. Order: the finishing-time version first (a specific
g(t), nearly free once the time warp exists), then the rate version,
then w(t), then T* variability if the hard ceiling bites. Unit
tests: w = 0 recovers the LBA, the first version reproduces LBA
choice probabilities exactly, the second with fixed w reproduces the
collapsing-threshold LBA.