Load Data

Load pilot data (from Julien Mayor’s lab).

d1 <- read_csv(here("pilot/data/first_session_babies.csv")) %>%
  rename(sd_LT_incongruent_trials = sd_LT_congruent_trials_1) 
## Warning: Duplicated column names deduplicated: 'sd_LT_congruent_trials' =>
## 'sd_LT_congruent_trials_1' [18]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_character(),
##   age_months = col_double(),
##   trial = col_double(),
##   test_trial_audio_length_ms = col_double(),
##   looking_time = col_double(),
##   mean_LT_congruent_trials = col_double(),
##   sd_LT_congruent_trials = col_double(),
##   mean_LT_incongruent_trials = col_double(),
##   sd_LT_congruent_trials_1 = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
d2 <- read_csv(here("pilot/data/second_session_babies.csv")) %>%
  rename(sd_LT_incongruent_trials = sd_LT_congruent_trials_1,
         looking_time = `looking time`) 
## Warning: Duplicated column names deduplicated: 'sd_LT_congruent_trials' =>
## 'sd_LT_congruent_trials_1' [18]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_character(),
##   age_months = col_double(),
##   trial = col_double(),
##   test_trial_audio_length_ms = col_double(),
##   `looking time` = col_double(),
##   mean_LT_congruent_trials = col_double(),
##   sd_LT_congruent_trials = col_double(),
##   mean_LT_incongruent_trials = col_double(),
##   sd_LT_congruent_trials_1 = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
cdat <- bind_rows(d1, d2) %>%
  filter(!is.na(id)) %>%
  mutate(labID = 1)

Loaded data from 19 participants. Their age distribution is shown below.

cdat %>% group_by(id, age_months) %>% 
  summarise(n=n()) %>%
  group_by(age_months) %>%
  summarise(Nch=n()) %>% kable()
## `summarise()` has grouped output by 'id'. You can override using the `.groups` argument.
age_months Nch
8 7
9 2
10 2
11 4
12 4

Plot

# test data
dat <- cdat %>% select(1:14) # remaining columns are averages and notes

fam_dat <- dat %>% filter(is.na(test_trial)) %>% 
  select(-test_trial, -test_trial_audio, -test_trial_audio_length_ms, 
         -test_trial_type, -test_version)

test_dat <- dat %>% filter(!is.na(test_trial)) %>% 
  select(-fam_trial_audio, -fam_version) %>%
  separate(test_trial_type, sep='_', 
           into = c("test_trial_type_ch", "test_trial_type_num")) %>%
  mutate(trial_type = ifelse(test_trial_type_ch == fam_condition, "same", "different"),
         log_looking_time = log(looking_time))

dag <- test_dat %>% group_by(id, fam_condition, age_months, trial_type) %>%
  summarise(looking_time = mean(log_looking_time)) %>% 
  group_by(trial_type, age_months) %>%
  tidyboot::tidyboot_mean(looking_time) # quite slow..
## `summarise()` has grouped output by 'id', 'fam_condition', 'age_months'. You can override using the `.groups` argument.
pos = position_dodge(width=.2)
ggplot(dag, aes(x=age_months, y=mean, group=trial_type, color=trial_type)) + 
  geom_point(aes(y=mean, x=age_months), pos=pos) + 
  ylab("log(looking time)") + xlab("Age (months)") + 
  geom_linerange(aes(ymin=ci_lower, ymax=ci_upper), pos=pos) + 
  theme_bw() + geom_smooth(method="lm")
## `geom_smooth()` using formula 'y ~ x'

Regression

Preregistered analysis, but without per-lab random effects and procedure fixed effects: Infants’ log(looking time) (DV) ~ 1 + familiarization order (ABB vs ABA) * trial_type + age * trial_type (same rule vs different rule at test) + trial_num * trial_type + (trial_num*trial_type | subject)

m1 <- lmer(log_looking_time ~ 1 + 
             fam_condition * trial_type +
             age_months * trial_type +
             trial * trial_type + (trial * trial_type | id), data=test_dat)
## boundary (singular) fit: see ?isSingular
# singular, pruning random effects

m2 <- lmer(log_looking_time ~ 1 + 
             fam_condition * trial_type +
             age_months * trial_type +
             trial * trial_type + (1 | id), data=test_dat)
tab_model(m2)
  log looking time
Predictors Estimates CI p
(Intercept) 8.30 7.14 – 9.45 <0.001
fam_condition [ABB] 0.33 -0.04 – 0.70 0.078
trial_type [same] -0.37 -1.29 – 0.54 0.424
age_months 0.00 -0.11 – 0.12 0.960
trial -0.04 -0.07 – -0.01 0.006
fam_condition [ABB] *
trial_type [same] 		-0.11 -0.39 – 0.18 0.459
trial_type [same] *
age_months 		0.03 -0.06 – 0.12 0.494
trial_type [same] * trial 0.04 -0.00 – 0.08 0.085
Random Effects
σ2 0.29
Ï„00 id 0.12
ICC 0.29
N id 19
Observations 228
Marginal R2 / Conditional R2 0.082 / 0.347 
plot_model(m2, sort.est = T) + theme_bw()

Bayesian Regression

bm1 <- brm(log_looking_time ~ 1 + 
             fam_condition * trial_type +
             age_months * trial_type +
             trial * trial_type + (1 | id), 
           data = test_dat, 
           warmup = 1000, iter = 2000, chains = 4, 
           inits = "random", cores = 4, seed=42) 
## 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(bm1)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: log_looking_time ~ 1 + fam_condition * trial_type + age_months * trial_type + trial * trial_type + (1 | id) 
##    Data: test_dat (Number of observations: 228) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~id (Number of levels: 19) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.37      0.09     0.23     0.58 1.00     1268     1880
## 
## Population-Level Effects: 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                           8.30      0.62     7.10     9.55 1.00
## fam_conditionABB                    0.33      0.20    -0.05     0.74 1.00
## trial_typesame                     -0.39      0.48    -1.34     0.52 1.00
## age_months                          0.00      0.06    -0.13     0.12 1.00
## trial                              -0.04      0.02    -0.07    -0.01 1.00
## fam_conditionABB:trial_typesame    -0.11      0.14    -0.39     0.18 1.00
## trial_typesame:age_months           0.03      0.05    -0.06     0.12 1.00
## trial_typesame:trial                0.04      0.02    -0.00     0.08 1.00
##                                 Bulk_ESS Tail_ESS
## Intercept                           1209     2027
## fam_conditionABB                    1129     1601
## trial_typesame                      1845     2246
## age_months                          1148     2165
## trial                               2472     2906
## fam_conditionABB:trial_typesame     3834     2877
## trial_typesame:age_months           1759     2562
## trial_typesame:trial                2535     2485
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.54      0.03     0.49     0.60 1.00     4248     2664
## 
## 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).

Generate posterior predictive values.

# posterior predictive
test_dat %>%
  data_grid(id, fam_condition, trial_type, age_months, trial) %>%
  add_fitted_draws(bm1) %>%
  ggplot(aes(x = .value, y = trial_type)) +
  facet_grid(. ~ fam_condition) +
  stat_pointinterval(.width = c(.66, .95)) + 
  theme_bw()