```{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(joeyr)
library(emmeans)
library(ggpubr)
```

# CBASS temperature logs

#### Import temperature logs for Arduino-controlled tanks
```{r}
readlog <- function(logfile, prefix) {
  read_csv(logfile) %>%
    # Remove internal header rows
    filter(PrintDate != "PrintDate") %>%
    # Format date and time
    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 to long format
    pivot_longer(starts_with("T"), names_to = "key", values_to = "temp") %>%
    # Remove rows where temp sensors gave error codes
    filter(temp > 0, temp < 50) %>%
    # Link arduino box and sensor positions (e.g., B1T2)
    mutate(tank = str_extract(key, "T[0-9]"),
           tank = paste0(prefix, tank),
           key = case_when(grepl("SP", key) ~ str_sub(key, 3, 4),
                           TRUE ~ str_sub(key, 1, 4)),
           key = tolower(key)) %>%
    # Drop rows that did not parse a date/time
    drop_na(dttm, temp) %>%
    # Create columns for set point and actual temperature
    pivot_wider(names_from = key, values_from = temp) %>%
    # Tidy column types
    mutate(tank = factor(tank), sp = as.numeric(sp), temp = as.numeric(temp)) %>%
    # Add maximum setpoint temperature as max_temp column
    group_by(tank) %>%
    mutate(max_temp = max(sp, na.rm = TRUE)) %>%
    ungroup()
}

# Read in arduino log files
b1log <- readlog("data/raw/20201014_UM_Acer/temperature/sdlog/20201014_BOX1_LOG.TXT", prefix = "B1")
b3log <- readlog("data/raw/20201014_UM_Acer/temperature/sdlog/20201014_BOX3_LOG.TXT", prefix = "B3")

# Tidy arduino log files for this run
arduinolog <- bind_rows(b1log, b3log) %>%
  # Filter to just the date of this run
  filter(date(dttm) == "2020-10-14") %>%
  # Specify which arduino boxes controlled port vs. starboard tank sets
  mutate(tank_set = if_else(str_detect(tank, "B1"), "starboard", "port")) %>%
  # Add tank_id column from tank design
  left_join(read_csv("data/tank_setup.csv")) %>%
  select(dttm, hm, sp, temp, max_temp, tank_set, tank_id)
```

#### Import HOBO temperature logs for Inkbird-controlled tanks
```{r}
# Read in HOBO temperature logs
hobologfiles <- list.files("data/raw/20201014_UM_Acer/temperature/hobo", pattern = ".xlsx", full.names = TRUE)
hobolog <- tibble(filename = hobologfiles) %>%
  mutate(log = map(filename, read_xlsx, skip = 1)) %>%
  unnest()

# Tidy HOBO temperature logs
hobolog <- hobolog %>%
  clean_names() %>%
  # Get tank_id from filename
  mutate(tank_id = str_sub(basename(filename), 1, 2)) %>%
  # Convert temps from F to C
  mutate(temp = (temp_f - 32) * 5/9) %>%
  # Select columns to keep
  select(dttm = date_time_gmt_0400, temp, tank_id) %>%
  mutate(hm = format(dttm, "%H:%M")) %>%
  # Join with tank_id information
  left_join(read_csv("data/tank_setup.csv")) %>%
  # Filter only target date
  filter(date(dttm) == "2020-10-14")

# Import inkbird setpoints
inkbird_sp <- read_csv("data/program_files/inkbird_settings.csv") %>%
  pivot_longer(2:5, names_to = "max_temp", values_to = "sp") %>%
  mutate(max_temp = as.numeric(str_extract(max_temp, "[0-9]+")),
         dttm = ymd_hms(paste("2020-10-14", time)), hm = format(dttm, "%H:%M"),
         temp = NA) %>%
  left_join(read_csv("data/tank_setup.csv")) %>%
  select(dttm, hm, sp, temp, max_temp, tank_set, tank_id)

hobolog <- bind_rows(inkbird_sp, hobolog)
```

### Plot CBASS temperature profiles and set points
```{r, fig.width = 10}
# Combine arduino and hobo temperature log data
templog <- bind_rows(arduinolog, hobolog)

# Import hobo offset data
hobocal <- read_csv("data/processed/hobo_calibrations.csv")
templog <- left_join(templog, hobocal) %>%
  mutate(temp.adj = case_when(!is.na(hobo) ~ temp - mean_offset, TRUE ~  temp)) %>%
  select(dttm, hm, sp, temp = temp.adj, max_temp, tank_set, tank_id)

# Plot
ggplot(templog, mapping = aes(x = dttm, y = temp, group = tank_id)) +
  geom_line(aes(color = tank_id), lwd = 0.2) +
  geom_step(aes(y = sp, color = tank_id), lwd = 0.5) +
  facet_wrap(~ tank_set, scales = "free") +
  scale_y_continuous(breaks = 30:38, limits = c(29, 39)) +
  scale_x_datetime(breaks = "hours", labels = label_date("%H:%M"),
                   limits = as.POSIXct(c("2020-10-14 12:00:00", "2020-10-14 21:00:00"), tz = "UTC")) +
  theme_hc() +
  theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "Time of day", y = "Temperature (°C)")

#write_csv(templog, "data/20201014_UM_Acer/temperature/UM_tempdata.csv")

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

# 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
```

#### Import data from DIVING PAM II
```{r}
# Read manually entered data from Diving PAM II
path <- "data/raw/20201014_UM_Acer/PAM/20201014dii_Acer.xlsx"
dpii <- path %>%
  excel_sheets() %>%
  set_names() %>%
  discard(.p = ~str_detect(.x, "Acer")) %>%
  map_df(~ read_excel(path = path, sheet = .x, range = "A1:H43", col_types = "text"), .id = "tank_id") %>%
  mutate(f = as.numeric(`F`), fm = as.numeric(`Fm'`), fvfm = as.numeric(`Y (II)`)) %>%
  select(tank_id, Geno, f, fm, fvfm) %>%
  clean_names() %>%
  mutate(geno = str_replace(geno, "-", ""),
         geno = str_replace(geno, " ", "")) %>%
  drop_na(geno)

# Combine with tank_id metadata
tanks <- read_csv("data/tank_setup.csv")
dpii <- full_join(dpii, tanks) %>%
  select(geno, f, fm, fvfm, target_temp = max_temp, tank_id, tank_set, cooler_side)

# Fix error in data -- genotype "8yy" should be "88y"
dpii$geno[dpii$geno == "8yy"] <- "88y"

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

# Initial data filtering and QC
```{r}
dpii <- dpii %>%
  mutate(nursery = "um") %>%
  select(nursery, geno, f, fm, fvfm, target_temp, max_temp, tank_id, tank_set, cooler_side)

# 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 > median(mahal) & 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))

# Filter fvfm values based on various outlier classifications
df <- dpii %>%
  # Save raw fvfm in new column
  mutate(fvfmraw = fvfm) %>%          
  # Merge with mahalanobis outlier info from above
  left_join(out) %>%
  # Identify problematic data points
  mutate(problem = case_when(
    is.na(fvfm) ~ "no signal",
    fvfm > 0.750 ~ "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") )
}
```

# Positional effects
```{r, fig.width = 10}
# Import maps of the position of each genotype in each tank
Pmap <- read_csv("data/raw/20201014_UM_Acer/P_map.csv") %>%
  pivot_longer(-Row, values_to = "geno", names_to = "Col") %>%
  drop_na() %>%
  unite(pos, Row, Col, sep = "", remove = FALSE)
Smap <- read_csv("data/raw/20201014_UM_Acer/S_map.csv") %>%
  pivot_longer(-Row, values_to = "geno", names_to = "Col") %>%
  drop_na() %>%
  unite(pos, Row, Col, sep = "", remove = FALSE)

df_pos <- bind_rows(Pmap, Smap) %>%
  full_join(df) %>%
  # Create new cooler_id column
  mutate(cooler = case_when(tank_id %in% c("P1", "P2") ~ "p_cooler_1",
                            tank_id %in% c("P3", "P4") ~ "p_cooler_2",
                            tank_id %in% c("P5", "P6") ~ "p_cooler_3",
                            tank_id %in% c("P7", "P8") ~ "p_cooler_4",
                            tank_id %in% c("S1", "S2") ~ "s_cooler_1",
                            tank_id %in% c("S3", "S4") ~ "s_cooler_2",
                            tank_id %in% c("S5", "S6") ~ "s_cooler_3",
                            tank_id %in% c("S7", "S8") ~ "s_cooler_4")) %>%
  # Recode cooler column position as number of columns from center of cooler
  mutate(cols_from_ctr = case_when(
    cooler_side == "right" ~ recode(Col, "5" = 4, "4" = 4, "3" = 3, "2" = 2, "1" = 1),
    cooler_side == "left" ~ recode(Col, "5" = 1, "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(tank_set~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, tank_set, cooler, 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-10") %>%
  write_csv("data/processed/um_fvfm_adj_clean.csv")
```

