suppressMessages(library(here))
suppressMessages(library(jpeg))
suppressMessages(library(grid))
suppressMessages(library(lmerTest))
suppressMessages(library(car))
suppressMessages(library(brms))
source(here::here("helper/common.R"))
##
## Attaching package: 'purrr'
## The following object is masked from 'package:car':
##
## some
##
## Attaching package: 'rlang'
## The following objects are masked from 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
## flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
## splice
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:rlang':
##
## set_names
## The following object is masked from 'package:purrr':
##
## set_names
##
## Attaching package: 'assertthat'
## The following object is masked from 'package:rlang':
##
## has_name
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3 ✓ readr 1.4.0
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x rlang::%@%() masks purrr::%@%()
## x rlang::as_function() masks purrr::as_function()
## x tidyr::expand() masks Matrix::expand()
## x tidyr::extract() masks magrittr::extract()
## x dplyr::filter() masks stats::filter()
## x rlang::flatten() masks purrr::flatten()
## x rlang::flatten_chr() masks purrr::flatten_chr()
## x rlang::flatten_dbl() masks purrr::flatten_dbl()
## x rlang::flatten_int() masks purrr::flatten_int()
## x rlang::flatten_lgl() masks purrr::flatten_lgl()
## x rlang::flatten_raw() masks purrr::flatten_raw()
## x tibble::has_name() masks assertthat::has_name(), rlang::has_name()
## x rlang::invoke() masks purrr::invoke()
## x dplyr::lag() masks stats::lag()
## x rlang::list_along() masks purrr::list_along()
## x rlang::modify() masks purrr::modify()
## x tidyr::pack() masks Matrix::pack()
## x rlang::prepend() masks purrr::prepend()
## x dplyr::recode() masks car::recode()
## x magrittr::set_names() masks rlang::set_names(), purrr::set_names()
## x purrr::some() masks car::some()
## x rlang::splice() masks purrr::splice()
## x tidyr::unpack() masks Matrix::unpack()
source(here("helper/preprocessing_helper.R"))
#knitr::opts_chunk$set(cache = TRUE, warn = FALSE, message = FALSE)
Intro
Pilot data analysis for MB2.
We can re-preprocess all of our data, this is set not to evaluate by default. In order to do this, you will need to register with Eyelink to get their binary package and then install edfR.
# labs <- dir(here::here("pilot_data"))
#
# for (lab in labs) {
# print(lab)
# source(here::here("pilot_data", lab, "import_scripts", "import.R"))
# }
File reading
labs <- dir(here::here("pilot_data"))
d <- labs %>%
map_df(function(lab) {
aoi_data <- read_csv(here(paste0("pilot_data/",
lab,"/processed_data/aoi_data.csv")))
subjects <- read_csv(here(paste0("pilot_data/",
lab,"/processed_data/subjects.csv")))
trials <- read_csv(here(paste0("pilot_data/",
lab,"/processed_data/trials.csv")))
datasets <- read_csv(here(paste0("pilot_data/",
lab,"/processed_data/datasets.csv")))
left_join(aoi_data, subjects) %>%
left_join(trials) %>%
left_join(datasets) %>%
select(lab_subject_id, lab_dataset_id, lab_trial_id, trial_id,
age, t, aoi, trial_num, error, experiment_num, sex) %>%
rename(subid = lab_subject_id,
lab = lab_dataset_id,
stimulus = lab_trial_id) %>%
mutate(aoi = case_when(grepl("target", aoi) ~ "target",
grepl("distractor", aoi) ~ "distractor",
is.na(aoi) == F ~ "other")
)
})
Exclusions
d$experiment = case_when(str_detect(d$experiment_num, "1b") ~ "pilot_1b",
d$experiment_num == "pilot_1a" ~ "pilot_1a_kid",
d$experiment_num == "pilot_1a_adult" ~ "pilot_1a_adult")
# exclude subject marked with any error
d <- d %>%
group_by(lab, subid, experiment_num) %>%
mutate(error_subj = any(error)) %>%
exclude_by(quo(error_subj), quiet=FALSE)
##
## filtering by error_subj
## [1] "This variable excludes NA trials, which is NA % of all trials."
## [1] "10 subjects, 7.4 %, have all trials where error_subj is true: exclude"
## [1] "Excluded: "
##
##
## |lab |subid |experiment_num |
## |:------------------|:-------|:--------------|
## |CEU_babylab |mb2_p14 |pilot_1a |
## |CEU_babylab |mb2_p8 |pilot_1a |
## |copenhagen_babylab |mb2p01 |pilot_1a |
## |copenhagen_babylab |mb2p05 |pilot_1a |
## |copenhagen_babylab |mb2p13 |pilot_1a |
# exclude trials under 32s (which are not complete trials)
# changed from 35s to 32 after pilot 1b because no_outcome
# trials are shorter
d <- ungroup(d) %>%
group_by(lab, trial_id, subid, experiment) %>%
mutate(time_range = (max(t) - min(t))/1000) %>%
exclude_by(quo(time_range <= 32), quiet=FALSE)
##
## filtering by time_range <= 32
## [1] "This variable excludes 25 trials, which is 2.5 % of all trials."
## [1] "0 subjects, 0 %, have all trials where time_range <= 32 is true: exclude"
## [1] "Excluded: "
##
##
## |lab | trial_id|subid |experiment |
## |:------------------|--------:|:-------|:------------|
## |CEU_babylab | 16|p1b_03 |pilot_1b |
## |CEU_babylab | 58|mb2_p17 |pilot_1a_kid |
## |CEU_babylab | 86|mb2_p5 |pilot_1a_kid |
## |copenhagen_babylab | 20|mb2p04 |pilot_1a_kid |
## |copenhagen_babylab | 22|mb2p04 |pilot_1a_kid |
# print trial time ranges by lab
ungroup(d) %>%
group_by(lab, experiment) %>%
summarise(shortest_trial=min(time_range),
longest_trial=max(time_range)) %>%
kable(digits=2)
| CEU_babylab |
pilot_1a_kid |
33.73 |
38.62 |
| CEU_babylab |
pilot_1b |
33.98 |
38.08 |
| copenhagen_babylab |
pilot_1a_kid |
32.67 |
38.62 |
| copenhagen_babylab |
pilot_1b |
35.05 |
38.08 |
| copenhagen_babylab_adults |
pilot_1a_adult |
38.23 |
38.60 |
| goettingen_babylab |
pilot_1a_kid |
38.30 |
38.73 |
| goettingen_babylab_adult |
pilot_1a_adult |
38.50 |
38.75 |
| lmu_babylab |
pilot_1a_kid |
35.62 |
38.65 |
| lmu_babylab_adult |
pilot_1a_adult |
35.90 |
38.67 |
| trento_babylab |
pilot_1a_kid |
34.70 |
38.75 |
| ubc_infantlab |
pilot_1a_kid |
36.45 |
38.70 |
| uchicago_babylab |
pilot_1a_kid |
32.42 |
38.62 |
| uchicago_babylab |
pilot_1b |
32.17 |
38.15 |
# exclude subjects who did not complete at least 6 trials (AST note: change it to 6 trials to make the criteria less strigent)
d <- ungroup(d) %>%
group_by(lab, subid, experiment) %>%
mutate(trials_completed = length(unique(trial_id))) %>%
exclude_by(quo(trials_completed < 6),quiet=FALSE) %>%
mutate(subid = paste(subid, lab, experiment, sep="_"))
##
## filtering by trials_completed < 6
## [1] "This variable excludes 21 trials, which is 2.2 % of all trials."
## [1] "6 subjects, 4.8 %, have all trials where trials_completed < 6 is true: exclude"
## [1] "Excluded: "
##
##
## |lab |subid |experiment |
## |:----------------|:------|:------------|
## |CEU_babylab |p1b_03 |pilot_1b |
## |trento_babylab |MB2_P1 |pilot_1a_kid |
## |ubc_infantlab |AC15 |pilot_1a_kid |
## |uchicago_babylab |MB2-1 |pilot_1a_kid |
## |uchicago_babylab |MB2-6 |pilot_1a_kid |
Descriptive Analysis
Descriptives
d %>%
group_by(experiment, lab, subid) %>%
summarise(age = mean(age)) %>%
summarise(n = n(),
age = mean(age)/30.25) %>%
kable(digits = 2)
| pilot_1a_adult |
copenhagen_babylab_adults |
3 |
1.06 |
| pilot_1a_adult |
goettingen_babylab_adult |
9 |
0.82 |
| pilot_1a_adult |
lmu_babylab_adult |
30 |
0.76 |
| pilot_1a_kid |
CEU_babylab |
14 |
24.83 |
| pilot_1a_kid |
copenhagen_babylab |
10 |
22.18 |
| pilot_1a_kid |
goettingen_babylab |
6 |
24.79 |
| pilot_1a_kid |
lmu_babylab |
13 |
22.99 |
| pilot_1a_kid |
trento_babylab |
4 |
20.98 |
| pilot_1a_kid |
ubc_infantlab |
10 |
21.22 |
| pilot_1a_kid |
uchicago_babylab |
8 |
23.89 |
| pilot_1b |
CEU_babylab |
3 |
25.16 |
| pilot_1b |
copenhagen_babylab |
4 |
25.43 |
| pilot_1b |
uchicago_babylab |
5 |
22.52 |
#mean age, min and max age
#pilot 1a and 1b children
d %>%
distinct(subid, .keep_all = TRUE) %>%
group_by(experiment) %>%
filter(experiment_num != "pilot_1a_adult") %>%
summarise(n = n(),
min_age = min(age, na.rm = TRUE)/30.25,
max_age = max(age, na.rm = TRUE)/30.25,
mean_age = mean(age, na.rm = TRUE)/30.25)
## # A tibble: 2 x 5
## experiment n min_age max_age mean_age
## <chr> <int> <dbl> <dbl> <dbl>
## 1 pilot_1a_kid 65 18.2 26.8 23.1
## 2 pilot_1b 12 19.1 26.0 24.1
#pilot 1 adult
d %>%
distinct(subid, .keep_all = TRUE) %>%
group_by(experiment) %>%
filter(experiment_num == "pilot_1a_adult") %>%
summarise(n = n(),
min_age = min(age, na.rm = TRUE),
max_age = max(age, na.rm = TRUE),
mean_age = mean(age, na.rm = TRUE))
## # A tibble: 1 x 5
## experiment n min_age max_age mean_age
## <chr> <int> <dbl> <dbl> <dbl>
## 1 pilot_1a_adult 42 19 53 24.1
# sex
d %>%
distinct(subid, .keep_all = TRUE) %>%
group_by(experiment, sex) %>%
tally()
## # A tibble: 10 x 3
## # Groups: experiment [3]
## experiment sex n
## <chr> <chr> <int>
## 1 pilot_1a_adult female 35
## 2 pilot_1a_adult male 5
## 3 pilot_1a_adult male/other 1
## 4 pilot_1a_adult NC- 1
## 5 pilot_1a_kid female 37
## 6 pilot_1a_kid male 28
## 7 pilot_1b F 3
## 8 pilot_1b female 3
## 9 pilot_1b M 2
## 10 pilot_1b male 4
Anticipation plot across all trials.
ms <- d %>%
group_by(t, trial_num, experiment_num) %>%
summarise(target = mean(aoi == "target", na.rm=TRUE),
distractor = mean(aoi == "distractor", na.rm=TRUE)) %>%
gather(region, looking, target, distractor)
ggplot(ms, aes(x = t, y = looking, col = region)) +
geom_line() +
geom_vline(xintercept = 120, col = "red", lty = 2) +
facet_grid(experiment_num ~ .) +
coord_cartesian(xlim = c(-4000+120, 4120)) #only consider the 4sec window before and after POD
In the primary time period of interest
ms <- d %>%
group_by(t, experiment_num) %>%
summarise(target = mean(aoi == "target", na.rm = TRUE),
distractor = mean(aoi == "distractor", na.rm = TRUE)) %>%
gather(region, looking, target, distractor)
ggplot(ms, aes(x = t, y = looking, col = region)) +
geom_point() +
xlim(-4000 + 120, 4000 + 120) +
geom_vline(xintercept = 120, col = "red", lty = 2) +
geom_text(x = -4000, y = .95, group = 1, col = "black",
label = "Anticipation", hjust = 0) +
geom_text(x = 200, y = .95, group = 1, col = "black",
label = "Reaction", hjust = 0) +
facet_grid(. ~ experiment_num)
Now, broken down by trial. Summary across anticipation window.
ms <- d %>%
filter(t > -4000, t < 120) %>%
group_by(lab, subid, trial_num, experiment_num) %>%
summarise(target = mean(aoi == "target", na.rm = TRUE)) %>%
group_by(trial_num, experiment_num) %>%
langcog::multi_boot_standard(col = "target", na.rm = TRUE)
ggplot(ms, aes(x = trial_num, y = mean)) +
geom_pointrange(aes(ymin = ci_lower, ymax = ci_upper)) +
geom_line() +
facet_grid(. ~ experiment_num)
Binned for cleaner curves
ms <- d %>%
mutate(block = ifelse(trial_num < 5, "Trials 1-4", "Trials 5-8")) %>%
group_by(t, block, experiment_num) %>%
summarise(target = mean(aoi == "target", na.rm = TRUE),
distractor = mean(aoi == "distractor", na.rm = TRUE)) %>%
gather(region, looking, target, distractor)
ggplot(ms, aes(x = t, y = looking, col = region)) +
geom_point() +
# geom_smooth(span = 2, se = FALSE) +
xlim(-4000 + 120, 4000 + 120) +
geom_vline(xintercept = 120, col = "black", lty = 3) +
annotate("text", x = -3800, y = 1, col = "black",
label = "Anticipation", hjust = 0) +
annotate("text", x = 200, y = 1, col = "black",
label = "Reaction", hjust = 0) +
ggthemes::scale_color_solarized(name = "Area of Interest") +
xlab("Time (msec)") +
ylab("Proportion gaze in AOI") +
theme(legend.position = "bottom") +
facet_wrap(experiment_num~block)
And by lab:
ms <- d %>%
mutate(block = ifelse(trial_num < 5, 1, 2)) %>%
group_by(t, lab, block, experiment_num) %>%
summarise(target = mean(aoi == "target", na.rm = TRUE),
distractor = mean(aoi == "distractor", na.rm = TRUE)) %>%
gather(region, looking, target, distractor)
ggplot(ms, aes(x = t, y = looking, col = region)) +
geom_point() +
# geom_smooth(span = 2, se = FALSE) +
xlim(-4000 + 120, 4000 + 120) +
geom_vline(xintercept = 120, col = "red", lty = 2) +
facet_grid(lab~block + experiment_num)
Set up priors for Analyses 1 and 2
For Bayesian analyses below. We are using a weakly informed prior for all analyses.
priors <-c(set_prior("normal(0, 2)", class = "Intercept"),
set_prior("normal(0, 2)", class = "b"),
set_prior("normal(0, .1)", class = "sd"),
set_prior("lkj(2)", class = "L"))
null_prior <- c(set_prior("normal(0, 2)", class = "b"),
set_prior("normal(0, .1)", class = "sd"),
set_prior("lkj(2)", class = "L"))
priors_diff <-c(set_prior("uniform(-.5, .5)", class = "Intercept"),
set_prior("normal(0, .1)", class = "b"),
set_prior("normal(0, .05)", class = "sd"),
set_prior("lkj(2)", class = "L"))
null_prior_diff <- c(set_prior("normal(0, .1)", class = "b"),
set_prior("normal(0, .05)", class = "sd"),
set_prior("lkj(2)", class = "L"))
Main analysis 1: First anticipatory look
First Look dataframe
Create a dataframe to identify the first look and see if the first look is >150ms (AT: change it to 150ms based on Dana’s comment on RR) Ref: Rayner, K., Smith, T. J., Malcolm, G. L., & Henderson, J. M. (2009). Eye Movements and Visual Encoding During Scene Perception. Psychological Science, 20(1), 6-10. doi:10.1111/j.1467-9280.2008.02243.x
# get time and location of first look (either target or distractor)
df_t_first_aoi <- d %>%
# for now only look at exp1, and at window between bear disappearing and point of disambiguation (+120 to account for time to react)
filter(., t >= -4000+120, t <= 120) %>%
# create variable for AOI in order to exploit rle, which handles NaN's neatly
mutate(aoi_recode = ifelse(aoi == "other", NA, aoi)) %>%
group_by(subid, trial_num) %>%
# 6 consecutive timepoints looking at target or distractor is the first look as we are using 150 ms
mutate(first_look_aoi = rle(aoi_recode)$values[rle(aoi_recode)$lengths >= 6][1],
first_look_t = t[sum(rle(aoi_recode)$lengths[0:(which.max(rle(aoi_recode)$lengths >= 6)-1)])+1]) %>%
ungroup()
# short df with first looks to check
first_looks_df <- df_t_first_aoi %>%
mutate(trial_num_factor = as.factor(paste0("Trial",sep = " ", trial_num))) %>%
#select(subid, trial_num, first_look_aoi, first_look_t, age, experiment_num) %>%
group_by(subid, trial_num, experiment_num,lab, trial_num_factor) %>%
summarize(aoi = unique(first_look_aoi)) %>%
mutate(aoi = case_when (is.na(aoi) ~ "other",
TRUE ~ as.character(aoi)))
Number of first look in differnet experiments, looks like pilot 1b outcome and no outcome conditions have similiar number of first look
Note that NA represents trials where the kid only looked at other places (not target/distractor) during the whole anticipatory window between -3880 to +120 ms
I checked online, didn’t find any concensus of how to analyze the “other” looking. Different people have different decisions on how to treat the “other” looking in the analyses.
I suggest removing “NA (other)” in the following analyses, as I think that looking at “distractor” and “other” are two different concepts. Theortically, we are interested in understanding whether kids looked at target more than distractor and we may not want to include the time when kids were bored and didn’t pay attention to trials/simply did not understand the task. In addition, the number of cases when kids only looked at “other” is low across experiments, so we are good.
first_looks_df %>%
group_by(aoi, experiment_num) %>%
arrange(experiment_num) %>%
count()
## # A tibble: 12 x 3
## # Groups: aoi, experiment_num [12]
## aoi experiment_num n
## <chr> <chr> <int>
## 1 distractor pilot_1a 131
## 2 distractor pilot_1a_adult 92
## 3 distractor pilot_1b_no_outcome 10
## 4 distractor pilot_1b_outcome 12
## 5 other pilot_1a 28
## 6 other pilot_1a_adult 13
## 7 other pilot_1b_no_outcome 4
## 8 other pilot_1b_outcome 1
## 9 target pilot_1a 340
## 10 target pilot_1a_adult 229
## 11 target pilot_1b_no_outcome 31
## 12 target pilot_1b_outcome 35
Data visualization across different test trials
### Pilot 1a kids
p1a_first_look <- ggplot(data = filter(first_looks_df, experiment_num == "pilot_1a")) +
geom_bar(mapping = aes(x = aoi, y = ..prop.., group = 1, stat = "count")) +
scale_y_continuous(labels = scales::percent_format()) +
facet_wrap(.~trial_num_factor, ncol = 4) +
labs(title = "Pilot 1 (Children)") +
ylab("Proportion")
ggsave(filename = "1a_first_look.jpg", p1a_first_look, width = 10, height = 7)
### Pilot 1a adult
p1adult_first_look <- ggplot(data = filter(first_looks_df, experiment_num == "pilot_1a_adult")) +
geom_bar(mapping = aes(x = aoi, y = ..prop.., group = 1, stat = "count")) +
scale_y_continuous(labels = scales::percent_format()) +
facet_wrap(.~trial_num_factor, ncol = 4)+
labs(title = "Pilot 1 adult") +
ylab("Proportion")
ggsave(filename = "1adult_first_look.jpg", p1adult_first_look, width = 10, height = 7)
### Pilot 1b
pilot_1b_bar <- first_looks_df %>%
filter(experiment_num %in% c("pilot_1b_outcome", "pilot_1b_no_outcome"))
#change labeller
exp_num_lab <- c("Pilot 2 outcome", "Pilot 2 no outcome")
names(exp_num_lab) <- c("pilot_1b_outcome", "pilot_1b_no_outcome")
trial_num_lab <- c("Trial 1", "Trial 2", "Trial 3", "Trial 4", "Trial 1", "Trial 2", "Trial 3", "Trial 4")
names(trial_num_lab) <- c("1", "2", "3", "4", "5", "6", "7", "8")
pilot_1b_first_look <- ggplot(data = pilot_1b_bar) +
geom_bar(mapping = aes(x = aoi, y = ..prop.., group = 1, stat = "count")) +
scale_y_continuous(labels = scales::percent_format()) +
facet_wrap(trial_num~experiment_num, ncol = 4,
labeller = labeller(trial_num = trial_num_lab,
experiment_num = exp_num_lab)) +
labs(title = "Pilot 2 (Children)")
ggsave(filename = "1b_first_look.jpg", pilot_1b_first_look, width = 10, height = 7)
Descriptive stats
Note that the following analyses removed NA (i.e., when kids looked at other all the time during the anticipatory window between -3880ms to +120ms)
Mean and sd of first look location in different trials. The means of all trials are above 50%
mean_sd_aoi <- first_looks_df %>%
ungroup() %>%
filter(!is.na(aoi)) %>%
mutate(aoi_dummy = ifelse(aoi == "target", 1, 0)) %>%
group_by(trial_num, experiment_num) %>%
summarise(number_of_kids_each_trial = n(),
mean(aoi_dummy),
sd(aoi_dummy))
Data visualization
ggplot(mean_sd_aoi, aes(x = trial_num, y = `mean(aoi_dummy)`)) +
geom_line()+
facet_grid(.~experiment_num)
Mixed-level analysis: First look location
Pilot 1a
first_looks_df_1a <- first_looks_df %>%
filter(experiment_num == "pilot_1a", !is.na(aoi)) %>%
mutate(aoi_dummy = ifelse(aoi == "target", 1, 0)) %>%
group_by(subid) %>%
mutate(trial_num_recoded = trial_num - 8) %>%
ungroup()
Pilot_1a_child: model_fit (Full model) (more info prior)
full_first_look_1a <- brms::brm(aoi_dummy ~ 1 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data = first_looks_df_1a,
family = bernoulli(link = "logit"),
prior = priors,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 123)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
summary(full_first_look_1a)
## Family: bernoulli
## Links: mu = logit
## Formula: aoi_dummy ~ 1 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: first_looks_df_1a (Number of observations: 499)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.23 1.00
## sd(trial_num_recoded) 0.05 0.04 0.00 0.14 1.00
## cor(Intercept,trial_num_recoded) 0.01 0.44 -0.78 0.80 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2938 2139
## sd(trial_num_recoded) 1762 2453
## cor(Intercept,trial_num_recoded) 3227 2540
##
## ~subid (Number of levels: 65)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.10 0.07 0.00 0.27 1.00
## sd(trial_num_recoded) 0.09 0.05 0.01 0.19 1.00
## cor(Intercept,trial_num_recoded) -0.13 0.43 -0.85 0.72 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2208 2038
## sd(trial_num_recoded) 1162 1838
## cor(Intercept,trial_num_recoded) 1559 2254
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.44 0.18 0.07 0.80 1.00 5324 3052
## trial_num_recoded -0.11 0.05 -0.22 -0.01 1.00 3791 2575
##
## Samples were drawn 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).
Pilot_1a_child: simpler model_fit (null model)
null_first_look_1a_full <- brm(aoi_dummy ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data = first_looks_df_1a,
family = bernoulli(link = "logit"),
prior = null_prior,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 789)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
summary(null_first_look_1a_full)
## Family: bernoulli
## Links: mu = logit
## Formula: aoi_dummy ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: first_looks_df_1a (Number of observations: 499)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.09 0.07 0.00 0.25 1.00
## sd(trial_num_recoded) 0.05 0.04 0.00 0.14 1.00
## cor(Intercept,trial_num_recoded) 0.03 0.44 -0.79 0.83 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2465 1880
## sd(trial_num_recoded) 1737 2006
## cor(Intercept,trial_num_recoded) 3538 2960
##
## ~subid (Number of levels: 65)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.10 0.07 0.00 0.27 1.00
## sd(trial_num_recoded) 0.10 0.05 0.01 0.20 1.00
## cor(Intercept,trial_num_recoded) -0.10 0.43 -0.84 0.74 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2062 2154
## sd(trial_num_recoded) 1007 1317
## cor(Intercept,trial_num_recoded) 1161 2082
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## trial_num_recoded -0.20 0.04 -0.28 -0.12 1.00 2691 2663
##
## Samples were drawn 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).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_pilot1a <- bridge_sampler(full_first_look_1a, silent = T)
marloglik_null_pilot1a <- bridge_sampler(null_first_look_1a_full, silent = T)
Bayes factor
BF_full_In_1a <- bayes_factor(marloglik_full_pilot1a, marloglik_null_pilot1a)
BF_full_In_1a
## Estimated Bayes factor in favor of x1 over x2: 1.51939
Pilot 1a adult
first_looks_df_1adult <- first_looks_df %>%
filter(experiment_num == "pilot_1a_adult", !is.na(aoi)) %>%
mutate(aoi_dummy = ifelse(aoi == "target", 1, 0)) %>%
group_by(subid) %>%
mutate(trial_num_recoded = trial_num - 8) %>%
ungroup()
Pilot_1a_adult: model_fit (Full model) (more info prior)
full_first_look_1adult <- brms::brm(aoi_dummy ~ 1 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data = first_looks_df_1adult,
family = bernoulli(link = "logit"),
prior = priors,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 123,
control = list(adapt_delta = 0.91))
summary(full_first_look_1adult)
## Family: bernoulli
## Links: mu = logit
## Formula: aoi_dummy ~ 1 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: first_looks_df_1adult (Number of observations: 334)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.22 1.00
## sd(trial_num_recoded) 0.05 0.04 0.00 0.16 1.00
## cor(Intercept,trial_num_recoded) 0.02 0.45 -0.81 0.83 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3021 1903
## sd(trial_num_recoded) 2433 2419
## cor(Intercept,trial_num_recoded) 4345 2397
##
## ~subid (Number of levels: 42)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.09 0.07 0.00 0.26 1.00
## sd(trial_num_recoded) 0.06 0.04 0.00 0.14 1.00
## cor(Intercept,trial_num_recoded) -0.02 0.45 -0.82 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2495 1594
## sd(trial_num_recoded) 1741 2184
## cor(Intercept,trial_num_recoded) 2676 2921
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.95 0.27 1.42 2.48 1.00 9514 2995
## trial_num_recoded 0.31 0.07 0.18 0.46 1.00 3889 2619
##
## Samples were drawn 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).
Pilot_1a_adult: simpler model_fit (null model)
null_first_look_1adult_full <- brm(aoi_dummy ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data = first_looks_df_1adult,
family = bernoulli(link = "logit"),
prior = null_prior,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 789,
control = list(adapt_delta = 0.99))
summary(null_first_look_1adult_full)
## Family: bernoulli
## Links: mu = logit
## Formula: aoi_dummy ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: first_looks_df_1adult (Number of observations: 334)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.37 0.07 0.24 0.51 1.00
## sd(trial_num_recoded) 0.06 0.05 0.00 0.17 1.00
## cor(Intercept,trial_num_recoded) 0.14 0.44 -0.73 0.86 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 5120 2771
## sd(trial_num_recoded) 2405 1911
## cor(Intercept,trial_num_recoded) 3984 3103
##
## ~subid (Number of levels: 42)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.09 0.07 0.00 0.25 1.00
## sd(trial_num_recoded) 0.05 0.04 0.00 0.14 1.00
## cor(Intercept,trial_num_recoded) -0.01 0.45 -0.81 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2820 1960
## sd(trial_num_recoded) 1529 1628
## cor(Intercept,trial_num_recoded) 2968 2444
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## trial_num_recoded 0.09 0.10 -0.11 0.28 1.00 2665 2634
##
## Samples were drawn 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).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_pilot1adult <- bridge_sampler(full_first_look_1adult, silent = T)
marloglik_null_pilot1adult <- bridge_sampler(null_first_look_1adult_full, silent = T)
Bayes factor
BF_full_In_1adult <- bayes_factor(marloglik_full_pilot1adult, marloglik_null_pilot1adult)
BF_full_In_1adult
## Estimated Bayes factor in favor of x1 over x2: 1077412032.60681
Pilot 1b outcome
first_looks_df_1b_o <- first_looks_df %>%
filter(experiment_num == "pilot_1b_outcome", !is.na(aoi)) %>%
mutate(aoi_dummy = ifelse(aoi == "target", 1, 0)) %>%
group_by(subid) %>%
mutate(trial_num_recoded = trial_num - 4) %>%
ungroup()
Pilot_1b_outcome: model_fit (Full model) (more info prior)
full_first_look_1b_o <- brms::brm(aoi_dummy ~ 1 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data = first_looks_df_1b_o,
family = bernoulli(link = "logit"),
prior = priors,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 123)
summary(full_first_look_1b_o)
## Family: bernoulli
## Links: mu = logit
## Formula: aoi_dummy ~ 1 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: first_looks_df_1b_o (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.22 1.00
## sd(trial_num_recoded) 0.08 0.06 0.00 0.22 1.00
## cor(Intercept,trial_num_recoded) -0.00 0.45 -0.82 0.81 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3151 1796
## sd(trial_num_recoded) 3891 2342
## cor(Intercept,trial_num_recoded) 7702 2792
##
## ~subid (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.23 1.00
## sd(trial_num_recoded) 0.08 0.06 0.00 0.23 1.00
## cor(Intercept,trial_num_recoded) -0.01 0.46 -0.83 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3906 1627
## sd(trial_num_recoded) 3306 2363
## cor(Intercept,trial_num_recoded) 7362 2854
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.69 0.52 -0.29 1.75 1.00 8633 2862
## trial_num_recoded -0.23 0.31 -0.84 0.36 1.00 8361 2809
##
## Samples were drawn 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).
Pilot_1b_outcome: simpler model_fit (null model)
null_first_look_1b_full_o <- brm(aoi_dummy ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data = first_looks_df_1b_o,
family = bernoulli(link = "logit"),
prior = null_prior,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 789)
summary(null_first_look_1b_full_o)
## Family: bernoulli
## Links: mu = logit
## Formula: aoi_dummy ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: first_looks_df_1b_o (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.22 1.00
## sd(trial_num_recoded) 0.08 0.06 0.00 0.23 1.00
## cor(Intercept,trial_num_recoded) -0.00 0.44 -0.80 0.80 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3150 2101
## sd(trial_num_recoded) 3436 2058
## cor(Intercept,trial_num_recoded) 7715 2955
##
## ~subid (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.23 1.00
## sd(trial_num_recoded) 0.08 0.06 0.00 0.22 1.00
## cor(Intercept,trial_num_recoded) -0.00 0.45 -0.81 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3808 2059
## sd(trial_num_recoded) 3119 2272
## cor(Intercept,trial_num_recoded) 7477 2770
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## trial_num_recoded -0.55 0.21 -1.00 -0.17 1.00 6536 2807
##
## Samples were drawn 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).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_pilot1b_o <- bridge_sampler(full_first_look_1b_o, silent = T)
marloglik_null_pilot1b_o <- bridge_sampler(null_first_look_1b_full_o, silent = T)
Bayes factor
BF_full_In_1b_o <- bayes_factor(marloglik_full_pilot1b_o, marloglik_null_pilot1b_o)
BF_full_In_1b_o
## Estimated Bayes factor in favor of x1 over x2: 0.54570
Pilot 1b no outcome
first_looks_df_1b_no <- first_looks_df %>%
filter(experiment_num == "pilot_1b_no_outcome", !is.na(aoi)) %>%
mutate(aoi_dummy = ifelse(aoi == "target", 1, 0)) %>%
group_by(subid) %>%
mutate(trial_num_recoded = trial_num - 4) %>%
ungroup()
Pilot_1b_no_outcome: model_fit (Full model) (more info prior)
full_first_look_1b_no <- brms::brm(aoi_dummy ~ 1 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data = first_looks_df_1b_no,
family = bernoulli(link = "logit"),
prior = priors,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 123)
summary(full_first_look_1b_no)
## Family: bernoulli
## Links: mu = logit
## Formula: aoi_dummy ~ 1 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: first_looks_df_1b_no (Number of observations: 45)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.22 1.00
## sd(trial_num_recoded) 0.08 0.06 0.00 0.22 1.00
## cor(Intercept,trial_num_recoded) -0.00 0.44 -0.81 0.81 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3226 1517
## sd(trial_num_recoded) 3055 1796
## cor(Intercept,trial_num_recoded) 6763 3025
##
## ~subid (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.23 1.00
## sd(trial_num_recoded) 0.08 0.06 0.00 0.22 1.00
## cor(Intercept,trial_num_recoded) 0.00 0.45 -0.81 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3037 1413
## sd(trial_num_recoded) 3230 2675
## cor(Intercept,trial_num_recoded) 6196 3241
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.03 0.74 -1.47 1.42 1.00 7991 2822
## trial_num_recoded 0.37 0.30 -0.20 0.98 1.00 5817 2636
##
## Samples were drawn 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).
Pilot_1b_outcome: simpler model_fit (null model)
null_first_look_1b_full_no <- brm(aoi_dummy ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data = first_looks_df_1b_no,
family = bernoulli(link = "logit"),
prior = null_prior,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 789)
summary(null_first_look_1b_full_no)
## Family: bernoulli
## Links: mu = logit
## Formula: aoi_dummy ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: first_looks_df_1b_no (Number of observations: 45)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.23 1.00
## sd(trial_num_recoded) 0.08 0.06 0.00 0.22 1.00
## cor(Intercept,trial_num_recoded) -0.02 0.45 -0.81 0.81 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3698 2247
## sd(trial_num_recoded) 3017 2071
## cor(Intercept,trial_num_recoded) 7011 2551
##
## ~subid (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.06 0.00 0.23 1.00
## sd(trial_num_recoded) 0.08 0.06 0.00 0.22 1.00
## cor(Intercept,trial_num_recoded) -0.01 0.45 -0.81 0.80 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3011 1500
## sd(trial_num_recoded) 2698 1893
## cor(Intercept,trial_num_recoded) 5762 3044
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## trial_num_recoded 0.37 0.15 0.08 0.66 1.00 5158 2917
##
## Samples were drawn 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).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
marloglik_full_pilot1b_no <- bridge_sampler(full_first_look_1b_no, silent = T)
marloglik_null_pilot1b_no <- bridge_sampler(null_first_look_1b_full_no, silent = T)
Bayes factor
BF_full_In_1b_no <- bayes_factor(marloglik_full_pilot1b_no, marloglik_null_pilot1b_no)
BF_full_In_1b_no
## Estimated Bayes factor in favor of x1 over x2: 0.34939
Main analysis 2: Differential time looking analysis
This analysis explores whether kids looked at the target more than distractor during anticipatory period.
Data visualization across different test trials
df_analysis_windowed <- d %>%
filter(t >= -4000+120 & t <= 120)
Pilot 1a kid
ggplot(data = filter(df_analysis_windowed, experiment_num == "pilot_1a")) +
geom_bar(mapping = aes(x = aoi, y = ..prop.., group = 1, stat = "count")) +
scale_y_continuous(labels = scales::percent_format()) +
facet_wrap(.~trial_num, ncol = 4)
Pilot 1a adult
ggplot(data = filter(df_analysis_windowed, experiment_num == "pilot_1a_adult")) +
geom_bar(mapping = aes(x = aoi, y = ..prop.., group = 1, stat = "count")) +
scale_y_continuous(labels = scales::percent_format()) +
facet_wrap(.~trial_num, ncol = 4)
Pilot 1b :
ggplot(data = filter(df_analysis_windowed, experiment_num %in% c("pilot_1b_outcome", "pilot_1b_no_outcome"))) +
geom_bar(mapping = aes(x = aoi, y = ..prop.., group = 1, stat = "count")) +
scale_y_continuous(labels = scales::percent_format()) +
facet_wrap(trial_num~experiment_num, ncol = 4)
For proportion score, focus on target and distractor
ms_proportion_trials <- df_analysis_windowed %>%
group_by(t, experiment_num) %>%
summarise(target = mean(aoi == "target", na.rm = TRUE),
distractor = mean(aoi == "distractor", na.rm = TRUE)) %>%
mutate(proportion_to_target = target / (target + distractor))
#change labeller
exp_num_lab <- c("Pilot 1 (Children)", "Pilot 1 (Adult)", "Pilot 2 outcome", "Pilot 2 no outcome")
names(exp_num_lab) <- c("pilot_1a", "pilot_1a_adult", "pilot_1b_outcome", "pilot_1b_no_outcome")
diff_plot <- ggplot(ms_proportion_trials,
aes(x = t, y = proportion_to_target)) +
geom_point(alpha = 0.3) +
ylim(0,1) +
# geom_smooth(span = 2, se = FALSE) +
xlim(-4000+120, 4120) +
geom_vline(xintercept = 120, col = "red", lty = 2) +
geom_text(x = -4000, y = .95, group = 1, col = "black",
label = "Anticipation", hjust = 0 ) +
geom_text(x = 200, y = .95, group = 1, col = "black",
label = "Reaction", hjust = 0 )+
facet_wrap(~experiment_num, labeller = labeller( experiment_num = exp_num_lab)) +
ylab("Proportion looking to the correct AOI (p)") +
xlab("Time (t)")
ggsave("diff_plot.jpg", diff_plot, width = 10, height = 7)
Traditional approach - Mixed level approach
ms_proportion_each_trial <- df_analysis_windowed %>%
group_by(lab, subid, experiment_num, trial_num) %>%
summarise(target = mean(aoi == "target", na.rm = TRUE),
distractor = mean(aoi == "distractor", na.rm = TRUE)) %>%
mutate(proportion_target = target/(target + distractor),
proportion_target_c = proportion_target - 0.5) %>%
filter(!is.na(proportion_target)) %>%
ungroup()
Pilot 1a
Mixed effect model for pilot 1a (kids) - full model (more info prior)
df_ms_proportion_1a <- ms_proportion_each_trial %>%
filter(experiment_num == "pilot_1a") %>%
mutate(trial_num_recoded = trial_num - 8)
m_lmer_full_1a <- brm(proportion_target_c ~ trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data=df_ms_proportion_1a,
family = gaussian(),
prior = priors_diff,
save_pars = save_pars(all = TRUE),
control = list(adapt_delta = 0.95), #get rid of divergent issues
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
seed = 1)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
m_lmer_full_1a
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: proportion_target_c ~ trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: df_ms_proportion_1a (Number of observations: 483)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.03 0.02 0.00 0.09 1.00
## sd(trial_num_recoded) 0.01 0.01 0.00 0.03 1.00
## cor(Intercept,trial_num_recoded) 0.09 0.43 -0.76 0.84 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1621 1492
## sd(trial_num_recoded) 1438 1879
## cor(Intercept,trial_num_recoded) 3140 3304
##
## ~subid (Number of levels: 65)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.02 0.00 0.09 1.00
## sd(trial_num_recoded) 0.01 0.01 0.00 0.02 1.00
## cor(Intercept,trial_num_recoded) 0.08 0.44 -0.76 0.83 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1273 1626
## sd(trial_num_recoded) 1249 1562
## cor(Intercept,trial_num_recoded) 2918 2736
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.16 0.03 0.10 0.23 1.00 5642 2704
## trial_num_recoded -0.01 0.01 -0.03 0.00 1.00 3067 2651
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.33 0.01 0.31 0.35 1.00 4850 2759
##
## Samples were drawn 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).
Null model pilot 1a (kid) - (more info prior)
m_lmer_null_1a <- brm(proportion_target_c ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data=df_ms_proportion_1a,
family = gaussian(),
prior = null_prior_diff,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.95),#need to adjust the delta value to a very high number in order to get rid of divergence messages
seed = 2)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
m_lmer_null_1a
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: proportion_target_c ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: df_ms_proportion_1a (Number of observations: 483)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.10 0.03 0.06 0.16 1.00
## sd(trial_num_recoded) 0.02 0.01 0.00 0.04 1.00
## cor(Intercept,trial_num_recoded) 0.28 0.42 -0.64 0.90 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2376 2404
## sd(trial_num_recoded) 1250 1527
## cor(Intercept,trial_num_recoded) 1975 2385
##
## ~subid (Number of levels: 65)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.02 0.00 0.09 1.00
## sd(trial_num_recoded) 0.01 0.01 0.00 0.02 1.00
## cor(Intercept,trial_num_recoded) 0.11 0.45 -0.78 0.86 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1341 1675
## sd(trial_num_recoded) 1349 1829
## cor(Intercept,trial_num_recoded) 2908 2576
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## trial_num_recoded -0.03 0.01 -0.06 -0.01 1.00 1914 2085
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.33 0.01 0.31 0.35 1.00 3894 2864
##
## Samples were drawn 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).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
prop_marloglik_full_pilot1a <- bridge_sampler(m_lmer_full_1a, silent = T)
prop_marloglik_null_pilot1a <- bridge_sampler(m_lmer_null_1a, silent = T)
Bayes factor
BF_prop_In_more_info_1a <- bayes_factor(prop_marloglik_full_pilot1a, prop_marloglik_null_pilot1a)
BF_prop_In_more_info_1a
## Estimated Bayes factor in favor of x1 over x2: 493.24825
Pilot 1a adult
Mixed effect model for pilot 1a (adult) - full model (more info prior)
df_ms_proportion_1adult <- ms_proportion_each_trial %>%
filter(experiment_num == "pilot_1a_adult") %>%
mutate(trial_num_recoded = trial_num - 8)
m_lmer_full_1adult <- brm(proportion_target_c ~ trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data=df_ms_proportion_1adult,
family = gaussian(),
prior = priors_diff,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.9999999), #failling to get rid of divergent issues
seed = 1)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
m_lmer_full_1adult
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: proportion_target_c ~ trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: df_ms_proportion_1adult (Number of observations: 324)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.03 0.02 0.00 0.09 1.00
## sd(trial_num_recoded) 0.02 0.02 0.00 0.06 1.00
## cor(Intercept,trial_num_recoded) 0.02 0.45 -0.79 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2531 1826
## sd(trial_num_recoded) 1525 2013
## cor(Intercept,trial_num_recoded) 3496 2855
##
## ~subid (Number of levels: 42)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.07 0.03 0.01 0.13 1.00
## sd(trial_num_recoded) 0.01 0.01 0.00 0.02 1.00
## cor(Intercept,trial_num_recoded) 0.00 0.43 -0.80 0.78 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1043 1069
## sd(trial_num_recoded) 1011 1346
## cor(Intercept,trial_num_recoded) 3220 2134
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.46 0.04 0.38 0.54 1.00 3525 3371
## trial_num_recoded 0.07 0.02 0.03 0.10 1.00 1162 841
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.30 0.01 0.27 0.32 1.00 2974 2514
##
## Samples were drawn 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).
Null model pilot 1a (adult) - (more info prior)
m_lmer_null_1adult <- brm(proportion_target_c ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data=df_ms_proportion_1adult,
family = gaussian(),
prior = null_prior_diff,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.9995),#need to adjust the delta value to a very high number in order to get rid of divergence messages
seed = 2)
m_lmer_null_1adult
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: proportion_target_c ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: df_ms_proportion_1adult (Number of observations: 324)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.18 0.03 0.14 0.24 1.00
## sd(trial_num_recoded) 0.02 0.02 0.00 0.06 1.00
## cor(Intercept,trial_num_recoded) 0.10 0.46 -0.75 0.86 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2387 2605
## sd(trial_num_recoded) 2057 2274
## cor(Intercept,trial_num_recoded) 3340 2792
##
## ~subid (Number of levels: 42)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.08 0.03 0.02 0.13 1.00
## sd(trial_num_recoded) 0.01 0.01 0.00 0.02 1.00
## cor(Intercept,trial_num_recoded) 0.03 0.44 -0.79 0.81 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1282 770
## sd(trial_num_recoded) 1252 1424
## cor(Intercept,trial_num_recoded) 4150 2943
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## trial_num_recoded 0.05 0.03 -0.02 0.10 1.00 2178 1492
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.30 0.01 0.27 0.32 1.01 4194 2272
##
## Samples were drawn 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).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
prop_marloglik_full_pilot1adult <- bridge_sampler(m_lmer_full_1adult, silent = T)
prop_marloglik_null_pilot1adult <- bridge_sampler(m_lmer_null_1adult, silent = T)
Bayes factor
BF_prop_In_more_info_adult <- bayes_factor(prop_marloglik_full_pilot1adult, prop_marloglik_null_pilot1adult)
BF_prop_In_more_info_adult
## Estimated Bayes factor in favor of x1 over x2: 4859188.92423
Pilot 1b outcome
Mixed effect model for pilot 1b (kids) - full model (more info prior)
df_ms_proportion_1b_o <- ms_proportion_each_trial %>%
filter(experiment_num == "pilot_1b_outcome") %>%
mutate(trial_num_recoded = trial_num - 4)
m_lmer_full_1b_o <- brm(proportion_target_c ~ trial_num_recoded +
(trial_num_recoded | lab) +
(trial_num_recoded | subid),
data=df_ms_proportion_1b_o,
family = gaussian(),
prior = priors_diff,
save_pars = save_pars(all = TRUE),
warmup = 20000,
iter = 40000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.999), #failing to get rid of divergent issues
seed = 1)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
m_lmer_full_1b_o
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: proportion_target_c ~ trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: df_ms_proportion_1b_o (Number of observations: 47)
## Samples: 4 chains, each with iter = 40000; warmup = 20000; thin = 1;
## total post-warmup samples = 80000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.03 0.03 0.00 0.10 1.00
## sd(trial_num_recoded) 0.03 0.02 0.00 0.09 1.00
## cor(Intercept,trial_num_recoded) 0.03 0.45 -0.80 0.83 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 62561 38625
## sd(trial_num_recoded) 55742 41248
## cor(Intercept,trial_num_recoded) 119343 57166
##
## ~subid (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.03 0.00 0.10 1.00
## sd(trial_num_recoded) 0.03 0.02 0.00 0.07 1.00
## cor(Intercept,trial_num_recoded) 0.04 0.45 -0.79 0.83 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 59937 39754
## sd(trial_num_recoded) 61221 39142
## cor(Intercept,trial_num_recoded) 122704 55170
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.13 0.09 -0.05 0.30 1.00 122668 61400
## trial_num_recoded -0.06 0.05 -0.15 0.04 1.00 106381 60582
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.36 0.04 0.30 0.45 1.00 124020 56495
##
## Samples were drawn 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).
Null model pilot 1b (kid) - (more info prior)
m_lmer_null_1b_o <- brm(proportion_target_c ~ 0 + trial_num_recoded +
(trial_num_recoded | lab) +
(trial_num_recoded | subid),
data=df_ms_proportion_1b_o,
family = gaussian(),
prior = null_prior_diff,
save_pars = save_pars(all = TRUE),
warmup = 20000,
iter = 40000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.9999),#need to adjust the delta value to a very high number in order to get rid of divergence messages
seed = 2)
m_lmer_null_1b_o
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: proportion_target_c ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: df_ms_proportion_1b_o (Number of observations: 47)
## Samples: 4 chains, each with iter = 40000; warmup = 20000; thin = 1;
## total post-warmup samples = 80000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.03 0.00 0.10 1.00
## sd(trial_num_recoded) 0.03 0.02 0.00 0.09 1.00
## cor(Intercept,trial_num_recoded) 0.03 0.45 -0.80 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 63132 41703
## sd(trial_num_recoded) 54172 39512
## cor(Intercept,trial_num_recoded) 104586 56116
##
## ~subid (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.03 0.00 0.10 1.00
## sd(trial_num_recoded) 0.03 0.02 0.00 0.07 1.00
## cor(Intercept,trial_num_recoded) 0.04 0.45 -0.80 0.83 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 56295 37568
## sd(trial_num_recoded) 60518 38750
## cor(Intercept,trial_num_recoded) 113682 56104
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## trial_num_recoded -0.10 0.04 -0.17 -0.02 1.00 78601 55684
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.37 0.04 0.30 0.46 1.00 121232 57643
##
## Samples were drawn 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).
Bayes factor
BF_prop_In_more_info_1b_o <- bayes_factor(m_lmer_full_1b_o, m_lmer_null_1b_o)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
BF_prop_In_more_info_1b_o
## Estimated Bayes factor in favor of m_lmer_full_1b_o over m_lmer_null_1b_o: 0.63439
Pilot 1b no outcome
Mixed effect model for pilot 1b (no outcome) - full model (more info prior)
df_ms_proportion_1b_no <- ms_proportion_each_trial %>%
filter(experiment_num == "pilot_1b_no_outcome") %>%
mutate(trial_num_recoded = trial_num - 4)
m_lmer_full_1b_no <- brm(proportion_target_c ~ trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data=df_ms_proportion_1b_no,
family = gaussian(),
prior = priors_diff,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.999),
seed = 1)
m_lmer_full_1b_no
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: proportion_target_c ~ trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: df_ms_proportion_1b_no (Number of observations: 43)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.03 0.00 0.10 1.00
## sd(trial_num_recoded) 0.03 0.02 0.00 0.09 1.00
## cor(Intercept,trial_num_recoded) -0.04 0.45 -0.84 0.83 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3566 2315
## sd(trial_num_recoded) 2587 1948
## cor(Intercept,trial_num_recoded) 4936 2934
##
## ~subid (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.03 0.00 0.11 1.00
## sd(trial_num_recoded) 0.03 0.02 0.00 0.07 1.00
## cor(Intercept,trial_num_recoded) -0.01 0.44 -0.80 0.79 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2949 1920
## sd(trial_num_recoded) 2571 2070
## cor(Intercept,trial_num_recoded) 4895 2746
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.10 0.12 -0.13 0.36 1.00 6422 2962
## trial_num_recoded 0.04 0.05 -0.06 0.13 1.00 5310 2922
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.36 0.04 0.29 0.46 1.00 3992 2659
##
## Samples were drawn 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).
Null model pilot 1b (no outcome) - (more info prior)
m_lmer_null_1b_no <- brm(proportion_target_c ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid),
data=df_ms_proportion_1b_no,
family = gaussian(),
prior = null_prior_diff,
save_pars = save_pars(all = TRUE),
warmup = 1000,
iter = 2000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.9999),#need to adjust the delta value to a very high number in order to get rid of divergence messages
seed = 2)
m_lmer_null_1b_no
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: proportion_target_c ~ 0 + trial_num_recoded + (trial_num_recoded | lab) + (trial_num_recoded | subid)
## Data: df_ms_proportion_1b_no (Number of observations: 43)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~lab (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.03 0.00 0.10 1.00
## sd(trial_num_recoded) 0.03 0.02 0.00 0.08 1.00
## cor(Intercept,trial_num_recoded) -0.04 0.45 -0.81 0.79 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2982 1527
## sd(trial_num_recoded) 2804 1924
## cor(Intercept,trial_num_recoded) 5045 2994
##
## ~subid (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.04 0.03 0.00 0.11 1.00
## sd(trial_num_recoded) 0.03 0.02 0.00 0.07 1.00
## cor(Intercept,trial_num_recoded) -0.01 0.45 -0.81 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2616 2129
## sd(trial_num_recoded) 2659 1984
## cor(Intercept,trial_num_recoded) 5263 2878
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## trial_num_recoded 0.07 0.03 0.01 0.13 1.00 3765 3129
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.36 0.04 0.28 0.45 1.00 6162 2754
##
## Samples were drawn 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).
Calculate the marginal log likelihoods for each hypothese using Bridge sampling
prop_marloglik_full_pilot1b_no <- bridge_sampler(m_lmer_full_1b_no, silent = T)
prop_marloglik_null_pilot1b_no <- bridge_sampler(m_lmer_null_1b_no, silent = T)
Bayes factor
BF_prop_In_more_info_1b_no <- bayes_factor(prop_marloglik_full_pilot1b_no, prop_marloglik_null_pilot1b_no)
BF_prop_In_more_info_1b_no
## Estimated Bayes factor in favor of x1 over x2: 0.47826
Mirman growth curve analysis
Create dataset with sub-sampled counts for growth curve analysis a la Mirman (2014).