# -----------------------------
# 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
parcelas
## # A tibble: 12 × 3
##    Etapa      Parcela ID           
##    <fct>        <int> <chr>        
##  1 Pionera          1 Pionera_P1   
##  2 Intermedia       1 Intermedia_P1
##  3 Tardia           1 Tardia_P1    
##  4 Madura           1 Madura_P1    
##  5 Pionera          2 Pionera_P2   
##  6 Intermedia       2 Intermedia_P2
##  7 Tardia           2 Tardia_P2    
##  8 Madura           2 Madura_P2    
##  9 Pionera          3 Pionera_P3   
## 10 Intermedia       3 Intermedia_P3
## 11 Tardia           3 Tardia_P3    
## 12 Madura           3 Madura_P3
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
# fila-por-fila
head(abund)
## # A tibble: 6 × 18
##   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…     8    22    57     2     6     1     7     3     4
## 2 Intermedia       1 Inte…     3    11    22    16     4     7     6    13     3
## 3 Tardia           1 Tard…     0     1     3     3    17     7    18     6     8
## 4 Madura           1 Madu…     1     5     1     1     9     1     1    23     3
## 5 Pionera          2 Pion…    35    29     9     2    12     8     0     0     4
## 6 Intermedia       2 Inte…    11    19    29    19     1     5    10     3    19
## # ℹ 6 more variables: Sp10 <int>, Sp11 <int>, Sp12 <int>, Sp13 <int>,
## #   Sp14 <int>, Sp15 <int>
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
head(mat)
##               Sp01 Sp02 Sp03 Sp04 Sp05 Sp06 Sp07 Sp08 Sp09 Sp10 Sp11 Sp12 Sp13
## Pionera_P1       8   22   57    2    6    1    7    3    4    7    0    3    0
## Intermedia_P1    3   11   22   16    4    7    6   13    3    8    3    6    5
## Tardia_P1        0    1    3    3   17    7   18    6    8    6    8   22    4
## Madura_P1        1    5    1    1    9    1    1   23    3    6    8   12    4
## Pionera_P2      35   29    9    2   12    8    0    0    4    6    0    1    6
## Intermedia_P2   11   19   29   19    1    5   10    3   19    2    5    0   12
##               Sp14 Sp15
## Pionera_P1       0    0
## Intermedia_P1    3    1
## Tardia_P1        2    3
## Madura_P1       10    7
## Pionera_P2       1    0
## Intermedia_P2    0    9
# -----------------------------
# FASE 2. Tema 2: Diversidad alfa (por parcela)
# -----------------------------

riqueza  <- vegan::specnumber(mat)                                              # Calcula riqueza S (número de especies con abundancia > 0)
head(riqueza)
##    Pionera_P1 Intermedia_P1     Tardia_P1     Madura_P1    Pionera_P2 
##            11            15            14            15            11 
## Intermedia_P2 
##            13
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      11    1.73   0.724  0.722
##  2 Intermedia_P1 Intermedia       1      15    2.45   0.895  0.905
##  3 Tardia_P1     Tardia           1      14    2.34   0.880  0.885
##  4 Madura_P1     Madura           1      15    2.32   0.875  0.856
##  5 Pionera_P2    Pionera          2      11    1.92   0.808  0.802
##  6 Intermedia_P2 Intermedia       2      13    2.29   0.883  0.895
##  7 Tardia_P2     Tardia           2      15    2.50   0.906  0.924
##  8 Madura_P2     Madura           2      14    2.29   0.879  0.869
##  9 Pionera_P3    Pionera          3      12    2.00   0.821  0.804
## 10 Intermedia_P3 Intermedia       3      15    2.27   0.859  0.839
## 11 Tardia_P3     Tardia           3      14    2.50   0.909  0.949
## 12 Madura_P3     Madura           3      14    2.39   0.890  0.906
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)     
## # 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…         11.3      0.577         1.88     0.137         0.785    0.0526 
## 2 Inter…         14.3      1.15          2.34     0.0982        0.879    0.0184 
## 3 Tardia         14.3      0.577         2.45     0.0971        0.899    0.0157 
## 4 Madura         14.3      0.577         2.33     0.0500        0.881    0.00761
## # ℹ 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