```{r setup}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

```{r}
library(janitor)
library(lubridate)
library(scales)
library(readxl)
library(ggthemes)
library(drc)
library(broom)
library(tidyverse)
library(emmeans)
library(ggpubr)
library(joeyr)
```

# CBASS temperature logs

## Import CBASS log files
```{r}
b1log <- read_csv("data/raw/20200820_KJSNSU_Acer/Log_Files/B1_LOG_0821.TXT")

b1log <- b1log %>%
  filter(PrintDate != "PrintDate") %>%
  mutate(date = as_date(Date, format = "%Y_%B_%d")) %>%
  unite(time, Th, Tm, Ts, sep = ":") %>%
  unite(dttm, date, time) %>%
  mutate(dttm = ymd_hms(dttm),
         hm = format(dttm, "%H:%M")) %>%
  select(dttm, hm, T1SP, TempT1, T2SP, TempT2, T3SP, TempT3, T4SP, TempT4) %>%
  pivot_longer(starts_with("T"), names_to = "key", values_to = "temp") %>%
  filter(temp > 0, temp < 50) %>%
  mutate(tank = str_extract(key, "T[0-9]"),
         tank = paste0("B1", tank),
         key = case_when(grepl("SP", key) ~ str_sub(key, 3, 4),
                         TRUE ~ str_sub(key, 1, 4))) %>%
  pivot_wider(names_from = key, values_from = temp)

b2log <- read_csv("data/raw/20200820_KJSNSU_Acer/Log_Files/B2_LOG_0821.TXT")

b2log <- b2log %>%
  filter(PrintDate != "PrintDate") %>%
  mutate(date = as_date(Date, format = "%Y_%B_%d")) %>%
  unite(time, Th, Tm, Ts, sep = ":") %>%
  unite(dttm, date, time) %>%
  mutate(dttm = ymd_hms(dttm),
         hm = format(dttm, "%H:%M")) %>%
  select(dttm, hm, T1SP, TempT1, T2SP, TempT2, T3SP, TempT3, T4SP, TempT4) %>%
  pivot_longer(starts_with("T"), names_to = "key", values_to = "temp") %>%
  filter(temp > 0, temp < 50) %>%
  mutate(tank = str_extract(key, "T[0-9]"),
         tank = paste0("B2", tank),
         key = case_when(grepl("SP", key) ~ str_sub(key, 3, 4),
                         TRUE ~ str_sub(key, 1, 4))) %>%
  pivot_wider(names_from = key, values_from = temp)

log <- bind_rows(b1log, b2log) %>%
  mutate(tank = factor(tank),
         SP = as.numeric(SP),
         Temp = as.numeric(Temp))
```

## Import CBASS set points (used for settings file)
```{r}
sp <- read_xlsx("data/raw/20200820_KJSNSU_Acer/Temp_Profiles/CRII_ACERV_2020.xlsx") %>%
  mutate(Time = format(Time, "%H:%M")) %>%
  filter(Time >= "13:00", Time <= "22:00") %>%
  crossing(date = as_date(c("2020-08-20", "2020-08-21"))) %>%
  mutate(dttm = ymd_hm(paste(date, Time))) %>%
  select(dttm, matches("B")) %>%
  pivot_longer(-dttm, names_to = "tank", values_to = "Temp") %>%
  arrange(dttm)
```

## Plot CBASS temperature profiles with set points
```{r, fig.width = 10}
ggplot(mapping = aes(x = dttm, y = Temp, group = tank)) +
  geom_line(data = filter(log, hm > "12:50", hm < "21:50"),
            aes(color = tank), lwd = 0.2) +
  geom_step(data = sp, lwd = 0.2) +
  facet_wrap(~date(dttm), scales = "free") +
  scale_y_continuous(breaks = 30:37) +
  scale_x_datetime(breaks = "hours", labels = label_date("%H:%M")) +
  theme_hc() +
  theme(legend.position = "none") +
  labs(x = "Time of day", y = "Temperature (°C)")

#  Calculate actual hold temperatures
actual_temps <- log %>%
  filter(hm > "16:15" & hm < "19:05") %>%
  group_by(tank_id = tank) %>%
  summarise(max_temp = mean(Temp, na.rm = TRUE)) %>%
  arrange(max_temp)

actual_temps
```

# PAM data

## Turn each filtering and adjustment step on or off
```{r, echo = TRUE}
# Initial pre-filtering (abnormally high, abnormally high w/ low Ft)
pre.filter <- T

# Post-model-fit filtering (high cook's distance outliers)
mod.filter <- T

# Adjust for positional effects
pos.adjust <- T

# Use actual max temperatures instead of target temperatures
temp.adjust <- T
```

# Data pre-filtering and QC
```{r}
dpii <- read_csv("data/raw/20200820_KJSNSU_Acer/PAM/KJS_pamdata_tidy.csv")

dpii <- dpii %>%
  mutate(nursery = "nsu",
         geno = paste0("kjs", Genotype)) %>%
  select(nursery, geno, f, fm, fvfm, target_temp = Final_Temp, tank_id = Tank) %>%
  drop_na(fvfm)

# Replace max_temp with actual max_temp calculated above
dpii <- dpii %>%
  left_join(actual_temps)

# Identify points at high temperatures where low background fluorescence (f) resulted in abnormally high values of fvfm
out <- dpii %>%
  filter(target_temp >= 36) %>%
  drop_na(fvfm) %>%
  mutate(mahal = tidy_mahalanobis(f, fvfm),
         is_out = mahal > quantile(mahal, 0.75) & f < median(f) & fvfm > median(fvfm, na.rm = T))
# Plot showing those outliers
ggplot(out, aes(x = f, y = fvfm)) +
  geom_point(aes(shape = is_out, color = factor(target_temp))) +
  scale_shape_manual(values = c(19, 25))

df <- dpii %>%
  # Save raw fvfm data in new column
  mutate(fvfmraw = fvfm) %>% 
  # Merge with mahalanobis outlier info from above
  left_join(out) %>%
  # Identify problematic data points
  mutate(problem = case_when(
    fvfm > 0.75 ~ "abnormally high",
    is_out ~ "abnormally high w/ low Ft",
    TRUE ~ "none")) %>%
  mutate(fvfm = case_when(
    problem == "abnormally high" ~ NA_real_,    # Change abnormally high values to NA
    problem == "abnormally high w/ low Ft" ~ NA_real_,
    problem == "none" ~ fvfm))

# Select pre-filtered or raw data to use moving forward based on set option
df <- if (pre.filter) { 
  df 
} else { 
  dpii %>% 
    mutate(fvfmraw = fvfm, 
           problem = case_when(is.na(fvfm) ~ "no signal",
                               TRUE ~ "none") )
}
```

# Adjust for positional effects
```{r}
# Import maps of the position of each genotype in each tank
map1 <- read_csv("data/raw/20200820_KJSNSU_Acer/20200820_map.csv") %>%
  mutate(across(.fns = as.character)) %>%
  pivot_longer(-Row, values_to = "geno", names_to = "Col") %>%
  drop_na() %>%
  unite(pos, Row, Col, sep = "", remove = FALSE)
map2 <- read_csv("data/raw/20200820_KJSNSU_Acer/20200821_map.csv") %>%
  pivot_longer(-Row, values_to = "geno", names_to = "Col") %>%
  mutate(across(.fns = as.character)) %>%
  drop_na() %>%
  unite(pos, Row, Col, sep = "", remove = FALSE)

df_pos <- bind_rows(map1, map2) %>%
  mutate(geno = paste0("kjs", geno)) %>%
  full_join(df) %>%
  # Create new cooler_side column
  mutate(cooler_side = case_when(tank_id %in% c("B2T2", "B2T4", "B1T2", "B1T4") ~ "left",
                                 tank_id %in% c("B2T1", "B2T3", "B1T1", "B1T3") ~ "right")) %>%
  # Recode cooler column position as number of columns from center of cooler
  mutate(cols_from_ctr = case_when(
    cooler_side == "right" ~ recode(Col, "4" = 4, "3" = 3, "2" = 2, "1" = 1),
    cooler_side == "left" ~ recode(Col, "4" = 1, "3" = 2, "2" = 3, "1" = 4))) %>%
  # Recode cooler row position as number of rows from center of cooler
  mutate(rows_from_ctr = case_when(Row == "A" ~ 3, Row == "B" ~ 2, Row == "C" ~ 1, 
                                   Row == "D" ~ 1, Row == "E" ~ 2, Row == "F" ~ 3))

# Visualize positional effects within each tank
ggplot(df_pos, aes(x = cols_from_ctr, y = fvfm, group = rows_from_ctr, 
                   color = factor(rows_from_ctr), shape = factor(rows_from_ctr))) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_grid(~target_temp)

# Fit model to adjust for positional effects
mod <- lm(fvfm ~ tank_id * cols_from_ctr * rows_from_ctr, data = df_pos, na.action = na.exclude)
anova(mod) %>% tidy() %>% knitr::kable()

# Calculate estimated marginal means for the "center" of each tank
emms <- emmeans(mod, specs  = c("tank_id", "cols_from_ctr", "rows_from_ctr")) %>%
  as_tibble() %>%
  select(tank_id, emmean)

# Adjust FvFm using partial residuals to account for positional effects within each tank
df_pos <- augment(mod, df_pos) %>%
  mutate(resid = .resid) %>%
  select(!starts_with("."))

df_pos <- left_join(emms, df_pos) %>%
  mutate(fvfm.adj = emmean + resid) %>%
  select(geno, f, fm, fvfm, emmean, resid, fvfm.adj, target_temp, max_temp, tank_id, cols_from_ctr, rows_from_ctr)

# Select whether to use adjusted or unadjusted fvfm data downstream
df <- if (pos.adjust) {
  # If position adjustment turned on, join with adjusted fvfm values
  left_join(df, df_pos) %>%
    # Select adjusted fvfm to use in modeling below
    select(nursery, geno, f, fm, fvfmraw, fvfm = fvfm.adj, target_temp, max_temp, problem)
} else { 
  # If position adjustment turned off, stick with original data
  df %>% select(nursery, geno, f, fm, fvfmraw, fvfm, target_temp, max_temp, problem)
}
```

# Dose-response curve model fitting
```{r model, results = 'hide'}
# Select whether to use target temp or actual recorded max temp for model fitting
df <- if (temp.adjust) {
  df %>% select(-target_temp)
} else {
  df %>% mutate(max_temp = target_temp) %>% select(-target_temp)
}

# Define function to fit 3-parameter LL model to data and return NULL if fitting error
ll3 <- function(data) {
  drm(fvfm ~ max_temp, data = data, 
      fct = LL.3(names = c("hill", "max", "ED50")),
      upperl = c(120, 0.72, 40),
      lowerl = c(10, 0.55, 30))}
tryll3 <- possibly(ll3, otherwise = NULL)

# Fit model to each coral, get parameters, fitted values, and residuals
initmods <- df %>%
  nest(data = c(max_temp, f, fm, fvfmraw, fvfm, problem)) %>%
  # Fit the model to each coral
  mutate(ll3 = map(data, tryll3)) %>%
  # Get model parameters and fitted values/residuals
  mutate(pars = map(ll3, tidy),
         pred = map2(ll3, data, ~augment(.x, drop_na(.y, fvfm))))

# Extract ed50 parameter values from model fits
ed50 <- initmods %>% 
  select(nursery, geno, pars) %>%
  unnest(pars) %>%
  filter(term == "ED50")

# Collect raw data, fitted values, and diagnostics
vals <- initmods %>%
  select(nursery, geno, pred) %>%
  unnest(pred) %>%
  full_join(ed50) %>%
  full_join(df) %>%
  rename(ed50 = estimate)

# Identify problematic data points based on cook's distance and residuals
counts <- vals %>% 
  group_by(nursery, geno) %>% 
  summarise(n = sum(!is.na(fvfm)))
dff <- vals %>%
  left_join(counts) %>%
  group_by(nursery, geno) %>%
  mutate(resid.thresh = 2*sd(.resid, na.rm = T)) %>%  # Calculate residual threshold as 2 standard deviations
  mutate(cooksd.thresh = 4/n) %>%   # Calculate cook's distance threshold as 4/n
  mutate(max_to_remove = floor(n * 0.2)) %>%
  ungroup() %>%
  mutate(problem = case_when(.cooksd > cooksd.thresh & .resid > 0 ~ "high cook's distance",
                             .resid > resid.thresh ~ "high residual", 
                             TRUE ~ problem)) %>%
  group_by(nursery, geno, outlier = problem %in% c("high cook's distance", "high residual")) %>%
  mutate(n.outliers = n(),
         rank.out = order(.cooksd, decreasing = TRUE)) %>%
  ungroup() %>%
  mutate(fvfm = case_when(outlier & rank.out <= max_to_remove ~ NA_real_, 
                          TRUE ~ fvfm)) 


# Refit models without problematic points
fmods <- dff %>%
  select(nursery, geno, max_temp, f, fm, fvfmraw, problem, fvfm) %>%
  nest(data = c(max_temp, f, fm, fvfmraw, fvfm, problem)) %>%
  # Fit the model to each coral
  mutate(ll3 = map(data, tryll3)) %>%
  # Get model parameters and fitted values/residuals
  mutate(pars = map(ll3, tidy),
         pred = map2(ll3, data, ~augment(.x, drop_na(.y, fvfm))))

# Extract ed50 parameter values from model fits
fed50 <- fmods %>% 
  select(nursery, geno, pars) %>%
  unnest(pars) %>%
  filter(term == "ED50")

# Collect raw data, fitted values, and ed50 estimates
fvals <- fmods %>%
  select(nursery, geno, pred) %>%
  unnest(pred) %>%
  full_join(fed50) %>%
  full_join(select(dff, nursery, geno, max_temp, f, fm, fvfmraw, problem, fvfm)) %>%
  rename(ed50 = estimate)
```

```{r plot, fig.width = 10, fig.height = 10, echo = FALSE}
# Define function to plot raw data, fitted values, and ed50 for each genotype
plotfits <- function(data) {
  ggplot(data = data, aes(x = max_temp)) + 
    geom_point(pch = 4, size = 1.25,
               aes(y = fvfmraw, color = factor(problem, levels = c("none", "no  signal",
                                              "abnormally high",  
                                              "abnormally high w/ low Ft",  
                                              "high residual", "high cook's distance")))) + 
    geom_point(aes(y = fvfm), pch = 1, size = 2) + 
    geom_line(data = drop_na(data, .fitted), 
              aes(y = .fitted)) + 
    geom_vline(data = distinct(data, nursery, geno, ed50),
               aes(xintercept = ed50), 
               lwd = 0.2, lty = 2) +
    geom_text(data = distinct(data, nursery, geno, ed50),
              aes(x = ed50, y = 0.05, label = round(ed50, 2)), 
              hjust = 1, nudge_x = -0.2, size = 3) +
    facet_wrap(~ nursery + geno, drop = TRUE) +
    scale_color_manual(name = "problem", drop = FALSE,
                       values = c("black", "yellow", "red", "orange", "blue", "turquoise"))
}
```

# Plot fits
```{r, fig.width = 10, fig.height = 10}
# Select whether to use model-based filtering or not
final <- if (mod.filter) { fvals } else { vals }

# Plot fits
plotfits(data = final)
```

# Summarize data filtering
```{r}
final %>%
  count(problem) %>%
  mutate(freq = n / sum(n)) %>%
  knitr::kable()

# Calculate pseudo R2 across dataset of the fitted models with given filtering options
finalmods <- if (mod.filter) { fmods } else { initmods }
calc.ss <- finalmods %>%
  mutate(aov = map(data, ~ aov(fvfm ~ factor(max_temp), data = .)),
         tss = map_dbl(aov, ~ sum(tidy(.)$sumsq)),
         drc.rss = map_dbl(ll3, ~ tidy(modelFit(.))$rss[2]),
         pseudo.R2 = 1 - (drc.rss / tss))

tibble(pseudo_R2 = 1 - (sum(calc.ss$drc.rss) / sum(calc.ss$tss))) %>%
  knitr::kable()
```

# Save processed data to file
```{r save}
# Write tidy data to file
final %>%
  select(nursery, geno, max_temp, fvfm) %>%
  drop_na() %>%
  mutate(date = "2020-08") %>%
  write_csv("data/processed/KJS_fvfm_adj_clean.csv")
```
