# Load required packages
library(dplyr)
library(lme4)
library(car)
library(multcomp)


## Experiment 1: Yellow-phase eel ## 

# Ensure factors are correctly specified
YEInitial <- YEInitial %>%
  mutate(
    Side = as.factor(Side),
    FS   = as.factor(FS)
  )

YEOverall <- YEOverall %>%
  mutate(
    Channel   = as.factor(Channel),
    FS        = as.factor(FS),
    Direction = as.factor(Direction),
    EelID     = as.factor(EelID)
  )

# Initial passage: Fisher's Exact test # 

dat <- matrix(c(9,7,4,7), nrow = 2,
              dimnames = list(Treatment = c("Control","Low"),
                              Channel   = c("E","N")))

# Check expected counts
chisq.test(dat)$expected

# Fisher's Exact Test (used if expected < 5)
fisher_res <- fisher.test(dat)
fisher_res

# Initial Avoidance # 
mod_initial_full<-glm(Any~Side +FS,family=binomial, data=YEInitial)
summary (mod_initial_full)

mod_initial_side <- glm(Any ~ Side, family = binomial, data = YEInitial)
summary(mod_initial_side)

mod_initial_fs <- glm(Any ~ FS, family = binomial, data = YEInitial)
summary(mod_initial_fs)


#Total Avoidance# 

mod_route <- glmer(RouteChange ~ Channel + FS + Direction + (1|EelID),
                   family = binomial, data = YEOverall)
summary(mod_route)
Anova(mod_route)

mod_reject <- glmer(Reject ~ Channel + FS + Direction + (1|EelID),
                    family = binomial, data = YEOverall)
summary(mod_reject)
Anova(mod_reject)

mod_react <- glmer(React ~ Channel + FS + Direction + (1|EelID),
                   family = binomial, data = YEOverall)
summary(mod_react)
Anova(mod_react)

mod_nochange <- glmer(NoChange ~ Channel + FS + Direction + (1|EelID),
                      family = binomial, data = YEOverall)
summary(mod_nochange)
Anova(mod_nochange)


# Split channels (EC/NEC) #

Electrified <- YEOverall %>% filter(Channel == "E")
NonElectric <- YEOverall %>% filter(Channel == "N")

# Electrified Channel (EC) #
mod_EC_route <- glmer(RouteChange ~ FS + (1|EelID),
                     family = binomial, data = Electrified)
Anova(mod_EC_route)

mod_EC_reject <- glmer(Reject ~ FS + (1|EelID),
                      family = binomial, data = Electrified)
Anova(mod_EC_reject)

mod_EC_react <- glmer(React ~ FS + (1|EelID),
                     family = binomial, data = Electrified)
Anova(mod_EC_react)

mod_EC_nochange <- glmer(NoChange ~ FS + (1|EelID),
                        family = binomial, data = Electrified)

Anova(mod_EC_nochange)

# Post hoc (Tukey)
posthoc_EC_reject <- glht(mod_EC_reject, linfct = mcp(FS = "Tukey"))
summary(posthoc_EC_reject)

posthoc_EC_react <- glht(mod_EC_react, linfct = mcp(FS = "Tukey"))
summary(posthoc_EC_react)

posthoc_EC_nochange <- glht(mod_EC_nochange, linfct = mcp(FS = "Tukey"))
summary(posthoc_EC_nochange)

# Non-Electrified Channel (NEC) # 

mod_NEC_route <- glmer(RouteChange ~ FS + (1|EelID),
                     family = binomial, data = NonElectric)
Anova(mod_NEC_route)

mod_NEC_reject <- glmer(Reject ~ FS + (1|EelID),
                      family = binomial, data = NonElectric)
Anova(mod_NEC_reject)

mod_NEC_react <- glmer(React ~ FS + (1|EelID),
                     family = binomial, data = NonElectric)
Anova(mod_NEC_react)

mod_NEC_nochange <- glmer(NoChange ~ FS + (1|EelID),
                        family = binomial, data = NonElectric)
Anova(mod_NEC_nochange)

# Direction of approach# 

fs_levels <- c("0.28", "0.37", "0.66")

Treat_EC <- Electrified %>% filter(FS %in% fs_levels)
Treat_NEC <- NonElectric %>% filter(FS %in% fs_levels)

### Electrified

mod_EC_dir_react <- glmer(React ~ Direction + (1|EelID),
                         family = binomial, data = Treat_EC)
summary(mod_EC_dir_react); Anova(mod_EC_dir_react)

mod_EC_dir_route <- glmer(RouteChange ~ Direction + (1|EelID),
                         family = binomial, data = Treat_EC)
summary(mod_EC_dir_route) 
Anova(mod_EC_dir_route)

mod_EC_dir_reject <- glmer(Reject ~ Direction + (1|EelID),
                          family = binomial, data = Treat_EC)
summary(mod_EC_dir_reject) 
Anova(mod_EC_dir_reject)

mod_EC_dir_nochange <- glmer(NoChange ~ Direction + (1|EelID),
                            family = binomial, data = Treat_EC)
summary(mod_EC_dir_nochange)
Anova(mod_EC_dir_nochange)

### Non-electrified

mod_NEC_dir_react <- glmer(React ~ Direction + (1|EelID),
                         family = binomial, data = Treat_NEC)
summary(mod_NEC_dir_react)
Anova(mod_NEC_dir_react)

mod_NEC_dir_route <- glmer(RouteChange ~ Direction + (1|EelID),
                         family = binomial, data = Treat_NEC)
summary(mod_NEC_dir_route)
Anova(mod_NEC_dir_route)

mod_NEC_dir_reject <- glmer(Reject ~ Direction + (1|EelID),
                          family = binomial, data = Treat_NEC)
summary(mod_NEC_dir_reject)
Anova(mod_NEC_dir_reject)

mod_NEC_dir_nochange <- glmer(NoChange ~ Direction + (1|EelID),
                            family = binomial, data = Treat_NEC)

summary(mod_NEC_dir_nochange)
Anova(mod_NEC_dir_nochange)


#Total channel passage# 

mod_pass_full <- glmer(Pass ~ Channel + FS + Direction + (1|EelID),
                       family = binomial, data = YEOverall)

summary(mod_pass_full)
Anova(mod_pass_full)

#######################################################

###### Experiment 2: Silver-phase eel #######
# Ensure correct data types
SEInitial <- SEInitial %>%
  mutate(
    Side      = as.factor(Side),
    FS        = as.factor(FS),
    Frequency = as.factor(Frequency)
  )

# Initial Passage: Goodness-of-fit Test #
#Observed counts
control <- c(8, 12)

# Expected probabilities must sum to 1
expected_prob <- c(0.35, 0.65)

chisq_pass <- chisq.test(control, p = expected_prob)
chisq_pass

# Initial Avoidance # 
mod_initial <- glm(Any ~ Side + FS + Frequency,
                family = binomial, data = SEInitial)

summary(mod_initial)
Anova(mod_initial)

# Split channel (Ec/NEC) #

Electrified <- SEInitial %>% filter(Side == "E")
NonElectric <- SEInitial %>% filter(Side == "N")

# Electrified Channel Analysis (EC) #

mod_EC_FS <- glm(Any ~ FS, family = binomial, data = Electrified)
summary(mod_EC_FS)
Anova(mod_EC_FS)

posthoc_EC_FS <- glht(mod_EC_FS, linfct = mcp(FS = "Tukey"))
summary(posthoc_EC_FS)


fs_levels <- c("0.18", "0.3")

Treat_EC <- Electrified %>% filter(FS %in% fs_levels)

mod_EC_freq <- glm(Any ~ Frequency, family = binomial, data = Treat_EC)
summary(mod_EC_freq)
Anova(mod_EC_freq)

# Non-Electrified Analysis (NEC) # 

mod_NEC_FS <- glm(Any ~ FS, family = binomial, data = NonElectric)
summary(mod_NEC_FS)
Anova(mod_NEC_FS)

Treat_NEC <- NonElectric %>% filter(FS %in% fs_levels)

mod_NEC_freq <- glm(Any ~ Frequency, family = binomial, data = Treat_NEC)
summary(mod_NEC_freq)
Anova(mod_NEC_freq)