Main visual precision descriptive analyses

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1          ✔ readr     2.2.0     
✔ forcats   1.0.1          ✔ stringr   1.6.0     
✔ ggplot2   4.0.3.9000     ✔ tibble    3.3.1     
✔ lubridate 1.9.5          ✔ tidyr     1.3.2     
✔ purrr     1.2.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
here() starts at /Users/tspri/Documents/visual-precision

Note: The package "relayer" is highly experimental. Use at your own risk.

Loading required package: viridisLite

Loading required package: lme4

Loading required package: Matrix


Attaching package: 'Matrix'


The following objects are masked from 'package:tidyr':

    expand, pack, unpack



Attaching package: 'lmerTest'


The following object is masked from 'package:lme4':

    lmer


The following object is masked from 'package:stats':

    step



Attaching package: 'cowplot'


The following object is masked from 'package:lubridate':

    stamp



Attaching package: 'scales'


The following object is masked from 'package:viridis':

    viridis_pal


The following object is masked from 'package:purrr':

    discard


The following object is masked from 'package:readr':

    col_factor



Attaching package: 'patchwork'


The following object is masked from 'package:cowplot':

    align_plots



Attaching package: 'rlang'


The following objects are masked from 'package:purrr':

    flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
    flatten_raw, invoke, splice

Warning in readLines(file): incomplete final line found on
'/Users/tspri/Documents/visual-precision/.env'

Load data

PROCESSED_DATA_PATH = here("data", DATA_FOLDER, "processed_data")
looking_time_resampled_clean <- read.csv(file.path(PROCESSED_DATA_PATH, "level-looks_data.csv"))
trial_metadata <- read.csv(here("data","metadata","level-trialtype_data.csv"))
trial_summary_data <- read.csv(file.path(PROCESSED_DATA_PATH,"level-trials_data.csv"))
all_similarities <- read.csv(here("data", "embeddings", "similarities-all_data.csv"))
openclip_similarities <- read.csv(here("data", "embeddings", "similarities-openclip_data.csv"))
layerwise <- read.csv(here("data/embeddings/similarities-layerwise_data.csv"))
all_sims <- trial_metadata |> left_join(layerwise, by=c("Trials.targetImage"="text1", "Trials.distractorImage"="text2")) |> left_join(all_similarities, by=c("Trials.targetImage"="text1", "Trials.distractorImage"="text2"))

# Filter by section if PROJECT_VERSION specifies a particular sample
if (!is.null(SECTION_FILTER)) {
  section_trials <- trial_metadata |> filter(section == SECTION_FILTER) |> distinct(Trials.trialID)
  trial_summary_data   <- trial_summary_data   |> semi_join(section_trials, by = "Trials.trialID")
  looking_time_resampled_clean <- looking_time_resampled_clean |> semi_join(section_trials, by = "Trials.trialID")
  all_looking_times    <- all_looking_times    |> semi_join(section_trials, by = "Trials.trialID")
}

looking_data_summarized <- trial_summary_data |>
  filter(trial_exclusion == 0 & exclude_participant == 0 & exclude_participant_insufficient_data == 0 & SubjectInfo.subjID != "VVI128") |> dplyr::select(-section) |>
  left_join(trial_metadata) |>
  arrange(AoA_Est_target) 

Joining with `by = join_by(Trials.trialID, Trials.targetImage,
Trials.distractorImage, Trials.imagePair)`

Warning in left_join(dplyr::select(filter(trial_summary_data, trial_exclusion == : Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 1 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

median_age <- looking_data_summarized %>%
  distinct(SubjectInfo.subjID, .keep_all = TRUE) %>%  # Keep only distinct subjects
  summarize(overall_median_age = median(SubjectInfo.testAge, na.rm = TRUE))  %>%
  pull(overall_median_age) / 30

median_age <- round(median_age, 2)

aoa_summary <- trial_metadata |>
  summarize(median_aoa = median(AoA_Est_target), 
            mean_aoa = mean(AoA_Est_target))
median_aoa <- aoa_summary$median_aoa

# Total counts per section
total_counts <- trial_summary_data %>%
  distinct(SubjectInfo.subjID, section) %>%
  count(section, name = "total")

# Filtered counts per section
filtered_counts <- looking_data_summarized %>%
  distinct(SubjectInfo.subjID, section) %>%
  count(section, name = "filtered") 

# Combine and compute percentage
summary_table <- total_counts %>%
  left_join(filtered_counts, by = "section") %>%
  mutate(
    filtered = coalesce(filtered, 0),
    percentage_filtered = round((filtered / total) * 100, 2)
  )

summary_table

  section total filtered percentage_filtered
1 sample1   124       90               72.58
2 sample2   106       80               75.47

util function to see what the order of trials were for individual participants

order_of_trials <- function(data, subjID) {
  curr_order <- data |>
  filter(SubjectInfo.subjID == subjID) |>
  distinct(Trials.trialID, .keep_all=TRUE) |>
  arrange(Trials.ordinal) |>
  dplyr::select(Trials.ordinal, Trials.trialID, Trials.targetImage, Trials.distractorImage, Trials.trialType, Trials.leftImage, Trials.rightImage, Trials.targetAudio)
  write_csv(curr_order, paste0(subjID,".csv"))
}

#order_of_trials(all_looking_times, "L7Y4Y6")

Overall timecourse plot of proportion target looking

#summarizing within subject for each time point
summarize_subj <- looking_time_resampled_clean %>%
  filter(trial_exclusion == 0 & exclude_participant ==0 & exclude_participant_insufficient_data == 0 & SubjectInfo.subjID != "VVI128") %>%
  filter(!is.na(accuracy_transformed)) %>%
  group_by(SubjectInfo.subjID, time_normalized_corrected, SubjectInfo.testAge) %>%
  summarized_data(., "SubjectInfo.subjID", "accuracy_transformed", c("SubjectInfo.testAge", "time_normalized_corrected")) %>%
  rename(mean_accuracy = mean_value)

Warning: There were 747 warnings in `summarize()`.
The first warning was:
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 1: `SubjectInfo.subjID = "VVI006"`, `SubjectInfo.testAge = 510`,
  `time_normalized_corrected = -3267`.
Caused by warning in `qt()`:
! NaNs produced
ℹ Run `dplyr::last_dplyr_warnings()` to see the 746 remaining warnings.

#summarizing across subjects for each time point
summarize_across_subj <- summarize_subj %>%
  group_by(time_normalized_corrected) %>%
  summarized_data(., "time_normalized_corrected", "mean_accuracy", c("time_normalized_corrected")) |>
  rename(accuracy = mean_value)

Warning: There were 8 warnings in `summarize()`.
The first warning was:
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 303: `time_normalized_corrected = 4667`.
Caused by warning in `qt()`:
! NaNs produced
ℹ Run `dplyr::last_dplyr_warnings()` to see the 7 remaining warnings.

looking_times <- ggplot(summarize_across_subj,aes(time_normalized_corrected,accuracy))+
  xlim(-2000,4000)+
  geom_errorbar(aes(ymin=accuracy-ci,ymax=accuracy+ci),width=0, alpha=0.2)+
  #geom_point(alpha=0.2)+
    geom_smooth(method="gam")+
  geom_vline(xintercept=0,size=1.5)+
  geom_hline(yintercept=0.5,size=1.2,linetype="dashed")+
  geom_vline(xintercept=300,linetype="dotted")+
  ylim(0,1)+
  xlab("Time (normalized to target word onset) in ms")+
  ylab("Proportion Target Looking") 

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

looking_times

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

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

ggsave(here("figures",PROJECT_VERSION,"prop_looking_across_time.png"),looking_times,width=9,height=6,bg = "white")

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

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

Descriptive plots

Usable trial stats

#Overall baseline-corrected proportion target looking by condition
looking_data_exclusions <- trial_summary_data |>
  mutate(trial_exclusion_reason = as.factor(ifelse(is.na(trial_exclusion_reason), "included", trial_exclusion_reason))) |>
  mutate(trial_exclusion_reason = factor(trial_exclusion_reason, levels=c(setdiff(unique(trial_exclusion_reason), "included"), "included")))

# Compute the number of included trials per subject
included_counts <- looking_data_exclusions |>
  filter(trial_exclusion_reason == "included") |>
  count(SubjectInfo.subjID, name = "included_trials")


excluded_subjects <- setdiff(unique(looking_data_exclusions$SubjectInfo.subjID), 
                             unique(looking_data_summarized$SubjectInfo.subjID))

# Ensure all subjects are included and replace NA with 0
exclusions_by_subject <- looking_data_exclusions |>
  summarize(N = n(), .by = c(SubjectInfo.subjID, trial_exclusion_reason)) |>
  left_join(included_counts, by = "SubjectInfo.subjID") |>
  mutate(included_trials = replace_na(included_trials, 0)) |>  # Replace NA values
  # Ordering subject IDs so that the 'red' colored label to mark excluded subjects is placed on the subjects who are actually excluded
  mutate(SubjectInfo.subjID = factor(SubjectInfo.subjID, 
                                     levels = unique(SubjectInfo.subjID[order(-included_trials)])),
         color = ifelse(SubjectInfo.subjID %in% excluded_subjects, "red", "black")) |>
  arrange(color)

looking_data_exclusions_combined <- looking_data_exclusions |>
  summarize(N = n(),
    .by=c(trial_exclusion_reason))
# Format labels with HTML <span> for colors

exclusions_color_labels <- exclusions_by_subject |> 
  distinct(SubjectInfo.subjID, .keep_all = TRUE) |>
  mutate(colored_labels = paste0("<span style='color:", color, "'>", SubjectInfo.subjID, "</span>")) |>
  arrange(desc(included_trials)) |>
  mutate(colored_labels = factor(colored_labels, levels = colored_labels))

ggplot(exclusions_by_subject, aes(x = SubjectInfo.subjID, y = N, fill = trial_exclusion_reason)) +
  geom_bar(stat = "identity", position = "stack") +  
  labs(x = "Subject ID", y = "Count (N)", title = "Trial Exclusions by Subject", fill = "Trial Exclusion Reason") +
  scale_x_discrete(labels = exclusions_color_labels$colored_labels) +  # Apply colored labels
    theme(
    axis.text.x = ggtext::element_markdown(angle = 45, hjust = 1),  # Enable HTML-styled labels
    strip.text = element_text(size = 5)
  ) +
  scale_fill_brewer(palette = "Set3")

#TODO: we have a participant who participated twice but their first session only had 2 trials. figure out how to deal with this/discard first session trials. Also have to figure out how to deal with participants who are excluded by default

Estimating item and subject-level noise

target_looking_item_subject_level <- summarized_data(looking_data_summarized, "Trials.targetImage", "corrected_target_looking", c("SubjectInfo.subjID", "SubjectInfo.testAge", "AoA_Est_target"))

Warning: There were 1256 warnings in `summarize()`.
The first warning was:
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 89: `Trials.targetImage = "airplane"`, `SubjectInfo.subjID =
  "VVI140"`, `SubjectInfo.testAge = 630`, `AoA_Est_target = 3.94`.
Caused by warning in `qt()`:
! NaNs produced
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1255 remaining warnings.

target_looking_item_level <- summarized_data(target_looking_item_subject_level |> rename(mean_target_looking = mean_value), "Trials.targetImage", "mean_target_looking", "AoA_Est_target")

target_looking_item_age_level <- summarized_data(target_looking_item_subject_level |> rename(mean_target_looking = mean_value) |> add_age_split(), "Trials.targetImage", "mean_target_looking", c("AoA_Est_target", "age_half"))

target_looking_trial_subject_level <- summarized_data(looking_data_summarized, "Trials.trialID", "corrected_target_looking", c("SubjectInfo.subjID", "SubjectInfo.testAge", "AoA_Est_distractor", "AoA_Est_target", "Trials.targetImage", "Trials.distractorImage", "Trials.imagePair"))

Warning: There were 2136 warnings in `summarize()`.
The first warning was:
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 157: `Trials.trialID = "easy-balloon-dresser"`, `SubjectInfo.subjID
  = "VVI139"`, `SubjectInfo.testAge = 630`, `AoA_Est_distractor = 4.28`,
  `AoA_Est_target = 4.37`, `Trials.targetImage = "balloon"`,
  `Trials.distractorImage = "dresser"`, `Trials.imagePair = "balloon-dresser"`.
Caused by warning in `qt()`:
! NaNs produced
ℹ Run `dplyr::last_dplyr_warnings()` to see the 2135 remaining warnings.

target_looking_trial_level <- summarized_data(target_looking_trial_subject_level |> rename(mean_target_looking = mean_value), "Trials.trialID", "mean_target_looking", c("AoA_Est_distractor", "AoA_Est_target", "Trials.targetImage", "Trials.distractorImage", "Trials.imagePair"))

target_looking_subject_level <- summarized_data(looking_data_summarized, "SubjectInfo.subjID", "corrected_target_looking", "SubjectInfo.testAge")

Item-level performance

item_performances <- ggplot(target_looking_item_level, aes(reorder(Trials.targetImage, mean_value), mean_value)) +
   geom_hline(yintercept=0,linetype="dashed")+
  geom_errorbar(aes(ymin=lower_ci,ymax=upper_ci),width=0, alpha=0.6)+
  (geom_point(aes(aoa_size=AoA_Est_target)) |> rename_geom_aes(new_aes = c("size" = "aoa_size")))+
  (geom_jitter(data=target_looking_item_subject_level, aes(x=Trials.targetImage, y=mean_value, color=SubjectInfo.subjID, age_size = SubjectInfo.testAge/30), alpha=0.3, width=0.2) |> rename_geom_aes(new_aes = c("size" = "age_size"))) +
  xlab("Target image")+
  ylab("Proportion of target looking")+
  ggtitle("Mean proportion of target looking across target items") +
  theme(axis.title.x = element_text(face="bold", size=15, vjust=-1),
        axis.text.x  = element_text(size=10,angle=45,vjust=0.5),
        axis.title.y = element_text(face="bold", size=15),
        axis.text.y  = element_text(size=10),
        strip.text.x = element_text(size = 10,face="bold"),
        plot.margin = margin(t = 10, r = 10, b = 30, l = 10)
        ) +
  scale_y_continuous(breaks = seq(-1, 1, by = 0.2)) +
  scale_size_c(aesthetics = "age_size",name = "Age of participant in months", range=c(2,4), guide = guide_legend(order = 2)) +
  scale_size_c( aesthetics = "aoa_size",name = "Est. AoA of target words in years", guide = guide_legend(order = 1)) +
  scale_color_viridis_d(name = "Participant IDs",option="D") +
  guides(
    color = "none"
  ) 

Warning in geom_point(aes(aoa_size = AoA_Est_target)): Ignoring unknown
aesthetics: aoa_size

Warning in geom_jitter(data = target_looking_item_subject_level, aes(x =
Trials.targetImage, : Ignoring unknown aesthetics: age_size

Warning: The `scale_name` argument of `continuous_scale()` is deprecated as of ggplot2
3.5.0.

Warning: The `trans` argument of `continuous_scale()` is deprecated as of ggplot2 3.5.0.
ℹ Please use the `transform` argument instead.

item_performances

ggsave(here("figures",PROJECT_VERSION,"item_performances.png"),item_performances,width=15,height=10,bg = "white")

trial_performances <- ggplot(summarized_data(looking_data_summarized, "Trials.trialID", "corrected_target_looking", c("Trials.trialID", "AoA_Est_target")), aes(reorder(Trials.trialID, mean_value), mean_value)) +
   geom_hline(yintercept=0,linetype="dashed")+
  geom_errorbar(aes(ymin=lower_ci,ymax=upper_ci),width=0, alpha=0.2)+
  (geom_point(aes(aoa_size=AoA_Est_target)) |> rename_geom_aes(new_aes = c("size" = "aoa_size")))+
  (geom_jitter(data=looking_data_summarized,aes(x=Trials.trialID, y=corrected_target_looking, color=SubjectInfo.subjID, age_size = SubjectInfo.testAge/30), alpha=0.3, width=0.2) |> rename_geom_aes(new_aes = c("size" = "age_size"))) +
  xlab("Trial type")+
  ylab("Proportion of target looking")+
  ggtitle("Mean proportion of target looking across trial types") +
  scale_size_c(aesthetics = "age_size",name = "Age of participant in months", range=c(2,4), guide = guide_legend(order = 2)) +
  scale_size_c( aesthetics = "aoa_size",name = "Est. AoA of target words in years", guide = guide_legend(order = 1)) +
  scale_color_viridis_d(name = "Participant IDs",option="D") +
  coord_cartesian(ylim = c(-1, 1)) +
    theme(axis.title.x = element_text(face="bold", size=15, vjust=-1),
        axis.text.x  = element_text(size=8,angle=45,hjust=1),
        axis.title.y = element_text(face="bold", size=15),
        axis.text.y  = element_text(size=15),
        strip.text.x = element_text(size = 10,face="bold"),
        aspect.ratio = 1,
        plot.margin = margin(t = 10, r = 10, b = 30, l = 10)
        ) +
  guides(
    color = "none"
  ) 

Warning in geom_point(aes(aoa_size = AoA_Est_target)): Ignoring unknown
aesthetics: aoa_size

Warning in geom_jitter(data = looking_data_summarized, aes(x = Trials.trialID,
: Ignoring unknown aesthetics: age_size

trial_performances

####Item-level timecourses

# First, let's make sure we have the correct data
summarize_target <- looking_time_resampled_clean %>%
  filter(trial_exclusion == 0 & exclude_participant == 0 & exclude_participant_insufficient_data == 0) %>%
  #dplyr::select(-section) |>
  left_join(trial_metadata) %>%
  filter(!is.na(accuracy_transformed)) %>%
  group_by(Trials.targetImage, time_normalized_corrected) %>%
  #filter(Trials.targetImage %in% c("acorn", "bulldozer", "cheese", "potato", "squirrel", "snail", "turkey", "turtle")) %>%
  summarized_data(., "Trials.targetImage", "accuracy_transformed", c("Trials.targetImage", "AoA_Est_target", "time_normalized_corrected")) %>%
  rename(mean_accuracy = mean_value)

Joining with `by = join_by(Trials.trialID, Trials.targetImage,
Trials.distractorImage, Trials.imagePair)`

Warning in left_join(., trial_metadata): Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 8 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

Warning: There were 418 warnings in `summarize()`.
The first warning was:
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 246: `Trials.targetImage = "airplane"`, `AoA_Est_target = 3.94`,
  `time_normalized_corrected = -4367`.
Caused by warning in `qt()`:
! NaNs produced
ℹ Run `dplyr::last_dplyr_warnings()` to see the 417 remaining warnings.

# Create a mapping of target images to their AoA values
target_aoa_mapping <- summarize_target %>%
  dplyr::select(Trials.targetImage, AoA_Est_target) %>%
  distinct() %>%
  arrange(AoA_Est_target)  # Sort by AoA

# Convert Trials.targetImage to a factor with levels ordered by AoA (descending)
summarize_target <- summarize_target %>%
  mutate(Trials.targetImage = factor(Trials.targetImage, 
                                     levels = target_aoa_mapping$Trials.targetImage))
aoa_values <- target_aoa_mapping$AoA_Est_target
# Normalize aoa_values to range between 0 and 1
aoa_scaled <- rescale(aoa_values, to = c(0, 1))

# Generate a high-resolution color gradient (e.g., 1000 colors)
num_colors <- 1000
color_palette <- viridis(num_colors, option = "E",alpha=0.8)

# Map scaled values to colors using interpolation
color_values <- color_palette[ceiling((1-aoa_scaled) * (num_colors - 1)) + 1]


# Create named color vector
target_colors <- setNames(color_values, target_aoa_mapping$Trials.targetImage)


# Create the plot with a discrete legend but gradient colors
looking_times_target <- ggplot(summarize_target, 
                              aes(time_normalized_corrected, 
                                  mean_accuracy, 
                                  color = AoA_Est_target, group=Trials.targetImage)) +  # Group by target image
  annotate("rect", xmin = 300, xmax = 3500, ymin = 0, ymax = 1, 
          fill = "gray95", alpha = 0.5) +  
 annotate("rect", xmin = -2000, xmax = 0, ymin = 0, ymax = 1, 
           fill = "gray95", alpha = 0.5) + 

  coord_cartesian(xlim = c(-2000, 3500), ylim = c(0, 1)) +
     geom_line (alpha=0.2, size=0.8) +
   geom_smooth(data=summarize_across_subj,
              aes(x=time_normalized_corrected,
                  y=accuracy,
                  ymin=lower_ci,  # Assuming 'se' is your standard error column
                  ymax=upper_ci,  # Adjust these to match your error metric
                  group=NA),
              stat="identity",  # Use identity instead of a smoothing method
              color="black",
              fill="gray40",
              size=1.5) +
   scale_color_viridis_c(option="B",direction=-1,name="Estimated target\nword AoA") +
  guides(color = guide_colorbar(reverse = TRUE))+
  geom_vline(xintercept = 0, size = 1.5) +
  geom_hline(yintercept = 0.5, size = 1.2, linetype = "dashed") +
  geom_vline(xintercept = 300, linetype = "dotted", size=1.2) +
  xlab("Time (normalized to target word onset) in ms") +
  ylab("Proportion\nTarget Looking") +
    # Add arrows
  annotate("segment", x = -2000, xend = -50, y = 0.85, yend = 0.85, 
           arrow = arrow(length = unit(0.2, "cm"), type="closed", ends="both"), size = 1, color = "black") + 
  annotate("segment", x = 340, xend = 3490, y = 0.85, yend = 0.85, 
           arrow = arrow(length = unit(0.2, "cm"),type="closed", ends="both"), size = 1, color = "black") + 

  # Add text above arrows
  annotate("text", x = -1000, y = 0.93, label = "Baseline window", 
           size = 6, fontface = "bold") +
  annotate("text", x = 2000, y = 0.93, label = "Critical window", 
           size = 6, fontface = "bold") +
  theme_minimal() +
  scale_x_continuous(breaks=seq(-2000,3500,by=1000)) +
  scale_y_continuous(breaks=seq(0,1,by=0.5))+
  theme(
    text = element_text(size=10,face = "bold"),              # All text bold
    # Increase space between x-axis and its title
    axis.title.x = element_text(
      face = "bold", 
      size = 18,
      margin = margin(t = 10, r = 0, b = 0, l = 0)
    ),
    
    # Increase space between y-axis and its title
    axis.title.y = element_text(
      face = "bold", 
      size= 18,
      margin = margin(t=0,r = 5, b = 0, l = 0)
    ),      
    axis.text = element_text(size=14,face = "bold"),         # Axis text bold
    legend.title = element_text(size = 14, face = "bold"),  # Small bold legend title
    legend.text = element_text(size = 14, face = "bold")    # Small bold legend text
  )

# Display the plot
looking_times_target

ggsave(here("figures",PROJECT_VERSION,"timecourse_plot.svg"),looking_times_target, device="pdf",width=8.5,height=3.3,bg = "white")
ggsave(here("figures",PROJECT_VERSION,"timecourse_plot.png"),looking_times_target,width=9.5,height=5.3,bg = "white")

Age plots

Age-split timecourses

summarize_across_age_halves <- summarize_subj |> add_age_split() |>
  group_by(time_normalized_corrected, age_half) |>
  dplyr::summarize(n=n(),
            non_na_n = sum(!is.na(mean_accuracy)),
            mean_value=mean(mean_accuracy,na.rm=TRUE),
            sd_accuracy=sd(mean_accuracy,na.rm=TRUE),
            se_accuracy=sd_accuracy/sqrt(n),
            ci = ifelse(n > 1, qt(0.975, n - 1) * sd_accuracy / sqrt(n), NA)) %>%
  filter(n > 5)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by time_normalized_corrected and age_half.
ℹ Output is grouped by time_normalized_corrected.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(time_normalized_corrected, age_half))` for
  per-operation grouping (`?dplyr::dplyr_by`) instead.

looking_times_age <- ggplot(summarize_across_age_halves,aes(time_normalized_corrected,mean_value,color=age_half))+
  xlim(-2000,4000)+
  geom_errorbar(aes(ymin=mean_value-ci,ymax=mean_value+ci),width=0, alpha=0.2)+
  #geom_point(alpha=0.2)+
    geom_smooth(method="gam")+
  geom_vline(xintercept=0,size=1.5)+
  geom_hline(yintercept=0.5,size=1.2,linetype="dashed")+
  geom_vline(xintercept=300,linetype="dotted")+
  ylim(0,1)+
  xlab("Time (normalized to target word onset) in ms")+
  ylab("Proportion Target Looking") +
  labs(title="Proportion of looking time across age and time", caption = (paste0("14-24 month olds age split at median=",median_age," months"))) +
  scale_color_brewer(palette = "Set1", name="Age half")

looking_times_age

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

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

ggsave(here("figures",PROJECT_VERSION,"prop_looking_across_time_age.png"),looking_times,width=9,height=6,bg = "white")

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

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

Age-split performance

# note: using shorter window just because there's clearly no age-related effect with longer window 
participant_age_halves = summarized_data(looking_data_summarized |> left_join(target_looking_item_level) |> rename(mean_looking = mean_value) |>  add_age_split() |> mutate(acc = corrected_target_looking), "age_half", "acc", 
                                         c("age_half", "SubjectInfo.subjID", "SubjectInfo.testAge")) 

Joining with `by = join_by(Trials.targetImage, AoA_Est_target)`

participant_age_halves$age_half <- factor(participant_age_halves$age_half, levels = rev(levels(factor(participant_age_halves$age_half))))

participant_age_halves_summarized = summarized_data(participant_age_halves |> rename(looking_time = mean_value), "age_half", "looking_time", c("age_half"))
accuracy_age <- ggplot(participant_age_halves,
       aes(x = age_half, y = mean_value, color = age_half)) +
    geom_violin(aes(color=age_half), alpha=0.3) +
  geom_jitter(size = 3, alpha = 0.3, width=0.2) +
  geom_hline(yintercept=0,size=1.2,linetype="dashed")+ 
  geom_point(data = participant_age_halves_summarized, size = 3) +
  geom_linerange(data = participant_age_halves_summarized, aes(ymin = mean_value - ci, ymax = mean_value + ci), alpha = 0.9) + 
  scale_color_brewer(palette = "Set1", name="Age half", direction=-1) +  # Using RColorBrewer for colors
  ylab("Baseline-corrected proportion of target looking") +
  xlab("Age half of child") +
  ggpubr::stat_cor(method = "pearson") +
  labs(caption = (paste0("14-24 month olds age split at median=",median_age, " months")), title="Word recognition acuracy by age")

Registered S3 method overwritten by 'car':
  method           from
  na.action.merMod lme4

accuracy_age

ggsave(here("figures",PROJECT_VERSION,"accuracy_age.png"),accuracy_age,width=9,height=6,bg = "white")

Age continuous

participant_age_jittered <- participant_age_halves %>%
  rowwise() %>%
  mutate(x_jittered = SubjectInfo.testAge/30 + runif(n(), -0.25, 0.25))
jitterer <- position_jitter(width = .5,seed=1)
age_continuous <- ggplot(participant_age_jittered, aes(x = x_jittered, y = mean_value)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(size = 4, alpha = 0.4, color="#2171B5") +
  geom_smooth(method="lm",color="#2171B5") +
  geom_linerange(aes(ymin = mean_value - ci, ymax = mean_value + ci), alpha = 0.1) + 
  scale_y_continuous(breaks = seq(-0.4, 0.3, by = 0.1), limits = c(-0.3, 0.3)) +
  ggpubr::stat_cor()+
  labs(x = "Age in months", y = "Baseline-corrected proportion target looking") +
  #ggtitle("Corrected proportion of target looking by age") +
  theme_minimal()
age_continuous

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 18 rows containing missing values or values outside the scale range
(`geom_segment()`).

participant_age_jittered <- participant_age_halves %>%
  rowwise() %>%
  mutate(x_jittered = SubjectInfo.testAge/30 + runif(n(), -0.25, 0.25))
jitterer <- position_jitter(width = .5,seed=1)
age_continuous_formatted <- ggplot(participant_age_jittered, aes(x = x_jittered, y = mean_value)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(size = 4, alpha = 0.4, color="#2171B5") +
  geom_smooth(method="lm",color="#2171B5") +
  geom_linerange(aes(ymin = mean_value - ci, ymax = mean_value + ci), alpha = 0.1) + 
  scale_y_continuous(breaks = seq(-0.4, 0.3, by = 0.1), limits = c(-0.3, 0.3)) +
  labs(x = "Age in months", y = "Baseline-corrected\nproportion target looking") +
  #ggtitle("Corrected proportion of target looking by age") +
  theme(
    text = element_text(size=10,face = "bold"),              # All text bold
    # Increase space between x-axis and its title
    axis.title.x = element_text(
      face = "bold", 
      size = 18,
      margin = margin(t = 10, r = 0, b = 0, l = 0)
    ),
    
    # Increase space between y-axis and its title
    axis.title.y = element_text(
      face = "bold", 
      size= 18,
      margin = margin(t=0,r = 5, b = 0, l = 0)
    ),      
    axis.text = element_text(size=14,face = "bold"),         # Axis text bold
    legend.title = element_text(size = 14, face = "bold"),  # Small bold legend title
    legend.text = element_text(size = 14, face = "bold")    # Small bold legend text
  )
age_continuous_formatted

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 18 rows containing missing values or values outside the scale range
(`geom_segment()`).

ggsave(here("figures",PROJECT_VERSION,"age_continuous_formatted.png"),age_continuous_formatted,width=9,height=6,bg = "white")

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 18 rows containing missing values or values outside the scale range
(`geom_segment()`).

younger group mean

younger_mean <- participant_age_halves |>
  filter(age_half == "younger") |>
  summarize(mean_value = mean(SubjectInfo.testAge)/30) |>
  pull(mean_value)

Order effects

First vs second instances for primary targets

first_instance_target <- looking_data_summarized |>
  group_by(SubjectInfo.subjID, Trials.targetImage) |>
  arrange(Trials.ordinal, .by_group = TRUE) |>
  filter(Trials.trialType %in% c("easy", "hard")) |>
  slice(1) |>
  ungroup() |>
  mutate(order_id = "one")
  
second_instance_target <- looking_data_summarized |>
  group_by(SubjectInfo.subjID, Trials.targetImage) |>
  arrange(Trials.ordinal, .by_group = TRUE) |>
  filter(Trials.trialType %in% c("easy", "hard")) |>
  slice(2) |>
  ungroup() |>
  mutate(order_id = "two")

ordered_conditions <- rbind(first_instance_target, second_instance_target)

order_plot <- half_violins_plot(ordered_conditions, "order_id", "corrected_target_looking", "SubjectInfo.subjID", c("one", "two"), input_xlab = "Target image order")

Warning: There was 1 warning in `summarize()`.
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 330: `order_id = "two"`, `SubjectInfo.subjID = "VVI232"`.
Caused by warning in `qt()`:
! NaNs produced

# A tibble: 2 × 8
  order_id mean_value sd_value     N     se     ci lower_ci upper_ci
  <chr>         <dbl>    <dbl> <int>  <dbl>  <dbl>    <dbl>    <dbl>
1 one          0.0634    0.141   170 0.0108 0.0213   0.0420   0.0847
2 two          0.0590    0.140   163 0.0109 0.0216   0.0374   0.0806

Warning in geom_path(aes(group = SubjectInfo.subjID), color = "black", fill =
NA, : Ignoring unknown parameters: `fill`

order_plot

Sections

violins

section_plot <- half_violins_plot(looking_data_summarized, "section", "corrected_target_looking", "SubjectInfo.subjID", c("sample1", "sample2"), input_xlab = "Data sample")

# A tibble: 2 × 8
  section mean_value sd_value     N      se     ci lower_ci upper_ci
  <chr>        <dbl>    <dbl> <int>   <dbl>  <dbl>    <dbl>    <dbl>
1 sample1     0.0500   0.0852    90 0.00898 0.0178   0.0322   0.0679
2 sample2     0.0566   0.0972    80 0.0109  0.0216   0.0350   0.0783

Warning in geom_path(aes(group = SubjectInfo.subjID), color = "black", fill =
NA, : Ignoring unknown parameters: `fill`

section_plot

`geom_path()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

timecourse

summarize_across_sections <- summarize_subj |>
  left_join(looking_data_summarized |> distinct(section, SubjectInfo.subjID)) |>
  group_by(time_normalized_corrected, section) |>
  dplyr::summarize(n=n(),
            non_na_n = sum(!is.na(mean_accuracy)),
            mean_value=mean(mean_accuracy,na.rm=TRUE),
            sd_accuracy=sd(mean_accuracy,na.rm=TRUE),
            se_accuracy=sd_accuracy/sqrt(n),
            ci = ifelse(n > 1, qt(0.975, n - 1) * sd_accuracy / sqrt(n), NA)) %>%
  filter(n > 5)

Joining with `by = join_by(SubjectInfo.subjID)`
`summarise()` has regrouped the output.

looking_times_section <- ggplot(summarize_across_sections,aes(time_normalized_corrected,mean_value,color=section))+
  xlim(-2000,4000)+
  geom_errorbar(aes(ymin=mean_value-ci,ymax=mean_value+ci),width=0, alpha=0.2)+
  #geom_point(alpha=0.2)+
    geom_smooth(method="gam")+
  geom_vline(xintercept=0,size=1.5)+
  geom_hline(yintercept=0.5,size=1.2,linetype="dashed")+
  geom_vline(xintercept=300,linetype="dotted")+
  ylim(0,1)+
  xlab("Time (normalized to target word onset) in ms")+
  ylab("Proportion Target Looking") +
  labs(title="Proportion of looking time by section and time") +
  scale_color_brewer(palette = "Set2", name="Section")

looking_times_section

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

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

Easy vs hard trial plots

Proportion of target looking for easy vs hard trials

avg_corrected_target_looking_by_condition <- summarized_data(looking_data_summarized,
                                                             "Trials.trialType", "corrected_target_looking", "SubjectInfo.subjID")
high_accuracy_subjects <- avg_corrected_target_looking_by_condition |>
  filter(lower_ci > 0 & N > 5) |>
  distinct(SubjectInfo.subjID)

overall_condition_plot <- half_violins_plot(looking_data_summarized |> group_by(SubjectInfo.subjID, Trials.targetImage) |>
  arrange(Trials.ordinal, .by_group = TRUE) |>
  #slice(2) |>
  ungroup() |> mutate(
  Trials.trialType = ifelse(Trials.trialType == "easy-distractor", "easy", Trials.trialType),
  Trials.trialType = ifelse(Trials.trialType == "hard-distractor", "hard", Trials.trialType)), "Trials.trialType", "corrected_target_looking", "SubjectInfo.subjID", c("easy", "hard"))

# A tibble: 2 × 8
  Trials.trialType mean_value sd_value     N      se     ci lower_ci upper_ci
  <chr>                 <dbl>    <dbl> <int>   <dbl>  <dbl>    <dbl>    <dbl>
1 easy                 0.0659    0.109   170 0.00837 0.0165   0.0493   0.0824
2 hard                 0.0408    0.115   170 0.00882 0.0174   0.0234   0.0582

Warning in geom_path(aes(group = SubjectInfo.subjID), color = "black", fill =
NA, : Ignoring unknown parameters: `fill`

overall_condition_plot

timecourse

summarize_trial_difficulty <- looking_time_resampled_clean %>%
  mutate(
  Trials.trialType = ifelse(Trials.trialType == "easy-distractor", "easy", Trials.trialType),
  Trials.trialType = ifelse(Trials.trialType == "hard-distractor", "hard", Trials.trialType)) |>
  filter(trial_exclusion == 0 & exclude_participant ==0 & exclude_participant_insufficient_data == 0) %>%
  filter(!is.na(accuracy_transformed)) %>%
  group_by(SubjectInfo.subjID, time_normalized_corrected, SubjectInfo.testAge) %>%
  summarized_data(., "Trials.trialID", "accuracy_transformed", c("Trials.trialType", "time_normalized_corrected")) %>%
  rename(mean_accuracy = mean_value) 

Warning: There were 929 warnings in `summarize()`.
The first warning was:
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 1: `Trials.trialID = "easy-acorn-key"`, `Trials.trialType = "easy"`,
  `time_normalized_corrected = -3800`.
Caused by warning in `qt()`:
! NaNs produced
ℹ Run `dplyr::last_dplyr_warnings()` to see the 928 remaining warnings.

summarize_across_trial_type <- summarize_trial_difficulty |>
  group_by(time_normalized_corrected, Trials.trialType) |>
  dplyr::summarize(n=n(),
            non_na_n = sum(!is.na(mean_accuracy)),
            mean_value=mean(mean_accuracy,na.rm=TRUE),
            sd_accuracy=sd(mean_accuracy,na.rm=TRUE),
            se_accuracy=sd_accuracy/sqrt(n),
            ci = ifelse(n > 1, qt(0.975, n - 1) * sd_accuracy / sqrt(n), NA)) 

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by time_normalized_corrected and
  Trials.trialType.
ℹ Output is grouped by time_normalized_corrected.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(time_normalized_corrected, Trials.trialType))` for
  per-operation grouping (`?dplyr::dplyr_by`) instead.

looking_times_trial_type <- ggplot(summarize_across_trial_type,aes(time_normalized_corrected,mean_value,color=Trials.trialType))+
  xlim(-2000,4000)+
  geom_errorbar(aes(ymin=mean_value-ci,ymax=mean_value+ci),width=0, alpha=0.2)+
  #geom_point(alpha=0.2)+
    geom_smooth(method="gam")+
  geom_vline(xintercept=0,size=1.5)+
  geom_hline(yintercept=0.5,size=1.2,linetype="dashed")+
  geom_vline(xintercept=300,linetype="dotted")+
  ylim(0,1)+
  xlab("Time (normalized to target word onset) in ms")+
  ylab("Proportion Target Looking") +
 # labs(title="Proportion of looking time across time, split at mean age-of-acquisition") +
  scale_color_manual(values=c("forestgreen", "#4292C6"), name="Trial type", labels=c( "easy", "hard"))+
   guides(color = guide_legend(reverse = TRUE))

looking_times_trial_type

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

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

Target image difficulty

# summarize average accuracy within participant (by word alone)
condition_based_looking <- looking_data_summarized  |>
  filter(!grepl("distractor", Trials.trialType)) |>
  distinct(SubjectInfo.subjID, Trials.trialID, Trials.targetImage, corrected_target_looking, Trials.trialType, AoA_Est_target)

target_looking_by_target_word <- looking_data_summarized  |>
  filter(!grepl("distractor", Trials.trialType)) |>
  group_by(SubjectInfo.subjID, Trials.targetImage) |>
  summarize(target_looking_diff = corrected_target_looking[match("easy", Trials.trialType)] - corrected_target_looking[match("hard", Trials.trialType)],
            baseline_window_looking_diff = mean_target_looking_baseline_window[match("easy", Trials.trialType)] - mean_target_looking_baseline_window[match("hard", Trials.trialType)], 
            critical_window_looking_diff = mean_target_looking_critical_window[match("easy", Trials.trialType)] - mean_target_looking_critical_window[match("hard", Trials.trialType)]) |>
  filter(!is.na(target_looking_diff))

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by SubjectInfo.subjID and Trials.targetImage.
ℹ Output is grouped by SubjectInfo.subjID.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(SubjectInfo.subjID, Trials.targetImage))` for
  per-operation grouping (`?dplyr::dplyr_by`) instead.

#clean names for individual images for plot
overall_target_looking_by_word <- target_looking_by_target_word %>%
  filter(!is.na(target_looking_diff)) %>%
  group_by(Trials.targetImage) %>%
  summarize(N=n(),
            corrected_target_looking=mean(target_looking_diff,na.rm=TRUE),
            ci=qt(0.975, N-1)*sd(target_looking_diff,na.rm=T)/sqrt(N),
            lower_ci=corrected_target_looking-ci,
            upper_ci=corrected_target_looking+ci
          ) 


word_prefs <- ggplot(overall_target_looking_by_word,aes(reorder(Trials.targetImage,corrected_target_looking),corrected_target_looking))+
  geom_hline(yintercept=0,linetype="dashed")+
  geom_errorbar(aes(ymin=lower_ci,ymax=upper_ci),width=0, alpha=0.3)+
  geom_jitter(data=target_looking_by_target_word, aes(x=Trials.targetImage, y=target_looking_diff, color=SubjectInfo.subjID), size = 2.5, alpha = 0.3, width=0.1) +
  geom_point(size = 2.4)+
  xlab("Target Word")+
  ylab("Easy trial prefered target looking")+
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=10,angle=90,vjust=0.5,hjust=1),
        axis.title.y = element_text(face="bold", size=15),
        axis.text.y  = element_text(size=10),
        strip.text.x = element_text(size = 10,face="bold")
        ) +
  scale_y_continuous(breaks = seq(-1, 1, by = 0.2)) +
  scale_size_continuous(name = "Number of participants") + 
  scale_color_viridis_d(name = "Participant IDs",option="D") +
  guides(color = "none")
  
word_prefs

Baseline image-pair preferences

saliency_effects <- looking_data_summarized |>
  # Calculating the proportion of time looking at the target word even if it isn't the target word in that particular study
 mutate(original_target_looking_baseline_window = ifelse(grepl("distractor", Trials.trialType), 1 - mean_target_looking_baseline_window,   mean_target_looking_baseline_window),
        original_target_looking_critical_window = ifelse(grepl("distractor", Trials.trialType), 1 - mean_target_looking_critical_window, mean_target_looking_critical_window)) |>
  add_age_split() |>
 summarize(mean_baseline_looking = mean(original_target_looking_baseline_window),
           mean_critical_looking = mean(original_target_looking_critical_window),
           .by = c(Trials.imagePair, age_half))
saliency_effects

    Trials.imagePair age_half mean_baseline_looking mean_critical_looking
1       bathtub-sock    older             0.3962419             0.4461523
2       bathtub-sock  younger             0.5151571             0.4193646
3     bathtub-shower    older             0.3918835             0.4166591
4     bathtub-shower  younger             0.4615050             0.4732180
5   bulldozer-orange  younger             0.4343250             0.3825479
6   bulldozer-orange    older             0.5748909             0.4443989
7       duck-chicken  younger             0.4577177             0.4573538
8       duck-chicken    older             0.4444205             0.4695415
9          truck-car  younger             0.4607519             0.4887008
10         truck-car    older             0.5455288             0.4914979
11    duck-butterfly    older             0.4835706             0.4586279
12    duck-butterfly  younger             0.4760166             0.4710987
13        cloud-tree  younger             0.4212441             0.4361528
14        cloud-tree    older             0.4978974             0.4666187
15         acorn-key  younger             0.4469638             0.4186809
16         acorn-key    older             0.4788331             0.4480183
17      cloud-zipper    older             0.3903522             0.3640776
18      cloud-zipper  younger             0.3629311             0.3796322
19    truck-airplane  younger             0.4555424             0.5407692
20    truck-airplane    older             0.3690387             0.5263529
21        snail-worm    older             0.3996053             0.4479406
22        snail-worm  younger             0.4054677             0.4260101
23         snail-cow    older             0.4532581             0.4402747
24         snail-cow  younger             0.3983266             0.3683577
25       turkey-swan    older             0.4403477             0.4877313
26       turkey-goat  younger             0.4106288             0.4254979
27       turkey-swan  younger             0.4929389             0.4957185
28       turkey-goat    older             0.4253035             0.4760688
29        cheese-mud  younger             0.5015125             0.5914350
30        cheese-mud    older             0.4545625             0.5276830
31      turtle-horse  younger             0.4527153             0.4231412
32      turtle-horse    older             0.4166157             0.3823349
33       turtle-frog    older             0.4551658             0.5246384
34       turtle-frog  younger             0.4927704             0.5100310
35   squirrel-monkey  younger             0.4322144             0.4753925
36   squirrel-monkey    older             0.3600953             0.3757038
37   balloon-dresser    older             0.5781278             0.5594834
38   balloon-dresser  younger             0.6761027             0.6824057
39     cheese-butter    older             0.5058167             0.5619523
40     cheese-butter  younger             0.6411010             0.5965848
41    potato-glasses    older             0.5362731             0.5233058
42    potato-glasses  younger             0.5746260             0.6033294
43      balloon-kite    older             0.5967721             0.5870634
44      balloon-kite  younger             0.6304199             0.6625484
45    squirrel-eagle    older             0.4370153             0.5333088
46    squirrel-eagle  younger             0.5074317             0.5251362
47    camel-elephant  younger             0.4301461             0.4033518
48    camel-elephant    older             0.3415470             0.3860543
49        potato-pot    older             0.5062483             0.5873245
50        potato-pot  younger             0.5675511             0.6441189
51     camel-ladybug    older             0.4186847             0.4559067
52     camel-ladybug  younger             0.4816635             0.5268343
53     lizard-walrus  younger             0.3960536             0.3964788
54     lizard-walrus    older             0.3716760             0.4093808
55      lizard-snake    older             0.3818615             0.3567609
56      lizard-snake  younger             0.3360660             0.3823818
57 bulldozer-tractor    older             0.4571464             0.4083064
58 bulldozer-tractor  younger             0.4087415             0.4220546
59      crow-penguin    older             0.4246100             0.4681515
60      crow-penguin  younger             0.4678887             0.5320632
61     crow-seahorse    older             0.3827343             0.4212359
62     crow-seahorse  younger             0.4822122             0.4780421
63     acorn-coconut  younger             0.5080370             0.5027286
64     acorn-coconut    older             0.5163894             0.5122532

CLIP similarity plots

Comparing proportion of target looking time to CLIP and other model cosine similarities

# CLIP
clip_data_summarized <- summarize_similarity_data(looking_data_summarized)
clip_data_summarized_low_aoa <- summarize_similarity_data(looking_data_summarized |> filter(AoA_Est_target < 4.78))
# |>filter(SubjectInfo.subjID %in% high_accuracy_subjects$SubjectInfo.subjID)))
clip_age_half_summarized <- looking_data_summarized |>
  add_age_split() |>
  summarize_similarity_data(extra_fields = c("age_half", "AoA_Est_target"))

clip_plots <- generate_multimodal_plots(clip_data_summarized, "CLIP")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

clip_plots_low_aoa <- generate_multimodal_plots(clip_data_summarized_low_aoa, "CLIP AoA < 4.78")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

input_median_age = median_age
clip_age_plots <- generate_multimodal_age_effect_plots(clip_age_half_summarized, model_type="CLIP", median_age=input_median_age)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

clip_plots

clip_age_plots

# CVCL
create_model_plots(all_similarities |> transmute(image_similarity=cvcl_image_similarity, text_similarity=cvcl_text_similarity, multimodal_similarity=cvcl_multimodal_similarity, text1, text2), name="CVCL", median_age=input_median_age)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

create_model_plots(all_similarities |> transmute(image_similarity=dino_imagenet100_vitb14_image_similarity, text1, text2), name="ImageNetVIT", median_age=input_median_age)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

create_model_plots(all_similarities |> transmute(image_similarity=dino_say_vitb14_image_similarity, text1, text2), name="SayCamVIT", median_age=input_median_age)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

create_model_plots(all_similarities |> transmute(image_similarity=dinov3.babyview_image_similarity, text1, text2), name="BV-DINO", median_age=input_median_age)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

create_model_plots(all_similarities |> transmute(image_similarity=clip_image_similarity, text_similarity=clip_text_similarity, multimodal_similarity=clip_multimodal_similarity,text1, text2), name="CLIP", median_age=input_median_age)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

combined

bv_dino_data <- looking_data_summarized |> left_join(all_similarities, by=c("Trials.targetImage"="text1", "Trials.distractorImage"="text2")) |> rename("BabyView Similarity"="dinov3.babyview_image_similarity")
group_vars = c("Trials.trialID", "Trials.targetImage", "Trials.distractorImage",  "image_similarity", "text_similarity", "multimodal_similarity", "BabyView Similarity", "section")
dino_data_summarized <- (summarized_data(
  bv_dino_data,
  "Trials.trialID", 
  "corrected_target_looking", 
  group_vars
))
dino_data_long <- dino_data_summarized |>
  rename(
    `Text Similarity` = text_similarity,
    `Image Similarity` = image_similarity,
    `Multimodal Similarity` = multimodal_similarity
  ) |>
  pivot_longer(
    cols = c("Text Similarity", "Image Similarity", "Multimodal Similarity", "BabyView Similarity"),
    names_to = "sim_type",
    values_to = "similarity"
  ) |>
  mutate(
    sim_type = factor(
      sim_type,
      levels = c("Text Similarity", "Image Similarity", "Multimodal Similarity", "BabyView Similarity")
    )
  )

all_sections_plot <- multiple_similarity_effects_plot_refactored(
  dino_data_long,
  x_var = "similarity",
  group_var = "sim_type",
  input_title = "Looking time and embedding correlations",
  label_filter = "Trials.targetImage == 'acorn'",
  facet_scales_x = list(
    `Text Similarity` = scale_x_continuous(
      breaks = seq(0.7, 0.9, by = 0.1),
      limits = c(0.62, 0.93)
    ),
    `Image Similarity` = scale_x_continuous(
      breaks = seq(0.4, 0.8, by = 0.1),
      limits = c(0.38, 0.86)
    ),
    `Multimodal Similarity` = scale_x_continuous(
      breaks = seq(0.1, 0.4, by=0.1),
      limits = c(0.05, 0.45)
    ),
        `Babyview Similarity` = scale_x_continuous(breaks = seq(0, 0.4, by = 0.2), limits = c(-0.1, 0.5))
  ),
  color_map = c("Text Similarity" = "#C76D48", "Multimodal Similarity" = "#749B69"),
  n_col = 4
)
all_sections_plot

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"clip_plot_all.svg"),all_sections_plot,width=20, height=10,bg = "white",device="pdf")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"clip_plot_all.png"),all_sections_plot,width=20, height=10,bg = "white")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

scaled axes.

dino_data_long_scaled <- dino_data_summarized |>
  transmute(
    `Text Similarity` = scale(text_similarity),
    `Image Similarity` = scale(image_similarity),
    `Multimodal Similarity` = scale(multimodal_similarity),
    `BabyView Similarity` = scale(`BabyView Similarity`),
    mean_value = mean_value,
    Trials.targetImage = Trials.targetImage,
    Trials.distractorImage = Trials.distractorImage
  ) |>
  pivot_longer(
    cols = c("Text Similarity", "Image Similarity", "Multimodal Similarity", "BabyView Similarity"),
    names_to = "sim_type",
    values_to = "similarity"
  ) |>
  mutate(
    sim_type = factor(
      sim_type,
      levels = c("Text Similarity", "Image Similarity", "Multimodal Similarity", "BabyView Similarity")
    )
  )

all_sections_plot <- multiple_similarity_effects_plot_refactored(
  dino_data_long_scaled,
  x_var = "similarity",
  group_var = "sim_type",
  input_title = "Looking time and embedding correlations",
  label_filter = "Trials.targetImage == 'acorn'",
  facet_scales_x = list(
    `Text Similarity` = scale_x_continuous(
      breaks = seq(-2, 2, by = 0.5),
      limits = c(-2, 2)
    ),
    `Image Similarity` = scale_x_continuous(
      breaks = seq(-2, 2, by = 0.5),
      limits = c(-2, 2)
    ),
    `Multimodal Similarity` = scale_x_continuous(
      breaks = seq(-2, 2, by = 0.5),
      limits = c(-2, 2)
    ),
        `Babyview Similarity` = scale_x_continuous(
      breaks = seq(-2, 2, by = 0.5),
      limits = c(-2, 2)
          )
  ),
  color_map = c("Text Similarity" = "#C76D48", "Multimodal Similarity" = "#749B69"),
  n_col = 4
)
all_sections_plot

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

section 1 plot

section_1_plot <- multiple_similarity_effects_plot_refactored(
  dino_data_long |> filter(section == "sample1" & sim_type %in% c("Image Similarity", "Text Similarity")) |> mutate(
    sim_type = factor(
      sim_type,
      levels = c("Text Similarity", "Image Similarity")
    )
  ),
  x_var = "similarity",
  group_var = "sim_type",
  input_title = "Looking time and embedding correlations",
  label_filter = "Trials.targetImage == 'acorn' | Trials.distractorImage == 'acorn'",
  facet_scales_x = list(
    `Text Similarity` = scale_x_continuous(
      breaks = seq(0.75, 0.9, by = 0.05),
      limits = c(0.7, 0.91)
    ),
    `Image Similarity` = scale_x_continuous(
      breaks = seq(0.5, 0.9, by = 0.1),
      limits = c(0.45, 0.9)
    )
  ),
  color_map = c("Text Similarity" = "#C76D48")
)
section_1_plot

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"clip_plot_section11.svg"),section_1_plot,width=15, height=10,bg = "white",device="pdf")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"clip_plot_section1.png"),section_1_plot,width=15, height=10,bg = "white")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

correlation between bv and CLIP models

clip_bv_df <- all_sims |>
  mutate(
    fitted = predict(lm(image_similarity ~ dinov3.babyview_image_similarity, data = all_sims)),
    residual = abs(image_similarity - fitted),
    is_outlier = residual > (mean(residual) + 0.2*sd(residual))
  ) 

clip_bv <- ggplot(clip_bv_df, aes(x=dinov3.babyview_image_similarity, y=image_similarity)) +
  geom_smooth(method="lm", color="#215D89", alpha=0.5, se=FALSE) +
  geom_smooth(method="lm", color="#215D89", se=TRUE, alpha=0.2) +
  geom_point(size=8, alpha=0.8, color="#215D89") +
  scale_y_continuous(breaks = seq(0.4, 0.8, by = 0.2)) +
  geom_label_repel(
    data = clip_bv_df |> filter(is_outlier),
    aes(label=Trials.imagePair),
    color="black",
    fill="white",
    alpha=0.8,
    segment.alpha=0.7,
    nudge_y=0.02,
    force=10,
    force_pull=0.1,
    segment.size=1.2,
    size=6
  ) +
  ylab("CLIP image similarity") +
  xlab("BabyView model similarity") +
  theme_classic() +
  theme(
    text = element_text(size=14, face="bold"),
    axis.title.x = element_text(face="bold", size=20, margin=margin(t=15, r=0, b=0, l=0)),
    axis.title.y = element_text(face="bold", size=20, margin=margin(t=0, r=10, b=0, l=0)),
    axis.text = element_text(size=16, face="bold"),
    legend.key = element_blank(),
    legend.title = element_text(size=22, face="bold"),
    legend.text = element_text(size=22, face="bold"),
    legend.position = "bottom",
    panel.spacing = unit(0.5, "cm")
  )
clip_bv

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures", PROJECT_VERSION,"clip_bv.svg"),clip_bv,width=6.67, height=6.67,units="in", bg = "white",device="pdf")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures", PROJECT_VERSION,"clip_bv.png"),clip_bv,width=6.67, height=6.67,units="in", bg = "white")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Some exploratory age effects with three age bins

clip_age_bins_summarized <- looking_data_summarized |>
  mutate(age_bin = case_when(
    SubjectInfo.testAge < 19*30 ~ "14-18",
    SubjectInfo.testAge < 22*30 ~ "18-21",
    TRUE ~ "21-24")) |>
  summarize_similarity_data(extra_fields = c("age_bin", "AoA_Est_target"))

similarity_age_plot(clip_age_bins_summarized, "image_similarity", age_grouping="age_bin", model_type="CLIP", median_age=median_age)

`geom_smooth()` using formula = 'y ~ x'

ggplot(clip_age_half_summarized |> left_join(all_similarities |> dplyr::select(text1, text2, dinov3.babyview_image_similarity), by=c("Trials.targetImage"="text1", "Trials.distractorImage"="text2")), aes(x = dinov3.babyview_image_similarity, y = mean_value, 
                   color = age_half)) +
    geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.7) +
    geom_point(size = 3, alpha = 0.7) +
    geom_smooth(method = "lm", se = TRUE, alpha = 0.3) +
    
    # Improved styling
    scale_color_brewer(palette = "Set1", name="Age half") +
    theme_minimal() +
    theme(
      legend.position = "bottom",
      axis.title = element_text(size = 11),
      legend.title = element_text(size = 10),
      plot.title = element_text(size = 12)
    ) +
    
    # Clear labels
    labs(
      title = "BabyView DINO correlation across age",
      x = paste("Target-Distractor", "Image Similarity"),
      y = "Baseline-corrected Proportion\nTarget Looking",
      caption = paste0("Median age: ", median_age, " months")
    ) 

`geom_smooth()` using formula = 'y ~ x'

AoA effects

input_median_age = median_age
facet_plot <- generate_aoa_facet_plots(
  data = clip_age_half_summarized, 
  model_type = "CLIP",
  median_age = input_median_age,  
  age_grouping = "age_half",
  include_multimodal = FALSE
)

`geom_smooth()` using formula = 'y ~ x'

# For separate plots by AoA group:
#low_aoa_plots <- generate_multimodal_age_effect_plots(
#  data = clip_age_half_summarized %>% filter(AoA_Est_target <= 4.44),
#  model_type = "CLIP",
#  median_age = input_median_age,  
#  age_grouping = "age_half",
#  include_multimodal = FALSE
#)
facet_plot

`geom_smooth()` using formula = 'y ~ x'

OpenCLIP training

looking_data_w_openclip <- looking_data_summarized |>
  dplyr::select(-text_similarity, -multimodal_similarity, -image_similarity) |>
  left_join(openclip_similarities, by = c("Trials.distractorImage"="text2", "Trials.targetImage"="text1"))

Warning in left_join(dplyr::select(looking_data_summarized, -text_similarity, : Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 51 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

openclip_data_summarized <- summarize_similarity_data(looking_data_w_openclip, extra_fields=c("epoch")) |> filter(!is.na(image_similarity))
openclip_age_half_summarized <- looking_data_w_openclip |>
  add_age_split() |>
  summarize_similarity_data(extra_fields = c("age_half", "epoch")) |> filter(!is.na(image_similarity))

# Function to calculate Pearson's correlation per epoch
calculate_correlations <- function(data, x_var, y_var, group_var = c("epoch"), conf_level = 0.95, method="pearson") {
  data |>
    group_by(across(all_of(group_var))) |>
    summarize(
      {
        cor_test <- cor.test(.data[[x_var]], .data[[y_var]], method = method, conf.level = conf_level)
        tibble(
          pearson_cor = cor_test$estimate,
          p_value = cor_test$p.value,
          ci_lower = cor_test$conf.int[1],
          ci_upper = cor_test$conf.int[2]
        )
      },
      .groups = "drop"
    )
}

image_correlation_results_age_split <- calculate_correlations(openclip_age_half_summarized, "image_similarity", "mean_value", group_var=c("epoch", "age_half"))

text_correlation_results <- calculate_correlations(openclip_data_summarized, "text_similarity", "mean_value")
image_correlation_results <- calculate_correlations(openclip_data_summarized, "image_similarity", "mean_value")
multimodal_correlation_results <- calculate_correlations(openclip_data_summarized, "multimodal_similarity", "mean_value")

ggplot(image_correlation_results, aes(x = log(epoch), y = pearson_cor)) +
  geom_point(aes(color = p_value < 0.05), size = 3) + 
  geom_smooth(span = 2) +
  labs(title = "Image similarity correlation across Open-CLIP training",
       x = "log(Epoch)",
       y = "Pearson Correlation") +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "gray")) +  # Set color for significance
  theme_minimal() +
  theme(legend.position = "none")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(text_correlation_results, aes(x = log(epoch), y = pearson_cor)) +
  geom_point(aes(color = p_value < 0.05), size = 3) + 
  geom_smooth(span = 2) +
  labs(title = "Word similarity correlation across Open-CLIP training",
       x = "log (Epoch)",
       y = "Pearson Correlation") +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "gray")) +  # Set color for significance
  theme_minimal() +
  theme(legend.position = "none")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(multimodal_correlation_results, aes(x = log(epoch), y = pearson_cor)) +
  geom_point(aes(color = p_value < 0.05), size = 3) + 
  geom_smooth(span = 2) +
  labs(title = "Multimodal similarity correlation across Open-CLIP training",
       x = "log (Epoch)",
       y = "Pearson Correlation") +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "gray")) +  # Set color for significance
  theme_minimal() +
  theme(legend.position = "none")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

image_openclip_age <- epoch_age_half_plot(calculate_correlations(openclip_age_half_summarized , "image_similarity", "mean_value", group_var=c("epoch", "age_half")), "image_similarity")
text_openclip_age <- epoch_age_half_plot(calculate_correlations(openclip_age_half_summarized, "text_similarity", "mean_value", group_var=c("epoch", "age_half")), "text_similarity")
image_openclip_age

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

text_openclip_age

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"openclip_image_embeddings_age.png"),image_openclip_age,width=9,height=6,bg = "white")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# --- set to FALSE to drop multimodal ---
include_multimodal <- TRUE

sim_colors <- c(
  "Image"      = "#4E79A7",  
  "Text"       = "#C0724A",  
  "Multimodal" = "#749B69"   
)

combined <- bind_rows(
  image_correlation_results %>% mutate(similarity = "Image"),
  text_correlation_results  %>% mutate(similarity = "Text"),
  if (include_multimodal)
    multimodal_correlation_results %>% mutate(similarity = "Multimodal")
) %>%
  mutate(similarity = factor(similarity, levels = names(sim_colors)))

ggplot(combined, aes(log(epoch), pearson_cor, color = similarity)) +
  geom_point(size = 3, alpha = 0.9) +
  geom_smooth(span = 2, se = TRUE) +
  scale_color_manual(values = sim_colors, drop = TRUE) +
  labs(x = "log(OpenCLIP epoch)",
       y = "Pearson's correlation (r)",
       color = NULL) +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(0.02, 0.98),
        legend.justification = c(0, 1),
        legend.background = element_rect(fill = "white", color = "grey80"))

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"openclip_all_measures.png"),width=9,height=6,bg = "white")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Examining linear correlation between similarity types

cor(data.frame(trial_metadata$multimodal_similarity, trial_metadata$image_similarity, trial_metadata$text_similarity, trial_metadata$cor))

                                     trial_metadata.multimodal_similarity
trial_metadata.multimodal_similarity                            1.0000000
trial_metadata.image_similarity                                 0.4565731
trial_metadata.text_similarity                                  0.6124357
trial_metadata.cor                                                     NA
                                     trial_metadata.image_similarity
trial_metadata.multimodal_similarity                       0.4565731
trial_metadata.image_similarity                            1.0000000
trial_metadata.text_similarity                             0.4915159
trial_metadata.cor                                                NA
                                     trial_metadata.text_similarity
trial_metadata.multimodal_similarity                      0.6124357
trial_metadata.image_similarity                           0.4915159
trial_metadata.text_similarity                            1.0000000
trial_metadata.cor                                               NA
                                     trial_metadata.cor
trial_metadata.multimodal_similarity                 NA
trial_metadata.image_similarity                      NA
trial_metadata.text_similarity                       NA
trial_metadata.cor                                    1

cor.test(trial_metadata$image_similarity, trial_metadata$text_similarity)

    Pearson's product-moment correlation

data:  trial_metadata$image_similarity and trial_metadata$text_similarity
t = 5.472, df = 94, p-value = 3.673e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3228437 0.6299275
sample estimates:
      cor 
0.4915159 

cor.test(all_sims$dino_imagenet100_vitb14_image_similarity, all_sims$image_similarity, method="spearman")

Warning in cor.test.default(all_sims$dino_imagenet100_vitb14_image_similarity,
: cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  all_sims$dino_imagenet100_vitb14_image_similarity and all_sims$image_similarity
S = 59783, p-value = 1.699e-10
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.5945241 

library(GGally)

similarity_cols <- all_sims |>
  dplyr::select(matches("similarity$")) |>
  dplyr::select(-matches("^layer\\d+_"))

ggpairs(
  similarity_cols,
  upper = list(continuous = wrap("cor", method = "spearman", size = 3)),
  lower = list(continuous = wrap("points", alpha = 0.3, size = 0.5)),
  diag = list(continuous = "densityDiag")
)

Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties
Warning in cor.test.default(x, y, method = method): cannot compute exact
p-value with ties

AoA plots

aoa_jittered <- target_looking_item_level %>%
  rowwise() %>%
  mutate(aoa_x_jittered = AoA_Est_target + runif(n(), -0.01, 0.01))

aoa_across_time <- ggplot(aoa_jittered, aes(x=aoa_x_jittered, y=mean_value))+
  geom_hline(yintercept=0,linetype="dashed")+
  geom_point(size=8, alpha=0.5, color="#CF784B")+
  #geom_linerange(aes(ymin = mean_value - ci, ymax = mean_value + ci), alpha = 0.1) + 
  geom_smooth(method="lm", color="#CF784B", alpha=0.5)+
  geom_label_repel(size=7,aes(label=Trials.targetImage),force_pull=10,nudge_y=-0.03) +
  xlab("Estimated target word age-of-acquisition (AoA)")+
  ylab("Baseline-corrected\nprop. target looking")+
  #ggtitle("Prop. of target looking by estimated AoA of target words")+
  theme_classic()+
 theme(
    text = element_text(size=10,face = "bold"),              # All text bold
    # Increase space between x-axis and its title
    axis.title.x = element_text(
      face = "bold", 
      size = 18,
      margin = margin(t = 10, r = 0, b = 0, l = 0)
    ),
    
    # Increase space between y-axis and its title
    axis.title.y = element_text(
      face = "bold", 
      size= 18,
      margin = margin(t = 10, r = 5, b = 0, l = 0)
    ),      
    axis.text = element_text(size=16,face = "bold"),         # Axis text bold
    legend.title = element_text(size = 12, face = "bold"),  # Small bold legend title
    legend.text = element_text(size = 12, face = "bold")    # Small bold legend text
  )+
 scale_x_continuous(breaks = seq(3, 6, by = 0.5)) +
   coord_cartesian(xlim = c(3.2,6.5), ylim = c(-0.12, 0.18))
aoa_across_time

`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"aoa_plot_new.svg"),aoa_across_time,width=7,height=7,bg = "white", device="pdf")

`geom_smooth()` using formula = 'y ~ x'

  #ggpubr::stat_cor(method = "spearman", label.y = max(target_looking_item_level$mean_value) + 0.05)


age_half_plot(target_looking_item_age_level, "AoA_Est_target", x_label="AoA of target word", title="Target looking by estimated AoA of target word across age half")

`geom_smooth()` using formula = 'y ~ x'

ggplot(target_looking_trial_level, aes(x=AoA_Est_distractor, y=mean_value)) +
  geom_point()+
  geom_label_repel(aes(label=Trials.imagePair)) +
  geom_smooth(method="gam")+
  geom_hline(yintercept=0,linetype="dashed")+
  xlab("Estimated AoA of distractor words")+
  ylab("Baseline-corrected prop. of target looking")+
  ggtitle("Prop. of target looking by estimated AoA of distractor words")+
  theme_minimal()+
  theme(axis.title.x = element_text(face="bold", size=15, vjust=-1),
        axis.text.x  = element_text(size=10,angle=0,vjust=0.5),
        axis.title.y = element_text(face="bold", size=15),
        axis.text.y  = element_text(size=10),
        strip.text.x = element_text(size = 10,face="bold"),
        
        ) +
  ggpubr::stat_cor(method = "spearman")

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

hist(target_looking_item_level$AoA_Est_target)

aoa_jittered <- target_looking_item_level %>%
  rowwise() %>%
  mutate(aoa_x_jittered = AoA_Est_target + runif(n(), -0.01, 0.01))

aoa_across_time <- ggplot(aoa_jittered, aes(x=aoa_x_jittered, y=mean_value))+
  geom_hline(yintercept=0,linetype="dashed")+
  geom_point(size=8, color="#E95454", alpha=0.6)+ #, aes(alpha=section)
  scale_alpha_manual(values = c(sample1 = 0.8, sample2 = 0.3)) +
  #geom_linerange(aes(ymin = mean_value - ci, ymax = mean_value + ci), alpha = 0.1) + 
  geom_smooth(method="lm", color="#E95454", alpha=0.5, se=FALSE) +
geom_smooth(method="lm", color="#E95454", se=TRUE, alpha=0.2) +
  geom_label_repel(data=aoa_jittered |> filter(AoA_Est_target > 5.7 | (AoA_Est_target < 5 & AoA_Est_target > 4.5) | Trials.targetImage %in% c("airplane", "mud", "worm", "car", "butterfly", "duck", "acorn", "coconut", "key", "truck") | (mean_value > 0.1 & AoA_Est_target > 4)), size=7,aes(label=Trials.targetImage),force_pull=15,nudge_y=-0.03) +
  xlab("Estimated target word age-of-acquisition (AoA)")+
  ylab("Baseline-corrected\nproportion target looking")+
  #ggtitle("Prop. of target looking by estimated AoA of target words")+
  theme_classic()+
  guides(shape = "none", alpha="none") +
 theme(
    text = element_text(size=10,face = "bold"),              # All text bold
    # Increase space between x-axis and its title
    axis.title.x = element_text(
      face = "bold", 
      size = 24,
      margin = margin(t = 10, r = 0, b = 0, l = 0)
    ),
    
    # Increase space between y-axis and its title
    axis.title.y = element_text(
      face = "bold", 
      size= 24,
      margin = margin(t = 10, r = 5, b = 0, l = 0)
    ),      
    axis.text = element_text(size=20),         # Axis text bold
    #legend.title = element_text(size = 12, face = "bold"),  # Small bold legend title
    #legend.text = element_text(size = 12, face = "bold")    # Small bold legend text
  )+
 scale_x_continuous(breaks = seq(3, 6, by = 0.5)) +
   coord_cartesian(xlim = c(3.2,6.5), ylim = c(-0.16, 0.18))
aoa_across_time

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"aoa_plot.svg"),aoa_across_time,height=6.936,width=9.437,units="in",bg = "white", device="pdf")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures",PROJECT_VERSION,"aoa_plot.png"),aoa_across_time,height=6.936,width=9.437,units="in",bg = "white")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

aoa_data <- looking_data_summarized |> add_aoa_split() 
aoa_data$aoa_half <- factor(aoa_data$aoa_half, levels = c("lower", "higher"))

aoa_plot <- half_violins_plot(aoa_data, "aoa_half", "corrected_target_looking", "SubjectInfo.subjID", c("lower", "higher"), input_xlab = "Age-of-acquisition of target word", input_caption=paste0("AoA split at median=", median_aoa))

# A tibble: 2 × 8
  aoa_half mean_value sd_value     N      se     ci lower_ci upper_ci
  <fct>         <dbl>    <dbl> <int>   <dbl>  <dbl>    <dbl>    <dbl>
1 lower        0.0746    0.121   170 0.00928 0.0183  0.0563    0.0929
2 higher       0.0284    0.124   170 0.00952 0.0188  0.00957   0.0471

Warning in geom_path(aes(group = SubjectInfo.subjID), color = "black", fill =
NA, : Ignoring unknown parameters: `fill`

aoa_plot

summarize_trial <- looking_time_resampled_clean %>%
  filter(trial_exclusion == 0 & exclude_participant ==0 & exclude_participant_insufficient_data == 0) %>%
  filter(!is.na(accuracy_transformed)) %>%
  group_by(SubjectInfo.subjID, time_normalized_corrected, SubjectInfo.testAge) %>%
  summarized_data(., "Trials.trialID", "accuracy_transformed", c("time_normalized_corrected")) %>%
  rename(mean_accuracy = mean_value)

Warning: There were 929 warnings in `summarize()`.
The first warning was:
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 1: `Trials.trialID = "easy-acorn-key"`, `time_normalized_corrected =
  -3800`.
Caused by warning in `qt()`:
! NaNs produced
ℹ Run `dplyr::last_dplyr_warnings()` to see the 928 remaining warnings.

 summarize_across_aoa <- summarize_trial |> left_join(trial_metadata) |> add_aoa_split() |>
  group_by(time_normalized_corrected, aoa_half) |>
  dplyr::summarize(n=n(),
            non_na_n = sum(!is.na(mean_accuracy)),
            mean_value=mean(mean_accuracy,na.rm=TRUE),
            sd_accuracy=sd(mean_accuracy,na.rm=TRUE),
            se_accuracy=sd_accuracy/sqrt(n),
            ci = ifelse(n > 1, qt(0.975, n - 1) * sd_accuracy / sqrt(n), NA)) 

Joining with `by = join_by(Trials.trialID)`

Warning in left_join(summarize_trial, trial_metadata): Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 19 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by time_normalized_corrected and aoa_half.
ℹ Output is grouped by time_normalized_corrected.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(time_normalized_corrected, aoa_half))` for
  per-operation grouping (`?dplyr::dplyr_by`) instead.

looking_times_aoa <- ggplot(summarize_across_aoa,aes(time_normalized_corrected,mean_value,color=aoa_half))+
  xlim(-2000,4000)+
  geom_errorbar(aes(ymin=mean_value-ci,ymax=mean_value+ci),width=0, alpha=0.2)+
  #geom_point(alpha=0.2)+
    geom_smooth(method="gam")+
  geom_vline(xintercept=0,size=1.5)+
  geom_hline(yintercept=0.5,size=1.2,linetype="dashed")+
  geom_vline(xintercept=300,linetype="dotted")+
  ylim(0,1)+
  xlab("Time (normalized to target word onset) in ms")+
  ylab("Proportion Target Looking") +
 # labs(title="Proportion of looking time across time, split at mean age-of-acquisition") +
  scale_color_manual(values=c("#6BAED6", "#4292C6"), name="Target word \nage-of-acquisition (AoA)", labels=c( "AoA >= 4.78", "AoA < 4.78"))+
   guides(color = guide_legend(reverse = TRUE))

looking_times_aoa

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

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

AoA correlations

aoa_wb_trials <- trial_metadata |> filter(!is.na(AoA_target_WB))
cor.test(aoa_wb_trials$AoA_target_WB, aoa_wb_trials$AoA_Est_target, method="spearman")

Warning in cor.test.default(aoa_wb_trials$AoA_target_WB,
aoa_wb_trials$AoA_Est_target, : cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  aoa_wb_trials$AoA_target_WB and aoa_wb_trials$AoA_Est_target
S = 4838.3, p-value = 0.1363
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.2607651 

cor.test(trial_metadata$AoA_Est_target, trial_metadata$image_similarity, method="pearson")

    Pearson's product-moment correlation

data:  trial_metadata$AoA_Est_target and trial_metadata$image_similarity
t = 1.4599, df = 94, p-value = 0.1477
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.05317716  0.33925485
sample estimates:
      cor 
0.1488962 

subjects_data <- read.csv(file.path(PROCESSED_DATA_PATH, "level-participant_added-trials_data.csv")) |> filter(exclude_participant_insufficient_data == 0)
trial_data <- looking_data_summarized |> distinct(SubjectInfo.subjID, SubjectInfo.testAge)
subjects_data <- subjects_data |> 
  left_join(trial_data)

Joining with `by = join_by(SubjectInfo.subjID)`

hist(subjects_data$SubjectInfo.testAge/30, breaks=10)

Saliency correlations

target_looking_trial_age_level <- summarized_data(target_looking_trial_subject_level |> rename(mean_target_looking = mean_value) |> add_age_split(), "Trials.trialID", "mean_target_looking", c("AoA_Est_target", "age_half"))
trialid_age_level_summary <- target_looking_trial_age_level |> left_join(trial_metadata)

Joining with `by = join_by(Trials.trialID, AoA_Est_target)`

Warning in left_join(target_looking_trial_age_level, trial_metadata): Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 19 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

trialid_level_summary <- target_looking_trial_level |> left_join(trial_metadata)

Joining with `by = join_by(Trials.trialID, AoA_Est_distractor, AoA_Est_target,
Trials.targetImage, Trials.distractorImage, Trials.imagePair)`

age_half_plot(trialid_age_level_summary, "MeanTargetSaliency", x_label="Mean target saliency")

`geom_smooth()` using formula = 'y ~ x'

age_half_plot(trialid_age_level_summary, "MeanSaliencyDiff", x_label="Mean saliency difference", title="Target looking by mean saliency difference across age and trial types")

`geom_smooth()` using formula = 'y ~ x'

target_looking_image_pair <- summarized_data(looking_data_summarized, "Trials.imagePair", "corrected_target_looking", c("SubjectInfo.subjID", "SubjectInfo.testAge"))

Warning: There were 327 warnings in `summarize()`.
The first warning was:
ℹ In argument: `ci = qt(0.975, N - 1) * sd_value/sqrt(N)`.
ℹ In group 180: `Trials.imagePair = "balloon-dresser"`, `SubjectInfo.subjID =
  "VVI141"`, `SubjectInfo.testAge = 750`.
Caused by warning in `qt()`:
! NaNs produced
ℹ Run `dplyr::last_dplyr_warnings()` to see the 326 remaining warnings.

target_looking_imagepair_level <- summarized_data(target_looking_image_pair |> rename("mean_looking"="mean_value") |> add_age_split(), "Trials.imagePair", "mean_looking", c("age_half")) |> left_join(trial_metadata |> distinct(Trials.imagePair, .keep_all=TRUE)) |> mutate(MeanSaliencyDiff = abs(MeanSaliencyDiff))

Joining with `by = join_by(Trials.imagePair)`

age_half_plot(target_looking_imagepair_level, "MeanSaliencyDiff", x_label="Mean saliency difference", title="Target looking by mean saliency difference across age half and image pairs")

`geom_smooth()` using formula = 'y ~ x'

Saliency vs baseline

baseline_looking_trial_age_level <- summarized_data(looking_data_summarized |> add_age_split(), "Trials.trialID", "mean_target_looking_baseline_window", c("age_half")) |> left_join(trial_metadata)

Joining with `by = join_by(Trials.trialID)`

Warning in left_join(summarized_data(add_age_split(looking_data_summarized), : Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 19 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

baseline_looking_trial_level <- summarized_data(looking_data_summarized |> rename("mean_looking"="mean_target_looking_critical_window"), "Trials.trialID", "mean_looking", c("Trials.trialID")) |> left_join(trial_metadata)

Joining with `by = join_by(Trials.trialID)`

age_half_plot(baseline_looking_trial_age_level, "MeanSaliencyDiff", x_label="Mean saliency difference", y_label = "Target looking in baseline window")

`geom_smooth()` using formula = 'y ~ x'

ggplot(baseline_looking_trial_level, aes(x=MeanSaliencyDiff, y=mean_value)) +
  geom_smooth(method="lm") +
  geom_point() +
  ylab("Baseline looking") +
  xlab("Mean target-distractor saliency difference") +
  ggpubr::stat_cor()

`geom_smooth()` using formula = 'y ~ x'

# Saliency plots
# Prepare the data by combining both datasets with a grouping variable
baseline_looking_trial_level <- summarized_data(looking_data_summarized |> rename("mean_looking"="mean_target_looking_baseline_window"), "Trials.trialID", "mean_looking", c("Trials.trialID", "section")) |> left_join(trial_metadata) 

Joining with `by = join_by(Trials.trialID, section)`

critical_looking_trial_level <- summarized_data(looking_data_summarized |> rename("mean_looking"="mean_target_looking_critical_window"), "Trials.trialID", "mean_looking", c("Trials.trialID", "section")) |> left_join(trial_metadata)

Joining with `by = join_by(Trials.trialID, section)`

trialid_level_summary$plot_type <- "Baseline-corrected"
baseline_looking_trial_level$plot_type <- "Baseline window"
critical_looking_trial_level$plot_type <- "Critical window"

combined_data <- bind_rows(trialid_level_summary, baseline_looking_trial_level, critical_looking_trial_level)

# Create labels for trials where target is acorn and distractor is coconut
combined_data$label_text <- ifelse(
  combined_data$Trials.targetImage == "acorn",
  paste("Target:", combined_data$Trials.targetImage, "\nDistractor:", combined_data$Trials.distractorImage),
  ""
)

label_data <- subset(combined_data, label_text != "")

# Per-facet dashed reference lines
hline_data <- data.frame(
  plot_type = c("Baseline window", "Critical window", "Baseline-corrected"),
  yint      = c(0.5, 0.5, 0.0)
)

combined_saliency_plot <- ggplot(combined_data, aes(x = MeanSaliencyDiff, y = mean_value)) +
  geom_smooth(method = "lm", color = "#182B49", fill = "#182B49", alpha = 0.15) +
  geom_hline(data = hline_data, aes(yintercept = yint),
             linetype = "dashed", color = "gray40", linewidth = 1) +
  geom_point(size = 8, color = "#182B49", alpha = 0.8) +
  #scale_alpha_manual(values = c(sample1 = 0.8, sample2 = 0.3)) +
  geom_label_repel(data = label_data,
                   aes(label = label_text),
                   color = "black",
                   fill = "white",
                   alpha = 0.8,
                   segment.alpha = 0.7,
                   nudge_y = 0.02,
                   force = 10,
                   force_pull = 0.1,
                   segment.size = 1.2,
                   max.overlaps = Inf,
                   size = 8) +
  facet_wrap(~ plot_type, scales = "free_y", ncol = 3) +
  xlab("Target-distractor GBVS mean saliency difference") +
  ylab("Proportion of target looking") +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),
    axis.title.x = element_text(face = "bold", size = 29,
                                margin = margin(t = 15, r = 0, b = 0, l = 0)),
    legend.key = element_blank(),
    axis.title.y = element_text(face = "bold", size = 29,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),
    axis.text = element_text(size = 24, face = "bold"),
    legend.title = element_text(size = 22, face = "bold"),
    legend.text = element_text(size = 22, face = "bold"),
    legend.position = "bottom",
    strip.text = element_text(size = 28, face = "bold"),
    strip.background = element_rect(fill = "gray90", color = NA),
    strip.text.x = element_text(margin = margin(t = 8, b = 8)),
    panel.spacing = unit(0.5, "cm")
  )

combined_saliency_plot

`geom_smooth()` using formula = 'y ~ x'

# Save the combined plot
ggsave(here("figures", PROJECT_VERSION, "saliency_plot_combined_full.svg"), 
       combined_saliency_plot, width = 15, height = 12, bg = "white", device = "pdf")

`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures", PROJECT_VERSION, "saliency_plot_combined_full.png"), 
       combined_saliency_plot, width = 15, height = 12, bg = "white")

`geom_smooth()` using formula = 'y ~ x'

combined_image_similarity_plot <- ggplot(combined_data, aes(x = image_similarity, y = mean_value)) +
  geom_smooth(method = "lm", color = "#182B49", fill = "#182B49", alpha = 0.15) +
  geom_hline(data = hline_data, aes(yintercept = yint),
             linetype = "dashed", color = "gray40", linewidth = 1) +
  geom_point(size = 8, color = "#182B49", alpha = 0.8) +
  #scale_alpha_manual(values = c(sample1 = 0.8, sample2 = 0.3)) +
  geom_label_repel(data = label_data,
                   aes(label = label_text),
                   color = "black",
                   fill = "white",
                   alpha = 0.8,
                   segment.alpha = 0.7,
                   nudge_y = 0.02,
                   force = 10,
                   force_pull = 0.1,
                   segment.size = 1.2,
                   max.overlaps = Inf,
                   size = 8) +
  facet_wrap(~ plot_type, scales = "free_y", ncol = 3) +
  xlab("Target-distractor image similarity difference") +
  ylab("Proportion of target looking") +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),
    axis.title.x = element_text(face = "bold", size = 29,
                                margin = margin(t = 15, r = 0, b = 0, l = 0)),
    legend.key = element_blank(),
    axis.title.y = element_text(face = "bold", size = 29,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),
    axis.text = element_text(size = 24, face = "bold"),
    legend.title = element_text(size = 22, face = "bold"),
    legend.text = element_text(size = 22, face = "bold"),
    legend.position = "bottom",
    strip.text = element_text(size = 28, face = "bold"),
    strip.background = element_rect(fill = "gray90", color = NA),
    strip.text.x = element_text(margin = margin(t = 8, b = 8)),
    panel.spacing = unit(0.5, "cm")
  )

combined_image_similarity_plot

`geom_smooth()` using formula = 'y ~ x'

# Save the combined plot
ggsave(here("figures", PROJECT_VERSION, "saliency_plot_combined_full.svg"), 
       combined_saliency_plot, width = 15, height = 12, bg = "white", device = "pdf")

`geom_smooth()` using formula = 'y ~ x'

ggsave(here("figures", PROJECT_VERSION, "image_sim_plot_combined_full.png"), 
       combined_image_similarity_plot, width = 15, height = 12, bg = "white")

`geom_smooth()` using formula = 'y ~ x'

Embedding type comparisons

cor.test(trial_metadata$text_similarity, trial_metadata$image_similarity)

    Pearson's product-moment correlation

data:  trial_metadata$text_similarity and trial_metadata$image_similarity
t = 5.472, df = 94, p-value = 3.673e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3228437 0.6299275
sample estimates:
      cor 
0.4915159 

cor.test(trial_metadata$multimodal_similarity, trial_metadata$image_similarity)

    Pearson's product-moment correlation

data:  trial_metadata$multimodal_similarity and trial_metadata$image_similarity
t = 4.9755, df = 94, p-value = 2.932e-06
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2818904 0.6019579
sample estimates:
      cor 
0.4565731 

cor.test(trial_metadata$AoA_Est_target, trial_metadata$text_similarity)

    Pearson's product-moment correlation

data:  trial_metadata$AoA_Est_target and trial_metadata$text_similarity
t = -0.098581, df = 94, p-value = 0.9217
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2102248  0.1907073
sample estimates:
        cor 
-0.01016738 

cor.test(trial_metadata$AoA_Est_target, trial_metadata$image_similarity)

    Pearson's product-moment correlation

data:  trial_metadata$AoA_Est_target and trial_metadata$image_similarity
t = 1.4599, df = 94, p-value = 0.1477
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.05317716  0.33925485
sample estimates:
      cor 
0.1488962 

cor.test(trial_metadata$MeanSaliencyDiff, trial_metadata$AoA_Est_target)

    Pearson's product-moment correlation

data:  trial_metadata$MeanSaliencyDiff and trial_metadata$AoA_Est_target
t = 0.15857, df = 94, p-value = 0.8743
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1847383  0.2161306
sample estimates:
      cor 
0.0163533 

layer-wise

layerwise_long <- layerwise |>
  pivot_longer(
    cols = starts_with("layer"),
    names_to = "layer",
    values_to = "image_similarity"
  ) |>
  mutate(layer = as.integer(str_extract(layer, "\\d+")))


looking_data_w_layerwise <- looking_data_summarized |>
  dplyr::select(-image_similarity) |>
  left_join(layerwise_long, by = c("Trials.distractorImage" = "text1", "Trials.targetImage" = "text2"), relationship="many-to-many")

layerwise_data_summarized <- summarize_similarity_data(looking_data_w_layerwise, extra_fields = c("layer"))

image_correlation_results_layerwise <- calculate_correlations(
  layerwise_data_summarized, "image_similarity", "mean_value", group_var = "layer", method="pearson"
)

ggplot(image_correlation_results_layerwise, aes(x = layer, y = pearson_cor)) +
  geom_point(aes(color = p_value < 0.05), size = 3) +
  geom_line() +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.15) +
  scale_x_continuous(breaks = 0:12) +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "gray")) +
  geom_hline(yintercept=0, linetype="dotted") +
  labs(
    title = "Image similarity correlation across CLIP ViT layers",
    x = "Layer",
    y = "Pearson Correlation"
  ) +
  theme_classic(base_size=16) +
  theme(legend.position = "none")

ggsave(here("figures", PROJECT_VERSION, "clip_layerwise.png"), width=12, height=9, bg="white")

layerwise_models <- layerwise_long |>
  group_by(layer) |>
  group_split() |>
  map(function(layer_data) {
    layer_num <- unique(layer_data$layer)
    
  d <- looking_data_summarized |> mutate(age_in_months = SubjectInfo.testAge) |>
  dplyr::select(-any_of("image_similarity")) |>
  left_join(layer_data, by = c("Trials.distractorImage" = "text2", "Trials.targetImage" = "text1"),
            relationship = "many-to-many")
    
    m <- lmer(scale(corrected_target_looking) ~ scale(image_similarity) * scale(age_in_months) 
              + scale(AoA_Est_target)
              + scale(MeanSaliencyDiff)
              + (scale(image_similarity) | SubjectInfo.subjID)
              + (1 | Trials.targetImage)
              + (1 | Trials.imagePair),
              data = d)
    
    coefs <- summary(m)$coefficients
    
    tibble(
      layer = layer_num,
      image_estimate = coefs["scale(image_similarity)", "Estimate"],
      image_se = coefs["scale(image_similarity)", "Std. Error"],
      image_p = coefs["scale(image_similarity)", "Pr(>|t|)"],
      interaction_estimate = coefs["scale(image_similarity):scale(age_in_months)", "Estimate"],
      interaction_se = coefs["scale(image_similarity):scale(age_in_months)", "Std. Error"],
      interaction_p = coefs["scale(image_similarity):scale(age_in_months)", "Pr(>|t|)"]
    )
  }) |>
  bind_rows()

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00472105 (tol = 0.002, component 1)
  See ?lme4::convergence and ?lme4::troubleshooting.

# Plot: image_similarity main effect
p1 <- ggplot(layerwise_models, aes(x = layer, y = image_estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_ribbon(aes(ymin = image_estimate - 1.96 * image_se,
                  ymax = image_estimate + 1.96 * image_se), alpha = 0.15) +
  geom_line() +
  geom_point(aes(color = image_p < 0.05), size = 3) +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "gray")) +
  scale_x_continuous(breaks = 0:12) +
  labs(title = "Image similarity effect across ViT layers",
       x = "Layer", y = "β (scaled)") +
  theme_minimal() + theme(legend.position = "none")

# Plot: interaction with age
p2 <- ggplot(layerwise_models, aes(x = layer, y = interaction_estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_ribbon(aes(ymin = interaction_estimate - 1.96 * interaction_se,
                  ymax = interaction_estimate + 1.96 * interaction_se), alpha = 0.15) +
  geom_line() +
  geom_point(aes(color = interaction_p < 0.05), size = 3) +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "gray")) +
  scale_x_continuous(breaks = 0:12) +
  labs(title = "Image similarity × age interaction across ViT layers",
       x = "Layer", y = "β (scaled)") +
  theme_minimal() + theme(legend.position = "none")

p1 / p2

cor.test(all_sims$layer12_image_similarity, all_sims$image_similarity)

    Pearson's product-moment correlation

data:  all_sims$layer12_image_similarity and all_sims$image_similarity
t = 6.154, df = 94, p-value = 1.838e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3757891 0.6649442
sample estimates:
      cor 
0.5359001 

cor.test(all_sims$clip.hf_image_similarity, all_sims$layer12_image_similarity)

    Pearson's product-moment correlation

data:  all_sims$clip.hf_image_similarity and all_sims$layer12_image_similarity
t = 6.149, df = 94, p-value = 1.88e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3754103 0.6646980
sample estimates:
      cor 
0.5355856 

cor.test(all_sims$clip.hf_image_similarity, all_sims$image_similarity)

    Pearson's product-moment correlation

data:  all_sims$clip.hf_image_similarity and all_sims$image_similarity
t = 5521.3, df = 94, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9999977 0.9999990
sample estimates:
      cor 
0.9999985 

CDI

cdi_data <- read.csv(here("data", "main", "data_to_analyze", "level-participant_source-cdi_added-percentile_data.csv"))
looking_data_with_cdi <- looking_data_summarized |> left_join(cdi_data, by=c("SubjectInfo.subjID"="local_id"))
# split-half percentile_rank and plot 
# Median split on CDI percentile_rank
pr_median <- median(looking_data_with_cdi$percentile_rank, na.rm = TRUE)

clip_cdi_bins_summarized <- looking_data_with_cdi |>
  filter(!is.na(percentile_rank)) |>
  mutate(cdi_percentile_bin = if_else(percentile_rank < pr_median,
                           "Low percentile", "High percentile")) |>
  summarize_similarity_data(extra_fields = c("cdi_percentile_bin", "AoA_Est_target"))

similarity_age_plot(clip_cdi_bins_summarized, "image_similarity",
                    age_grouping = "cdi_percentile_bin", model_type = "CLIP",
                    median_age = median_age)

`geom_smooth()` using formula = 'y ~ x'

vocab_score_median <- median(looking_data_with_cdi$vocabulary_score, na.rm = TRUE)
clip_cdi_bins_summarized <- looking_data_with_cdi |>
  filter(!is.na(vocabulary_score)) |>
  mutate(cdi_vocab_bin = if_else(percentile_rank < vocabulary_score,
                           "Low percentile", "High percentile")) |>
  summarize_similarity_data(extra_fields = c("cdi_vocab_bin", "AoA_Est_target"))

similarity_age_plot(clip_cdi_bins_summarized, "image_similarity",
                    age_grouping = "cdi_vocab_bin", model_type = "CLIP",
                    median_age = median_age)

`geom_smooth()` using formula = 'y ~ x'

target_looking_subject_level_cdi <- target_looking_subject_level |> left_join(cdi_data, by=c("SubjectInfo.subjID"="local_id"))
ggplot(data=target_looking_subject_level_cdi, aes(x=mean_value, y=percentile_rank)) +
  geom_point(size=4, alpha=0.6) +
  geom_smooth(method="lm") + theme_classic() + xlab("Baseline-corrected proportion target looking") +
  ylab("Vocabulary percentile rank")

`geom_smooth()` using formula = 'y ~ x'

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

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

cor.test(target_looking_subject_level_cdi$percentile_rank,target_looking_subject_level_cdi$vocabulary_score)

    Pearson's product-moment correlation

data:  target_looking_subject_level_cdi$percentile_rank and target_looking_subject_level_cdi$vocabulary_score
t = 6.4431, df = 36, p-value = 1.785e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5380067 0.8521453
sample estimates:
      cor 
0.7318226