# -----------------------------
# 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
mat
##               Sp01 Sp02 Sp03 Sp04 Sp05 Sp06 Sp07 Sp08 Sp09 Sp10 Sp11 Sp12 Sp13
## Pionera_P1      22   17   17   17    8    6    8    1    1    1    2    4    3
## Intermedia_P1   14    4   24   10    1    7    4    2    4    6    5    4    9
## Tardia_P1       12   23    2   15   16    3   11   10    9   12    5    9    7
## Madura_P1        5   10    4    7    3    2   11   16    4   12    1    6   12
## Pionera_P2      24   34   15    5   15    0   17    1    1   13    0    3    3
## Intermedia_P2    5   18    7    9    8    4    1   10    2    0    1    8    0
## Tardia_P2        2    2    4    3   11   18    0    6    7   11   22    4    7
## Madura_P2        1    5   14    0   19    0    5    1   12   14   17   21    8
## Pionera_P3      19   37   22    3   10    1    4    1    2    0    4    5    5
## Intermedia_P3    7   17   13    9    5    0   10    9    4    4    7    4    5
## Tardia_P3        8    4    5   23   23   12    7   12    2    7    5   11    3
## Madura_P3        8    0    4    3    6    4    1    3    5    5   12   22    6
##               Sp14 Sp15
## Pionera_P1       0    0
## Intermedia_P1    3    0
## Tardia_P1        3   12
## Madura_P1       13   12
## Pionera_P2       3    0
## Intermedia_P2    1   17
## Tardia_P2        7    4
## Madura_P2       19    5
## Pionera_P3      11    0
## Intermedia_P3    2    5
## Tardia_P3        5   16
## Madura_P3        5    7
# -----------------------------
# 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 × 22
##    Etapa     Parcela ID     Sp01  Sp02  Sp03  Sp04  Sp05  Sp06  Sp07  Sp08  Sp09
##    <fct>       <int> <chr> <int> <int> <int> <int> <int> <int> <int> <int> <int>
##  1 Pionera         1 Pion…    22    17    17    17     8     6     8     1     1
##  2 Intermed…       1 Inte…    14     4    24    10     1     7     4     2     4
##  3 Tardia          1 Tard…    12    23     2    15    16     3    11    10     9
##  4 Madura          1 Madu…     5    10     4     7     3     2    11    16     4
##  5 Pionera         2 Pion…    24    34    15     5    15     0    17     1     1
##  6 Intermed…       2 Inte…     5    18     7     9     8     4     1    10     2
##  7 Tardia          2 Tard…     2     2     4     3    11    18     0     6     7
##  8 Madura          2 Madu…     1     5    14     0    19     0     5     1    12
##  9 Pionera         3 Pion…    19    37    22     3    10     1     4     1     2
## 10 Intermed…       3 Inte…     7    17    13     9     5     0    10     9     4
## 11 Tardia          3 Tard…     8     4     5    23    23    12     7    12     2
## 12 Madura          3 Madu…     8     0     4     3     6     4     1     3     5
## # ℹ 10 more variables: Sp10 <int>, Sp11 <int>, Sp12 <int>, Sp13 <int>,
## #   Sp14 <int>, Sp15 <int>, Riqueza <int>, Shannon <dbl>, Simpson <dbl>,
## #   Pielou <dbl>
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.7      0.577         2.12     0.0538        0.850     0.0148
## 2 Inter…         13.7      0.577         2.37     0.122         0.888     0.0171
## 3 Tardia         14.7      0.577         2.48     0.0805        0.902     0.0124
## 4 Madura         14        1             2.42     0.0967        0.897     0.0129
## # ℹ 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