############################################
# PBZ × LIGHT EXPERIMENT IN LIME
# COMPLETE FRBD ANALYSIS
############################################

#--------------------------------------------------
# 1. LOAD PACKAGES
#--------------------------------------------------
library(tidyverse)
library(lme4)
library(lmerTest)
library(emmeans)
library(multcomp)
library(openxlsx)
library(car)
library(ggplot2)

#--------------------------------------------------
# 2. IMPORT DATA
#--------------------------------------------------
year1 <- read.csv("Lime_Year1_Data.csv")
year2 <- read.csv("Lime_Year2_Data.csv")

#--------------------------------------------------
# 3. DATA PREPARATION
#--------------------------------------------------
prep_fun <- function(df, yr){
  df %>%
    mutate(
      Replication = factor(Replication),
      PBZ = factor(PBZ),
      Light = factor(Light),
      Year = factor(yr)
    )
}

year1 <- prep_fun(year1, "Y1")
year2 <- prep_fun(year2, "Y2")

combined <- bind_rows(year1, year2)

#--------------------------------------------------
# 4. RESPONSE VARIABLES
#--------------------------------------------------
responses <- c(
  "NoFruits","FruitWt_g","FruitSet_pct","JuiceContent_pct",
  "ShootLen_cm","NoShoots","NoLeaves","LeafArea_cm2",
  "SPAD","DaysToFlower","NoFlowers","TSS_pct","Acidity_pct"
)

#--------------------------------------------------
# 5. OUTPUT FOLDERS
#--------------------------------------------------
dir.create("Plots", showWarnings = FALSE)
dir.create("Diagnostics", showWarnings = FALSE)

#--------------------------------------------------
# 6. EXCEL WORKBOOK
#--------------------------------------------------
wb <- createWorkbook()

#--------------------------------------------------
# 7. LOOP THROUGH TRAITS
#--------------------------------------------------
for (response in responses) {
  
  cat("\nAnalyzing:", response, "\n")
  
  #--------------------------------------------
  # YEAR-WISE ANOVA
  #--------------------------------------------
  mod_y1 <- aov(as.formula(paste(response, "~ Replication + PBZ*Light")), data = year1)
  mod_y2 <- aov(as.formula(paste(response, "~ Replication + PBZ*Light")), data = year2)
  
  tab_y1 <- anova(mod_y1)
  tab_y2 <- anova(mod_y2)
  
  #--------------------------------------------
  # MIXED MODEL (Year + Replication random)
  #--------------------------------------------
  model <- lmer(
    as.formula(paste(response,
                     "~ PBZ*Light + (1|Year) + (1|Year:Replication)")),
    data = combined
  )
  
  tab_pool <- Anova(model, type = 3)
  
  #--------------------------------------------
  # ESTIMATED MEANS
  #--------------------------------------------
  emm <- emmeans(model, ~ PBZ * Light)
  
  #--------------------------------------------
  # MEAN SEPARATION (DESCENDING ORDER)
  #--------------------------------------------
  tukey <- cld(
    emm,
    adjust = "tukey",
    Letters = letters,
    sort = TRUE,
    decreasing = TRUE
  )
  
  lsd <- cld(
    emm,
    adjust = "none",
    Letters = letters,
    sort = TRUE,
    decreasing = TRUE
  )
  
  tukey_df <- as.data.frame(tukey) %>% arrange(desc(emmean))
  lsd_df   <- as.data.frame(lsd)   %>% arrange(desc(emmean))
  
  #--------------------------------------------
  # INTERACTION PLOT
  #--------------------------------------------
  plot_df <- as.data.frame(emm)
  
  p <- ggplot(plot_df,
              aes(PBZ, emmean, color = Light, group = Light)) +
    geom_point(size = 3) +
    geom_line() +
    labs(title = paste("PBZ × Light Interaction -", response),
         y = response) +
    theme_bw()
  
  ggsave(paste0("Plots/", response, "_interaction.png"),
         p, width = 6, height = 5)
  
  #--------------------------------------------
  # DIAGNOSTICS
  #--------------------------------------------
  png(paste0("Diagnostics/", response, "_diagnostics.png"), 800, 800)
  par(mfrow=c(2,2))
  plot(model)
  qqnorm(residuals(model)); qqline(residuals(model))
  hist(residuals(model))
  dev.off()
  
  #--------------------------------------------
  # EXPORT TO EXCEL
  #--------------------------------------------
  addWorksheet(wb, paste0(response, "_Y1_ANOVA"))
  writeData(wb, paste0(response, "_Y1_ANOVA"), tab_y1)
  
  addWorksheet(wb, paste0(response, "_Y2_ANOVA"))
  writeData(wb, paste0(response, "_Y2_ANOVA"), tab_y2)
  
  addWorksheet(wb, paste0(response, "_Pooled_ANOVA"))
  writeData(wb, paste0(response, "_Pooled_ANOVA"), as.data.frame(tab_pool))
  
  addWorksheet(wb, paste0(response, "_Tukey"))
  writeData(wb, paste0(response, "_Tukey"), tukey_df)
  
  addWorksheet(wb, paste0(response, "_LSD"))
  writeData(wb, paste0(response, "_LSD"), lsd_df)
}

#--------------------------------------------------
# 8. SAVE EXCEL FILE
#--------------------------------------------------
saveWorkbook(wb, "PBZ_Light_All_Results.xlsx", overwrite = TRUE)

cat("\nAnalysis complete ✔ Highest mean now labeled 'a'\n")