In [ ]:
In [25]:
library(brms)
library(ggplot2)
In [22]:
rt_dat <- read.csv("./rt_dat.csv")
In [23]:
head(rt_dat)
A data.frame: 6 × 4
XRunMotionCongruencyRT
<int><chr><chr><dbl>
11run-03congruent 1.594
22run-03congruent 0.669
33run-03incongruent0.429
44run-03congruent 0.459
55run-03congruent 0.801
66run-03incongruent0.477
In [26]:
ggplot(
    rt_dat, aes(x=MotionCongruency, y = RT, fill=MotionCongruency)) + 
    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", size=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(),
        legend.title = element_blank(),
        legend.text = element_text(size = 18),
        panel.background = element_rect(color = "#eeeeee", fill = "#eeeeee"),
        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_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")
)
In [ ]:
In [37]:
mf_ <- bf(
    RT ~ MotionCongruency, sigma ~ MotionCongruency
)

priors_ = c(
    set_prior('normal(0, .1)', class='b'),
    set_prior('normal(0, .5)', class='b', dpar = "sigma")
)

fit_exg <- brm(
    mf_,
    data = rt_dat,
    prior = priors_,
    sample_prior = TRUE,
    family = exgaussian(link = "identity"),
    chains = 4,
    cores = 4,
    save_pars = save_pars(all=TRUE),
    control = list(adapt_delta = 0.9, max_treedepth = 15)
)

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

Compiling Stan program...

Start sampling

Compiling Stan program...

Start sampling
In [38]:
summary(fit_exg)

 Family: exgaussian 
  Links: mu = identity; sigma = log; beta = identity 
Formula: RT ~ MotionCongruency 
         sigma ~ MotionCongruency
   Data: rt_dat (Number of observations: 261) 
  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
Intercept                             0.63      0.01     0.60     0.65 1.00
sigma_Intercept                      -2.84      0.16    -3.16    -2.52 1.00
MotionCongruencyincongruent          -0.08      0.01    -0.10    -0.05 1.00
sigma_MotionCongruencyincongruent    -0.86      0.25    -1.36    -0.39 1.00
                                  Bulk_ESS Tail_ESS
Intercept                             3185     2329
sigma_Intercept                       3238     2723
MotionCongruencyincongruent           3099     2545
sigma_MotionCongruencyincongruent     3081     2773

Family Specific Parameters: 
     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
beta     0.15      0.01     0.13     0.18 1.00     2407     2154

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).
In [39]:
summary(fit_sln)

 Family: shifted_lognormal 
  Links: mu = identity; sigma = log; ndt = identity 
Formula: RT ~ MotionCongruency 
         sigma ~ MotionCongruency
   Data: rt_dat (Number of observations: 261) 
  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
Intercept                            -0.99      0.08    -1.13    -0.83 1.00
sigma_Intercept                      -0.92      0.10    -1.13    -0.73 1.00
MotionCongruencyincongruent          -0.21      0.05    -0.30    -0.12 1.00
sigma_MotionCongruencyincongruent     0.04      0.09    -0.13     0.21 1.00
                                  Bulk_ESS Tail_ESS
Intercept                             1207     1479
sigma_Intercept                       1377     1543
MotionCongruencyincongruent           2454     2417
sigma_MotionCongruencyincongruent     2710     2261

Family Specific Parameters: 
    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
ndt     0.22      0.02     0.16     0.25 1.00     1098     1217

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).

Difference between congruent and incongruent for each run

In [40]:
for (run_ in unique(rt_dat$Run)) {
    print(
        median(filter(rt_dat, MotionCongruency=="congruent" & Run==run_)$RT) - median(filter(rt_dat, MotionCongruency=="incongruent" & Run==run_)$RT)
    )
}

[1] 0.177
[1] 0.049
[1] 0.0985
[1] 0.1085
[1] 0.069
In [45]:
median(filter(rt_dat, MotionCongruency=="congruent")$RT) - median(filter(rt_dat, MotionCongruency=="incongruent")$RT)
0.111
In [41]:
pp_check(fit_exg, ndraws=100)
In [42]:
pp_check(fit_sln, ndraws=100)
In [ ]:
In [43]:
hyp_exg <- hypothesis(fit_exg, 'MotionCongruencyincongruent = 0')
hyp_exg[1]
$hypothesis =
A data.frame: 1 × 8
HypothesisEstimateEst.ErrorCI.LowerCI.UpperEvid.RatioPost.ProbStar
<chr><dbl><dbl><dbl><dbl><dbl><dbl><chr>
(MotionCongruencyincongruent) = 0-0.07707670.0127186-0.1017294-0.05214908-2.383316e-44-2.383316e-44*
In [44]:
hyp_sln <- hypothesis(fit_sln, 'MotionCongruencyincongruent = 0')
hyp_sln[1]
$hypothesis =
A data.frame: 1 × 8
HypothesisEstimateEst.ErrorCI.LowerCI.UpperEvid.RatioPost.ProbStar
<chr><dbl><dbl><dbl><dbl><dbl><dbl><chr>
(MotionCongruencyincongruent) = 0-0.20850420.0460133-0.299113-0.11923154.139881e-114.139881e-11*
In [ ]: