# pulled from 06_main-analysis-mixed-model -- was there other preprocessing in 07_main-analysis_results we need??
merged_data <- read_csv("data/processed_data/03_merged_data_sans_exclusions.csv")
## Rows: 8043 Columns: 61
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (46): participant_id, lab_id, unique_participant_id, trial_type, stimulu...
## dbl (12): trial_num, looking_time_seconds, trial_time_seconds, participant_a...
## lgl  (3): session_exclusion_reason_2, session_exclusion_reason_3, session_ex...
## 
## ℹ 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.
# add necessary variables
proc_data <- merged_data %>%
  rename(trial_type_char = trial_type) |>
  mutate(
    log_looking_time = log(looking_time_seconds),
    multilingual_exposure = as.vector(scale((100 - lang1_exposure) / 100)), # proportion of not L1
    age = as.vector(scale(participant_age_days / 30.42)), # convert age to months
    repetition_char = ifelse(trial_type_char != "fam" & grepl("abb", stimulus_filename), "ABB", "ABA"),
    repetition = ifelse(repetition_char=="ABB", .5, -.5),
    trial_type = ifelse(trial_type_char=="consistent", .5, -.5),
    trial_num_sc = scale(trial_num),
  ) %>% # add repetition variable based on test trials (ABB and ABA)
  rename(experimental_method = method)

# remove familiarisation trials
proc_data <- proc_data %>%
  filter(trial_type_char != "fam") %>%
  filter(!is.infinite(log_looking_time))

models <- list()

Original analysis (full model)

Pruning order:

# Define the function to fit models in the specified pruning order
prune_lmer_model <- function(data) {
  # Shared part of the formula
  shared_formula <- "log_looking_time ~ 1 + 
                   familiarization_rule * trial_type + 
                   age * trial_type + 
                   experimental_method * trial_type + 
                   multilingual_exposure * trial_type + 
                   trial_num_sc * trial_type + 
                   trial_num_sc * age + 
                   repetition * trial_type"
  
  # Define the list of formulas in pruning order, with only the random effects changing
  formulas <- list(
    full_model = paste(shared_formula, "+ (trial_num_sc * trial_type | unique_participant_id) + (trial_type + test_order | lab_id)"),
    prun1_model = paste(shared_formula, "+ (trial_num_sc + trial_type | unique_participant_id) + (trial_type + test_order | lab_id)"),
    prun2_model = paste(shared_formula, "+ (trial_num_sc + trial_type | unique_participant_id) + (test_order | lab_id)"),
    prun3_model = paste(shared_formula, "+ (trial_num_sc + trial_type | unique_participant_id) + (1 | lab_id)"),
    prun4_model = paste(shared_formula, "+ (trial_type | unique_participant_id) + (1 | lab_id)"),
    prun5_model = paste(shared_formula, "+ (1 | unique_participant_id) + (1 | lab_id)"),
    prun6_model = paste(shared_formula, "+ (1 | unique_participant_id)")
  )
  
  # Iterate over the formulas and fit the model
  for (i in seq_along(formulas)) {
    formula <- formulas[[i]]
    
    # Print the formula being tried
    cat(paste0("Trying model: ", names(formulas)[i], "\n"))
    
    # Fit the model
    model <- tryCatch({
      lmer(as.formula(formula), data = data, REML = FALSE,
           control=lmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e5)))
    }, warning = function(w) {
      if(grepl('failed to converge', w$message)) {
        cat(paste0("Model ", names(formulas)[i], " failed to converge with warning: ", w$message, "\n"))
      }
      invokeRestart("muffleWarning")  # Muffle the warning so it doesn't stop execution
      return(NULL)
    }, error = function(e) {
      cat(paste0("Model ", names(formulas)[i], " failed with error: ", e$message, "\n"))
      return(NULL)
    })
    
    # Check for convergence only if model fitting was successful
    if (!is.null(model)) {
      if (!is.null(model@optinfo$derivs)) {
        if(converged(model)) { 
          cat(paste0("Model ", names(formulas)[i], " converged.\n"))
          return(model)
        }
      } else {
        cat(paste0("Model ", names(formulas)[i], " did not converge properly (no derivatives found).\n"))
      }
    }
  }
  
  # Return NULL if no model converges
  cat("No model converged.\n")
  return(NULL)
}

models[["original"]] <- prune_lmer_model(proc_data)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Model full_model failed to converge with warning: Model failed to converge with 2 negative eigenvalues: -1.4e-01 -2.7e-01
## Model full_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Model prun1_model failed to converge with warning: Model failed to converge with 2 negative eigenvalues: -2.2e-02 -2.1e+01
## Model prun1_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Model prun2_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -1.5e+01
## Model prun2_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
broom.mixed::tidy(models[["original"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.627 0.084 19.373 44.578 0.000 1.458 1.796
fixed NA trial_num_sc -0.153 0.007 -22.435 6362.083 0.000 -0.166 -0.140
fixed NA age:trial_num_sc -0.018 0.007 -2.673 6386.438 0.008 -0.032 -0.005

Robustness analyses

The goal here is to run the robustness analysis as defined in the RR, focussing on the following decisions. For each of the proposed robustness analyses, we will: 1. Re-define the data to be used 2. Re-run the main analysis 3. Re-run the pruning 4. Check outcome of main results (comparing the fixed-effect coefficients and significance)

1. Minimal valid looking time per trial of 2s instead of 1s as in the main analyses;

# for each, let's report dataset size before and after change, and change in significance (+coefficient size?)
proc_data_lt2 <- proc_data |> filter(looking_time_seconds >= 2)
models[["min_2s_LT"]] <- prune_lmer_model(proc_data_lt2)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Model full_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -2.6e+01
## Model full_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Model prun1_model failed to converge with warning: Model failed to converge with 2 negative eigenvalues: -3.6e-01 -4.7e+01
## Model prun1_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Model prun2_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -4.6e-01
## Model prun2_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
broom.mixed::tidy(models[["min_2s_LT"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.738 0.070 24.895 43.858 0.000 1.597 1.879
fixed NA trial_num_sc -0.116 0.006 -18.272 5731.104 0.000 -0.129 -0.104
fixed NA age:trial_num_sc -0.017 0.006 -2.673 5750.576 0.008 -0.029 -0.005

This keeps 85.6% of the data (6328 of 7390 trials).

2. Minimum number of trials: include thresholds of 3, 4, 6, 9, 12 respectively to test the robustness of the results

(in all cases also assuming a minimum of 1 same/1 different trial);

get_participants_with_minimum_n_trials <- function(data, minimum) {
  keepers <- data |> group_by(participant_id) |>
    count() |> filter(n >= minimum) |> pull(participant_id)
  data |> filter(participant_id %in% keepers)
}

proc_data_3trials <- get_participants_with_minimum_n_trials(proc_data, minimum=3)
proc_data_4trials <- get_participants_with_minimum_n_trials(proc_data, minimum=4)
proc_data_6trials <- get_participants_with_minimum_n_trials(proc_data, minimum=6)
proc_data_9trials <- get_participants_with_minimum_n_trials(proc_data, minimum=9)
proc_data_12trials <- get_participants_with_minimum_n_trials(proc_data, minimum=12)

models[["min_3_trials"]] <- prune_lmer_model(proc_data_3trials)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Model full_model failed to converge with warning: Model failed to converge with 2 negative eigenvalues: -2.3e+00 -4.1e+00
## Model full_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Model prun1_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -1.4e+00
## Model prun1_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Model prun2_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -2.4e+01
## Model prun2_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
models[["min_4_trials"]] <- prune_lmer_model(proc_data_4trials) 
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Model full_model failed to converge with warning: Model failed to converge with 2 negative eigenvalues: -6.3e+00 -1.1e+01
## Model full_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Model prun1_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -1.4e-01
## Model prun1_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Model prun2_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -1.4e+01
## Model prun2_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
models[["min_6_trials"]] <- prune_lmer_model(proc_data_6trials)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Model prun1_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -9.3e+00
## Model prun1_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
models[["min_9_trials"]] <- prune_lmer_model(proc_data_9trials)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Model full_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -2.8e-01
## Model full_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Model prun1_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -3.0e-01
## Model prun1_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
models[["min_12_trials"]] <- prune_lmer_model(proc_data_12trials)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Model prun1_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -2.4e+00
## Model prun1_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Model prun2_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -2.8e+00
## Model prun2_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
broom.mixed::tidy(models[["min_3_trials"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.627 0.084 19.372 44.188 0.000 1.457 1.796
fixed NA trial_num_sc -0.153 0.007 -22.456 6351.588 0.000 -0.167 -0.140
fixed NA age:trial_num_sc -0.018 0.007 -2.669 6374.087 0.008 -0.032 -0.005 
broom.mixed::tidy(models[["min_4_trials"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.627 0.084 19.424 44.131 0.000 1.458 1.796
fixed NA trial_num_sc -0.153 0.007 -22.488 6352.530 0.000 -0.167 -0.140
fixed NA age:trial_num_sc -0.018 0.007 -2.695 6375.108 0.007 -0.032 -0.005 
broom.mixed::tidy(models[["min_6_trials"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.620 0.085 19.168 41.284 0.000 1.449 1.790
fixed NA trial_num_sc -0.154 0.007 -22.515 6289.480 0.000 -0.168 -0.141
fixed NA age:trial_num_sc -0.019 0.007 -2.781 6297.587 0.005 -0.032 -0.006 
broom.mixed::tidy(models[["min_9_trials"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.638 0.087 18.858 39.515 0.000 1.463 1.814
fixed NA trial_num_sc -0.152 0.007 -22.022 6113.077 0.000 -0.166 -0.139
fixed NA age:trial_num_sc -0.018 0.007 -2.612 6117.428 0.009 -0.032 -0.005 
broom.mixed::tidy(models[["min_12_trials"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.670 0.092 18.230 33.551 0.000 1.484 1.857
fixed NA trial_num_sc -0.148 0.007 -20.184 5406.535 0.000 -0.162 -0.133
fixed NA age:trial_num_sc -0.015 0.007 -2.087 5410.631 0.037 -0.030 -0.001

Keeping participants with at least 3, 4, 6, 9, or 12 trials, this keeps (respectively) 99.9%, 99.8%, 98.8%, 96.2%, or 84.7% of the data.

3. Test-trials within a block: only include infants if they had a pair (one ‘same’, one ‘different’) of valid test-trials within a block of 4 trials.

# iterates over participant data and excludes if any sliding window of 4 doesn't contain variation in trial_type
# (trial_type is 0.5 or -.5, so there is no variation if sum across windowsize 4 is 2 or -2)
get_participants_with_varying_trial_types_in_window <- function(data, window_size=4) {
  ids_to_remove <- c()
  
  for (p in unique(data$unique_participant_id)) {
    this_dat <- data |> filter(unique_participant_id == p) |> arrange(trial_num)
    trial_types <- this_dat$trial_type
    
    # Skip participants with fewer trials than the window size
    if (length(trial_types) < window_size) {
      ids_to_remove <- c(ids_to_remove, p)
      next
    }
    
    # Check for any sliding window with no variation
    no_variation_found <- FALSE
    for (i in 1:(nrow(this_dat) - window_size + 1)) {
      window_sum <- sum(trial_types[i:(i + window_size - 1)])
      if (window_sum == 2 || window_sum == -2) {
        no_variation_found <- TRUE
        break
      }
    }
    if (no_variation_found) ids_to_remove <- c(ids_to_remove, p)
  }
  data |> filter(!unique_participant_id %in% ids_to_remove)
}


proc_data_variation_in_win4 <- get_participants_with_varying_trial_types_in_window(proc_data)
models[["varying_trials"]] <- prune_lmer_model(proc_data_variation_in_win4)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Model full_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -1.6e+01
## Model full_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Model prun2_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -9.7e+00
## Model prun2_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
broom.mixed::tidy(models[["varying_trials"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.622 0.083 19.430 40.920 0 1.454 1.791
fixed NA trial_num_sc -0.157 0.007 -22.281 5950.078 0 -0.171 -0.143

4. Second-session infants: in the main model we include both first- and second-session infants; here we exclude second-session infants;

proc_data_session1 <- proc_data |> filter(second_session == "no")

models["session1"] <- prune_lmer_model(proc_data_session1)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
## Warning in `[<-`(`*tmp*`, "session1", value = new("lmerModLmerTest",
## vcov_varpar = structure(c(0.0006506983799297, : implicit list embedding of S4
## objects is deprecated
broom.mixed::tidy(models[["session1"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.622 0.088 18.461 43.319 0.000 1.445 1.799
fixed NA trial_num_sc -0.158 0.007 -21.768 5548.394 0.000 -0.173 -0.144
fixed NA age:trial_num_sc -0.016 0.007 -2.242 5570.679 0.025 -0.030 -0.002

This keeps 86.2% of the data (6373 of 7390 trials).

5. Indistinguishable syllables: Should there be labs that report that some of the syllables are non-distinguishable in the dominant language in some infants we will exclude those infants here.

# waiting to hear if there are any

6. Trial spill-over effects: exclude all trials that follow an initial missing trial due to infant factors as some of the effects may spill-over from the missing trial.

# ToDo: test! look at output for specific participants. George wrote it and it worked the first time, so he doesn't trust it. ;P
exclude_all_trials_after_a_missing_trial <- function(dat) {
  for(s in unique(dat$unique_participant_id)) {
    s_ind <- which(dat$unique_participant_id==s)
    missing_subject_trial_inds = sort(setdiff(1:12, dat[s_ind,]$trial_num))
    if(length(missing_subject_trial_inds)>0) { 
      trial_indices_to_exclude = missing_subject_trial_inds[1]:12 
      dat <- dat[!(dat$unique_participant_id == s & dat$trial_num %in% trial_indices_to_exclude), ]
    }
  }
  return(dat)
}

no_spillover_data <- exclude_all_trials_after_a_missing_trial(proc_data)
models[["no_spillover_trials"]] <- prune_lmer_model(no_spillover_data)
## Trying model: full_model
## boundary (singular) fit: see help('isSingular')
## Model full_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -7.1e+01
## Model full_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun1_model
## boundary (singular) fit: see help('isSingular')
## Model prun1_model failed to converge with warning: Model failed to converge with 2 negative eigenvalues: -9.2e-03 -4.2e+01
## Model prun1_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun2_model
## boundary (singular) fit: see help('isSingular')
## Model prun2_model failed to converge with warning: Model failed to converge with 1 negative eigenvalue: -1.9e+00
## Model prun2_model failed with error: no 'restart' 'muffleWarning' found
## Trying model: prun3_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun4_model
## boundary (singular) fit: see help('isSingular')
## Trying model: prun5_model
## Model prun5_model converged.
broom.mixed::tidy(models[["no_spillover_trials"]], conf.int = TRUE) |> 
  filter(effect=="fixed", p.value<.05) |>
  kableExtra::kable(digits=3)
effect group term estimate std.error statistic df p.value conf.low conf.high
fixed NA (Intercept) 1.616 0.082 19.740 42.563 0.000 1.450 1.781
fixed NA trial_num_sc -0.147 0.007 -20.642 5995.250 0.000 -0.161 -0.133
fixed NA age:trial_num_sc -0.016 0.007 -2.247 6027.070 0.025 -0.030 -0.002

This keeps 92.6% of the data (6846 of 7390 trials).

7. Missing eye-tracker data: In the main analysis, trials are included when at least 50% of trial duration has recorded data. Here we will set this percentage to 70, and 90% respectively.

# what is the variable that indicates amount of missing eye-tracker data per trial? (perhaps missing from proc_data)
# only experimental_method=="singlescreen_eyetracker" ?

Save models

saveRDS(models, "08_robustness_models.rds")

Comparing coefficients across robustness analyses

# extract fixed effects from all our models and collate
collate_model_coefficients <- function(models) {
  model_coeffs <- list()
  
  for (i in seq_along(models)) {
    # extract fixed effects
    if (!is.null(models[[i]])) {
      coeffs <- fixef(models[[i]])
      model_coeffs[[i]] <- coeffs
    } else {
      model_coeffs[[i]] <- NA
    }
  }
  # Combine fixed effects into a data frame
  coeffs_df <- do.call(rbind, lapply(model_coeffs, function(x) t(as.data.frame(x))))
  coeffs_df <- as.data.frame(coeffs_df)
  rownames(coeffs_df) <- names(models)
  
  return(coeffs_df)
}

coefficients_table <- collate_model_coefficients(models)

coefficients_table |> kableExtra::kable(digits=3)
(Intercept) familiarization_rule trial_type age experimental_methodsinglescreen_cf experimental_methodsinglescreen_eyetracker multilingual_exposure trial_num_sc repetition familiarization_rule:trial_type trial_type:age trial_type:experimental_methodsinglescreen_cf trial_type:experimental_methodsinglescreen_eyetracker trial_type:multilingual_exposure trial_type:trial_num_sc age:trial_num_sc trial_type:repetition
original 1.627 -0.028 0.011 -0.009 0.085 -0.183 0.023 -0.153 -0.017 0.018 0.009 -0.033 -0.010 0.008 0.003 -0.018 -0.114
min_2s_LT 1.738 -0.007 0.009 -0.008 0.045 -0.144 0.019 -0.116 -0.065 -0.038 -0.005 -0.034 -0.009 0.008 0.002 -0.017 -0.017
min_3_trials 1.627 -0.027 0.012 -0.009 0.086 -0.180 0.023 -0.153 -0.017 0.018 0.009 -0.034 -0.010 0.009 0.001 -0.018 -0.112
min_4_trials 1.627 -0.028 0.012 -0.008 0.085 -0.179 0.023 -0.153 -0.017 0.019 0.008 -0.034 -0.011 0.009 0.002 -0.018 -0.113
min_6_trials 1.620 -0.030 0.015 -0.005 0.093 -0.158 0.026 -0.154 -0.017 0.019 0.008 -0.036 -0.015 0.010 0.001 -0.019 -0.111
min_9_trials 1.638 -0.028 0.013 -0.006 0.071 -0.173 0.026 -0.152 0.010 0.047 0.009 -0.035 -0.013 0.010 0.000 -0.018 -0.108
min_12_trials 1.670 -0.030 0.010 -0.001 0.063 -0.190 0.030 -0.148 0.011 0.046 0.011 -0.031 -0.008 0.012 0.010 -0.015 -0.100
varying_trials 1.622 -0.031 0.006 -0.012 0.097 -0.166 0.019 -0.157 -0.038 -0.003 0.010 -0.027 -0.010 0.008 0.002 -0.014 -0.142
session1 1.622 -0.015 0.010 -0.005 0.088 -0.171 0.021 -0.158 -0.022 0.023 0.006 -0.025 -0.012 0.009 -0.006 -0.016 -0.112
no_spillover_trials 1.616 -0.032 0.003 -0.009 0.118 -0.140 0.027 -0.147 -0.011 0.024 0.009 -0.025 -0.006 0.016 0.010 -0.016 -0.134 
# Reshape the data for heatmap plotting
coefs_long <- coefficients_table |>
  rownames_to_column(var = "Model") |>
  pivot_longer(cols = -Model, names_to = "Coefficient", values_to = "Value")

coefs_long <- coefs_long |> filter(Coefficient!="(Intercept)")

coefs_long |> filter(Coefficient!="(Intercept)") |>
  ggplot(aes(x = Coefficient, y = Model, fill = Value)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0, 
                       limit = c(min(coefs_long$Value), max(coefs_long$Value))) +
  labs(x = "Coefficient", y = "Model", fill = "Coefficient Value") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))