## ----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, results = 'hide'---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# libs
library(bayes4psy)
library(ggplot2)

# load data
data <- stroop_simple

# reading incongruent vs reading neutral
ri_vs_rn <- data$reading_incongruent - data$reading_neutral

# fit
fit_ri_vs_rn <- b_ttest(ri_vs_rn,
                        iter=200, warmup=100, chains=1)

# traceplot
#plot_trace(fit_ri_vs_rn)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# libs
# compare the fit with 0
comparison <- compare_means(fit_ri_vs_rn, mu=0)

## ----message=FALSE, warning=FALSE, results = 'hide'---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# priors
mu_prior <- b_prior(family="normal", pars=c(60, 30))
sigma_prior <- b_prior(family="uniform", pars=c(0, 120))
priors <- list(c("mu", mu_prior),
               c("sigma", sigma_prior))

# fit
fit_reading_neutral <- b_ttest(data$reading_neutral,
                               priors=priors,
                               iter=200, warmup=100, chains=1)

fit_reading_incongruent <- b_ttest(data$reading_incongruent,
                                   priors=priors,
                                   iter=200, warmup=100, chains=1)

fit_naming_neutral <- b_ttest(data$naming_neutral,
                              priors=priors,
                              iter=200, warmup=100, chains=1)

fit_naming_incongruent <- b_ttest(data$naming_incongruent,
                                  priors=priors,
                                  iter=200, warmup=100, chains=1)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# trace plots
#plot_trace(fit_reading_neutral)
#plot_trace(fit_reading_incongruent)
#plot_trace(fit_naming_neutral)
plot_trace(fit_naming_incongruent)

# check fit
# the commands below are commented out for the sake of brevity
#print(fit_reading_neutral)
#print(fit_reading_incongruent)
#print(fit_naming_neutral)
#print(fit_naming_incongruent)

# visual inspection
#plot(fit_reading_neutral)
#plot(fit_reading_incongruent)
#plot(fit_naming_neutral)
plot(fit_naming_incongruent)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# join all fits but the first in a list
fit_list <- c(fit_reading_incongruent,
              fit_naming_neutral,
              fit_naming_incongruent)

# compare all fits simultaneously
multiple_comparison <- compare_means(fit_reading_neutral,
                                     fits=fit_list)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
plot_means(fit_reading_neutral, fits=fit_list) +
  scale_fill_hue(labels=c("Reading neutral",
                          "Reading incongruent",
                          "Naming neutral",
                          "Naming incongruent")) +
  theme(legend.title=element_blank())

## ----message=FALSE, warning=FALSE, fig.height=9, fig.width=8------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
plot_means_difference(fit_reading_neutral, fits=fit_list)