Ingreso de información

# -----------------------------
# 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
view(abund)
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)
view(riqueza)
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.14   0.850  0.859
##  2 Intermedia_P1 Intermedia       1      14    2.30   0.868  0.871
##  3 Tardia_P1     Tardia           1      15    2.50   0.910  0.924
##  4 Madura_P1     Madura           1      14    2.46   0.900  0.932
##  5 Pionera_P2    Pionera          2      12    2.02   0.835  0.812
##  6 Intermedia_P2 Intermedia       2      15    2.25   0.855  0.831
##  7 Tardia_P2     Tardia           2      15    2.40   0.888  0.885
##  8 Madura_P2     Madura           2      14    2.39   0.896  0.907
##  9 Pionera_P3    Pionera          3      14    2.22   0.856  0.842
## 10 Intermedia_P3 Intermedia       3      15    2.50   0.908  0.924
## 11 Tardia_P3     Tardia           3      15    2.36   0.883  0.871
## 12 Madura_P3     Madura           3      15    2.56   0.914  0.946
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      1.15          2.13     0.102         0.847    0.0108 
## 2 Inter…         14.7      0.577         2.35     0.134         0.877    0.0278 
## 3 Tardia         15        0             2.42     0.0746        0.894    0.0143 
## 4 Madura         14.3      0.577         2.47     0.0842        0.903    0.00906
## # ℹ 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

`

cat("\n==============================\n")
## 
## ==============================
cat("FASE 3: Diversidad beta (Bray-Curtis) y NMDS\n")
## FASE 3: Diversidad beta (Bray-Curtis) y NMDS
cat("==============================\n")
## ==============================
d_bray <- vegan::vegdist(mat, method = "bray")
nmds <- vegan::metaMDS(mat, distance = "bray", k = 2, trymax = 100)
## Wisconsin double standardization
## Run 0 stress 0.09696105 
## Run 1 stress 0.09696125 
## ... Procrustes: rmse 0.001071754  max resid 0.002428227 
## ... Similar to previous best
## Run 2 stress 0.09696122 
## ... Procrustes: rmse 0.001093041  max resid 0.002466439 
## ... Similar to previous best
## Run 3 stress 0.1164435 
## Run 4 stress 0.1142937 
## Run 5 stress 0.09696113 
## ... Procrustes: rmse 0.001012058  max resid 0.002287718 
## ... Similar to previous best
## Run 6 stress 0.118626 
## Run 7 stress 0.1186262 
## Run 8 stress 0.09696092 
## ... New best solution
## ... Procrustes: rmse 0.00076363  max resid 0.001716072 
## ... Similar to previous best
## Run 9 stress 0.09696081 
## ... New best solution
## ... Procrustes: rmse 0.0003446332  max resid 0.0007691888 
## ... Similar to previous best
## Run 10 stress 0.09696104 
## ... Procrustes: rmse 0.0005143398  max resid 0.001155989 
## ... Similar to previous best
## Run 11 stress 0.096961 
## ... Procrustes: rmse 0.0004751445  max resid 0.001066309 
## ... Similar to previous best
## Run 12 stress 0.1147603 
## Run 13 stress 0.1142937 
## Run 14 stress 0.09696081 
## ... Procrustes: rmse 3.399584e-05  max resid 7.615779e-05 
## ... Similar to previous best
## Run 15 stress 0.1142938 
## Run 16 stress 0.1164437 
## Run 17 stress 0.1175782 
## Run 18 stress 0.118626 
## Run 19 stress 0.1164435 
## Run 20 stress 0.1147603 
## *** Best solution repeated 4 times
cat("OK: NMDS ejecutado. Stress = ", round(nmds$stress, 3), "\n", sep = "")
## OK: NMDS ejecutado. Stress = 0.097
# Importante: usar SOLO "sites" para evitar el error de filas (parcelas vs especies)
scores_nmds <- as.data.frame(vegan::scores(nmds, display = "sites")) %>%
  tibble::rownames_to_column("ID") %>%
  left_join(res_alfa %>% select(ID, Etapa), by = "ID")

ggplot(scores_nmds, aes(x = NMDS1, y = NMDS2, shape = Etapa)) +
  geom_point(size = 3) +
  theme_minimal() +
  labs(title = paste0("NMDS (Bray-Curtis) - stress = ", round(nmds$stress, 3)))

gamma_n <- sum(colSums(mat) > 0)
cat("OK: Diversidad gamma (especies presentes): ", gamma_n, "\n", sep = "")
## OK: Diversidad gamma (especies presentes): 15
cat("\nResumen FASE 3:\n")
## 
## Resumen FASE 3:
cat("- Se calculo disimilitud Bray-Curtis.\n")
## - Se calculo disimilitud Bray-Curtis.
cat("- Se ejecuto NMDS 2D y se graficaron coordenadas de parcelas (sites).\n")
## - Se ejecuto NMDS 2D y se graficaron coordenadas de parcelas (sites).
cat("- Se calculo diversidad gamma.\n")
## - Se calculo diversidad gamma.
cat("Comandos clave: vegan::vegdist(), vegan::metaMDS(), vegan::scores(display='sites'), colSums().\n")
## Comandos clave: vegan::vegdist(), vegan::metaMDS(), vegan::scores(display='sites'), colSums().