# Load packages -----------------------------------------------------------
library(BayesFactor) # needed to calculate Bayes factors
library(here) # needed to assign working directory for relative paths
library(tidyverse) # needed for importing, processing, and plotting data
library(ggplot2)
library(ggpubr)
library(brms)
library(janitor)
library(tidybayes)

Loading required package: coda

Loading required package: Matrix

************
Welcome to BayesFactor 0.9.12-4.4. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).

Type BFManual() to open the manual.
************

here() starts at H:/projects/scientific/beams

── Attaching packages ──────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.4.1     ✔ purrr   1.0.1
✔ tibble  3.1.8     ✔ dplyr   1.1.0
✔ tidyr   1.3.0     ✔ stringr 1.5.0
✔ readr   2.1.4     ✔ forcats 1.0.0
── Conflicts ─────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand() masks Matrix::expand()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ tidyr::pack()   masks Matrix::pack()
✖ tidyr::unpack() masks Matrix::unpack()

Attaching package: 'janitor'


The following objects are masked from 'package:stats':

    chisq.test, fisher.test



Attaching package: 'tidybayes'


The following objects are masked from 'package:brms':

    dstudent_t, pstudent_t, qstudent_t, rstudent_t

rt_dat <- read.csv(file = 'C:\\Users\\mariu\\Documents\\rt_dat.csv', sep=',')

head(rt_dat)

A data.frame: 6 × 4

	X	Run	MotionCongruency	RT
	<int>	<chr>	<chr>	<dbl>
1	1	run-02	congruent	0.775
2	2	run-02	incongruent	0.618
3	3	run-02	congruent	0.666
4	4	run-02	congruent	0.780
5	5	run-02	incongruent	0.500
6	6	run-02	incongruent	0.648

ggplot(rt_dat, aes(x=MotionCongruency, y = RT, fill=MotionCongruency)) + 
        geom_boxplot(alpha = 0.7, notch = TRUE, width=0.3, outlier.shape = NA, na.rm=TRUE) +
        # geom_violin(col=NA, alpha = 0.7, width=0.5) +
        geom_point(position=position_jitterdodge(0.15), alpha=0.3, aes(group=MotionCongruency)) + #, col=Run)) +
        stat_summary(fun=median, geom="line", linewidth=1, aes(group = Run, col=Run)) +
        labs(y = "RT (sec)") +
        guides(col = 'none') +       
            theme(
                plot.title = element_text(size = (15), face = "bold"),
                axis.text.y=element_text(size = 16, color = "black"),
                axis.text.x=element_blank(),
                # axis.ticks.x =  element_line(size = 1, color = "black"),
                legend.title = element_blank(),
                legend.text = element_text(size = 18),
                panel.background = element_rect(color = "#eeeeee", fill = "#eeeeee"),
                # panel.background = element_blank(),
                panel.grid.major.y = element_line(color = "#eeeeee", linewidth = 0.5),
                panel.grid.minor.y = element_line(color = "#eeeeee", linewidth = 0.5),
                panel.grid.major.x = element_line(color = "#eeeeee", linewidth = 0.5),
                panel.grid.minor.x = element_line(color = "#eeeeee", linewidth = 0.5),
                axis.title.y = element_text(hjust = 0.5, size = 18),
                axis.title.x = element_text(hjust = 0.5, size = 18),
                # axis.line = element_blank(),
                # axis.line = element_line(),
                axis.line = element_line(color = "white"),
                strip.background = element_rect(fill = "#444444", color = "white"),
                strip.text = element_text(
                    size = 18, face = "bold", hjust = 0.5, vjust = 0.5, color="white")
                )

shifted_lognormal (simplest model)

priors_ = c(
    set_prior('normal(log(.4), 0.05)', class='b', coef='Intercept'),
    # set_prior('normal(-1, 0.2)', class='Intercept', dpar = "sigma"),
    set_prior('normal(0, 0.05)', class='b', coef='MotionCongruencyincongruent'),
    # set_prior('normal(-0.01, .1)', class='b', dpar = "sigma"),
    set_prior('normal(.3, 0.15)', class='sigma', lb=0)#,
    # set_prior('normal(0.3, .2)', class='sd', dpar = "sigma", lb = 0)
)

mf_ <- bf(
    RT ~ 0 + Intercept + MotionCongruency # + (MotionCongruency || Run)#, sigma ~ MotionCongruency + (MotionCongruency | Run)
) 

hyp_string <- 'MotionCongruencyincongruent = 0'

brms_mm = brm(
    mf_,
    data = rt_dat,
    prior = priors_,
    sample_prior = TRUE,
    family = shifted_lognormal(link = "identity"),
    chains = 4,
    # iter = 10,
    cores = 4,
    save_pars = save_pars(all=TRUE),
    control = list(adapt_delta = 0.9, max_treedepth = 15))

Compiling Stan program...

Start sampling

Prior predictive check (sample_prior = 'only')

pp_check(brms_mm, ndraws=50)

Posterior predictive check

pp_check(brms_mm, ndraws=50)

summary(brms_mm)

 Family: shifted_lognormal 
  Links: mu = identity; sigma = identity; ndt = identity 
Formula: RT ~ 0 + Intercept + MotionCongruency 
   Data: rt_dat (Number of observations: 593) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                      -1.05      0.05    -1.13    -0.96 1.00     1153
MotionCongruencyincongruent    -0.04      0.02    -0.08    -0.00 1.00     1560
                            Tail_ESS
Intercept                       1413
MotionCongruencyincongruent     1706

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.28      0.01     0.25     0.31 1.00     1207     1786
ndt       0.10      0.02     0.07     0.13 1.00     1113     1471

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

hyp <- hypothesis(brms_mm, "MotionCongruencyincongruent = 0")
hyp

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (MotionCongruency... = 0    -0.04      0.02    -0.08        0       0.38
  Post.Prob Star
1      0.28    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

plot(hyp)

exgaussian (simplest model)

priors_ = c(
    set_prior('normal(400, 50)', class='b', coef='Intercept'),
    # set_prior('normal(-1, 0.2)', class='Intercept', dpar = "sigma"),
    set_prior('normal(0, 50)', class='b', coef='MotionCongruencyincongruent'),
    # set_prior('normal(-0.01, .1)', class='b', dpar = "sigma"),
    set_prior('normal(300, 100)', class='sigma', lb=0)#,
    # set_prior('normal(0.3, .2)', class='sd', dpar = "sigma", lb = 0)
)

mf_ <- bf(
    RT ~ 0 + Intercept + MotionCongruency# + (MotionCongruency || Run)#, sigma ~ MotionCongruency + (MotionCongruency | Run)
) 

hyp_string <- 'MotionCongruencyincongruent = 0'
    
rt_dat_ms <- rt_dat
rt_dat_ms$RT <- rt_dat_ms$RT * 1000

brms_mm = brm(
    mf_,
    data = rt_dat_ms,
    prior = priors_,
    sample_prior = TRUE,
    family = exgaussian(link = "identity"),
    chains = 4,
    # iter = 10,
    cores = 4,
    save_pars = save_pars(all=TRUE),
    control = list(adapt_delta = 0.9, max_treedepth = 15))

Compiling Stan program...

Start sampling

Prior predictive check (sample_prior = 'only')

pp_check(brms_mm, ndraws=50)

Posterior predictive check

pp_check(brms_mm, ndraws=50)

summary(brms_mm)

 Family: exgaussian 
  Links: mu = identity; sigma = identity; beta = identity 
Formula: RT ~ 0 + Intercept + MotionCongruency 
   Data: rt_dat_ms (Number of observations: 593) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                     457.60      5.10   447.38   467.56 1.00     2462
MotionCongruencyincongruent    -2.02      6.02   -13.86     9.82 1.00     2709
                            Tail_ESS
Intercept                       2571
MotionCongruencyincongruent     2730

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    43.15      3.68    35.95    50.41 1.00     2251     2167
beta     95.64      5.74    84.73   107.34 1.00     2137     2350

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

hyp <- hypothesis(brms_mm, "MotionCongruencyincongruent = 0")
hyp

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (MotionCongruency... = 0    -2.02      6.02   -13.86     9.82       7.74
  Post.Prob Star
1      0.89     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

plot(hyp)