library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.8
## ✓ tidyr   1.2.0     ✓ stringr 1.4.0
## ✓ readr   2.1.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(here)
## here() starts at /Users/caoanjie/Desktop/projects/MB5_AP/data_analysis
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
#library(lmerTest)
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
lt_df <- read_csv(here("data/02_processed_lt_data.csv"))
## Rows: 228843 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (10): subject, stimulus_processed, complexity, gaze_location, phase, sti...
## dbl  (8): window_width, window_height, block_number, exposure_time, x, y, t,...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
proportion_df <- lt_df %>% 
  filter(phase == "pref") %>% 
  group_by(block_number, exposure_time,
           complexity, gaze_location_type, subject) %>% 
  summarise(sum_dwell_time = sum(dwell_time, na.rm = TRUE)) %>% 
  pivot_wider(names_from = gaze_location_type, 
              values_from = sum_dwell_time) %>% 
  mutate(familiar = ifelse(is.na(familiar), 0, familiar), 
         novel = ifelse(is.na(novel), 0, novel)) %>% 
  mutate(novelty_looking_proportion = novel / (familiar+novel))
## `summarise()` has grouped output by 'block_number', 'exposure_time',
## 'complexity', 'gaze_location_type'. You can override using the `.groups`
## argument.

first preregistered model

lmer(log(novel_looking_time) ~ familiarization_time * stimulus_complexity + trial_number +log(familiar_looking_time)+ (1+familiarization_time*stimulus_complexity+trial_number+log(familiar_looking_time)|participant) + (1+familiarization_time+trial_number+log(familiar_looking_time)|novel_item)+ (1+familiarization_time+trial_number+log(familiar_looking_time)|familiar_item))

data wrangling

stimulus_info_df <- lt_df %>% 
  filter(phase == "pref") %>% 
  mutate(familiar_item = case_when(
    left_stimulus == "familiar" ~ stimulus_processed_left, 
    right_stimulus == "familiar" ~ stimulus_processed_right
  ), 
        novel_item = case_when(
    left_stimulus == "novel" ~ stimulus_processed_left, 
    right_stimulus == "novel" ~ stimulus_processed_right
  )) %>% 
  select(subject, block_number, familiar_item, novel_item) %>% 
  distinct(subject, block_number, .keep_all = TRUE)
  

lt_comparison_df <- lt_df %>% 
  filter(phase == "pref") %>% 
  group_by(subject, exposure_time, block_number, complexity, gaze_location_type) %>% 
  summarise(sum_dwell_time = sum(dwell_time)) %>%
  pivot_wider(names_from = gaze_location_type, 
              values_from = sum_dwell_time) %>% 
  ungroup() %>% 
  left_join(stimulus_info_df, by = c("subject", "block_number")) %>% 
  rename(familiarization_time = exposure_time, 
         stimulus_complexity = complexity, 
         trial_number = block_number, 
        familiar_looking_time = familiar, 
         novel_looking_time = novel, 
         participant = subject) %>% 
  # currently excluding all the empty trials 
  filter(!is.na(familiar_looking_time) & !is.na(novel_looking_time))
## `summarise()` has grouped output by 'subject', 'exposure_time', 'block_number',
## 'complexity'. You can override using the `.groups` argument.
lt_comparison_df <- lt_comparison_df %>% 
  mutate(familiarization_time_scaled_centered =
           log10(familiarization_time) - log10(2000), 
         contrast_coded_complexity = ifelse(stimulus_complexity == "complex", 0.5, -0.5))

model

# failed to converge
# boundary (singular) fit: see help('isSingular')

# m_lt <- lmer(log(novel_looking_time) ~
# familiarization_time_scaled_centered * contrast_coded_complexity + trial_number +log(familiar_looking_time) +
# (1+familiarization_time_scaled_centered*contrast_coded_complexity+trial_number+log(familiar_looking_time)|participant) + (1+familiarization_time_scaled_centered+trial_number+log(familiar_looking_time)|novel_item)+
# (1+familiarization_time_scaled_centered+trial_number+log(familiar_looking_time)|familiar_item),
# data = lt_comparison_df)
# 
# # failed to converge
# # boundary (singular) fit: see help('isSingular')
# # Warning: Model failed to converge with 1 negative eigenvalue: -1.0e+00
# m2_lt <- lmer(log(novel_looking_time) ~
# familiarization_time_scaled_centered * contrast_coded_complexity + trial_number +log(familiar_looking_time) +
# (1+familiarization_time_scaled_centered*contrast_coded_complexity+trial_number+log(familiar_looking_time)|participant) + (1|novel_item)+
# (1|familiar_item),
# data = lt_comparison_df)
# 
# # Warning: Model failed to converge with max|grad| = 0.00203738 (tol = 0.002, component 1)
# 
# m3_lt <- lmer(log(novel_looking_time) ~
# familiarization_time_scaled_centered * contrast_coded_complexity + trial_number +log(familiar_looking_time) +
# (1+trial_number+log(familiar_looking_time)|participant) + (1|novel_item)+
# (1|familiar_item),
# data = lt_comparison_df)
# 
# 
# #boundary (singular) fit: see help('isSingular')
# m4_lt <- lmer(log(novel_looking_time) ~
# familiarization_time_scaled_centered * contrast_coded_complexity + trial_number +log(familiar_looking_time) +
# (1|participant) + (1|novel_item)+
# (1|familiar_item),
# data = lt_comparison_df)


m5_lt <- lmer(log(novel_looking_time) ~
familiarization_time_scaled_centered * contrast_coded_complexity + trial_number +log(familiar_looking_time) +
(1|participant),
data = lt_comparison_df)


summary(m5_lt)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(novel_looking_time) ~ familiarization_time_scaled_centered *  
##     contrast_coded_complexity + trial_number + log(familiar_looking_time) +  
##     (1 | participant)
##    Data: lt_comparison_df
## 
## REML criterion at convergence: 1073.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5270 -0.2430  0.1966  0.4778  1.9263 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept) 0.1926   0.4389  
##  Residual                0.3591   0.5992  
## Number of obs: 518, groups:  participant, 77
## 
## Fixed effects:
##                                                                 Estimate
## (Intercept)                                                     9.965399
## familiarization_time_scaled_centered                            0.118884
## contrast_coded_complexity                                       0.070776
## trial_number                                                   -0.025719
## log(familiar_looking_time)                                     -0.344000
## familiarization_time_scaled_centered:contrast_coded_complexity -0.088548
##                                                                Std. Error
## (Intercept)                                                      0.267682
## familiarization_time_scaled_centered                             0.051460
## contrast_coded_complexity                                        0.054025
## trial_number                                                     0.007959
## log(familiar_looking_time)                                       0.036369
## familiarization_time_scaled_centered:contrast_coded_complexity   0.102381
##                                                                t value
## (Intercept)                                                     37.228
## familiarization_time_scaled_centered                             2.310
## contrast_coded_complexity                                        1.310
## trial_number                                                    -3.231
## log(familiar_looking_time)                                      -9.459
## familiarization_time_scaled_centered:contrast_coded_complexity  -0.865
## 
## Correlation of Fixed Effects:
##             (Intr) fml___ cntr__ trl_nm lg(__)
## fmlrztn_t__ -0.161                            
## cntrst_cdd_ -0.037  0.004                     
## trial_numbr -0.247 -0.003 -0.046              
## lg(fmlr_l_) -0.960  0.167  0.043  0.074       
## fmlrz___:__ -0.081  0.036 -0.030 -0.021  0.088

second preregistered mode

data wrangling

p_df <- lt_comparison_df %>% 
  mutate(novelty_preference = familiar_looking_time / (familiar_looking_time + novel_looking_time), 
         item_pairs = paste0(familiar_item, "_", novel_item)) %>% 
    mutate(familiarization_time_scaled_centered =
           log10(familiarization_time) - log10(2000), 
         contrast_coded_complexity = ifelse(stimulus_complexity == "complex", 0.5, -0.5))

model

# Error: number of observations (=518) <= number of random effects (=766) for term (1 + familiarization_time_scaled_centered | item_pairs); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable

# pm_1 <- lmer(novelty_preference ~ familiarization_time_scaled_centered * contrast_coded_complexity + (1+familiarization_time_scaled_centered * contrast_coded_complexity|participant) + (1+familiarization_time_scaled_centered|item_pairs),
#      data = p_df)


#boundary (singular) fit: see help('isSingular')

# pm_2 <- lmer(novelty_preference ~ familiarization_time_scaled_centered * contrast_coded_complexity + (1+familiarization_time_scaled_centered * contrast_coded_complexity|participant) + (1|item_pairs),
#      data = p_df)


pm_3 <- lmer(novelty_preference ~ familiarization_time_scaled_centered * contrast_coded_complexity + (1|participant) + (1|item_pairs),
     data = p_df)


summary(pm_3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## novelty_preference ~ familiarization_time_scaled_centered * contrast_coded_complexity +  
##     (1 | participant) + (1 | item_pairs)
##    Data: p_df
## 
## REML criterion at convergence: -83.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.91859 -0.73039 -0.07448  0.62163  2.50285 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  item_pairs  (Intercept) 0.0044095 0.06640 
##  participant (Intercept) 0.0007197 0.02683 
##  Residual                0.0429268 0.20719 
## Number of obs: 518, groups:  item_pairs, 383; participant, 77
## 
## Fixed effects:
##                                                                Estimate
## (Intercept)                                                     0.41170
## familiarization_time_scaled_centered                           -0.08018
## contrast_coded_complexity                                      -0.02478
## familiarization_time_scaled_centered:contrast_coded_complexity -0.03382
##                                                                Std. Error
## (Intercept)                                                       0.01037
## familiarization_time_scaled_centered                              0.01805
## contrast_coded_complexity                                         0.01972
## familiarization_time_scaled_centered:contrast_coded_complexity    0.03611
##                                                                t value
## (Intercept)                                                     39.690
## familiarization_time_scaled_centered                            -4.442
## contrast_coded_complexity                                       -1.257
## familiarization_time_scaled_centered:contrast_coded_complexity  -0.937
## 
## Correlation of Fixed Effects:
##             (Intr) fml___ cntr__
## fmlrztn_t__ -0.032              
## cntrst_cdd_ -0.032 -0.014       
## fmlrz___:__ -0.014  0.015 -0.034

Model check

predictions vs mean

m5_lt <- lmer(log(novel_looking_time) ~ familiarization_time_scaled_centered * contrast_coded_complexity + trial_number +log(familiar_looking_time) + (1|participant), data = lt_comparison_df)

LT model

lt_comparison_df$pred_value <- predict(m5_lt)

lt_comparison_df %>% 
  ggplot(aes(x = log(novel_looking_time), y = pred_value)) + 
  geom_point() + 
  geom_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

### proportion model

pm_3 <- lmer(novelty_preference ~ familiarization_time_scaled_centered * contrast_coded_complexity + (1|participant) + (1|item_pairs), data = p_df)

p_df$pred_value <- predict(pm_3) 
p_df %>% 
  ggplot(aes(x = novelty_preference, y = pred_value)) + 
  geom_point() + 
  geom_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

variance

p_df %>% 
  mutate(complexity = ifelse(contrast_coded_complexity == -0.5, "complex", "simple")) %>% 
  group_by(trial_number, complexity, familiarization_time) %>% 
  summarise(var = sd(novelty_preference) ^ 2) %>% 
  ggplot(aes(x = trial_number, y = var)) + 
  geom_point() + 
  geom_smooth(method = "lm") + 
  facet_grid(complexity  ~ familiarization_time) + 
  scale_x_continuous(breaks = seq(1, 12, 1))
## `summarise()` has grouped output by 'trial_number', 'complexity'. You can
## override using the `.groups` argument.
## `geom_smooth()` using formula 'y ~ x'