################################################################################
#  Daniela Milano, Sonia A. Crichigno, Pablo E. Reggi, Miguel A. Battini, 
#  Leonardo M. Buria & Juan P. Barriga. Temperature-driven differences in 
#  feeding performance of native and non-native fish of Patagonia: a functional
#  response analysis. Freshwater Biology
################################################################################

# Following datasets are required:

# Galaxias platei simple treatment data (Gplatei_simple_GP1.csv)

# Oncorhynchus mykiss simple treatment data (Omykiss_simple_OM1.csv)

# Galaxias platei-Galaxias platei in conspecific treatment data 
# (Gplatei_cospecific_GP2.csv)

# Galaxias platei-Oncorhynchus mykiss in heterospecific treatment data 
#(Gplatei_Omykiss_heterospecific_GPOM.csv)

# Galaxias platei in heterospecific treatment data 
# (Gplatei_heterospecific_GP_GPOM.csv)

# Oncorhynchus mykiss in heterospecific treatment data 
#(Omykiss_heterospecific_OM_GPOM.csv)


###############################################################
# Variable description:
# den: prey density offered
# T13: eaten prey at 13 °C 
# T16= eaten prey at 16 °C 
# T19= eaten prey at 19 °C 
# T21= eaten prey at 21 °C

##############################################################
# Functional Response Analysis Workflow
# Steps: (1) model selection; (2) model fitting; 
#        (3) bootstrap confidence intervals; 
#        (4) graphical comparison of FR curves; 
#        (5) comparison of FR parameters and FRR between curves
##############################################################

# --- Packages ---
library(frair)
library(RColorBrewer)
library(dplyr)
library(tidyr)
library(purrr)

##############################################################
# (1) Model selection
##############################################################
# Example: simple treatment (Gp1, Om1) at 13–21 °C

# G. platei
frair_test(T13 ~ den, data = Gp1)
frair_test(T16 ~ den, data = Gp1)
frair_test(T19 ~ den, data = Gp1)
frair_test(T21 ~ den, data = Gp1)

# O. mykiss
frair_test(T13 ~ den, data = Om1)
frair_test(T16 ~ den, data = Om1)
frair_test(T19 ~ den, data = Om1)
frair_test(T21 ~ den, data = Om1)

##############################################################
# (2) Model fitting (Rogers’ type II)
##############################################################
# G. platei
GP13 <- frair_fit(T13 ~ den, data=Gp1, response="rogersII",
                  start=list(a=1.2, h=0.015), fixed=list(T=1))
GP16 <- frair_fit(T16 ~ den, data=Gp1, response="rogersII",
                  start=list(a=1.2, h=0.015), fixed=list(T=1))
GP19 <- frair_fit(T19 ~ den, data=Gp1, response="rogersII",
                  start=list(a=1.2, h=0.015), fixed=list(T=1))
GP21 <- frair_fit(T21 ~ den, data=Gp1, response="rogersII",
                  start=list(a=1.2, h=0.015), fixed=list(T=1))

# O. mykiss
OM13 <- frair_fit(T13 ~ den, data=Om1, response="rogersII",
                  start=list(a=1.2, h=0.015), fixed=list(T=1))
OM16 <- frair_fit(T16 ~ den, data=Om1, response="rogersII",
                  start=list(a=1.2, h=0.015), fixed=list(T=1))
OM19 <- frair_fit(T19 ~ den, data=Om1, response="rogersII",
                  start=list(a=1.2, h=0.015), fixed=list(T=1))
OM21 <- frair_fit(T21 ~ den, data=Om1, response="rogersII",
                  start=list(a=1.2, h=0.015), fixed=list(T=1))

##############################################################
# (3) Bootstrap confidence intervals (3000 iterations)
##############################################################
GP13boot <- frair_boot(GP13, nboot=3000)
GP16boot <- frair_boot(GP16, nboot=3000)
GP19boot <- frair_boot(GP19, nboot=3000)
GP21boot <- frair_boot(GP21, nboot=3000)

OM13boot <- frair_boot(OM13, nboot=3000)
OM16boot <- frair_boot(OM16, nboot=3000)
OM19boot <- frair_boot(OM19, nboot=3000)
OM21boot <- frair_boot(OM21, nboot=3000)

##############################################################
# (4) Graphical comparison of FR curves
##############################################################
# Example for GP13:
plot(GP13boot, xlim=c(0,65), type="n", main="Empirical 95% CI")
drawpoly(GP13boot, col=adjustcolor("blue", alpha=0.1))
points(GP13boot, pch=20, col=rgb(0,0,0,0.2))
lines(GP13)

# Repeat as needed for other datasets

##############################################################
# (4.5) Create unified bootstrap dataset for parameter contrasts
##############################################################
# Function to extract a/h from bootcoefs
extract_boot <- function(obj, sp, tem){
  df <- as.data.frame(obj$bootcoefs)
  df$FRR <- df$a / df$h
  df$sp <- sp
  df$tem <- tem
  return(df)
}

# Combine all bootstraps
Simple <- bind_rows(
  extract_boot(GP13boot,"Gp",13),
  extract_boot(GP16boot,"Gp",16),
  extract_boot(GP19boot,"Gp",19),
  extract_boot(GP21boot,"Gp",21),
  extract_boot(OM13boot,"Om",13),
  extract_boot(OM16boot,"Om",16),
  extract_boot(OM19boot,"Om",19),
  extract_boot(OM21boot,"Om",21)
)

# Export interations
write.csv(Simple, "Bootstraps_all.csv", row.names = FALSE)

# Generate a new dataset with the following variables:
# iter: number of iteration
# a: attack rate
# h: handling time
# FRR: functional response ratio (calculated as a/h)
# tem: experimental temperature
# sp: species (Gp or Om)

# Note that in conspecific or heterospecific treatments instead of “sp” should
# be used type of curve as “Observed” or “Expected”.

##############################################################
# (5) Comparisons of parameters and FRR
##############################################################
set.seed(123)

Simple_long <- Simple %>%
  mutate(
    sp  = factor(sp, levels = c("Gp","Om")),
    tem = factor(tem, levels = sort(unique(tem)))
  ) %>%
  pivot_longer(c(a,h,FRR), names_to="param", values_to="value")

chk <- Simple_long %>% count(sp, tem, param)
stopifnot(all(chk$n > 0))

n_resample <- 3000

# --- Species contrasts (Om − Gp) per temperature ---
species_res <- map_dfr(1:n_resample, function(i){
  samp <- Simple_long %>%
    group_by(sp, tem, param) %>%
    slice_sample(n=1) %>% ungroup()
  
  samp %>%
    pivot_wider(names_from=sp, values_from=value) %>%
    mutate(diff = Om - Gp, iter=i) %>%
    select(param, tem, diff, iter)
})

species_results <- species_res %>%
  group_by(param, tem) %>%
  summarise(
    est = mean(diff), med = median(diff),
    lwr = quantile(diff,0.025), upr = quantile(diff,0.975),
    p_gt0 = mean(diff > 0),
    p_two = 2 * pmin(p_gt0, 1 - p_gt0),
    .groups="drop"
  ) %>%
  arrange(param, tem)

print(species_results, n=Inf)

# --- Temperature contrasts (t2 − t1) within each species ---
temp_pairs <- combn(levels(Simple_long$tem), 2, simplify=FALSE)

temp_res <- map_dfr(1:n_resample, function(i){
  samp <- Simple_long %>%
    group_by(sp, tem, param) %>%
    slice_sample(n=1) %>% ungroup()
  
  map_dfr(unique(samp$param), function(p){
    map_dfr(unique(samp$sp), function(s){
      v <- filter(samp, param==p, sp==s) %>%
        pull(value) %>% setNames(as.character(unique(samp$tem)))
      map_dfr(temp_pairs, function(pr){
        tibble(param=p, sp=s, t1=pr[1], t2=pr[2],
               diff=v[pr[2]] - v[pr[1]], iter=i)
      })
    })
  })
})

temp_contrasts <- temp_res %>%
  group_by(param, sp, t1, t2) %>%
  summarise(
    est = mean(diff), med = median(diff),
    lwr = quantile(diff,0.025), upr = quantile(diff,0.975),
    p_gt0 = mean(diff > 0),
    p_two = 2 * pmin(p_gt0, 1 - p_gt0),
    .groups="drop"
  ) %>%
  mutate(contrast = paste0(t2," - ",t1)) %>%
  select(param, sp, contrast, est, med, lwr, upr, p_gt0, p_two) %>%
  arrange(param, sp, contrast)

print(temp_contrasts, n=Inf)