pkbb_ABBA_analysis

rm(list = ls())

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(performance)

source("process_pairs.R")
source("process_ABBA.R")

1 load data

exp2_adult <- read_csv("adult_stimuli_type.csv")

Rows: 44775 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (8): prolific_id, trial_type, block_type, background_stimulus, deviant_s...
dbl (4): total_rt, block_number, trial_number, total_trial_number

ℹ 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.

# log LT
exp2_adult$total_rt_log <- log(exp2_adult$total_rt)

# create bg_image, dv_image (contain all four info: pose, identity id, number, animacy)
exp2_adult <- exp2_adult %>%
  mutate(
    bg_image = paste(background_type, str_extract(background_stimulus, "\\d+"), sep = "_"),
    dv_image = paste(deviant_type, str_extract(deviant_stimulus, "\\d+"), sep = "_")
  ) %>%
  mutate(dv_image = ifelse(violation_type == "background", bg_image, dv_image))

# calculate dishab
exp2_test_only <- exp2_adult %>% filter(trial_number == total_trial_number) %>% mutate(my_trial_type = "test")
exp2_last_fam_only <- exp2_adult %>% filter(trial_number == total_trial_number - 1) %>% mutate(my_trial_type = "last_fam")
exp2_first_fam_only <- exp2_adult %>% filter(trial_number == 1) %>% mutate(my_trial_type = "first_fam")

exp2_core <- pivot_wider(bind_rows(exp2_first_fam_only, exp2_last_fam_only, exp2_test_only),
                         id_cols = c(prolific_id, block_number, violation_type, total_trial_number, bg_image, dv_image),
                         names_from = my_trial_type,
                         values_from = total_rt_log) %>%
  mutate(dishab = test - last_fam) %>%
  mutate(onefam = ifelse(total_trial_number == 2, 1, 0))

1.1 data quality check

1.1.1 raw data

# DataExplorer::create_report(exp2_adult)
DataExplorer::plot_histogram(exp2_adult)

DataExplorer::plot_missing(exp2_adult)

DataExplorer::plot_qq(exp2_adult)

skimr::skim(exp2_adult)

Data summary

Name	exp2_adult
Number of rows	44775
Number of columns	15
_______________________
Column type frequency:
character	10
numeric	5
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
prolific_id	0	1.00	24	24	0	468	0
trial_type	0	1.00	7	10	0	2	0
block_type	0	1.00	13	16	0	2	0
background_stimulus	0	1.00	39	42	0	64	0
deviant_stimulus	14931	0.67	39	42	0	64	0
background_type	0	1.00	17	22	0	8	0
deviant_type	14931	0.67	17	22	0	8	0
violation_type	0	1.00	4	10	0	5	0
bg_image	0	1.00	21	26	0	128	0
dv_image	0	1.00	21	26	0	128	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
total_rt	0	1	2347.60	2161.48	500.00	949.00	1568.00	2981.5	20526.00	▇▁▁▁▁
block_number	0	1	12.49	6.92	1.00	6.00	13.00	18.0	24.00	▇▇▆▇▇
trial_number	0	1	2.84	1.57	1.00	2.00	3.00	4.0	6.00	▇▂▂▂▂
total_trial_number	0	1	4.67	1.49	2.00	4.00	6.00	6.0	6.00	▃▁▅▁▇
total_rt_log	0	1	7.46	0.74	6.21	6.86	7.36	8.0	9.93	▇▇▆▂▁

janitor::get_dupes(exp2_adult)

No variable names specified - using all columns.

No duplicate combinations found of: prolific_id, trial_type, total_rt, block_type, block_number, background_stimulus, deviant_stimulus, background_type, deviant_type, ... and 6 other variables

# A tibble: 0 × 16
# ℹ 16 variables: prolific_id <chr>, trial_type <chr>, total_rt <dbl>,
#   block_type <chr>, block_number <dbl>, background_stimulus <chr>,
#   deviant_stimulus <chr>, background_type <chr>, deviant_type <chr>,
#   violation_type <chr>, trial_number <dbl>, total_trial_number <dbl>,
#   total_rt_log <dbl>, bg_image <chr>, dv_image <chr>, dupe_count <int>

1.1.2 dishab data

output:
- no duplicates
- missing data: due to missing trials in the raw data
- first_fam (0.55%)
- last_fam (0.35%)
- test (0.22%)
- dishab (0.55%) * a coincidence that number of missing trial for first_fam and dishab is the same

DataExplorer::plot_histogram(exp2_core)

DataExplorer::plot_missing(exp2_core)

DataExplorer::plot_qq(exp2_core)

Warning: Removed 188 rows containing non-finite outside the scale range
(`stat_qq()`).

Warning: Removed 188 rows containing non-finite outside the scale range
(`stat_qq_line()`).

skimr::skim(exp2_core)

Data summary

Name	exp2_core
Number of rows	11232
Number of columns	11
_______________________
Column type frequency:
character	4
numeric	7
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
prolific_id	0	1	24	24	0	468	0
violation_type	0	1	4	10	0	5	0
bg_image	0	1	21	26	0	128	0
dv_image	0	1	21	26	0	128	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
block_number	0	1.00	12.50	6.92	1.00	6.75	12.50	18.25	24.00	▇▇▆▇▇
total_trial_number	0	1.00	4.00	1.64	2.00	2.00	4.00	6.00	6.00	▇▁▇▁▇
first_fam	62	0.99	7.72	0.73	6.22	7.15	7.68	8.24	9.92	▃▇▇▃▁
last_fam	39	1.00	7.46	0.74	6.22	6.84	7.35	8.00	9.92	▇▇▆▂▁
test	25	1.00	7.44	0.72	6.21	6.85	7.33	7.96	9.93	▇▇▆▂▁
dishab	62	0.99	-0.02	0.50	-3.20	-0.23	0.00	0.21	2.95	▁▁▇▁▁
onefam	0	1.00	0.34	0.47	0.00	0.00	0.00	1.00	1.00	▇▁▁▁▅

2 unique pairs

## unique combinations of stimuli pairs in test trials
exp2_unique_stimuli_pairs <- exp2_test_only %>%
  group_by(bg_image, dv_image, violation_type) %>%
  summarise(count = n(), .groups = "drop") %>%
  mutate(pairID = row_number())

### number of unique combinations of deviant trials used in study:
nrow(filter(exp2_unique_stimuli_pairs, violation_type != "background"))

[1] 2674

In theory there are 1282+1282+12815+1282=2688 unique pairs

### number of background trials used in study
nrow(filter(exp2_unique_stimuli_pairs, violation_type == "background"))

[1] 128

In theory there are 128 background trial stimuli.

2.1 pose

within_POP_ID: 1 (left to right), 2 (right to left)

exp2_pairs_pose <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "pose")) 

2.2 number

within_POP_ID: 1 (pair to single), 2 (single to pair)

exp2_pairs_number <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "number")) 

2.3 animacy

within_POP_ID: 1 (animate to inanimate), 2 (inanimate to animate)

## no ignored types
exp2_pairs_animacy <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "animacy"))
check_single_within_PoP(exp2_pairs_animacy) # 484

[1] 484

## combing across pose
exp2_pairs_animacy_nopose <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "animacy"), ignore_violation_types = c("pose"))
check_single_within_PoP(exp2_pairs_animacy_nopose) # 126

[1] 126

## combing across identity
exp2_pairs_animacy_noidentity <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "animacy"), ignore_violation_types = c("identity"))
check_single_within_PoP(exp2_pairs_animacy_noidentity) # 0, but there are only 4 pair_of_pair_IDs. 

[1] 0

## combining across pose and number
exp2_pairs_animacy_nopose_nonumber <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "animacy"), ignore_violation_types = c("pose", "number"))
check_single_within_PoP(exp2_pairs_animacy_nopose_nonumber) # 9

[1] 9

2.4 identity (archive)

## no ignored types
exp2_pairs_identity <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "identity"))

Warning in process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == :
NAs introduced by coercion

check_single_within_PoP(exp2_pairs_identity) # 440

[1] 440

## combing across pose
exp2_pairs_identity_nopose <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "identity"), ignore_violation_types = c("pose"))

Warning in process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == :
NAs introduced by coercion

check_single_within_PoP(exp2_pairs_identity_nopose) # 116

[1] 116

## combing across pose and number
exp2_pairs_identity_nopose_nonumber <- process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == "identity"), ignore_violation_types = c("pose", "number"))

Warning in process_pairs(filter(exp2_unique_stimuli_pairs, violation_type == :
NAs introduced by coercion

check_single_within_PoP(exp2_pairs_identity_nopose_nonumber) # 6

[1] 6

3 ABBA: POSE

ABBA_pose <- exp2_core %>% filter(violation_type == "pose") %>%
  left_join(exp2_pairs_pose, by = c("violation_type", "bg_image", "dv_image"))

3.1 DA check

# data quality check
# DataExplorer::create_report(ABBA_pose)
DataExplorer::plot_histogram(ABBA_pose)

DataExplorer::plot_missing(ABBA_pose)

DataExplorer::plot_qq(ABBA_pose)

Warning: Removed 26 rows containing non-finite outside the scale range
(`stat_qq()`).

Warning: Removed 26 rows containing non-finite outside the scale range
(`stat_qq_line()`).

ABBA_pose <- ABBA_pose %>% filter(!is.na(first_fam) & !is.na(dishab))

3.1.1 summary of mean dishab

ABBA_pose %>% group_by(change_direction) %>% rstatix::get_summary_stats(dishab, type = "common")

# A tibble: 2 × 11
  change_direction variable     n   min   max median   iqr   mean    sd    se
  <chr>            <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl> <dbl>
1 left to right    dishab     939 -2.30  2.92 -0.029 0.394 -0.074 0.453 0.015
2 right to left    dishab     917 -2.54  2.52 -0.028 0.433 -0.078 0.526 0.017
# ℹ 1 more variable: ci <dbl>

ABBA_pose %>% group_by(change_direction, total_trial_number) %>% rstatix::get_summary_stats(dishab, type = "common")

# A tibble: 6 × 12
  total_trial_number change_direction variable     n   min   max median   iqr
               <dbl> <chr>            <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl>
1                  2 left to right    dishab     320 -2.30 0.962 -0.238 0.588
2                  4 left to right    dishab     307 -1.97 1.56   0.008 0.345
3                  6 left to right    dishab     312 -1.79 2.92   0.02  0.281
4                  2 right to left    dishab     324 -2.54 1.32  -0.228 0.663
5                  4 right to left    dishab     287 -1.89 1.78   0.018 0.329
6                  6 right to left    dishab     306 -1.61 2.52   0.043 0.355
# ℹ 4 more variables: mean <dbl>, sd <dbl>, se <dbl>, ci <dbl>

ABBA_pose %>%
    ggplot(aes(x = as.factor(change_direction), y = dishab)) +
    geom_boxplot(outlier.color = "red", outlier.shape = 1, outlier.alpha = 0.2) +
    geom_point(alpha = 0.1) +
  stat_summary(fun = "mean", geom = "point", position = position_dodge(.9), color = "red") +
    stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.2, position = position_dodge(.9), color = "red") +
    facet_grid(~total_trial_number) +
    labs(x = "Direction of Change",
         y = "Dishabituation (log(test) - log(last fam))",
         title = "Dishabituation by different pairs of stimuli (pose change)") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme_bw()

ABBA_pose$change_direction_dc <- as.factor(ABBA_pose$change_direction)
contrasts(ABBA_pose$change_direction_dc) <- contr.sum(2)/2

ABBA_pose$onefam_dc <- as.factor(ABBA_pose$onefam)

ABBA_pose$total_trial_number_c <- scale(ABBA_pose$total_trial_number, center = TRUE, scale = FALSE)
ABBA_pose$block_number_c <- scale(ABBA_pose$block_number, center = TRUE, scale = FALSE)

ABBA_pose$total_trial_number_log <- log(ABBA_pose$total_trial_number)
ABBA_pose$block_number_log <- log(ABBA_pose$block_number)

ABBA_pose$total_trial_number_log_c <- scale(ABBA_pose$total_trial_number_log, center = TRUE, scale = FALSE)
ABBA_pose$block_number_log_c <- scale(ABBA_pose$block_number_log, center = TRUE, scale = FALSE)

3.2 plots

3.2.1 raw dishab plots

plot_raw_dishab(dataset = ABBA_pose)

plot_raw_dishab(dataset = ABBA_pose, by_trialnumber = TRUE)

3.2.2 diffemean plot

sumstats_ABBA_pose_alltrials <- calculate_sumstats(ABBA_pose)
plot_diffmean(sumstats = sumstats_ABBA_pose_alltrials, title_suffix = "(all trials)")

get_stimuli_for_diffmean(sumstats_ABBA_pose_alltrials, exp2_pairs_pose, "pose")

# A tibble: 64 × 14
   pair_of_pair_ID pair_of_pair.x       mean_diff pooled_sd    n1    n2 ci_lower
             <int> <chr>                    <dbl>     <dbl> <int> <int>    <dbl>
 1              39 inanimate_pair_left…     0.354     0.449    18    13  0.196  
 2              21 animate_single_left…     0.305     0.508    13    16  0.120  
 3              17 animate_single_left…     0.284     0.761    10    17 -0.00359
 4              64 inanimate_single_le…     0.272     0.739    11    16 -0.00634
 5              26 animate_single_left…     0.225     0.603    12    11 -0.0213 
 6              36 inanimate_pair_left…     0.207     0.451    17    18  0.0580 
 7              53 inanimate_single_le…     0.190     0.508    21     8  0.00504
 8              43 inanimate_pair_left…     0.186     0.556     6    23 -0.0163 
 9              58 inanimate_single_le…     0.165     0.448    16    16  0.0102 
10              27 animate_single_left…     0.161     0.378    11    14  0.0130 
# ℹ 54 more rows
# ℹ 7 more variables: ci_upper <dbl>, bg_image <chr>, dv_image <chr>,
#   violation_type <chr>, pair_of_pair.y <chr>, within_PoP_ID <dbl>,
#   change_direction <chr>

pose_change_plot_with_stimuli <- plot_diffmean_with_images(sumstats_ABBA_pose_alltrials, title_suffix = "")
ggsave("pose_change_diffmean_withimage.pdf", pose_change_plot_with_stimuli, width = 26, height = 7)  ## CAUTION: LONG RUNNING TIME

sumstats_ABBA_pose_1fam <- calculate_sumstats(ABBA_pose, total_trial_number = 2)

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
  always returns an ungrouped data frame and adjust accordingly.

plot1 <- plot_diffmean(sumstats = sumstats_ABBA_pose_1fam, title_suffix = "(1 fam)", print_plot = FALSE, ylim = 1.5)

sumstats_ABBA_pose_3fam <- calculate_sumstats(ABBA_pose, total_trial_number = 4)
plot2 <- plot_diffmean(sumstats = sumstats_ABBA_pose_3fam, title_suffix = "(3 fams)", print_plot = FALSE, ylim = 1.5)

sumstats_ABBA_pose_5fam <- calculate_sumstats(ABBA_pose, total_trial_number = 6)

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
  always returns an ungrouped data frame and adjust accordingly.

plot3 <- plot_diffmean(sumstats = sumstats_ABBA_pose_5fam, title_suffix = "(5 fams)", print_plot = FALSE, ylim = 1.5)

cowplot::plot_grid(plot1, plot2, plot3,ncol = 3)

3.3 models

3.3.1 effect from change direction

m_pose_0 <- lmerTest::lmer(dishab ~ 1 + (1|prolific_id), data = ABBA_pose)
m_pose_1 <- lmerTest::lmer(dishab ~ change_direction_dc + (1|prolific_id), data = ABBA_pose)

anova(m_pose_0, m_pose_1)

refitting model(s) with ML (instead of REML)

Data: ABBA_pose
Models:
m_pose_0: dishab ~ 1 + (1 | prolific_id)
m_pose_1: dishab ~ change_direction_dc + (1 | prolific_id)
         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
m_pose_0    3 2622.7 2639.3 -1308.3   2616.7                    
m_pose_1    4 2624.7 2646.8 -1308.3   2616.7 0.036  1     0.8496

3.3.2 model building

# adding total trial number and block number
m_pose_2 <- lmerTest::lmer(dishab ~ total_trial_number_c + block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_pose)
summary(m_pose_2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
dishab ~ total_trial_number_c + block_number_c + change_direction_dc +  
    (1 | prolific_id)
   Data: ABBA_pose

REML criterion at convergence: 2410.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.9078 -0.4514  0.0321  0.5045  6.1170 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.004404 0.06637 
 Residual                0.207004 0.45498 
Number of obs: 1856, groups:  prolific_id, 468

Fixed effects:
                       Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)          -7.577e-02  1.100e-02  4.583e+02  -6.889 1.86e-11 ***
total_trial_number_c  1.027e-01  6.470e-03  1.851e+03  15.865  < 2e-16 ***
block_number_c       -2.778e-03  1.524e-03  1.422e+03  -1.823   0.0685 .  
change_direction_dc1  1.285e-03  2.131e-02  1.829e+03   0.060   0.9519    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ttl___ blck__
ttl_trl_nm_  0.000              
blck_nmbr_c  0.000  0.000       
chng_drct_1 -0.011 -0.007  0.012

## adding interaction of total_trial_number, block_number, and change_direction
m_pose_3 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_pose)
m_pose_4 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c + change_direction_dc * total_trial_number_c + (1|prolific_id), data = ABBA_pose) # rank deficient

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

m_pose_5 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c + change_direction_dc * block_number_c + (1|prolific_id), data = ABBA_pose)
m_pose_6 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c * change_direction_dc + (1|prolific_id), data = ABBA_pose)

anova(m_pose_0, m_pose_1, m_pose_2, m_pose_3, m_pose_5, m_pose_6)

refitting model(s) with ML (instead of REML)

Data: ABBA_pose
Models:
m_pose_0: dishab ~ 1 + (1 | prolific_id)
m_pose_1: dishab ~ change_direction_dc + (1 | prolific_id)
m_pose_2: dishab ~ total_trial_number_c + block_number_c + change_direction_dc + (1 | prolific_id)
m_pose_3: dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1 | prolific_id)
m_pose_5: dishab ~ total_trial_number * block_number_c + change_direction_dc * block_number_c + (1 | prolific_id)
m_pose_6: dishab ~ total_trial_number * block_number_c * change_direction_dc + (1 | prolific_id)
         npar    AIC    BIC  logLik deviance    Chisq Df Pr(>Chisq)    
m_pose_0    3 2622.7 2639.3 -1308.3   2616.7                           
m_pose_1    4 2624.7 2646.8 -1308.3   2616.7   0.0360  1    0.84956    
m_pose_2    6 2389.8 2422.9 -1188.9   2377.8 238.8971  2    < 2e-16 ***
m_pose_3    7 2386.1 2424.8 -1186.0   2372.1   5.6867  1    0.01709 *  
m_pose_5    8 2387.4 2431.7 -1185.7   2371.4   0.6475  1    0.42100    
m_pose_6   10 2389.5 2444.8 -1184.8   2369.5   1.8966  2    0.38739    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m_pose_3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
dishab ~ total_trial_number_c * block_number_c + change_direction_dc +  
    (1 | prolific_id)
   Data: ABBA_pose

REML criterion at convergence: 2416.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.9616 -0.4359  0.0183  0.5047  6.0615 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.004402 0.06635 
 Residual                0.206475 0.45439 
Number of obs: 1856, groups:  prolific_id, 468

Fixed effects:
                                      Estimate Std. Error         df t value
(Intercept)                         -7.577e-02  1.099e-02  4.582e+02  -6.897
total_trial_number_c                 1.030e-01  6.463e-03  1.850e+03  15.929
block_number_c                      -2.780e-03  1.522e-03  1.421e+03  -1.827
change_direction_dc1                 1.310e-03  2.128e-02  1.828e+03   0.062
total_trial_number_c:block_number_c  2.219e-03  9.311e-04  1.850e+03   2.383
                                    Pr(>|t|)    
(Intercept)                         1.77e-11 ***
total_trial_number_c                 < 2e-16 ***
block_number_c                        0.0679 .  
change_direction_dc1                  0.9509    
total_trial_number_c:block_number_c   0.0173 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ttl___ blck__ chn__1
ttl_trl_nm_  0.000                     
blck_nmbr_c  0.000  0.000              
chng_drct_1 -0.011 -0.007  0.012       
ttl_tr__:__  0.000  0.019 -0.001  0.000

## add onefam
# m_pose_3 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_pose)
m_pose_7 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_pose)
m_pose_8 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc * change_direction_dc + (1|prolific_id), data = ABBA_pose)

anova(m_pose_3, m_pose_7, m_pose_8)

refitting model(s) with ML (instead of REML)

Data: ABBA_pose
Models:
m_pose_3: dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1 | prolific_id)
m_pose_7: dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1 | prolific_id)
m_pose_8: dishab ~ total_trial_number_c * block_number_c + onefam_dc * change_direction_dc + (1 | prolific_id)
         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
m_pose_3    7 2386.1 2424.8 -1186.0   2372.1                          
m_pose_7    8 2347.8 2392.1 -1165.9   2331.8 40.2489  1  2.236e-10 ***
m_pose_8    9 2348.7 2398.4 -1165.3   2330.7  1.1332  1     0.2871    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m_pose_7)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: dishab ~ total_trial_number_c * block_number_c + onefam_dc +  
    change_direction_dc + (1 | prolific_id)
   Data: ABBA_pose

REML criterion at convergence: 2380.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.9187 -0.4133 -0.0248  0.4865  6.2361 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.004549 0.06745 
 Residual                0.201929 0.44937 
Number of obs: 1856, groups:  prolific_id, 468

Fixed effects:
                                      Estimate Std. Error         df t value
(Intercept)                          2.415e-02  1.910e-02  1.436e+03   1.265
total_trial_number_c                 3.050e-02  1.305e-02  1.845e+03   2.338
block_number_c                      -3.066e-03  1.505e-03  1.421e+03  -2.036
onefam_dc1                          -2.879e-01  4.521e-02  1.847e+03  -6.369
change_direction_dc1                -9.115e-04  2.106e-02  1.827e+03  -0.043
total_trial_number_c:block_number_c  2.224e-03  9.213e-04  1.849e+03   2.414
                                    Pr(>|t|)    
(Intercept)                           0.2062    
total_trial_number_c                  0.0195 *  
block_number_c                        0.0419 *  
onefam_dc1                           2.4e-10 ***
change_direction_dc1                  0.9655    
total_trial_number_c:block_number_c   0.0159 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ttl___ blck__ onfm_1 chn__1
ttl_trl_nm_ -0.716                            
blck_nmbr_c -0.025  0.026                     
onefam_dc1  -0.822  0.872  0.030              
chng_drct_1 -0.020  0.011  0.012  0.017       
ttl_tr__:__  0.001  0.009 -0.001 -0.001  0.000

## add first_fam
# m_pose_7 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_pose)
m_pose_9 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1|prolific_id), data = ABBA_pose)
m_pose_10 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * total_trial_number_c + change_direction_dc + (1|prolific_id), data = ABBA_pose)
m_pose_11 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_pose)
m_pose_12 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam + onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_pose)
m_pose_13 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam * onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_pose) # rank deficient

fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients

m_pose_14 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam * onefam_dc * change_direction_dc + (1|prolific_id), data = ABBA_pose) # rank deficient

fixed-effect model matrix is rank deficient so dropping 8 columns / coefficients

anova(m_pose_7, m_pose_9, m_pose_10, m_pose_11, m_pose_12)

refitting model(s) with ML (instead of REML)

Data: ABBA_pose
Models:
m_pose_7: dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1 | prolific_id)
m_pose_9: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1 | prolific_id)
m_pose_10: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * total_trial_number_c + change_direction_dc + (1 | prolific_id)
m_pose_11: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * block_number_c + change_direction_dc + (1 | prolific_id)
m_pose_12: dishab ~ total_trial_number_c * block_number_c * first_fam + onefam_dc + change_direction_dc + (1 | prolific_id)
          npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
m_pose_7     8 2347.8 2392.1 -1165.9   2331.8                          
m_pose_9     9 2285.1 2334.8 -1133.5   2267.1  64.785  1  8.354e-16 ***
m_pose_10   10 2188.2 2243.5 -1084.1   2168.2  98.838  1  < 2.2e-16 ***
m_pose_11   10 2284.2 2339.5 -1132.1   2264.2   0.000  0               
m_pose_12   12 2173.7 2240.1 -1074.9   2149.7 114.473  2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# FIXME
# global_model_pose <- lmerTest::lmer(dishab ~ total_trial_number_log_c * block_number_log_c * first_fam * change_direction_dc + (1|prolific_id), data = ABBA_pose, REML = TRUE, na.action = na.fail)
# 
# best_model <- cAIC4::stepcAIC(global_model_pose, direction = "backward", trace = TRUE)
# summary(best_model)
# best_model

3.3.3 best fitting model

check_model(m_pose_12)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_pose_11)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_pose_10)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_pose_9)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_pose_7)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

best_model_pose <- m_pose_9
# dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1|prolific_id)

summary(best_model_pose)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: dishab ~ total_trial_number_c * block_number_c + onefam_dc +  
    first_fam + change_direction_dc + (1 | prolific_id)
   Data: ABBA_pose

REML criterion at convergence: 2322.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.5142 -0.4796 -0.0151  0.5043  6.3931 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.006699 0.08185 
 Residual                0.192985 0.43930 
Number of obs: 1856, groups:  prolific_id, 468

Fixed effects:
                                      Estimate Std. Error         df t value
(Intercept)                          9.597e-01  1.167e-01  9.793e+02   8.226
total_trial_number_c                 3.277e-02  1.282e-02  1.839e+03   2.556
block_number_c                      -5.885e-03  1.513e-03  1.500e+03  -3.890
onefam_dc1                          -2.774e-01  4.444e-02  1.842e+03  -6.242
first_fam                           -1.216e-01  1.497e-02  9.788e+02  -8.126
change_direction_dc1                -4.863e-04  2.067e-02  1.818e+03  -0.024
total_trial_number_c:block_number_c  2.019e-03  9.061e-04  1.849e+03   2.228
                                    Pr(>|t|)    
(Intercept)                         6.11e-16 ***
total_trial_number_c                0.010679 *  
block_number_c                      0.000104 ***
onefam_dc1                          5.35e-10 ***
first_fam                           1.33e-15 ***
change_direction_dc1                0.981237    
total_trial_number_c:block_number_c 0.026020 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ttl___ blck__ onfm_1 frst_f chn__1
ttl_trl_nm_ -0.091                                   
blck_nmbr_c -0.231  0.020                            
onefam_dc1  -0.103  0.872  0.022                     
first_fam   -0.987 -0.024  0.230 -0.030              
chng_drct_1 -0.002  0.012  0.012  0.018 -0.002       
ttl_tr__:__ -0.028  0.007  0.006 -0.003  0.029  0.000

model_performance(best_model_pose)

# Indices of model performance

AIC      |     AICc |      BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
--------------------------------------------------------------------------------
2340.833 | 2340.931 | 2390.569 |      0.200 |      0.172 | 0.034 | 0.432 | 0.439

plot(effects::allEffects(best_model_pose))

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

4 ABBA: NUMBER

ABBA_number <- exp2_core %>% filter(violation_type == "number") %>%
  left_join(exp2_pairs_number, by = c("violation_type", "bg_image", "dv_image"))

4.1 DA check

# data quality check
# DataExplorer::create_report(ABBA_number)
DataExplorer::plot_histogram(ABBA_number)

DataExplorer::plot_missing(ABBA_number)

DataExplorer::plot_qq(ABBA_number)

Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_qq()`).

Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_qq_line()`).

# clean up and prep dataset
ABBA_number <- ABBA_number %>% filter(!is.na(first_fam) & !is.na(dishab))

4.1.1 summary of mean dishab

ABBA_number %>% group_by(change_direction) %>% rstatix::get_summary_stats(dishab, type = "common")

# A tibble: 2 × 11
  change_direction variable     n   min   max median   iqr   mean    sd    se
  <chr>            <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl> <dbl>
1 pair to single   dishab     943 -3.20  1.28 -0.048 0.416 -0.117 0.453 0.015
2 single to pair   dishab     909 -2.46  2.39  0.003 0.404 -0.035 0.464 0.015
# ℹ 1 more variable: ci <dbl>

ABBA_number %>% group_by(change_direction, total_trial_number) %>% rstatix::get_summary_stats(dishab, type = "common")

# A tibble: 6 × 12
  total_trial_number change_direction variable     n    min   max median   iqr
               <dbl> <chr>            <fct>    <dbl>  <dbl> <dbl>  <dbl> <dbl>
1                  2 pair to single   dishab     311 -3.20   0.94 -0.306 0.622
2                  4 pair to single   dishab     312 -1.10   1.28 -0.016 0.289
3                  6 pair to single   dishab     320 -1.89   1.28  0.026 0.314
4                  2 single to pair   dishab     314 -2.46   1.15 -0.223 0.588
5                  4 single to pair   dishab     286 -1.95   2.39  0.032 0.298
6                  6 single to pair   dishab     309 -0.927  1.99  0.074 0.327
# ℹ 4 more variables: mean <dbl>, sd <dbl>, se <dbl>, ci <dbl>

ABBA_number %>%
    ggplot(aes(x = as.factor(change_direction), y = dishab)) +
    geom_boxplot(outlier.color = "red", outlier.shape = 1, outlier.alpha = 0.2) +
    geom_point(alpha = 0.1) +
  stat_summary(fun = "mean", geom = "point", position = position_dodge(.9), color = "red") +
    stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.2, position = position_dodge(.9), color = "red") +
    facet_grid(~total_trial_number) +
    labs(x = "Direction of Change",
         y = "Dishabituation (log(test) - log(last fam))",
         title = "Dishabituation by different pairs of stimuli (number change)") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme_bw()

ABBA_number$change_direction_dc <- as.factor(ABBA_number$change_direction)
contrasts(ABBA_number$change_direction_dc) <- contr.sum(2)/2

ABBA_number$onefam_dc <- as.factor(ABBA_number$onefam)

ABBA_number$total_trial_number_c <- scale(ABBA_number$total_trial_number, center = TRUE, scale = FALSE)
ABBA_number$block_number_c <- scale(ABBA_number$block_number, center = TRUE, scale = FALSE)

ABBA_number$total_trial_number_log <- log(ABBA_number$total_trial_number)
ABBA_number$block_number_log <- log(ABBA_number$block_number)

ABBA_number$total_trial_number_log_c <- scale(ABBA_number$total_trial_number_log, center = TRUE, scale = FALSE)
ABBA_number$block_number_log_c <- scale(ABBA_number$block_number_log, center = TRUE, scale = FALSE)

4.2 plots

4.2.1 raw dishab plots

plot_raw_dishab(dataset = ABBA_number)

plot_raw_dishab(dataset = ABBA_number, by_trialnumber = TRUE)

4.2.2 diffemean plot

sumstats_ABBA_number_alltrials <- calculate_sumstats(ABBA_number)
plot_diffmean(sumstats = sumstats_ABBA_number_alltrials, title_suffix = "(all trials)")

get_stimuli_for_diffmean(sumstats_ABBA_number_alltrials, exp2_pairs_number, "number")

# A tibble: 64 × 14
   pair_of_pair_ID pair_of_pair.x       mean_diff pooled_sd    n1    n2 ci_lower
             <int> <chr>                    <dbl>     <dbl> <int> <int>    <dbl>
 1              11 animate_pair_left_0…     0.427     0.556    16    19   0.243 
 2              41 inanimate_pair_left…     0.425     0.417    19    13   0.280 
 3               6 animate_pair_left_0…     0.418     0.582    12    11   0.181 
 4              32 animate_pair_right_…     0.395     0.372    19    14   0.268 
 5              15 animate_pair_left_0…     0.318     0.881    11    15  -0.0203
 6              10 animate_pair_left_0…     0.298     0.457    21    10   0.137 
 7              31 animate_pair_right_…     0.288     0.658    21    15   0.0733
 8              46 inanimate_pair_left…     0.278     0.390    11    11   0.115 
 9              55 inanimate_pair_righ…     0.272     0.473    17    12   0.0994
10              50 inanimate_pair_righ…     0.256     0.405    17    14   0.113 
# ℹ 54 more rows
# ℹ 7 more variables: ci_upper <dbl>, bg_image <chr>, dv_image <chr>,
#   violation_type <chr>, pair_of_pair.y <chr>, within_PoP_ID <dbl>,
#   change_direction <chr>

number_change_plot_with_stimuli <- plot_diffmean_with_images(sumstats_ABBA_number_alltrials, title_suffix = "")
ggsave("number_change_diffmean_withimage.pdf", number_change_plot_with_stimuli, width = 26, height = 7)  ## CAUTION: LONG RUNNING TIME

sumstats_ABBA_number_1fam <- calculate_sumstats(ABBA_number, total_trial_number = 2)

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
  always returns an ungrouped data frame and adjust accordingly.

plot1 <- plot_diffmean(sumstats = sumstats_ABBA_number_1fam, title_suffix = "(1 fam)", print_plot = FALSE, ylim = 1.5)

sumstats_ABBA_number_3fam <- calculate_sumstats(ABBA_number, total_trial_number = 4)

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
  always returns an ungrouped data frame and adjust accordingly.

plot2 <- plot_diffmean(sumstats = sumstats_ABBA_number_3fam, title_suffix = "(3 fams)", print_plot = FALSE, ylim = 1.5)

sumstats_ABBA_number_5fam <- calculate_sumstats(ABBA_number, total_trial_number = 6)
plot3 <- plot_diffmean(sumstats = sumstats_ABBA_number_5fam, title_suffix = "(5 fams)", print_plot = FALSE, ylim = 1.5)

cowplot::plot_grid(plot1, plot2, plot3,ncol = 3)

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).

4.3 models

4.3.1 effect from change direction

m_number_0 <- lmerTest::lmer(dishab ~ 1 + (1|prolific_id), data = ABBA_number)
m_number_1 <- lmerTest::lmer(dishab ~ change_direction_dc + (1|prolific_id), data = ABBA_number)

anova(m_number_0, m_number_1)

refitting model(s) with ML (instead of REML)

Data: ABBA_number
Models:
m_number_0: dishab ~ 1 + (1 | prolific_id)
m_number_1: dishab ~ change_direction_dc + (1 | prolific_id)
           npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
m_number_0    3 2385.6 2402.2 -1189.8   2379.6                        
m_number_1    4 2372.9 2395.0 -1182.5   2364.9 14.63  1  0.0001308 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.3.2 model building

# adding total trial number and block number
m_number_2 <- lmerTest::lmer(dishab ~ total_trial_number_c + block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_number)
summary(m_number_2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
dishab ~ total_trial_number_c + block_number_c + change_direction_dc +  
    (1 | prolific_id)
   Data: ABBA_number

REML criterion at convergence: 2140.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.6706 -0.4531  0.0213  0.5112  5.7229 

Random effects:
 Groups      Name        Variance  Std.Dev.
 prolific_id (Intercept) 0.0001384 0.01176 
 Residual                0.1829286 0.42770 
Number of obs: 1852, groups:  prolific_id, 468

Fixed effects:
                       Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)          -7.617e-02  9.955e-03  4.501e+02  -7.651 1.22e-13 ***
total_trial_number_c  9.964e-02  6.041e-03  1.844e+03  16.493  < 2e-16 ***
block_number_c       -2.718e-03  1.432e-03  1.416e+03  -1.898   0.0579 .  
change_direction_dc1 -8.503e-02  1.989e-02  1.834e+03  -4.275 2.01e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ttl___ blck__
ttl_trl_nm_  0.000              
blck_nmbr_c  0.000  0.005       
chng_drct_1 -0.018 -0.009  0.011

## adding interaction of total_trial_number, block_number, and change_direction
m_number_3 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_number)
m_number_4 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c + change_direction_dc * total_trial_number_c + (1|prolific_id), data = ABBA_number) # rank deficient

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

m_number_5 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c + change_direction_dc * block_number_c + (1|prolific_id), data = ABBA_number)
m_number_6 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c * change_direction_dc + (1|prolific_id), data = ABBA_number)

anova(m_number_0, m_number_1, m_number_2, m_number_3, m_number_5, m_number_6)

refitting model(s) with ML (instead of REML)

Data: ABBA_number
Models:
m_number_0: dishab ~ 1 + (1 | prolific_id)
m_number_1: dishab ~ change_direction_dc + (1 | prolific_id)
m_number_2: dishab ~ total_trial_number_c + block_number_c + change_direction_dc + (1 | prolific_id)
m_number_3: dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1 | prolific_id)
m_number_5: dishab ~ total_trial_number * block_number_c + change_direction_dc * block_number_c + (1 | prolific_id)
m_number_6: dishab ~ total_trial_number * block_number_c * change_direction_dc + (1 | prolific_id)
           npar    AIC    BIC  logLik deviance    Chisq Df Pr(>Chisq)    
m_number_0    3 2385.6 2402.2 -1189.8   2379.6                           
m_number_1    4 2372.9 2395.0 -1182.5   2364.9  14.6300  1  0.0001308 ***
m_number_2    6 2119.2 2152.4 -1053.6   2107.2 257.7189  2  < 2.2e-16 ***
m_number_3    7 2114.4 2153.1 -1050.2   2100.4   6.8337  1  0.0089454 ** 
m_number_5    8 2114.3 2158.5 -1049.2   2098.3   2.0444  1  0.1527695    
m_number_6   10 2116.3 2171.6 -1048.2   2096.3   1.9989  2  0.3680801    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m_number_3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
dishab ~ total_trial_number_c * block_number_c + change_direction_dc +  
    (1 | prolific_id)
   Data: ABBA_number

REML criterion at convergence: 2145.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.5769 -0.4587  0.0210  0.5122  5.7215 

Random effects:
 Groups      Name        Variance  Std.Dev.
 prolific_id (Intercept) 0.0004752 0.0218  
 Residual                0.1820179 0.4266  
Number of obs: 1852, groups:  prolific_id, 468

Fixed effects:
                                      Estimate Std. Error         df t value
(Intercept)                         -7.605e-02  9.967e-03  4.502e+02  -7.631
total_trial_number_c                 9.952e-02  6.032e-03  1.843e+03  16.500
block_number_c                      -2.632e-03  1.429e-03  1.415e+03  -1.843
change_direction_dc1                -8.480e-02  1.985e-02  1.832e+03  -4.271
total_trial_number_c:block_number_c  2.268e-03  8.672e-04  1.846e+03   2.615
                                    Pr(>|t|)    
(Intercept)                         1.41e-13 ***
total_trial_number_c                 < 2e-16 ***
block_number_c                       0.06558 .  
change_direction_dc1                2.04e-05 ***
total_trial_number_c:block_number_c  0.00899 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ttl___ blck__ chn__1
ttl_trl_nm_  0.000                     
blck_nmbr_c  0.000  0.005              
chng_drct_1 -0.018 -0.009  0.011       
ttl_tr__:__  0.005 -0.007  0.022  0.004

## add onefam
# m_number_3 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_number)
m_number_7 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_number) # singular

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

m_number_8 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc * change_direction_dc + (1|prolific_id), data = ABBA_number) # singular

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

anova(m_number_3, m_number_7, m_number_8)

refitting model(s) with ML (instead of REML)

Data: ABBA_number
Models:
m_number_3: dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1 | prolific_id)
m_number_7: dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1 | prolific_id)
m_number_8: dishab ~ total_trial_number_c * block_number_c + onefam_dc * change_direction_dc + (1 | prolific_id)
           npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
m_number_3    7 2114.4 2153.1 -1050.2   2100.4                          
m_number_7    8 2071.5 2115.7 -1027.7   2055.5 44.9047  1  2.069e-11 ***
m_number_8    9 2073.4 2123.1 -1027.7   2055.4  0.0778  1     0.7803    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## add first fam
m_number_9 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1|prolific_id), data = ABBA_number)
m_number_10 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * total_trial_number_c + change_direction_dc + (1|prolific_id), data = ABBA_number)
m_number_11 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_number)
m_number_12 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam + onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_number)
m_number_13 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam * onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_number) # rank deficient

fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients

m_number_14 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam * onefam_dc * change_direction_dc + (1|prolific_id), data = ABBA_number) # rank deficient

fixed-effect model matrix is rank deficient so dropping 8 columns / coefficients

anova(m_number_3, m_number_9, m_number_10, m_number_11, m_number_12)

refitting model(s) with ML (instead of REML)

Data: ABBA_number
Models:
m_number_3: dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1 | prolific_id)
m_number_9: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1 | prolific_id)
m_number_10: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * total_trial_number_c + change_direction_dc + (1 | prolific_id)
m_number_11: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * block_number_c + change_direction_dc + (1 | prolific_id)
m_number_12: dishab ~ total_trial_number_c * block_number_c * first_fam + onefam_dc + change_direction_dc + (1 | prolific_id)
            npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
m_number_3     7 2114.4 2153.1 -1050.20   2100.4                         
m_number_9     9 2018.5 2068.2 -1000.26   2000.5 99.865  2  < 2.2e-16 ***
m_number_10   10 1956.0 2011.3  -968.02   1936.0 64.485  1  9.726e-16 ***
m_number_11   10 2001.9 2057.2  -990.95   1981.9  0.000  0               
m_number_12   12 1935.6 2001.9  -955.81   1911.6 70.294  2  5.444e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.3.3 best fitting model

check_model(m_number_12)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_number_11)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_number_10)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_number_9)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_number_3)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

best_model_number <- m_number_9
# dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1|prolific_id)

summary(best_model_number)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: dishab ~ total_trial_number_c * block_number_c + onefam_dc +  
    first_fam + change_direction_dc + (1 | prolific_id)
   Data: ABBA_number

REML criterion at convergence: 2057.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.9858 -0.4775 -0.0190  0.5051  5.4463 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.004104 0.06406 
 Residual                0.169135 0.41126 
Number of obs: 1852, groups:  prolific_id, 468

Fixed effects:
                                      Estimate Std. Error         df t value
(Intercept)                          8.238e-01  1.079e-01  1.001e+03   7.633
total_trial_number_c                 2.960e-02  1.188e-02  1.841e+03   2.493
block_number_c                      -5.166e-03  1.423e-03  1.507e+03  -3.632
onefam_dc1                          -2.779e-01  4.133e-02  1.838e+03  -6.723
first_fam                           -1.044e-01  1.382e-02  9.923e+02  -7.556
change_direction_dc1                -7.907e-02  1.933e-02  1.817e+03  -4.090
total_trial_number_c:block_number_c  2.104e-03  8.442e-04  1.838e+03   2.492
                                    Pr(>|t|)    
(Intercept)                         5.33e-14 ***
total_trial_number_c                0.012770 *  
block_number_c                      0.000291 ***
onefam_dc1                          2.37e-11 ***
first_fam                           9.45e-14 ***
change_direction_dc1                4.50e-05 ***
total_trial_number_c:block_number_c 0.012789 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ttl___ blck__ onfm_1 frst_f chn__1
ttl_trl_nm_ -0.107                                   
blck_nmbr_c -0.244 -0.007                            
onefam_dc1  -0.118  0.869 -0.012                     
first_fam   -0.987 -0.006  0.248 -0.011              
chng_drct_1  0.049  0.012 -0.002  0.019 -0.054       
ttl_tr__:__ -0.027  0.000  0.028  0.004  0.027  0.003

model_performance(best_model_number)

# Indices of model performance

AIC      |     AICc |      BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
--------------------------------------------------------------------------------
2075.323 | 2075.421 | 2125.039 |      0.204 |      0.185 | 0.024 | 0.406 | 0.411

plot(effects::allEffects(best_model_number))

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

5 ABBA: ANIMACY

# clean up and prep dataset
ABBA_animacy_noidentity <- exp2_core %>% filter(violation_type == "animacy") %>%
  mutate(
    bg_image = str_replace_all(bg_image, "[0-9]+", ""),
    dv_image = str_replace_all(dv_image, "[0-9]+", "")
  ) %>%
  left_join(exp2_pairs_animacy_noidentity, by = c("violation_type", "bg_image", "dv_image"))

5.1 DA check

# data quality check
# DataExplorer::create_report(ABBA_animacy)
DataExplorer::plot_histogram(ABBA_animacy_noidentity)

DataExplorer::plot_missing(ABBA_animacy_noidentity)

DataExplorer::plot_qq(ABBA_animacy_noidentity)

Warning: Removed 28 rows containing non-finite outside the scale range
(`stat_qq()`).

Warning: Removed 28 rows containing non-finite outside the scale range
(`stat_qq_line()`).

# clean up and prep dataset
ABBA_animacy_noidentity <- ABBA_animacy_noidentity %>%
  filter(!is.na(first_fam) & !is.na(dishab))

5.1.1 summary of mean dishab

ABBA_animacy_noidentity %>% group_by(change_direction) %>% rstatix::get_summary_stats(dishab, type = "common")

# A tibble: 2 × 11
  change_direction     variable     n   min   max median   iqr  mean    sd    se
  <chr>                <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
1 animate to inanimate dishab     931 -2.74  2.53  0.103 0.474 0.124 0.544 0.018
2 inanimate to animate dishab     925 -2.10  2.37  0.165 0.521 0.2   0.508 0.017
# ℹ 1 more variable: ci <dbl>

ABBA_animacy_noidentity %>% group_by(change_direction, total_trial_number) %>% rstatix::get_summary_stats(dishab, type = "common")

# A tibble: 6 × 12
  total_trial_number change_direction    variable     n   min   max median   iqr
               <dbl> <chr>               <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl>
1                  2 animate to inanima… dishab     299 -2.74  1.3  -0.125 0.469
2                  4 animate to inanima… dishab     318 -1.79  2.53  0.178 0.431
3                  6 animate to inanima… dishab     314 -1.64  2.48  0.248 0.508
4                  2 inanimate to anima… dishab     279 -2.10  1.91 -0.025 0.469
5                  4 inanimate to anima… dishab     312 -1.77  2.20  0.262 0.536
6                  6 inanimate to anima… dishab     334 -1.09  2.37  0.235 0.467
# ℹ 4 more variables: mean <dbl>, sd <dbl>, se <dbl>, ci <dbl>

ABBA_animacy_noidentity %>%
    ggplot(aes(x = as.factor(change_direction), y = dishab)) +
    geom_boxplot(outlier.color = "red", outlier.shape = 1, outlier.alpha = 0.2) +
    geom_point(alpha = 0.1) +
  stat_summary(fun = "mean", geom = "point", position = position_dodge(.9), color = "red") +
    stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.2, position = position_dodge(.9), color = "red") +
    facet_grid(~total_trial_number) +
    labs(x = "Direction of Change",
         y = "Dishabituation (log(test) - log(last fam))",
         title = "Dishabituation by different pairs of stimuli (animacy change)") +
  theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

ABBA_animacy_noidentity$change_direction_dc <- as.factor(ABBA_animacy_noidentity$change_direction)
contrasts(ABBA_animacy_noidentity$change_direction_dc) <- contr.sum(2)/2

ABBA_animacy_noidentity$onefam_dc <- as.factor(ABBA_animacy_noidentity$onefam)

ABBA_animacy_noidentity$total_trial_number_c <- scale(ABBA_animacy_noidentity$total_trial_number, center = TRUE, scale = FALSE)
ABBA_animacy_noidentity$block_number_c <- scale(ABBA_animacy_noidentity$block_number, center = TRUE, scale = FALSE)

ABBA_animacy_noidentity$total_trial_number_log <- log(ABBA_animacy_noidentity$total_trial_number)
ABBA_animacy_noidentity$block_number_log <- log(ABBA_animacy_noidentity$block_number)

ABBA_animacy_noidentity$total_trial_number_log_c <- scale(ABBA_animacy_noidentity$total_trial_number_log, center = TRUE, scale = FALSE)
ABBA_animacy_noidentity$block_number_log_c <- scale(ABBA_animacy_noidentity$block_number_log, center = TRUE, scale = FALSE)

5.2 plots

5.2.1 raw dishab plots

plot_raw_dishab(dataset = ABBA_animacy_noidentity)

plot_raw_dishab(dataset = ABBA_animacy_noidentity, by_trialnumber = TRUE)

5.2.2 diffemean plot

sumstats_ABBA_animacy_alltrials <- calculate_sumstats(ABBA_animacy_noidentity)
 plot_diffmean_with_images(sumstats = sumstats_ABBA_animacy_alltrials, title_suffix = "(all trials)")

get_stimuli_for_diffmean(sumstats_ABBA_animacy_alltrials, exp2_pairs_animacy_noidentity, "animacy")

# A tibble: 4 × 14
  pair_of_pair_ID pair_of_pair.x        mean_diff pooled_sd    n1    n2 ci_lower
            <int> <chr>                     <dbl>     <dbl> <int> <int>    <dbl>
1               1 animate_pair_left_<-…    0.100      0.570   235   234  0.0485 
2               2 animate_pair_right_<…    0.0754     0.506   225   231  0.0290 
3               3 animate_single_left_…    0.0732     0.522   236   240  0.0264 
4               4 animate_single_right…    0.0527     0.505   235   220  0.00622
# ℹ 7 more variables: ci_upper <dbl>, bg_image <chr>, dv_image <chr>,
#   violation_type <chr>, pair_of_pair.y <chr>, within_PoP_ID <dbl>,
#   change_direction <chr>

animacy_change_plot_with_stimuli <- plot_diffmean_with_images(sumstats_ABBA_animacy_alltrials, title_suffix = "")
ggsave("animacy_change_diffmean_withimage.pdf", animacy_change_plot_with_stimuli, width = 5, height = 4)

sumstats_ABBA_animacy_1fam <- calculate_sumstats(ABBA_animacy_noidentity, total_trial_number = 2)
plot1 <- plot_diffmean_with_images(sumstats = sumstats_ABBA_animacy_1fam, title_suffix = "(1 fam)", print_plot = FALSE, ylim = 0.5)

sumstats_ABBA_animacy_3fam <- calculate_sumstats(ABBA_animacy_noidentity, total_trial_number = 4)
plot2 <- plot_diffmean_with_images(sumstats = sumstats_ABBA_animacy_3fam, title_suffix = "(3 fams)", print_plot = FALSE, ylim = 0.5)

sumstats_ABBA_animacy_5fam <- calculate_sumstats(ABBA_animacy_noidentity, total_trial_number = 6)
plot3 <- plot_diffmean_with_images(sumstats = sumstats_ABBA_animacy_5fam, title_suffix = "(5 fams)", print_plot = FALSE, ylim = 0.5)

cowplot::plot_grid(plot1, plot2, plot3,ncol = 3)

5.3 models

5.3.1 effect from change direction

m_animacy_0 <- lmerTest::lmer(dishab ~ 1 + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_1 <- lmerTest::lmer(dishab ~ change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)

anova(m_animacy_0, m_animacy_1)

refitting model(s) with ML (instead of REML)

Data: ABBA_animacy_noidentity
Models:
m_animacy_0: dishab ~ 1 + (1 | prolific_id)
m_animacy_1: dishab ~ change_direction_dc + (1 | prolific_id)
            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
m_animacy_0    3 2888.8 2905.3 -1441.4   2882.8                        
m_animacy_1    4 2880.4 2902.5 -1436.2   2872.4 10.384  1   0.001271 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.3.2 model building

# adding total trial number and block number
m_animacy_2 <- lmerTest::lmer(dishab ~ total_trial_number_c + block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)

## adding interaction of total_trial_number, block_number, and change_direction
m_animacy_3 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_4 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c + change_direction_dc * total_trial_number_c + (1|prolific_id), data = ABBA_animacy_noidentity) # rank deficient

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

m_animacy_5 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c + change_direction_dc * block_number_c + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_6 <- lmerTest::lmer(dishab ~ total_trial_number * block_number_c * change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)

anova(m_animacy_0, m_animacy_1, m_animacy_2, m_animacy_3, m_animacy_5, m_animacy_6)

refitting model(s) with ML (instead of REML)

Data: ABBA_animacy_noidentity
Models:
m_animacy_0: dishab ~ 1 + (1 | prolific_id)
m_animacy_1: dishab ~ change_direction_dc + (1 | prolific_id)
m_animacy_2: dishab ~ total_trial_number_c + block_number_c + change_direction_dc + (1 | prolific_id)
m_animacy_3: dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1 | prolific_id)
m_animacy_5: dishab ~ total_trial_number * block_number_c + change_direction_dc * block_number_c + (1 | prolific_id)
m_animacy_6: dishab ~ total_trial_number * block_number_c * change_direction_dc + (1 | prolific_id)
            npar    AIC    BIC  logLik deviance    Chisq Df Pr(>Chisq)    
m_animacy_0    3 2888.8 2905.3 -1441.4   2882.8                           
m_animacy_1    4 2880.4 2902.5 -1436.2   2872.4  10.3843  1   0.001271 ** 
m_animacy_2    6 2685.2 2718.4 -1336.6   2673.2 199.1765  2  < 2.2e-16 ***
m_animacy_3    7 2671.2 2709.9 -1328.6   2657.2  16.0308  1  6.232e-05 ***
m_animacy_5    8 2672.3 2716.6 -1328.2   2656.3   0.8452  1   0.357921    
m_animacy_6   10 2673.1 2728.4 -1326.6   2653.1   3.2092  2   0.200966    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## add onefam
# m_animacy_3 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_7 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_8 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc * change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)

anova(m_animacy_3, m_animacy_7, m_animacy_8)

refitting model(s) with ML (instead of REML)

Data: ABBA_animacy_noidentity
Models:
m_animacy_3: dishab ~ total_trial_number_c * block_number_c + change_direction_dc + (1 | prolific_id)
m_animacy_7: dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1 | prolific_id)
m_animacy_8: dishab ~ total_trial_number_c * block_number_c + onefam_dc * change_direction_dc + (1 | prolific_id)
            npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
m_animacy_3    7 2671.2 2709.9 -1328.6   2657.2                          
m_animacy_7    8 2615.7 2659.9 -1299.8   2599.7 57.5239  1  3.339e-14 ***
m_animacy_8    9 2615.7 2665.4 -1298.8   2597.7  1.9574  1     0.1618    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## add first fam
# m_animacy_7 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_9 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_10 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * total_trial_number_c + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_11 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * block_number_c + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_12 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam + onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)
m_animacy_13 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam * onefam_dc + change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity) # rank deficient

fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients

m_animacy_14 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c * first_fam * onefam_dc * change_direction_dc + (1|prolific_id), data = ABBA_animacy_noidentity)  # rank deficient

fixed-effect model matrix is rank deficient so dropping 8 columns / coefficients

anova(m_animacy_7, m_animacy_9, m_animacy_10, m_animacy_11, m_animacy_12)

refitting model(s) with ML (instead of REML)

Data: ABBA_animacy_noidentity
Models:
m_animacy_7: dishab ~ total_trial_number_c * block_number_c + onefam_dc + change_direction_dc + (1 | prolific_id)
m_animacy_9: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1 | prolific_id)
m_animacy_10: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * total_trial_number_c + change_direction_dc + (1 | prolific_id)
m_animacy_11: dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam * block_number_c + change_direction_dc + (1 | prolific_id)
m_animacy_12: dishab ~ total_trial_number_c * block_number_c * first_fam + onefam_dc + change_direction_dc + (1 | prolific_id)
             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
m_animacy_7     8 2615.7 2659.9 -1299.8   2599.7                         
m_animacy_9     9 2591.5 2641.3 -1286.8   2573.5 26.119  1  3.211e-07 ***
m_animacy_10   10 2530.9 2586.1 -1255.4   2510.9 62.661  1  2.455e-15 ***
m_animacy_11   10 2591.0 2646.3 -1285.5   2571.0  0.000  0               
m_animacy_12   12 2521.0 2587.3 -1248.5   2497.0 74.025  2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.3.3 best fitting model

check_model(m_number_12)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_number_11)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_number_10)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_number_9)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

check_model(m_number_7)

Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.
Some of the variables were in matrix-format - probably you used
  `scale()` on your data?
  If so, and you get an error, please try `datawizard::standardize()` to
  standardize your data.

best_model_animacy <- m_animacy_9
# dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1|prolific_id)

summary(best_model_animacy)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: dishab ~ total_trial_number_c * block_number_c + onefam_dc +  
    first_fam + change_direction_dc + (1 | prolific_id)
   Data: ABBA_animacy_noidentity

REML criterion at convergence: 2626.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.6306 -0.5079 -0.0620  0.4633  4.6888 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.02183  0.1477  
 Residual                0.21601  0.4648  
Number of obs: 1856, groups:  prolific_id, 468

Fixed effects:
                                      Estimate Std. Error         df t value
(Intercept)                          9.687e-01  1.333e-01  1.164e+03   7.269
total_trial_number_c                 1.221e-02  1.355e-02  1.826e+03   0.901
block_number_c                       1.753e-03  3.240e-03  1.547e+03   0.541
onefam_dc1                          -3.669e-01  4.755e-02  1.834e+03  -7.716
first_fam                           -8.943e-02  1.706e-02  1.161e+03  -5.242
change_direction_dc1                -1.164e-01  4.347e-02  1.455e+03  -2.679
total_trial_number_c:block_number_c  4.019e-03  1.020e-03  1.812e+03   3.939
                                    Pr(>|t|)    
(Intercept)                         6.63e-13 ***
total_trial_number_c                 0.36778    
block_number_c                       0.58849    
onefam_dc1                          1.96e-14 ***
first_fam                           1.89e-07 ***
change_direction_dc1                 0.00748 ** 
total_trial_number_c:block_number_c 8.48e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ttl___ blck__ onfm_1 frst_f chn__1
ttl_trl_nm_ -0.099                                   
blck_nmbr_c -0.219  0.045                            
onefam_dc1  -0.116  0.861  0.034                     
first_fam   -0.989  0.004  0.218  0.005              
chng_drct_1  0.083 -0.031 -0.861 -0.028 -0.081       
ttl_tr__:__  0.000  0.028 -0.046  0.027  0.000  0.032

model_performance(best_model_animacy)

# Indices of model performance

AIC      |     AICc |      BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
--------------------------------------------------------------------------------
2644.563 | 2644.660 | 2694.298 |      0.230 |      0.152 | 0.092 | 0.447 | 0.465

plot(effects::allEffects(best_model_animacy))

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors total_trial_number_c, block_number_c are one-column matrices that
were converted to vectors

6 ABBA: across violation types

# is it valid to separate by different violation type
all_vio <- bind_rows(ABBA_pose, ABBA_number, ABBA_animacy_noidentity)
m_all_0 <- lmerTest::lmer(dishab ~ 1 + (1|prolific_id), data = all_vio)
m_all_1 <- lmerTest::lmer(dishab ~ violation_type + (1|prolific_id), data = all_vio)
anova(m_all_0, m_all_1)

refitting model(s) with ML (instead of REML)

Data: all_vio
Models:
m_all_0: dishab ~ 1 + (1 | prolific_id)
m_all_1: dishab ~ violation_type + (1 | prolific_id)
        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
m_all_0    3 8207.2 8227.0 -4100.6   8201.2                         
m_all_1    5 7925.6 7958.8 -3957.8   7915.6 285.52  2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# double check with best fitting model
m_all_2 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + (1|prolific_id), data = all_vio)
m_all_3 <- lmerTest::lmer(dishab ~ total_trial_number_c * block_number_c + onefam_dc + first_fam + change_direction_dc + violation_type + (1|prolific_id), data = all_vio) # rank deficient, became the same model as m_all_2

fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients

7 check baseline looking

7.1 pose

ABBA_pose <- ABBA_pose %>%
  mutate(bg_image_type = ifelse(str_detect(bg_image, "left") == TRUE, "left", "right"))

ABBA_pose$bg_image_type <- as.factor(ABBA_pose$bg_image_type)
ABBA_pose %>% group_by(bg_image_type) %>% rstatix::get_summary_stats(first_fam, type = "common")

# A tibble: 2 × 11
  bg_image_type variable      n   min   max median   iqr  mean    sd    se    ci
  <fct>         <fct>     <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 left          first_fam   939  6.25  9.88   7.69  1.12  7.73 0.735 0.024 0.047
2 right         first_fam   917  6.24  9.89   7.67  1.02  7.72 0.726 0.024 0.047

contrasts(ABBA_pose$bg_image_type) <- contr.sum(2)/2
m1 <- lmerTest::lmer(first_fam ~ block_number_log_c + bg_image_type + (1|prolific_id), data = ABBA_pose)
m2 <- lmerTest::lmer(first_fam ~ block_number_log_c * total_trial_number_log_c + bg_image_type + (1|prolific_id), data = ABBA_pose)
m3 <- lmerTest::lmer(first_fam ~ block_number_log_c * total_trial_number_log_c * bg_image_type + (1|prolific_id), data = ABBA_pose)
anova(m1, m2, m3)

refitting model(s) with ML (instead of REML)

Data: ABBA_pose
Models:
m1: first_fam ~ block_number_log_c + bg_image_type + (1 | prolific_id)
m2: first_fam ~ block_number_log_c * total_trial_number_log_c + bg_image_type + (1 | prolific_id)
m3: first_fam ~ block_number_log_c * total_trial_number_log_c * bg_image_type + (1 | prolific_id)
   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
m1    5 3189.4 3217.1 -1589.7   3179.4                     
m2    7 3193.0 3231.6 -1589.5   3179.0 0.4713  2     0.7901
m3   10 3196.4 3251.6 -1588.2   3176.4 2.5886  3     0.4595

summary(m1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: first_fam ~ block_number_log_c + bg_image_type + (1 | prolific_id)
   Data: ABBA_pose

REML criterion at convergence: 3197.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2739 -0.5554 -0.0920  0.4780  4.1045 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.2908   0.5392  
 Residual                0.2011   0.4485  
Number of obs: 1856, groups:  prolific_id, 468

Fixed effects:
                     Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)         7.726e+00  2.702e-02  4.666e+02 285.978   <2e-16 ***
block_number_log_c -2.449e-01  1.300e-02  1.404e+03 -18.843   <2e-16 ***
bg_image_type1     -8.928e-03  2.306e-02  1.495e+03  -0.387    0.699    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) blc___
blck_nmbr__  0.000       
bg_img_typ1 -0.005  0.015

plot(effects::allEffects(m1))

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictor block_number_log_c is a one-column matrix that was converted to a
vector

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictor block_number_log_c is a one-column matrix that was converted to a
vector

7.2 number

ABBA_number <- ABBA_number %>%
  mutate(bg_image_type = ifelse(str_detect(bg_image, "pair") == TRUE, "pair", "single"))

ABBA_number$bg_image_type <- as.factor(ABBA_number$bg_image_type)
ABBA_number %>% group_by(bg_image_type) %>% rstatix::get_summary_stats(first_fam, type = "common")

# A tibble: 2 × 11
  bg_image_type variable      n   min   max median   iqr  mean    sd    se    ci
  <fct>         <fct>     <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 pair          first_fam   943  6.26  9.90   7.71  1.11  7.76 0.751 0.024 0.048
2 single        first_fam   909  6.27  9.84   7.65  1.02  7.68 0.717 0.024 0.047

contrasts(ABBA_number$bg_image_type) <- contr.sum(2)/2
m1 <- lmerTest::lmer(first_fam ~ block_number_log_c + bg_image_type + (1|prolific_id), data = ABBA_number)
m2 <- lmerTest::lmer(first_fam ~ block_number_log_c * total_trial_number_log_c + bg_image_type + (1|prolific_id), data = ABBA_number)
m3 <- lmerTest::lmer(first_fam ~ block_number_log_c * total_trial_number_log_c * bg_image_type + (1|prolific_id), data = ABBA_number)
anova(m1, m2, m3)

refitting model(s) with ML (instead of REML)

Data: ABBA_number
Models:
m1: first_fam ~ block_number_log_c + bg_image_type + (1 | prolific_id)
m2: first_fam ~ block_number_log_c * total_trial_number_log_c + bg_image_type + (1 | prolific_id)
m3: first_fam ~ block_number_log_c * total_trial_number_log_c * bg_image_type + (1 | prolific_id)
   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
m1    5 3270.4 3298.0 -1630.2   3260.4                     
m2    7 3271.4 3310.0 -1628.7   3257.4 3.0577  2     0.2168
m3   10 3276.7 3331.9 -1628.3   3256.7 0.6569  3     0.8833

summary(m1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: first_fam ~ block_number_log_c + bg_image_type + (1 | prolific_id)
   Data: ABBA_number

REML criterion at convergence: 3278.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1361 -0.5751 -0.0983  0.4734  3.9713 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.2797   0.5289  
 Residual                0.2157   0.4645  
Number of obs: 1852, groups:  prolific_id, 468

Fixed effects:
                     Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)           7.72608    0.02673  465.20617 289.002   <2e-16 ***
block_number_log_c   -0.26953    0.01353 1400.80598 -19.916   <2e-16 ***
bg_image_type1        0.05171    0.02386 1499.25323   2.167   0.0304 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) blc___
blck_nmbr__ -0.001       
bg_img_typ1 -0.008  0.029

plot(effects::allEffects(m1))

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictor block_number_log_c is a one-column matrix that was converted to a
vector

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictor block_number_log_c is a one-column matrix that was converted to a
vector

7.3 animacy

ABBA_animacy_noidentity <- ABBA_animacy_noidentity %>%
  mutate(bg_image_type = ifelse(str_detect(bg_image, "inanimate") == TRUE, "inanimate", "animate"))

ABBA_animacy_noidentity$bg_image_type_c <- as.factor(ABBA_animacy_noidentity$bg_image_type)
ABBA_animacy_noidentity %>% group_by(bg_image_type) %>% rstatix::get_summary_stats(first_fam, type = "common")

# A tibble: 2 × 11
  bg_image_type variable      n   min   max median   iqr  mean    sd    se    ci
  <chr>         <fct>     <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 animate       first_fam   931  6.24  9.92   7.48  1.01  7.58 0.71  0.023 0.046
2 inanimate     first_fam   925  6.22  9.82   7.91  1.10  7.86 0.722 0.024 0.047

contrasts(ABBA_animacy_noidentity$bg_image_type_c) <- contr.sum(2)/2
m1 <- lmerTest::lmer(first_fam ~ block_number_log_c + bg_image_type + (1|prolific_id), data = ABBA_animacy_noidentity)
m2 <- lmerTest::lmer(first_fam ~ block_number_log_c * total_trial_number_log_c + bg_image_type + (1|prolific_id), data = ABBA_animacy_noidentity)
m3 <- lmerTest::lmer(first_fam ~ block_number_log_c * total_trial_number_log_c * bg_image_type + (1|prolific_id), data = ABBA_animacy_noidentity)
anova(m1, m2, m3)

refitting model(s) with ML (instead of REML)

Data: ABBA_animacy_noidentity
Models:
m1: first_fam ~ block_number_log_c + bg_image_type + (1 | prolific_id)
m2: first_fam ~ block_number_log_c * total_trial_number_log_c + bg_image_type + (1 | prolific_id)
m3: first_fam ~ block_number_log_c * total_trial_number_log_c * bg_image_type + (1 | prolific_id)
   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
m1    5 3253.1 3280.8 -1621.6   3243.1                     
m2    7 3255.6 3294.3 -1620.8   3241.6 1.4949  2     0.4736
m3   10 3257.2 3312.5 -1618.6   3237.2 4.4014  3     0.2213

summary(m1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: first_fam ~ block_number_log_c + bg_image_type + (1 | prolific_id)
   Data: ABBA_animacy_noidentity

REML criterion at convergence: 3260.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8180 -0.5850 -0.0789  0.5108  3.7168 

Random effects:
 Groups      Name        Variance Std.Dev.
 prolific_id (Intercept) 0.2654   0.5152  
 Residual                0.2153   0.4640  
Number of obs: 1856, groups:  prolific_id, 468

Fixed effects:
                         Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)               7.78930    0.03103  869.88914 251.050  < 2e-16 ***
block_number_log_c       -0.33821    0.02079 1435.54090 -16.265  < 2e-16 ***
bg_image_typeinanimate   -0.12986    0.03352 1415.01675  -3.874 0.000112 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) blc___
blck_nmbr__ -0.413       
bg_mg_typnn -0.539  0.765

plot(effects::allEffects(m1))

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictor block_number_log_c is a one-column matrix that was converted to a
vector

Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictor block_number_log_c is a one-column matrix that was converted to a
vector

7.4 plot across violation types

exp2_baseline_check <- bind_rows(ABBA_pose, ABBA_number, ABBA_animacy_noidentity)

exp2_baseline_check$bg_image_type <- as.factor(exp2_baseline_check$bg_image_type)
exp2_baseline_check$bg_image_type <- fct_relevel(exp2_baseline_check$bg_image_type, "left", "right", "single", "pair", "inanimate", "animate")

exp2_baseline_check %>%
    ggplot(aes(x = as.factor(bg_image_type), y = first_fam, fill = violation_type)) + 
    geom_boxplot(alpha = 0.5, width = 0.5, outlier.color = "red", outlier.shape = 1, outlier.alpha = 0.05) + 
    geom_point(alpha = 0.02) +
  stat_summary(fun = "mean", geom = "point", position = position_dodge(.9), color = "red") +
    stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.2, position = position_dodge(.9), color = "red") +
    labs(x = "Feature of Familiarization Image", 
         y = "Looking Time at First Familiarization (s)", 
         title = "Baseline Looking for Different Stimuli") +
    theme_classic()