library(dplyr)
#install.packages(c("purrr", "tibble"))
library(purrr)
library(tidyr)
library(tibble)
#install.packages(c("FSA", "rstatix"))
library(FSA) 
library(rstatix)
library(car)
setwd("C:/Users/PAULA/Desktop/Datos doctorado/inercia/Archivos para paper")

datos <- read.csv("datosANOVAtipoIII.csv")
ggplot(datos, aes(x = Diffusivity)) +
  geom_histogram(bins = 30, fill = "skyblue", color = "black") +
  theme_bw()

ggplot(datos, aes(x = Season, y = Diffusivity, fill = Sp)) +
  geom_boxplot() +
  facet_wrap(~Depth) +
  theme_bw()

# --- 4. Ajustar modelo GLM (Gamma con enlace log)
modelo_glm <- glm(Diffusivity ~ Sp * Depth * Season,
                  family = Gamma(link = "log"),
                  data = datos)

anova_glm <- Anova(modelo_glm, type = 3)   # guarda el resultado del ANOVA
summary_glm <- summary(modelo_glm)         # guarda el resumen del modelo

# Convertir ambos a data.frame
anova_table <- as.data.frame(anova_glm)
summary_table <- as.data.frame(summary_glm$coefficients)

# Exportar como CSV
write.csv(anova_table, "ANOVA_TipoIII_GLM.csv", row.names = TRUE)
write.csv(summary_table, "Resumen_GLM_Coefficients.csv", row.names = TRUE)
# --- 7. Diagnóstico del modelo
par(mfrow = c(2, 2))
plot(modelo_glm)

# Parámetro de dispersión (cercano a 1 = buen ajuste)
dispersion <- sum(residuals(modelo_glm, type = "pearson")^2) / modelo_glm$df.residual
print(paste("Parámetro de dispersión:", round(dispersion, 3)))

# --- 8. Comparaciones post-hoc (si hay efectos significativos)
emmeans(modelo_glm, pairwise ~ Sp | Season * Depth)
# --- 11. Exportar comparaciones post-hoc (EMMEANS)

# Calcular comparaciones múltiples
posthoc <- emmeans(modelo_glm, pairwise ~ Sp | Season * Depth)

# Convertir resultados a data frame
# Las comparaciones están en posthoc$contrasts
posthoc_table <- as.data.frame(posthoc$contrasts)

# Exportar CSV
write.csv(posthoc_table, "PostHoc_Comparisons.csv", row.names = FALSE)

# Si querés también redondear los valores:
posthoc_table_round <- posthoc_table
posthoc_table_round[, sapply(posthoc_table_round, is.numeric)] <- 
  round(posthoc_table_round[, sapply(posthoc_table_round, is.numeric)], 3)

write.csv(posthoc_table_round, "PostHoc_Comparisons_redondeado.csv", row.names = FALSE)



# --- 9. Visualización de resultados (opcional)
ggplot(datos, aes(x = Season, y = Diffusivity, fill = Sp)) +
  geom_boxplot() +
  facet_wrap(~Depth) +
  theme_bw() +
  labs(title = "Comparación de Diffusivity por especie, estación y profundidad",
       y = "Diffusivity", x = "Season")