Nothing
## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## -----------------------------------------------------------------------------
### uncomment to install packages:
# install.packages("devtools", "coda", "RcppArmadillo", "RcppProgress", "Rglpk", "quadprog")
# devtools::install_github("danheck/multinomineq")
library("multinomineq")
library("coda")
set.seed(1234)
## -----------------------------------------------------------------------------
# total number of observations for each paired comparison:
n <- 25
# Frequencies of choosing "Option H" in Description condition:
HER_desc <- c(9, 16, 16, 7, 12, 16)
# Frequencies of choosing "Option H" in Experience condition:
HER_exp <- c(22, 11, 7, 14, 5, 3)
## -----------------------------------------------------------------------------
# Underweighting of small probabilities CPT:
# (1 = choose "Option H"; 0 = do not choose "Option H")
V <- matrix(
c(
0, 0, 0, 0, 0, 0, # first predicted pattern
0, 0, 0, 1, 0, 0, # second predicted pattern
0, 0, 1, 1, 0, 0, # ...
1, 0, 0, 0, 0, 0,
1, 0, 0, 1, 0, 0,
1, 0, 1, 1, 0, 0,
1, 1, 0, 0, 0, 0,
1, 1, 0, 0, 1, 0,
1, 1, 0, 1, 1, 0,
1, 1, 0, 0, 1, 1,
1, 1, 0, 1, 0, 0,
1, 1, 0, 1, 1, 1,
1, 1, 1, 1, 0, 0,
1, 1, 1, 1, 1, 0,
1, 1, 1, 1, 1, 1
),
ncol = 6, byrow = TRUE
)
## -----------------------------------------------------------------------------
# define a specific vector of choice probabilities:
p_observed <- c(.25, .48, .93, .10, .32, .50)
# check whether p is inside the convex hull defined by V:
inside(p_observed, V = V)
# to find a (random) probability that is inside the convex hull:
find_inside(V = V, random = TRUE)
## -----------------------------------------------------------------------------
A <- matrix(
c(
0, 0, -1, 0, 0, 0, # first inequality
0, 0, 0, 0, 0, -1, # second inequality
-1, 1, 0, 0, 0, 0, # ...
0, -1, 0, 0, 1, 0,
0, 0, 0, 0, -1, 1,
0, 0, 1, -1, 0, 0,
0, 0, 0, 1, 0, 0,
1, 0, 0, 0, 0, 0
),
ncol = 6, byrow = TRUE
)
b <- c(0, 0, 0, 0, 0, 0, 1, 1)
## -----------------------------------------------------------------------------
p_observed <- c(.25, .48, .93, .10, .32, .50)
all(A %*% p_observed <= b)
# corresponding function in multinomineq:
inside(p_observed, A, b)
# find a point that satisfies the constraints:
find_inside(A, b, random = TRUE)
## ---- eval = FALSE------------------------------------------------------------
# ### (not functional for R >= 4.0.0):
# devtools::install_github("TasCL/rPorta")
# Ab <- V_to_Ab(V)
# Ab
## ---- eval = FALSE------------------------------------------------------------
# # fit data from Description and Experience condition:
# p_D <- sampling_binom(k = HER_desc, n = n, V = V, M = 2000, progress = FALSE)
# p_E <- sampling_binom(k = HER_exp, n = n, V = V, M = 2000, progress = FALSE)
## -----------------------------------------------------------------------------
# fit data from Description and Experience condition:
p_D <- sampling_binom(
k = HER_desc, n = n, A = A, b = b,
M = 2000, progress = FALSE
)
p_E <- sampling_binom(
k = HER_exp, n = n, A = A, b = b,
M = 2000, progress = FALSE
)
## ---- fig.width=7, fig.height=5.5---------------------------------------------
# check convergence of MCMC sampler (should look like hairy caterpillars):
plot(p_D)
# summarize posterior samples:
summary(p_D)
## -----------------------------------------------------------------------------
# posterior-predictive p-value for Description condition:
ppp_binom(prob = p_D, k = HER_desc, n = n)
ppp_binom(prob = p_D, k = HER_exp, n = n)
## ----eval = FALSE-------------------------------------------------------------
# # compute Bayes factor using the V-representation
# bf_D <- bf_binom(k = HER_desc, n = n, V = V, M = 10000, progress = FALSE)
# bf_E <- bf_binom(k = HER_exp, n = n, V = V, M = 10000, progress = FALSE)
## -----------------------------------------------------------------------------
# compute Bayes factor using the Ab-representation
bf_D <- bf_binom(HER_desc, n = n, A = A, b = b, M = 10000, progress = FALSE)
bf_E <- bf_binom(HER_exp, n = n, A = A, b = b, M = 10000, progress = FALSE)
bf_D
bf_E
## -----------------------------------------------------------------------------
# proportion of *prior* samples satisfying inequality constraints
cnt <- count_binom(k = 0, n = 0, A = A, b = b, M = 50000, progress = FALSE)
cnt
# proportion of *prior* samples satisfying inequality constraints
# (dataset: Experience)
cnt_E <- count_binom(k = HER_exp, n = n, A = A, b = b, M = 10000, progress = FALSE)
cnt_E
# obtain Bayes factor (including error of approximation)
count_to_bf(posterior = cnt_E, prior = cnt)
## -----------------------------------------------------------------------------
# use a stepwise counting approach for the "Description" condition:
# ("cmin" = minimum of samples that satisfy constraints before sampling stops)
cnt_D <- count_binom(
k = HER_desc, n = n, A = A, b = b,
M = 5000, cmin = 1, steps = c(3, 5, 7), progress = FALSE
)
cnt_D
# obtain Bayes factor (including error of approximation)
count_to_bf(posterior = cnt_D, prior = cnt)
## -----------------------------------------------------------------------------
# binomial data format
n <- 25
HER_desc <- c(9, 16, 16, 7, 12, 16)
# multinomial data format ("k" must be a vector!)
k_multinom <- c(
9, 16, # first binary gamble
16, 9, # second binary gamble
16, 9, # ...
7, 17,
12, 13,
16, 9
)
options <- c(2, 2, 2, 2, 2, 2) # 2 options for each type
## -----------------------------------------------------------------------------
mn <- binom_to_multinom(HER_desc, n)
mn
## -----------------------------------------------------------------------------
posterior <- sampling_multinom(k = mn$k, options = mn$options, A = A, b = b)
bayesfactor <- bf_multinom(k = k_multinom, options = options, A = A, b = b)
bayesfactor
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.