## ----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)
library(ggplot2)

# load data
data <- adaptation_level

# separate groups and parts
group1_part2 <- data %>% filter(group == 1 & part == 2)
group2_part2 <- data %>% filter(group == 2 & part == 2)

## ----message=FALSE, warning=FALSE, results = 'hide'---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# fit
fit1 <- b_linear(x=group1_part2$sequence,
                 y=group1_part2$response,
                 s=group1_part2$subject,
                 iter=200, warmup=100, chains=1)

fit2 <- b_linear(x=group2_part2$sequence,
                 y=group2_part2$response,
                 s=group2_part2$subject,
                 iter=200, warmup=100, chains=1)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# trace plots
plot_trace(fit1)
plot_trace(fit2)

# the commands below are commented out for the sake of brevity
#print(fit1)
#print(fit2)

# visual inspection
plot(fit1)
plot(fit2)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
summary(fit1)
summary(fit2)

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
comparison_results <- compare_means(fit1, fit2=fit2)

plot_means_difference(fit1, fit2=fit2, par="intercept")

## ----message=FALSE, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
plot_distributions(fit1, fit2) +
  labs(title="Part II", x="measurement number", y="") +
  theme(legend.position="none") +
  scale_x_continuous(limits=c(1, 10), breaks=seq(1:10)) +
  ylim(0, 10)