bhsdtr: Creating and fitting ordinal models with ordered thresholds

In boryspaulewicz/bhsdtr2: Bayesian (non-)hierarchical ordinal models with ordered thresholds

Creating and fitting ordinal models with ordered thresholds

Description

This is the main method of fitting ordinal models in the bhsdtr2 package. The model type (e.g., SDT, UV SDT, DPSDT, meta-d') is inferred automatically based on the supplied model formulae. An EV SDT model is fitted if the user provides model formulae for dprim and thr, a meta-d' model is fitted if there is a formula for the metad parameter vector, an UV version is fitted if there is formula for the sd_ratio parameter. Fitting an DPSDT model is a bit more complicated because there are two different versions of this model and because the probability of recall may be non-zero only for one stimulus class (e.g., "single-interval" tasks) or for both stimulus classes (2AFC tasks). In order to fit an DPSDT model the user has to introduce a recall indicator in the response_formula, e.g., resp ~ stim + possible_recall, where possible_recall is 1 iff recall is possible (a correct response with maximum confidence) on that trial. If d' and R may be correlated, then the user may model the dprimr parameter vector (the first element is the d' parameter, the second element is R), otherwise three formulas have to be provided: one for dprim, one for thr, and one for R. At present, only the logit link function is implemented for the R parameter.

Usage

bhsdtr(
  model_formulae,
  response_formula,
  data,
  links = list(gamma = "log_distance"),
  method = "jmap",
  prior = list(),
  sample.prior = F,
  thresholds_scale = 2,
  force.id_log = F,
  fit = T,
  ...
)

Arguments

model_formulae	a vector of two or more lmer-style model formulae for the parameters of the model. Every such vector must contain a model formula for the thr (=thresholds) parameter, and a model formula for either the dprim, the metad, or the mean (general ordinal model) parameter. This vector may have three elements iff it also contains a formula for the sdratio (Unequal Variance models) parameter.
response_formula	a formula specifying the names of the response variable and the stimulus variable. If this is an SDT-like model this must be of the form resp ~ stim [+ modifier], where resp is the name of the combined response variable or the binary response variable, stim is the name of the (binary) stimulus variable, and the optional modifier indicates something model-specific (e.g., the possibility of recall in DPSDT).
links	a list specyfing the link functions, the default is list(gamma = 'log_distance'). This list will be filled with the default link function specifications for the remaining parameters if these are not specified. Currently you can set gamma to 'log_distance' (the default), 'log_ratio', 'twoparameter', 'parsimonious', and 'softmax'. The log_distance and log_ratio link functions are described in the Readme file in the github repository, the parsimonious link function represents the thresholds using two parameters, the unconstrained main criterion and the logarithm of the scaling factor which is used to stretch the "unbiased" thresholds so that the distribution of responses fits. The softmax link function is described in the bhsdtr preprint. You can also choose between delta = 'log' (the default, which covers both the d' and the meta-d' parameters), and delta = 'identity' (i.e., no conversion between delta and d' or meta-d').
method	[= 'jmap'] a string which specifies how to fit the model. The default is 'jmap', which results in an attempt to quickly maximize the joint posterior by using the rstan::optimizing function, but you can also choose 'stan', or, if you do not want the model to be fitted right away, you can use anything else.
prior	[= NULL] a list with non-default prior settings which will be passed to the set.prior function. See the set.prior function documentation for details.
thresholds_scale	thresholds scaling factor for the softmax gamma link function (see the bhsdtr preprint).
force.id_log	If FALSE (the default) separate intercepts parametrization will be required for the delta / d' / meta-d' model matrices when using the id_log link function.
fit	If TRUE (the default) then also fit the model
...	arguments to be passed to the stan function.

Value

a bhsdtr_model object, which is an S3 class object. If it was fitted using the jmap (stan) method, $jmapfit ($stanfit) contains the result. You can refit the model using either method (fit(model, ...)) and the $jmapfit or the $stanfit slot will be updated.

Examples

    gabor$r = combined.response(gabor$stim, gabor$rating, gabor$acc)
    ## A hierarchical SDT model
    m.jmap = bhsdtr(c(dprim ~ duration * order + (duration | id), thr ~ order + (1 | id)),
                  r ~ stim,
                  gabor)

    ## Posterior samples (here single samples = maximum joint
    ## posterior point estimates) of d' for every unique combination
    ## of the predictors, here duration and order, because
    ## dprim ~ duration * order.
    samples(m.jmap, 'dprim')

    ## Thresholds
    samples(m.jmap, 'thr')

    ## Fitting the model using the Stan sampler
    m.stan = bhsdtr(c(dprim ~ duration * order + (duration | id), thr ~ order + (1 | id)),
                    r ~ stim,
                    gabor, method = 'stan')
    (smp = samples(m.stan, 'thr'))

    ## This is how you can access the stanfit object
    print(m.stan$stanfit, probs = c(.025, .975), pars = c('delta_fixed', 'gamma_fixed'))

    ## The plot method lets you see if the model fits. Here we
    ## specify the variables (duration and order) to see only
    ## the duration x order panels.
    library(ggplot2)
    plot(m.stan, vs = c('duration', 'order'), verbose = F)

    ## A simple contrast calculated on the posterior threshold samples
    round(t(apply(smp, 2, function(x)quantile(x[,'DECISION-RATING'] - x[,'RATING-DECISION'],
                                              c(.025, .975)))), 2)

boryspaulewicz/bhsdtr2 documentation built on July 17, 2024, 8:22 p.m.