002b-fix-et-data.qmd

---
title: "MB2 Fix up eye-tracking data"
format: html
---

```{r}
library(tidyverse)
library(here)
library(assertthat)
library(patchwork) # MBP2; to create plots with panels.

source(here('helper','ensure_repo_structure.R'))
FIRST_TIME = FALSE
```

# Standardize media names

Remove file media extensions

```{r reload_resamps}
load(file = here(INTERMEDIATE_FOLDER, INTERMEDIATE_002A))

xy <- xy %>% 
  mutate(media_name = tools::file_path_sans_ext(media_name))
```

Begin standardizing data. The main thing we want to do here is validate the media names to make sure that we can use them for merge later. 

When there are invalid media names, you need to put them in the right `txt` files below. 

```{r standardize and validate data}
vec_renaming <- read_csv(here('metadata', 
                              'media_renaming.csv')) %>%
  {setNames(as.character(.$target), .$original)}

media_deletion <- readLines(here('metadata', 
                                 'media_names_to_remove.txt'))

media_names_valid <- readLines(here('metadata', 
                                    'media_names_validate.txt'))

data <- xy |>
  filter(!is.na(media_name) & !(media_name %in% media_deletion)) %>% 
  mutate(media_name = ifelse(media_name %in% names(vec_renaming), vec_renaming[as.character(media_name)], media_name)) %>% 
  group_by(lab_id, participant_id) %>% 
  mutate(event_num = cumsum(c(1, na.omit(media_name != lag(media_name))))) %>% 
  ungroup()
```

Check that all media names fit the appropriate schema. 

```{r checking_names}

# this object helps to determine where invalid media names are coming from in cases where the fix is non-obvious
invalid_media_names <- data %>% 
  distinct(lab_id, participant_id, media_name) %>% 
  filter(!media_name %in% media_names_valid)

invalid_media_names$media_name
unique(invalid_media_names$media_name)

assert_that(nrow(invalid_media_names) == 0)
rm(invalid_media_names)
```


# Column-wise validation

Let's go column by column. 

## lab_id

First, let's unify the lab_ids, then check that they match demographics.

```{r}
lab_ids <- read_csv(here("metadata","labids.csv")) |>
  rename(lab_id = LabID) 

#harmonize lab ids
data <- data %>%
  mutate(
    lab_id = case_when(
      lab_id == "PLUS" ~ "ToMcdlSalzburg",
      lab_id == "Corbitlab" ~ "CorbitLab",
      lab_id == "BabylabNijmegen" ~ "babylabNijmegen",
      lab_id == "MINIDundee" ~ "MiniDundee",
      lab_id == "ccluniri" ~ "cclUNIRI",
      lab_id == "oxfordBabylab" ~ "babylabOxford",
      TRUE ~ lab_id)
  )

assert_that(!any(is.na(data$lab_id)))
```

Now check for match to demos. 


```{r}
load(file = here(INTERMEDIATE_FOLDER, INTERMEDIATE_001_ADULT))
load(file = here(INTERMEDIATE_FOLDER, INTERMEDIATE_001_TODDLER))

```

```{r}
et_adult_labids <- unique(filter(data, 
                                 age_cohort == "adults")$lab_id)
et_toddler_labids <- unique(filter(data, 
                                   age_cohort == "toddlers")$lab_id)

demo_adult_labids <- unique(adult_demo$lab_id)
demo_toddler_labids <- unique(toddler_demo$lab_id)
```

Check overlap/non-overlap in IDs. First adults.

```{r}
print(paste("demo IDs not in ET:", 
            dplyr::setdiff(demo_adult_labids, et_adult_labids)))

print(paste("ET IDs not in demo:", 
            dplyr::setdiff(et_adult_labids, demo_adult_labids)))
```
Now toddlers. 

```{r}
print(paste("demo IDs not in ET:", 
            dplyr::setdiff(demo_toddler_labids, et_toddler_labids)))

print(paste("ET IDs not in demo:", 
            dplyr::setdiff(et_toddler_labids, demo_toddler_labids)))
```

Now the test.

```{r}
assert_that(is_empty(dplyr::setdiff(demo_adult_labids, et_adult_labids)))
assert_that(is_empty(dplyr::setdiff(et_adult_labids, demo_adult_labids)))
assert_that(is_empty(dplyr::setdiff(demo_toddler_labids, et_toddler_labids)))
assert_that(is_empty(dplyr::setdiff(et_toddler_labids, demo_toddler_labids)))

```

## participant_id

```{r}
correct_ids <- function(df) {
  renaming <- read.csv2(here('metadata','participant_id_renaming.csv')) %>% 
    mutate(
      sample = paste0(sample,'s'),
      old = paste0(old,'-',sample),
      new = paste0(new,'-',sample)
      )
  
  # Merge the input dataframe with the renaming dataframe
  df <- df %>%
    left_join(renaming, by = c("lab_id" = "lab", "age_cohort" = "sample", "participant_id" = "old"))
  
  # Replace the participant_id with the correct id
  df <- df %>%
    mutate(participant_id = ifelse(is.na(new), participant_id, new)) %>%
    select(-new)
  
  # Additional cleanup
  df <- df %>%
    mutate(participant_id = case_when(
      lab_id == 'gertlabLancaster' ~ paste0("LU_", participant_id),
      age_cohort == 'toddlers' & lab_id == 'irlConcordia' ~ sub("MB2_", "", participant_id),
      age_cohort == 'adults' & lab_id == 'careylabHarvard' ~ paste0('careylabHarvard_', participant_id),
      age_cohort == 'adults' & lab_id == 'MEyeLab' ~ sub("\\.edf", "", participant_id),
      age_cohort == 'toddlers' & lab_id == 'careylabHarvard' ~ sub("careylabHarvard", "careylabharard", participant_id),
      lab_id == 'UIUCinfantlab' ~ paste0("0", substr(participant_id, 1, 2), "-", age_cohort),
      TRUE ~ participant_id
    ))

  return(df)
}

data <- correct_ids(data)
```

## Clean demographics

Create combined demographics files and prep for joining with eye-tracking data.

```{r}
#bind together adult and toddler demographics
#necessary in order to avoid issues with the merge
combined_demo <- adult_demo |>
  mutate(test_order = as.character(test_order)) |>
  mutate(age_cohort = "adults") |> #just to make sure that nothing goes wrong with the merge
  bind_rows(
    toddler_demo |>
      mutate(test_order = as.character(test_order)) |>
      mutate(age_cohort = "toddlers")
    )
```

## Remove the participants with session errors from the et data 

```{r}
# extract pilot ids
pilot_ids_all <- combined_demo %>%
  filter(pilot == "yes") %>%
  select(participant_id, lab_id,age_cohort)

#extract ids for participants with session errors
error_ids_all <- combined_demo %>%
  filter(session_error == 'error' & pilot != "yes") %>%
  select(participant_id, lab_id, session_error, session_error_info, session_error_notes,age_cohort)

#filter session error and pilot participants from the data
data <- data %>%
  anti_join(error_ids_all, by=join_by(participant_id, lab_id, age_cohort)) %>%
  anti_join(pilot_ids_all, by=join_by(participant_id, lab_id, age_cohort))

# remove the id-cohort suffixes for matching with metadatafiles and better readability in the output
error_ids_all <- error_ids_all %>%
  mutate(participant_id = str_replace(participant_id, "-[^-]*$", ""))
pilot_ids_all <- pilot_ids_all %>%
  mutate(participant_id = str_replace(participant_id, "-[^-]*$", ""))

known_missings <- read.csv(here("metadata","known_missing.csv"))
still_missings <- known_missings %>%
  anti_join(error_ids_all, by=join_by(participant_id, lab_id, age_cohort)) %>% 
  anti_join(pilot_ids_all, by=join_by(participant_id, lab_id, age_cohort))

write_csv(error_ids_all,
          here(INTERMEDIATE_FOLDER, INTERMEDIATE_000_EXCL_SESSION))

if(nrow(still_missings) > 0){
  write.csv(still_missings, here(INTERMEDIATE_FOLDER, INTERMEDIATE_000_EXCL_TEMP))
}

```


```{r}
et_adult_pids <- filter(data, age_cohort == "adults") |>
  select(lab_id, participant_id) |>
  distinct() 

et_toddler_pids <- filter(data, age_cohort == "toddlers") |>
  select(lab_id, participant_id) |>
  distinct() 

demo_adult_pids <- adult_demo |>
  filter(session_error == "noerror" & pilot != "yes") |>
  select(lab_id, participant_id) |>
  distinct() 

demo_toddler_pids <- toddler_demo |>
  filter(session_error == "noerror" & pilot != "yes") |>
  select(lab_id, participant_id) |>
  distinct()
```

Check overlap/non-overlap in IDs. First adults.

```{r}
metadata <- read.csv(here("metadata","known_missing.csv"))

# Filter for adults
metadata_adults <- metadata %>% filter(age_cohort == "adults") %>% mutate(participant_id = paste0(participant_id,'-adults'))

# Remove matches from et_adult_pids
et_ids_not_in_demo_adults <- anti_join(et_adult_pids, demo_adult_pids) %>%
  anti_join(metadata_adults, by = c("lab_id", "participant_id"))

# Remove matches from demo_adult_pids
demo_ids_not_in_et_adults <- anti_join(demo_adult_pids, et_adult_pids) %>%
  anti_join(metadata_adults, by = c("lab_id", "participant_id"))

unique(c(et_ids_not_in_demo_adults$lab_id, demo_ids_not_in_et_adults$lab_id))
```

Now toddlers. 

```{r}
metadata <- read.csv(here("metadata","known_missing.csv"))

metadata_toddlers <- metadata %>% filter(age_cohort == "toddlers") %>% mutate(participant_id = paste0(participant_id,'-toddlers'))

# Remove matches from et_toddler_pids
et_ids_not_in_demo_toddlers <- anti_join(et_toddler_pids, demo_toddler_pids) %>%
  anti_join(metadata_toddlers, by = c("lab_id", "participant_id"))

# Remove matches from demo_toddler_pids
demo_ids_not_in_et_toddlers <- anti_join(demo_toddler_pids, et_toddler_pids) %>%
  anti_join(metadata_toddlers, by = c("lab_id", "participant_id"))

unique(c(et_ids_not_in_demo_toddlers$lab_id, demo_ids_not_in_et_toddlers$lab_id))
```

Now the test.

```{r}
assert_that(nrow(anti_join(et_adult_pids, demo_adult_pids))==0)
#assert_that(is_empty(anti_join(demo_adult_pids, et_adult_pids))) # unfixable, as data is missing
assert_that(nrow(anti_join(et_toddler_pids, demo_toddler_pids))==0)
#assert_that(is_empty(anti_join(demo_toddler_pids, et_toddler_pids))) # unfixable, as data is missing
```

Test for uniqueness: first, uniqueness across labs.
```{r}
#joining participant and lab id to create a unique participant id
data <- data |> 
  unite(participant_lab_id, lab_id, participant_id, remove = F)

ppt_id_per_lab <- data |> 
  select(participant_lab_id, lab_id) |> 
  distinct() |> 
  group_by(participant_lab_id) |> 
  summarise(n = n())
assert_that(all(ppt_id_per_lab$n == 1)) 
```

Then, uniqueness within labs. 
Our strategy is to detect instances at which participant_id changes; given the assumptions about the data structure, these should uniquely capture the rows at which the data switch to a new participant.
```{r}
ppt_id_unique <- data |> 
  mutate(t_lag = lag(t),
         ppt_lag = lag(participant_lab_id)) |> 
  filter(participant_lab_id != ppt_lag)

assert_that(nrow(ppt_id_unique) + 1 == n_distinct(data$participant_lab_id))
```


## media_name
Check for uniqueness within participants.
```{r}
media_unique <- data |> 
  mutate(med_lag = lag(media_name)) |> 
  filter(media_name != med_lag) |> 
  group_by(participant_lab_id, media_name) |> 
  summarise(n_media = n())
assert_that(all(media_unique |> 
                  filter(media_name != "star_calib") |> 
                  pull(n_media) == 1))
# some participants viewed star_calib twice, which is okay
```


## x

```{r}
data |>
  group_by(lab_id) |>
  summarise(min_x = min(x, na.rm=TRUE), 
            max_x = max(x, na.rm=TRUE),
            median_x=median(x,na.rm=TRUE),
            zero_count = sum(x==0,na.rm=T)) 

#visualize distribution of x coordinates for all labs
ggplot(data,aes(x))+
  geom_histogram()+
  facet_wrap(~lab_id, scales="free")

#Follow up on x distribution for outlier labs. 
data |>
  filter(lab_id == "babylabINCC") |>
  ggplot(aes(x = x)) + 
  geom_histogram()

```

Tests.

```{r, eval=F}
# doesn't work---there are negative values and values on the order of 1e5
assert_that(all(data$x > 0))
assert_that(all(data$x < 5000))
```


## y

Lots of very negative data. Not sure what to do with this. 

```{r}
data |>
  group_by(lab_id) |>
  summarise(min_y = min(y, na.rm=TRUE), 
            max_y = max(y, na.rm=TRUE),
            median_x=median(x,na.rm=TRUE),
            zero_count = sum(x==0,na.rm=T))

#visualize distribution of y coordinates for all labs
ggplot(data,aes(y))+
  geom_histogram()+
  facet_wrap(~lab_id, scales="free")
```

Tests.

```{r, eval=F}
# doesn't work---there are negative values and values on the order of 1e5
assert_that(all(data$y > 0))
assert_that(all(data$y < 5000))
```

## t

t looks good (difference between web-based and in-person is due to not running the long trials online).

```{r}
data |>
  group_by(lab_id, participant_trial_id) |>
  mutate(t_norm = t - t[1]) |>
  group_by(lab_id) |>
  summarise(min_t = min(t_norm, na.rm=TRUE), 
            max_t = max(t_norm, na.rm=TRUE), 
            mean_diff = median(diff(t_norm), na.rm=TRUE))

```

## pupil_left & pupil_right

```{r}
data |>
  group_by(lab_id) |>
  summarise(min_pupil_left = min(pupil_left, na.rm=TRUE),
            max_pupil_left = max(pupil_left, na.rm=TRUE), 
            min_pupil_right = min(pupil_right, na.rm=TRUE),
            max_pupil_right = max(pupil_right, na.rm=TRUE)) 

# filter out web-based labs
data.pupil <- data |> filter(lab_id!="affcogUTSC" & lab_id!="beinghumanWroclaw" & lab_id!="CogConcordia" & lab_id!="collabUIOWA" & lab_id!="saxelab" & lab_id!="socallabUCSD")

# look at summaries again
data.pupil |>
  group_by(lab_id) |>
  summarise(min_pupil_left = min(pupil_left, na.rm=TRUE),
            max_pupil_left = max(pupil_left, na.rm=TRUE), 
            min_pupil_right = min(pupil_right, na.rm=TRUE),
            max_pupil_right = max(pupil_right, na.rm=TRUE)) 

# warnings from babylabMPIB, babylingOslo, tauccd, kidsdevUniofNewcastle due to Inf. Remove for now, then follow up.
# looking at second test trial only.
data.pupil.overview <- data.pupil |>
  filter(lab_id!="babylabMPIB" & lab_id!="babylingOslo" & lab_id!="tauccd" & lab_id!="kidsdevUniofNewcastle") |>
  filter(grepl("CON", media_name) | grepl("INC", media_name)) |>
  mutate(pupil_left = na_if(pupil_left, 0)) |>
  mutate(pupil_right = na_if(pupil_right, 0)) |>
  group_by(lab_id) |>
  summarise(mean_pupil_left = mean(pupil_left, na.rm=TRUE),
            mean_pupil_right = mean(pupil_right, na.rm=TRUE),
            sd_pupil_left = sd(pupil_left, na.rm=TRUE),
            sd_pupil_right = sd(pupil_right, na.rm=TRUE),
            min_pupil_left = min(pupil_left, na.rm=TRUE),
            max_pupil_left = max(pupil_left, na.rm=TRUE), 
            min_pupil_right = min(pupil_right, na.rm=TRUE),
            max_pupil_right = max(pupil_right, na.rm=TRUE)) 

# Here shorter version:
data.pupil.overview.short <- data.pupil |>
  filter(lab_id!="babylabMPIB" & lab_id!="babylingOslo" & lab_id!="tauccd" & lab_id!="kidsdevUniofNewcastle") |>
  filter(grepl("CON", media_name) | grepl("INC", media_name)) |>
  mutate(pupil_left = na_if(pupil_left, 0)) |>
  mutate(pupil_right = na_if(pupil_right, 0)) |>
  group_by(lab_id, age_cohort) |>
  filter(mean(pupil_left, na.rm=T) < 10) |>
  ungroup() |>
  group_by(age_cohort) |>
  summarise(mean_pupil_left = mean(pupil_left, na.rm=TRUE),
            mean_pupil_right = mean(pupil_right, na.rm=TRUE),
            sd_pupil_left = sd(pupil_left, na.rm=TRUE),
            sd_pupil_right = sd(pupil_right, na.rm=TRUE),
            min_pupil_left = min(pupil_left, na.rm=TRUE),
            max_pupil_left = max(pupil_left, na.rm=TRUE), 
            min_pupil_right = min(pupil_right, na.rm=TRUE),
            max_pupil_right = max(pupil_right, na.rm=TRUE),
            N = length(unique(participant_lab_id)),
            mean_X = mean(x, na.rm=TRUE),
            sd_X = sd(x, na.rm=TRUE),
            min_X = min(x, na.rm=TRUE),
            max_X = max(x, na.rm=TRUE),
            mean_Y = mean(y, na.rm=TRUE),
            sd_Y = sd(y, na.rm=TRUE),
            min_Y = min(y, na.rm=TRUE),
            max_Y = max(y, na.rm=TRUE),
            )

write.table(data.pupil.overview.short, file="./MB2P/summaries4simulation.txt", col.names=T, row.names=F, sep="\t")

to.plot <- 0 # Set this to 1 if you want 1000+ plots to be created.
if(to.plot == 1){
#uniq.labs <- unique(data.pupil$lab_id)
uniq.labs <- unique(data.pupil$participant_lab_id)
for(a in 1:length(uniq.labs)){
  #lab.data <- data.pupil |> filter(lab_id == uniq.labs[a]) |> filter(grepl("CON", media_name) | grepl("INC", media_name))
  # ggplot(lab.data, aes(t, pupil_left)) + geom_point(alpha = 0.1, size=0.1) +   facet_wrap(~age_cohort, scales="free")
  # ggplot(lab.data, aes(t, pupil_left)) + geom_line(alpha = 0.1, size=0.1) +   facet_wrap(~age_cohort, scales="free")
  # ggplot(lab.data, aes(pupil_right, pupil_left)) + geom_point(alpha = 0.1, size=0.1) +   facet_wrap(~age_cohort, scales="free")

  lab.data <- data.pupil |> filter(participant_lab_id == uniq.labs[a]) |> filter(grepl("CON", media_name) | grepl("INC", media_name))
  suppressMessages({
    suppressWarnings({
  plot1 <- ggplot(lab.data, aes(pupil_left)) + geom_histogram() 
  plot2 <- ggplot(lab.data, aes(t, pupil_left)) + geom_point(size=0.2) + ggtitle(as.character(uniq.labs[a])) + theme(plot.title = element_text(size = 10))
  plot3 <- ggplot(lab.data, aes(pupil_right, pupil_left)) + geom_point(size=0.1)
  plot4 <- ggplot(lab.data, aes(pupil_right)) + geom_histogram() 
  plot5 <- ggplot(lab.data, aes(t, pupil_right)) + geom_point(size=0.2)
  plot6 <- ggplot(lab.data, aes(x, y)) + geom_point(size=0.1)
  combined_plot <- (plot1 | plot2 | plot3) / (plot4 | plot5 | plot6)
  ggsave(plot = combined_plot, filename=paste("./MB2P/MB2P-diagnosticPlots/", uniq.labs[a],".png"), bg="white",width=9, height=6, device="png", dpi=100)
  })
  })
  #print(combined_plot)
}
}

# visualize distribution of pupil left for all labs
ggplot(data.pupil,aes(pupil_left))+
  geom_histogram()+
  facet_wrap(~lab_id, scales="free")+
  ylab("Left Pupil Size")

# visualize distribution of pupil right for all labs
ggplot(data.pupil,aes(pupil_right))+
  geom_histogram()+
  facet_wrap(~lab_id, scales="free")+
  ylab("Right Pupil Size")

```

## age_cohort

```{r}
assert_that(all(data$age_cohort %in% c("adults","toddlers")))
```

## event_num

Why are there some event numbers that go to 17?

Only for babylabOxford. Other labs go up to 9 event numbers which is due to calibration taking up an event number.
```{r}
unique(data$event_num)
```

```{r}
hist(data$event_num, breaks = 0:17)
```

Event 7 is overrepresented as it refers to calibration

```{r}
data |>
  group_by(lab_id, age_cohort, event_num) |>
  count()

```

# Exclusions

## Age Exclusions

First, process age information in the demographic file
```{r}
hist(combined_demo$age_days_num)

combined_demo <- combined_demo |>
  mutate(age_mo = case_when(
    #check age conversion we want to use here
    age_cohort=="toddlers" ~ age_days_num / 365.25*12,
    TRUE ~ NA_real_))

#quick checks
hist(combined_demo$age_mo)
hist(combined_demo$age_years_n)
table(filter(combined_demo,age_cohort=="adults")$age_years_n)
```

Now add columns tracking exclusions based on age. 

```{r}
# applying exclusions from the registered report and lab manual
# toddlers: 548 - 822 days
# adults: 18 - 55 years
combined_demo <- combined_demo |>
  mutate(
    age_exclusion = case_when(
      age_cohort == "adults" & age_years_n < 18 ~ "yes",
      age_cohort == "adults" & age_years_n >= 56 ~ "yes",
      age_cohort == "adults" & is.na(age_years_n) ~ "yes",
      age_cohort == "toddlers" & age_days_num < 548 ~ "yes",
      age_cohort == "toddlers" & age_days_num > 822 ~ "yes",
      age_cohort == "toddlers" & is.na(age_days_num) ~ "yes",
      TRUE ~ "no"
    )
  )
View(select(combined_demo,age_cohort,lab_id,participant_id,age_years_n,age_mo,age_days_num,age_exclusion))
```

# Saving 

Now save the merged xy data locally.

```{r saving}
save(data, file = here(INTERMEDIATE_FOLDER, INTERMEDIATE_002))
save(combined_demo, file = here(INTERMEDIATE_FOLDER, INTERMEDIATE_002B))
```