# -----------------------------
# FASE 1. Generar dataset simulado (abundancias por parcela)
# -----------------------------

etapas <- c("Pionera", "Intermedia", "Tardia", "Madura") # Define 4 etapas sucesionales (categorías)

parcelas <- expand.grid(Etapa = etapas, Parcela = 1:3) %>% # Crea combinaciones: 4 etapas x 3 parcelas
 as_tibble() %>% # Convierte a tibble
 mutate(ID = paste0(Etapa, "_P", Parcela))  # Crea un ID único por parcela (ej. Pionera_P1)

especies <- paste0("Sp", sprintf("%02d", 1:15))  # Crea nombres Sp01...Sp15 (15 especies)

dirichlet_probs <- function(alpha){  # Define función para generar probabilidades tipo Dirichlet
 w <- rgamma(length(alpha), shape = alpha, rate = 1) # Genera pesos positivos con distribución gamma
 w / sum(w) # Normaliza para que sumen 1 (probabilidades)
} # Fin de función

alpha_por_etapa <- list( # Lista con "alpha" por etapa (controla dominancia)
 Pionera  = c(6,6,5, 2,2, 1,1,1,1,1, 0.5,0.5,0.5,0.5,0.5), # Pionera: pocas especies dominan mucho
 Intermedia = c(3,3,3, 3,3, 2,2,2,1.5,1.5, 1,1,1,1,1), # Intermedia: dominancia más equilibrada
 Tardia = c(1.5,1.5,1.5, 2,2, 2.5,2.5,2.5,2,2, 2,1.8,1.8,1.6,1.6), # Tardía: muchas especies con peso medio
Madura = c(1.2,1.2,1.2, 1.5,1.5, 1.8,1.8,1.8,1.8,1.8, 2,2,2,2,2)  # Madura: más equidad, varias especies importantes
)

simular_parcela <- function(etapa, total_ind = sample(90:150, 1)){ # Función que simula una parcela (conteos por especie)
 alpha <- alpha_por_etapa[[etapa]] # Selecciona el alpha que corresponde a la etapa
 p <- dirichlet_probs(alpha) # Genera probabilidades para las 15 especies (suman 1)
 as.integer(rmultinom(1, size = total_ind, prob = p))  # Genera conteos multinomiales (enteros) según p
} # Fin de función

abund <- parcelas %>% # Toma la tabla parcelas (Etapa, Parcela, ID)
 rowwise() %>% # Indica que operaciones se harán fila por fila
 mutate(conteos = list(simular_parcela(Etapa))) %>%  # Para cada fila, simula conteos de 15 especies
 unnest_wider(conteos, names_sep = "_") %>% # Expande la lista a 15 columnas (conteos_1...conteos_15)
 rename_with(~ especies, starts_with("conteos_")) %>%# Renombra esas columnas a Sp01...Sp15
 ungroup() # Quita el modo fila-por-fila

mat <- abund %>% select(all_of(especies)) %>% as.matrix()  # Extrae solo especies y las convierte a matriz
rownames(mat) <- abund$ID # Pone como nombres de fila los ID de parcelas
view(mat)
# -----------------------------
# FASE 2. Tema 2: Diversidad alfa (por parcela)
# -----------------------------

riqueza <- vegan::specnumber(mat)  # Calcula riqueza S (número de especies con abundancia > 0)
shannon <- vegan::diversity(mat, index = "shannon") # Calcula índice de Shannon H'
simpson <- vegan::diversity(mat, index = "simpson")  # Calcula índice de Simpson (diversidad)
pielou  <- shannon / log(pmax(riqueza, 1)) # Calcula equitatividad de Pielou (H'/ln(S))

res_alfa <- abund %>% # Parte del dataset con abundancias
 select(ID, Etapa, Parcela) %>%  # Se queda con variables de identificación
 mutate(Riqueza = riqueza, # Agrega riqueza por parcela
Shannon = shannon,  #Agrega Shannon por parcela
Simpson = simpson,  #Agrega Simpson por parcela
Pielou = pielou)  # Agrega Pielou por parcela

print(res_alfa)# Imprime tabla de diversidad alfa por parcela
## # A tibble: 12 × 7
##    ID            Etapa      Parcela Riqueza Shannon Simpson Pielou
##    <chr>         <fct>        <int>   <int>   <dbl>   <dbl>  <dbl>
##  1 Pionera_P1    Pionera          1      12    2.04   0.842  0.822
##  2 Intermedia_P1 Intermedia       1      13    2.43   0.900  0.948
##  3 Tardia_P1     Tardia           1      14    2.38   0.887  0.900
##  4 Madura_P1     Madura           1      13    2.36   0.892  0.920
##  5 Pionera_P2    Pionera          2      12    1.93   0.816  0.778
##  6 Intermedia_P2 Intermedia       2      15    2.49   0.906  0.918
##  7 Tardia_P2     Tardia           2      14    2.55   0.915  0.964
##  8 Madura_P2     Madura           2      15    2.50   0.902  0.922
##  9 Pionera_P3    Pionera          3      13    2.21   0.862  0.861
## 10 Intermedia_P3 Intermedia       3      13    2.30   0.874  0.896
## 11 Tardia_P3     Tardia           3      13    2.26   0.871  0.880
## 12 Madura_P3     Madura           3      15    2.36   0.879  0.871
res_etapa <- res_alfa %>%  #Usa resultados por parcela
group_by(Etapa) %>% # Agrupa por etapa sucesional
 summarise(  # Resume con promedios y desviaciones estándar
  Riqueza_prom = mean(Riqueza), Riqueza_sd = sd(Riqueza), # Media y sd de riqueza
  Shannon_prom = mean(Shannon), Shannon_sd = sd(Shannon), # Media y sd de Shannon
  Simpson_prom = mean(Simpson), Simpson_sd = sd(Simpson), # Media y sd de Simpson
  Pielou_prom = mean(Pielou), Pielou_sd = sd(Pielou), # Media y sd de Pielou
 .groups = "drop"  # Quita agrupamiento para devolver tibble normal
 )

print(res_etapa) # Imprime tabla resumen por etapa
## # A tibble: 4 × 9
##   Etapa  Riqueza_prom Riqueza_sd Shannon_prom Shannon_sd Simpson_prom Simpson_sd
##   <fct>         <dbl>      <dbl>        <dbl>      <dbl>        <dbl>      <dbl>
## 1 Pione…         12.3      0.577         2.06     0.138         0.840     0.0230
## 2 Inter…         13.7      1.15          2.40     0.0970        0.893     0.0170
## 3 Tardia         13.7      0.577         2.39     0.144         0.891     0.0222
## 4 Madura         14.3      1.15          2.41     0.0792        0.891     0.0117
## # ℹ 2 more variables: Pielou_prom <dbl>, Pielou_sd <dbl>
ggplot(res_alfa, aes(x = Etapa, y = Shannon)) + # Inicia gráfico: etapa vs Shannon
 geom_boxplot() +  # Dibuja caja y bigotes (distribución)
 geom_jitter(width = 0.12) +  #Pone puntos con leve dispersión horizontal
 theme_minimal() +  #Estilo simple del gráfico
 labs(title = "Diversidad (Shannon) por etapa de sucesión")  # Título del gráfico