## ----message=FALSE, warning=FALSE, echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# knitr options
knitr::opts_chunk$set(fig.width=6, fig.height=4.5)
options(width=800)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# libs
library(bayes4psy)
library(dplyr)

# load data
data <- flanker

# analyse only correct answers, load test and control data
control_rt <- data %>% filter(result == "correct" &
                                 group == "control")

test_rt <- data %>% filter(result == "correct" &
                                 group == "test")

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# control group subject indexes range is 22..45 cast it to 1..23
# test group indexes are OK
control_rt$subject <- control_rt$subject - 21

## ----message=FALSE, warning=FALSE, results = 'hide'---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# prior
uniform_prior <- b_prior(family="uniform", pars=c(0, 3))

# attach priors to relevant parameters
priors <- list(c("mu_m", uniform_prior))

# fit
rt_control_fit <- b_reaction_time(t=control_rt$rt,
                                  s=control_rt$subject,
                                  priors=priors,
                                  chains=1, iter=200, warmup=100)

rt_test_fit <- b_reaction_time(t=test_rt$rt,
                               s=test_rt$subject,
                               priors=priors,
                               chains=1, iter=200, warmup=100)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# plot trace
plot_trace(rt_control_fit)
plot_trace(rt_test_fit)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# the commands below are commented out for the sake of brevity
#print(rt_control_fit)
#print(rt_test_fit)

## ----eval=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  rt_control_fit <- b_reaction_time(t=control_rt$rt,
#                                    s=control_rt$subject,
#                                    iter=4000)
#  
#  rt_test_fit <- b_reaction_time(t=test_rt$rt,
#                                 s=test_rt$subject,
#                                 iter=4000)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# check fits
plot(rt_control_fit)

plot(rt_test_fit)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# set rope (region of practical equivalence) interval to +/- 10ms
rope <- 0.01

# which mean is smaller/larger
rt_control_test <- compare_means(rt_control_fit, fit2=rt_test_fit, rope=rope)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# difference plot
plot_means_difference(rt_control_fit, fit2=rt_test_fit, rope=rope)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# visual comparsion of mean difference
plot_means(rt_control_fit, fit2=rt_test_fit)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# map correct/incorrect/timeout to 1 (success) or 0 (fail)
data$result_numeric <- 0
data[data$result == "correct", ]$result_numeric <- 1

# split to control and test groups
control_sr <- data %>% filter(group == "control")
test_sr <- data %>% filter(group == "test")

# control group subject indexes range is 22..45 cast it to 1..23
# test group indexes are OK
control_sr$subject <- control_sr$subject - 21

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# priors
p_prior <- b_prior(family="beta", pars=c(1, 1))

# attach priors to relevant parameters
priors <- list(c("p", p_prior))

# fit
sr_control_fit <- b_success_rate(r=control_sr$result_numeric,
                                 s=control_sr$subject,
                                 priors=priors,
                                 chains=1, iter=200, warmup=100)

sr_test_fit <- b_success_rate(r=test_sr$result_numeric,
                              s=test_sr$subject,
                              priors=priors,
                              chains=1, iter=200, warmup=100)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# inspect control group fit
plot_trace(sr_control_fit)
plot(sr_control_fit, subjects=FALSE)
print(sr_control_fit)

# inspect test group fit
plot_trace(sr_test_fit)
plot(sr_test_fit, subjects=FALSE)
#print(sr_test_fit)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# which one is higher
sr_control_test <- compare_means(sr_control_fit, fit2=sr_test_fit)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# difference plot
plot_means_difference(sr_control_fit, fit2=sr_test_fit)