Nothing
## ---- 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)
# load data
data <- flanker
# analyse only correct answers, load test and control data
control_rt <- data %>% filter(result == "correct" &
group == "control")
test_rt <- data %>% filter(result == "correct" &
group == "test")
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# control group subject indexes range is 22..45 cast it to 1..23
# test group indexes are OK
control_rt$subject <- control_rt$subject - 21
## ---- message=FALSE, warning=FALSE, results = 'hide'--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# prior
uniform_prior <- b_prior(family="uniform", pars=c(0, 3))
# attach priors to relevant parameters
priors <- list(c("mu_m", uniform_prior))
# fit
rt_control_fit <- b_reaction_time(t=control_rt$rt,
s=control_rt$subject,
priors=priors,
chains=1, iter=200, warmup=100)
rt_test_fit <- b_reaction_time(t=test_rt$rt,
s=test_rt$subject,
priors=priors,
chains=1, iter=200, warmup=100)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# plot trace
plot_trace(rt_control_fit)
plot_trace(rt_test_fit)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# the commands below are commented out for the sake of brevity
#print(rt_control_fit)
#print(rt_test_fit)
## ---- eval=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# rt_control_fit <- b_reaction_time(t=control_rt$rt,
# s=control_rt$subject,
# iter=4000)
#
# rt_test_fit <- b_reaction_time(t=test_rt$rt,
# s=test_rt$subject,
# iter=4000)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# check fits
plot(rt_control_fit)
plot(rt_test_fit)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# set rope (region of practical equivalence) interval to +/- 10ms
rope <- 0.01
# which mean is smaller/larger
rt_control_test <- compare_means(rt_control_fit, fit2=rt_test_fit, rope=rope)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# difference plot
plot_means_difference(rt_control_fit, fit2=rt_test_fit, rope=rope)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# visual comparsion of mean difference
plot_means(rt_control_fit, fit2=rt_test_fit)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# map correct/incorrect/timeout to 1 (success) or 0 (fail)
data$result_numeric <- 0
data[data$result == "correct", ]$result_numeric <- 1
# split to control and test groups
control_sr <- data %>% filter(group == "control")
test_sr <- data %>% filter(group == "test")
# control group subject indexes range is 22..45 cast it to 1..23
# test group indexes are OK
control_sr$subject <- control_sr$subject - 21
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# priors
p_prior <- b_prior(family="beta", pars=c(1, 1))
# attach priors to relevant parameters
priors <- list(c("p", p_prior))
# fit
sr_control_fit <- b_success_rate(r=control_sr$result_numeric,
s=control_sr$subject,
priors=priors,
chains=1, iter=200, warmup=100)
sr_test_fit <- b_success_rate(r=test_sr$result_numeric,
s=test_sr$subject,
priors=priors,
chains=1, iter=200, warmup=100)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# inspect control group fit
plot_trace(sr_control_fit)
plot(sr_control_fit, subjects=FALSE)
print(sr_control_fit)
# inspect test group fit
plot_trace(sr_test_fit)
plot(sr_test_fit, subjects=FALSE)
#print(sr_test_fit)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# which one is higher
sr_control_test <- compare_means(sr_control_fit, fit2=sr_test_fit)
## ---- message=FALSE, warning=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# difference plot
plot_means_difference(sr_control_fit, fit2=sr_test_fit)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.