Bayesian Inference Cognitive Model

bayes() fits a Bayesian cognitive model, updating beliefs about the probability of discrete event outcomes based on the frequencies of outcomes.

bayes_beta_c() fits a model for 2 outcomes (beta-binomial) for continuous responses
bayes_beta_d() fits a model for 2 outcomes (beta-binomial) for discrete responses
bayes_dirichlet_c() fits a model for n > 2 outcomes (dirichlet-categorical/multinomial) for continuous responses
bayes_dirichlet_d() fits a model for n > 2 outcomes (dirichlet-categorical/multinomial) for discrete responses

Usage

bayes_beta_c(
  formula,
  data,
  fix = NULL,
  format = c("raw", "count", "cumulative"),
  prior_sum = NULL,
  ...
)

bayes_beta_d(formula, data, fix = NULL, format = NULL, prior_sum = NULL, ...)

bayes_dirichlet_d(
  formula,
  data,
  fix = NULL,
  format = NULL,
  prior_sum = NULL,
  ...
)

bayes_dirichlet_c(
  formula,
  data,
  fix = NULL,
  format = NULL,
  prior_sum = NULL,
  ...
)

bayes(
  formula,
  data = data.frame(),
  fix = list(),
  format = c("raw", "count", "cumulative"),
  type = NULL,
  discount = 0L,
  options = list(),
  prior_sum = NULL,
  ...
)

Arguments

formula

A formula, the variables in data to be modeled. For example, y ~ x1 + x2 models response y as function of one stimulus with features x1, x2.

data

A data frame, the data to be modeled.

fix

(optional) A list with parameter-value pairs of fixed parameters. If missing all free parameters are estimated. If set to "start" all parameters are fixed to their start values. Model parameter names depend on formula and can be delta, x1, x2 (see details - model parameters).

list(x1 = 3.09) sets parameter x1 equal to 3.09.
list(x1 = "x2") sets parameter x1 equal to parameter x2 (estimates x2).
list(x2 = "x1", x1 = 3.09) sets parameter x2 equal to parameter x1 and sets x1 equal to 3.09 (estimates none of the two).
list(x1 = NA) omits the parameter x1, if possible.
"start" sets all parameters equal to their initial values (estimates none). Useful for building a first test model.

format

(optional) A string, the format the data to be modeled, can be abbreviated, default is "raw"; allowed values:

"raw" means that the data are trial-by-trial binary occurrence indicators: 1, 0, 1, ... means the event happened in trial with a value of 1.
"cumulative" means the data are trial-by-trial cumulative counts of events: 0, 1, 1, 2, ... counts how often the event happened up to the trial.
"count" means the data are total events counts, ignoring the trial-by-trial order of events: 2, 10, ... means the event happened 2 times, then (starting from zero!) it happened 10 times.

prior_sum

(optional) A number; the prior hyperparameter will be constrained to sum to this number; defaults to the number of prior parameters; if prior_sum = NA no sum constraint is placed.

...

other arguments, ignored.

type

(optional) A string, the type of inference, "beta-binomial" or "dirichlet-multinomial". Can be abbreviated. Will be inferred, if missing.

discount

A number, how many initial trials to not use during parameter fitting.

options

(optional) A list, list entries change the modeling procedure. For example, list(lb = c(k=0)) changes the lower bound of parameter k to 0, or list(fit_measure = "mse") changes the goodness of fit measure in parameter estimation to mean-squared error, for all options, see cm_options.

Value

Returns a cognitive model object, which is an object of class cm. A model, that has been assigned to m, can be summarized with summary(m) or anova(m). The parameter space can be viewed using pa. rspace(m), constraints can be viewed using constraints(m).

Details

The model models -- as response -- the belief about the occurrence of the first event in the formula as follows:

y ~ x1 models the beliefe about event x1 occurring versus it not occurring.
y ~ x1 + x2 models beliefs about x1 versus x2 occurring.
y ~ x1 + x2 + x3 models beliefs about x1, x2, and x3 occurring.

Model Parameters

The model has n + 1 (n = number of events) free parameters, which are:

delta is the learning rate, it weights the observation during learning, value < 1 causes conservatism, > 1 causes liberal learning, and 1 is optimal Bayesian.
x1, x2 (dynamic names) are the prior parameter, their names correspond to the right side of formula. Also known as the hyperparameter of the prior belief distribution before trial 1. If they are constrainted to sum to n and n - 1 parameter are estimated.
In bayes_beta_d() or bayes_dirichlet_d(): If choicerule = "softmax": tau is the temperature or choice softness, higher values cause more equiprobable choices. If choicerule = "epsilon": eps is the error proportion, higher values cause more errors from maximizing.

References

Griffiths, T. L., & Yuille, A. (2008). Technical Introduction: A primer on probabilistic inference. In N. Chater & M. Oaksford (Eds.), The Probabilistic Mind: Prospects for Bayesian Cognitive Science (pp. 1 - 2). Oxford University Press. https://doi.org/10.1093/acprof:oso/9780199216093.003.0002

Tauber, S., Navarro, D. J., Perfors, A., & Steyvers, M. (2017). Bayesian models of cognition revisited: Setting optimality aside and letting data drive psychological theory. Psychological Review, 124(4), 410 - 441. http://dx.doi.org/10.1037/rev0000052

Author

Markus Steiner

Jana B. Jarecki, jj@janajarecki.com

Examples

D <- data.frame(
  a = c(0,0,1,1,1),              # event A, e.g. coin toss "heads"
  b = c(1,1,0,0,0),              # event B, complement of A
  y = c(0.5,0.3,0.2,0.3,0.5))    # participants' beliefs about A

M <- bayes_beta_c(
     formula = y ~ a + b,
     data = D)   # fit all parameters
#> Fitting free parameters [delta, a, sigma] by maximizing loglikelihood (truncnorm pdf) with solnp.
predict(M)                        # predict posterior means
#> [1] 0.4879909 0.2804911 0.1968065 0.3813738 0.4969667
summary(M)                        # summarize model
#> 
#> Model:
#>   with no choice rule
#> Call:
#> y ~ a + b
#> 
#> (Constrained) Parameters:
#>       Estimate
#> delta 1.48    
#> a     0.976   
#> b     1.024   
#> sigma 0.03786 
#> 
#> Fit Measures:
#> MSE: 0.0014, LL: 9.3, AIC: -9.5, BIC: -14
#> 
parspace(M)                       # view parameter space
#> 
#> Parameter space of the cognitive model 'Bayesian model':
#>               lb         ub      start         na
#> delta  0.0000000 10.0000000  1.0000000  1.0000000
#> a      0.0000100  2.0000000  1.0000000  1.0000000
#> b      0.0000100  2.0000000  1.0000000  1.0000000
#> sigma  0.0000001  0.3000000  0.1500000         NA
#> ---
#> Note. lb = lower bound, ub = upper bound, start = start value.
anova(M)                          # anova-like table
#> Sum Sq. Table
#>  N Par    Sum Sq   Mean Sq
#>      3 0.0071659 0.0014332
logLik(M)                         # loglikelihood
#> 'log Lik.' 9.274954 (df=3)
MSE(M)                            # mean-squared error   
#> [1] 0.001433182


# Predictions ----------------------------------------------
predict(M, type = "mean")                  # posterior mean
#> [1] 0.4879909 0.2804911 0.1968065 0.3813738 0.4969667
predict(M, type = "max")                   # maximum posterior
#> [1]         -Inf -0.016233531 -0.008116766  0.327922156  0.495941617
predict(M, type = "sd")                    # posterior SD
#> [1] 0.2885919 0.2122563 0.1628695 0.1780914 0.1674263
predict(M, type = "posteriorpar")          # posterior hyper-par.
#>           pr_a     pr_s
#> [1,] 0.9759817 1.024018
#> [2,] 0.9759817 2.503564
#> [3,] 0.9759817 3.983111
#> [4,] 2.4555279 3.983111
#> [5,] 3.9350741 3.983111
predict(M, type = "draws", ndraws = 3)     #  --"--  3 draws
#> Warning: These arguments to 'predict.cogscimodel' are invalid and dropped: list(ndraws = 3)
#>          pr_a1     pr_a2       pr_a3      pr_s1     pr_s2      pr_s3
#> [1,] 0.9555926 0.7922289 0.914759423 0.04440736 0.2077711 0.08524058
#> [2,] 0.1652925 0.6759333 0.421011213 0.83470749 0.3240667 0.57898879
#> [3,] 0.2483440 0.1477724 0.001690943 0.75165602 0.8522276 0.99830906
#> [4,] 0.6732781 0.3551378 0.430624132 0.32672191 0.6448622 0.56937587
#> [5,] 0.7924639 0.6221137 0.530650827 0.20753607 0.3778863 0.46934917


# Fix parameter ---------------------------------------------
bayes_beta_c(~a+b, D, list(delta=1, priorpar=c(1, 1)))  # delta=1, uniform prior
#> Fitting free parameters [sigma] by maximizing loglikelihood (truncnorm pdf) with solnp.
#> Error: Can't compute goodness of fit, because the model has no observed resonses..
#>   * Did you forget a left side in 'formula'? (such as 'y' in y ~ x1 + x2)
bayes_beta_c(~a+b, D, list(delta=1, a=1, b=1))          # -- (same) --
#> Fitting free parameters [sigma] by maximizing loglikelihood (truncnorm pdf) with solnp.
#> Error: Can't compute goodness of fit, because the model has no observed resonses..
#>   * Did you forget a left side in 'formula'? (such as 'y' in y ~ x1 + x2)
bayes_beta_c(~a+b, D, fix = "start")                    # fix to start values
#> Bayesian model | choice rule: none
#> Call:
#> bayes_beta_c(formula = ~a + b, data = D, fix = "start")
#> 
#> Constrained and fixed parameters:
#> delta      a      b  sigma  
#>   1.0    1.0    1.0    0.5  
#> 
#> ---
#> Note:  No free parameters. View constraints by constraints(.), view parameter space by parspace(.)


# Parameter fitting ----------------------------------------
# Use a response variable, y, to which we fit parameter
bayes(y ~ a + b, D, fix = "start")              # "start" fixes all par., fit none 
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = "start")
#> 
#> Constrained and fixed parameters:
#> delta      a      b  sigma  
#>  1.00   1.00   1.00   0.15  
#> 
#> ---
#> Note:  No free parameters. View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, fix = list(delta=1))         # fix delta, fit priors 
#> Fitting free parameters [a, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = list(delta = 1))
#> 
#> Free parameters: estimates 
#>     a  sigma  
#> 0.911  0.044  
#> 
#> Constrained and fixed parameters:
#> delta      b  
#>   1.0    1.1  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, fix = list(a=1, b=1))        # fix priors, fit delta 
#> Fitting free parameters [delta, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = list(a = 1, b = 1))
#> 
#> Free parameters: estimates 
#> delta  sigma  
#> 1.560  0.038  
#> 
#> Constrained and fixed parameters:
#>    a     b  
#>    1     1  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, fix = list(delta=1, a=1))    # fix delta & prior on "a"
#> Fitting free parameters [sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = list(delta = 1, a = 1))
#> 
#> Free parameters: estimates 
#> sigma  
#> 0.052  
#> 
#> Constrained and fixed parameters:
#> delta      a      b  
#>     1      1      1  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, list(delta=1, b=1))          # fix delta & prior on "b"
#> Fitting free parameters [sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = list(delta = 1, b = 1))
#> 
#> Free parameters: estimates 
#> sigma  
#> 0.052  
#> 
#> Constrained and fixed parameters:
#> delta      a      b  
#>     1      1      1  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)


### Parameter meanings
# ---------------------------------------
# delta parameter: the learning rate or evidence weight
bayes(y ~ a + b, D, c(delta = 0))             # 0   -> no learning
#> Fitting free parameters [a, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = c(delta = 0))
#> 
#> Free parameters: estimates 
#>     a  sigma  
#>  0.72   0.12  
#> 
#> Constrained and fixed parameters:
#> delta      b  
#>   0.0    1.3  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, c(delta = 0.1))           # 0.1 -> slow learning
#> Fitting free parameters [a, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = c(delta = 0.1))
#> 
#> Free parameters: estimates 
#>     a  sigma  
#>  0.74   0.10  
#> 
#> Constrained and fixed parameters:
#> delta      b  
#>   0.1    1.3  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, c(delta = 9))             # 9   -> fast learning
#> Fitting free parameters [a, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = c(delta = 9))
#> 
#> Free parameters: estimates 
#>     a  sigma  
#>  1.10   0.11  
#> 
#> Constrained and fixed parameters:
#> delta      b  
#>   9.0    0.9  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, c(a=1.5, b=0.5))                # prior: a more likely
#> Fitting free parameters [delta, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = c(a = 1.5, b = 0.5))
#> 
#> Free parameters: estimates 
#> delta  sigma  
#>  3.33   0.12  
#> 
#> Constrained and fixed parameters:
#>    a     b  
#>  1.5   0.5  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, list(priorpar=c(1.5, 0.5)))     # -- (same) --
#> Fitting free parameters [delta, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = list(priorpar = c(1.5,  ...
#> 
#> Free parameters: estimates 
#> delta  sigma  
#>  3.33   0.12  
#> 
#> Constrained and fixed parameters:
#>    a     b  
#>  1.5   0.5  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, c(a = 0.1, b=1.9))              # prior: b more likely
#> Fitting free parameters [delta, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = c(a = 0.1, b = 1.9))
#> 
#> Free parameters: estimates 
#> delta  sigma  
#>  2.68   0.25  
#> 
#> Constrained and fixed parameters:
#>    a     b  
#>  0.1   1.9  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)
bayes(y ~ a + b, D, list(priorpar = c(0.1, 1.9)))   # -- (same) --
#> Fitting free parameters [delta, sigma] by maximizing loglikelihood (normal pdf) with solnp.
#> Bayesian model | choice rule: none
#> Call:
#> bayes(formula = y ~ a + b, data = D, fix = list(priorpar = c(0.1,  ...
#> 
#> Free parameters: estimates 
#> delta  sigma  
#>  2.68   0.25  
#> 
#> Constrained and fixed parameters:
#>    a     b  
#>  0.1   1.9  
#> 
#> ---
#> Note:  View constraints by constraints(.), view parameter space by parspace(.)