# 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)
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")
)
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)
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)