---
title: "R3 final: Parietal cortex's role in tracking adaptive memory feature"
author: "Xi Yang, Yufei Zhao"
date: "`r format(Sys.time(), '%b-%d-%Y %H:%M')`"
output: 
  flexdashboard::flex_dashboard:
    orientation: rows
    social: menu
    source_code: embed
---

```{r pkg, message = FALSE, warning = FALSE}
# to run the codes, needed to install {pacman} package in addition to {ggpubr}; the {magrittr} pkg is also needed
pacman::p_load(tidyverse, here, fs, glue, afex, emmeans, knitr, psych, flexdashboard)
# afex is used for stats modeling
```

```{r default,  message = FALSE, warning = FALSE}
# separate in different lines to keep things clear
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
knitr::opts_chunk$set(dpi = 300)
# knitr::opts_chunk$set(message = FALSE, warning = FALSE, dpi = 300)
# The revision in this line is based on suggestions from Cameron and Krista. 
afex_options(emmeans_model = "multivariate")
theme_set(ggpubr::theme_pubclean())
```

```{r dir-setting}
# Directories to bids form data
# through the whole scripts we are using path in case someone are going to run the codes over
# different operation system

bids_dir <- path(here())
deriv_dir <- path(bids_dir, "derivatives")
data_dir <- path(deriv_dir, "pattern_similarity")
# data_dir <- here("derivatives", "pattern_similarity")
# This last line is suggested by Cameron and Krista. We included it here to show a different way to code.
```

```{r sub-list}
# Subject list
subj_list <- 
  read_tsv(path(bids_dir, "participants.tsv"), col_types = "cicc") %>% 
  select(participant_id) %>% 
  separate(participant_id, c("prefix", "id"), sep = "-") %>% 
  pull(id)

# subj_list <- read_tsv(here("derivatives",  "pattern_similarity",  "participants.tsv"), col_types = "cicc")
# This last line is suggested by Cameron, Krista and Daniel. We included it here to show a different way to code
```

```{r preprocess-flavor}
# fMRI data preprocessing flavor
# pattern similarity data were processed with python
preproc_id <-  "hp0p01-smooth1p7"
estimate_id <- "tmaps"
flavor_id <- ""
```

```{r helper-function}
# read in codes with char to avoid possible padding zeros to be droped.
# it is not an issue here, but this works for all my different projects.

# function 1
# read in pattern similarity label map (item-wise)
read_labels <-
  function(label_file) {
    read_tsv(label_file, na = "999") %>%
    mutate_if(is.double, as.integer)
  }
# Changes made based on Cameron's suggestions

# function 2
# read in pattern similarity data map
read_simi <-
  function(simi_file, label_file) {
    # load data
    simi <- read_tsv(simi_file, col_types = "ccccdd")
    # load labels
    labels <- read_labels(label_file)
    # merge labels into data table
    left_join(simi, labels, by = c("subj_id", "pair_index"))
  }
# Changes made based on Cameron's suggestions
```

```{r read-data}
# map 1
simi_list_rtv <- 
  map_chr(
  subj_list, ~ path(
    data_dir,
    glue("sub-", .x),
    "roi_pearson_similarity",
    glue(
      "sub-",
      .x,
      "_space-T1w_task-retrieval_desc-{preproc_id}-{estimate_id}-zr{flavor_id}_similarity.tsv"
    )
  )
)

# map 2
labels_list <-
  map_chr(
    subj_list, ~ path(
      data_dir,
      glue("sub-", .x),
      "beh_labels",
      glue("sub-", .x, "_labels_pairwise.tsv")
    )
  )

# map2_*
# get the combined data
simi_rtv_raw <- map2_df(simi_list_rtv, labels_list, ~ read_simi(.x, .y))
```

```{r reduce-unrelated-pairs}
# the raw data have full pairwise pattern similarity results
# reduce the two items that are not in the same run
simi_rtv_reduced <- simi_rtv_raw %>% 
  filter(!i1_run_type == i2_run_type,
         i1_group == i2_group) %>% 
  mutate(group = (i1_group + i2_group)/2) %>% 
  select(ends_with("id"), 
         roi_hemi, 
         similarity, 
         ends_with("item"),
         group,
         ends_with("obj"))
# more concise style changes based on Cameron's suggestions

# atlas_id: the brain altas name
# roi_id: the region of interest
# roi_hemi: left/right/bi
```

```{r separate-within-item}
# get within-item data and non within-item data
simi_rtv_within <- simi_rtv_reduced %>% 
  filter(i1_item == i2_item)
simi_rtv_non_within <- simi_rtv_reduced %>% 
  filter(!i1_item == i2_item)
```

```{r label-same-pair}
# for non within-item data label info that can help with lablling category info
simi_rtv_non_within_pair <- simi_rtv_non_within %>% 
  mutate(
    pair_match = case_when(
      (i1_item < i2_item) ~ sprintf("%02d%02d", i1_item,i2_item),
      (i1_item > i2_item) ~ sprintf("%02d%02d", i2_item,i1_item)
    ),
    pair_match_obj = case_when(
      (i1_obj < i2_obj) ~ sprintf("%02d%02d", i1_obj,i2_obj),
      (i1_obj > i2_obj) ~ sprintf("%02d%02d", i2_obj,i1_obj)
    )
  ) %>% 
  group_by(subj_id, atlas_id, roi_id, roi_hemi, pair_match, pair_match_obj) %>% 
  summarise(similarity = mean(similarity),
            group = mean(group)) %>% 
  separate(col = pair_match, into = c("i1_item", "i2_item"), sep = 2) %>% 
  separate(col = pair_match_obj, into = c("i1_obj", "i2_obj"), sep = 2) %>%
  ungroup()

simi_rtv_non_within_pair <- simi_rtv_non_within_pair %>%
  mutate_at(vars(starts_with("i")), as.integer)
# changes made based on Cameron's suggestions
simi_rtv_unlabel <- rbind(simi_rtv_within, simi_rtv_non_within_pair)
```

```{r label-category}
# lable all categories
# we dont choose any "better" way to shorten this part
# because the category labeling is under conceptual discussion
# in my actual data set we make even more categories to examine the 
# research question closely. so this part can't be simplified

simi_rtv <-  simi_rtv_unlabel %>% 
  mutate(
    simi_cond = case_when(
      (i1_item == i2_item) ~ 'within_item',
      (i1_item == 1 & i2_item == 3) ~ 'within_pair',
      (i1_item == 2 & i2_item == 4) ~ 'within_pair',
      (i1_item == 5 & i2_item == 7) ~ 'within_pair',
      (i1_item == 6 & i2_item == 8) ~ 'within_pair',
      (i1_item == 9 & i2_item == 11) ~ 'within_pair',
      (i1_item == 10 & i2_item == 12) ~ 'within_pair',
      (i1_item == 13 & i2_item == 15) ~ 'within_pair',
      (i1_item == 14 & i2_item == 16) ~ 'within_pair',
      (i1_item == 17 & i2_item == 19) ~ 'within_pair',
      (i1_item == 18 & i2_item == 20) ~ 'within_pair',
      (i1_item == 21 & i2_item == 23) ~ 'within_pair',
      (i1_item == 22 & i2_item == 24) ~ 'within_pair',
      (i1_item == 1 & i2_item == 2) ~ 'same_color',
      (i1_item == 3 & i2_item == 4) ~ 'same_color',
      (i1_item == 5 & i2_item == 6) ~ 'same_color',
      (i1_item == 7 & i2_item == 8) ~ 'same_color',
      (i1_item == 9 & i2_item == 10) ~ 'same_color',
      (i1_item == 11 & i2_item == 12) ~ 'same_color',
      (i1_item == 13 & i2_item == 14) ~ 'same_color',
      (i1_item == 15 & i2_item == 16) ~ 'same_color',
      (i1_item == 17 & i2_item == 18) ~ 'same_color',
      (i1_item == 19 & i2_item == 20) ~ 'same_color',
      (i1_item == 21 & i2_item == 22) ~ 'same_color',
      (i1_item == 23 & i2_item == 24) ~ 'same_color',
      (i1_item == 3 & i2_item == 5) ~ 'between_pair',
      (i1_item == 3 & i2_item == 6) ~ 'between_pair',
      (i1_item == 4 & i2_item == 5) ~ 'between_pair',
      (i1_item == 4 & i2_item == 6) ~ 'between_pair', 
      (i1_item == 7 & i2_item == 9) ~ 'between_pair',
      (i1_item == 7 & i2_item == 10) ~ 'between_pair',
      (i1_item == 8 & i2_item == 9) ~ 'between_pair',
      (i1_item == 8 & i2_item == 10) ~ 'between_pair',
      (i1_item == 15 & i2_item == 17) ~ 'between_pair',
      (i1_item == 15 & i2_item == 18) ~ 'between_pair',
      (i1_item == 16 & i2_item == 17) ~ 'between_pair',
      (i1_item == 16 & i2_item == 18) ~ 'between_pair',
      (i1_item == 19 & i2_item == 21) ~ 'between_pair',
      (i1_item == 19 & i2_item == 22) ~ 'between_pair',
      (i1_item == 20 & i2_item == 21) ~ 'between_pair',
      (i1_item == 20 & i2_item == 22) ~ 'between_pair',
      TRUE ~ "other"
    )
  )
```


```{r}
# Average similarity across different condition (within-item, within-pair, between-pair) within subject
# select regions in parietal cortex, visual cortex, and hippocampus
dat_rtv <- simi_rtv %>%
  group_by(group, subj_id, atlas_id, roi_id, roi_hemi, simi_cond) %>%
  summarise(similarity = mean(similarity)) %>%
  ungroup() %>% 
  filter(!simi_cond == 'other') %>% 
  filter(roi_id %in% c("ANG", "IPS","SPL","HPC","OTC"))
dat_rtv$group <- recode_factor(dat_rtv$group, `1` = "Paired", `2` = "Control")
dat_rtv$simi_cond <- factor(dat_rtv$simi_cond, levels = c("within_item", "same_color", "between_pair", "within_pair"))
```

### Pattern stimilarity results for each condition by group
```{r}
# average data across subject
# select only bilateral roi
dat_rtv_all <- dat_rtv %>% 
  group_by(group, atlas_id, roi_id, roi_hemi, simi_cond) %>% 
  summarise(similarity = mean(similarity)) %>% 
  spread(key = "simi_cond", value = "similarity") %>% 
  filter(roi_hemi == "bilateral")

kable(dat_rtv_all, digits = 2)
# Specifying digits for readability. This change is based on Krista's suggestions. 
```

### Pure color representation: same_color - between_pair
```{r}
# pure color: same_color - between_pair
# nest
tmp <- dat_rtv %>% 
  group_by(atlas_id, roi_id, roi_hemi, subj_id, group) %>% 
  nest() %>% 
  mutate(
    data = map(data, 
               ~ data.frame(spread(data = ., key = simi_cond, value = similarity))),
    group_contrast = map(data, 
                         ~ mutate(., a = (same_color - between_pair)) %>% # substract between_pair similarity from the same_color similarity
                           pull()),
    group_contrast = as.numeric(group_contrast)
  ) %>% 
  select(-data) %>% 
  unnest() %>% 
  group_by(atlas_id, roi_id, roi_hemi) %>% 
  nest() %>% 
  mutate(
    one_way = map(data, 
                  ~ aov_ez(data = ., id = "subj_id", dv = "group_contrast", within = "group")) # one-way within-subject anova
  )

# get model based summary table for plotting
# error bar indicates within-subject errors
dat_plot <-
  tmp %>%
  mutate(
    model_based_summary = map(
      one_way,
      ~ afex_plot(.x, x = "group", error = "within", error_ci = FALSE, return = "data")
    ),
    summary_table = map(
      model_based_summary,
      ~ .x[[1]])) %>% 
  unnest(summary_table) %>% 
  select(atlas_id, roi_id, roi_hemi, x, y, error, lower, upper) %>%
  rename(group = x, similarity = y)

# removed magrittr::extract2 based on Cameron's suggestions
```

```{r }
dat_plot %>% 
  filter(atlas_id == "Destrieux",
         roi_hemi == "bilateral") %>% 
  ggplot(aes(x = group, y = similarity, fill = group)) +
  facet_wrap(~ roi_id) +
  geom_bar(stat = "identity", width = 0.5) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, size = 0.7)
```