############################################################
# ODER~SO | Oder River Disaster 2022 | Data | SP08a
# Daphnia toxicity analysis: strain × concentration × time

# This script:
# 1) Imports the raw dataset (CSV)
# 2) Reshapes it into long format (24/48/72 h)
# 3) Fits a quasibinomial GLM (logit link)
# 4) Performs Type II ANOVA
# 5) Computes Tukey-adjusted post hoc comparisons of strains
#    to assess differences in mortality among strains
#    at each time point (24, 48, and 72 h)
############################################################
# Clear workspace
rm(list = ls())
# Required packages
library(car)       # Type II ANOVA
library(emmeans)   # Estimated marginal means and post hoc tests
############################################################
# 1) Load data
############################################################
# The CSV file should be located in the same directory as this script
Mortality24 <- read.csv(
  "Daphnia_Toxicity_Strains_Density_20241212.csv",
  stringsAsFactors = FALSE
)
# Inspect data structure
str(Mortality24)
############################################################
# 2) Convert to long format (24, 48, 72 h)
############################################################
Mortality24_long <- data.frame(
  percent_dead = c(
    Mortality24$percent_dead_24h,
    Mortality24$percent_dead_48h,
    Mortality24$percent_dead_72h
  ),
  concentration = rep(Mortality24$concentration, 3),
  strain = rep(Mortality24$strain, 3),
  time = factor(rep(c(24, 48, 72), each = nrow(Mortality24)))
)
# Inspect reshaped data
str(Mortality24_long)

############################################################
# 3) Fit quasibinomial GLM (logit scale)
############################################################
glm_model <- glm(
  percent_dead ~ strain * concentration * time,
  family = quasibinomial(link = "logit"),
  data = Mortality24_long)
summary(glm_model)
############################################################
# 4) Type II ANOVA
############################################################
Anova(glm_model, type = "II")
############################################################
# 5) Post hoc comparisons:
#    Pairwise strain comparisons within each time point (24, 48, and 72 h),
#    Tukey-adjusted
############################################################
emm_strain_time <- emmeans(glm_model, ~ strain | time)
strain_comparisons <- pairs(
  emm_strain_time,
  adjust = "tukey")
summary(strain_comparisons)
############################################################
# End of script
############################################################
