# -----------------------------
# 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
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… 37 23 29 9 8 9 0 0 1
## 2 Intermedia 1 Inte… 10 5 18 10 11 5 5 20 3
## 3 Tardia 1 Tard… 0 5 3 6 8 8 10 8 21
## 4 Madura 1 Madu… 7 5 8 2 10 1 0 11 4
## 5 Pionera 2 Pion… 27 17 25 12 3 7 2 6 2
## 6 Intermedia 2 Inte… 17 13 8 26 9 8 20 7 5
## # ℹ 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 37 23 29 9 8 9 0 0 1 3 0 3 0
## Intermedia_P1 10 5 18 10 11 5 5 20 3 4 3 1 0
## Tardia_P1 0 5 3 6 8 8 10 8 21 13 15 9 5
## Madura_P1 7 5 8 2 10 1 0 11 4 1 17 9 18
## Pionera_P2 27 17 25 12 3 7 2 6 2 2 0 0 0
## Intermedia_P2 17 13 8 26 9 8 20 7 5 13 0 3 6
## Sp14 Sp15
## Pionera_P1 0 1
## Intermedia_P1 0 2
## Tardia_P1 2 10
## Madura_P1 6 3
## Pionera_P2 0 8
## Intermedia_P2 9 1
# ---------------------------------------------
# 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
## 10 13 14 14 11
## Intermedia_P2
## 14
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 10 1.84 0.803 0.797
## 2 Intermedia_P1 Intermedia 1 13 2.29 0.877 0.891
## 3 Tardia_P1 Tardia 1 14 2.49 0.907 0.945
## 4 Madura_P1 Madura 1 14 2.38 0.892 0.903
## 5 Pionera_P2 Pionera 2 11 2.04 0.841 0.852
## 6 Intermedia_P2 Intermedia 2 14 2.44 0.900 0.924
## 7 Tardia_P2 Tardia 2 13 2.34 0.887 0.913
## 8 Madura_P2 Madura 2 15 2.53 0.911 0.936
## 9 Pionera_P3 Pionera 3 13 2.16 0.861 0.844
## 10 Intermedia_P3 Intermedia 3 15 2.56 0.915 0.947
## 11 Tardia_P3 Tardia 3 14 2.47 0.904 0.936
## 12 Madura_P3 Madura 3 13 2.38 0.896 0.928
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… 11.3 1.53 2.01 0.166 0.835 0.0295
## 2 Inter… 14 1 2.43 0.139 0.897 0.0194
## 3 Tardia 13.7 0.577 2.44 0.0819 0.899 0.0110
## 4 Madura 14 1 2.43 0.0881 0.900 0.00969
## # ℹ 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
# 4. FASE 3 - Diversidad beta y NMDS
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.07584372
## Run 1 stress 0.07584372
## ... Procrustes: rmse 6.079355e-05 max resid 0.0001679119
## ... Similar to previous best
## Run 2 stress 0.07584378
## ... Procrustes: rmse 0.0001056972 max resid 0.0002928782
## ... Similar to previous best
## Run 3 stress 0.0758438
## ... Procrustes: rmse 0.0001211416 max resid 0.0003356847
## ... Similar to previous best
## Run 4 stress 0.07584373
## ... Procrustes: rmse 1.533565e-05 max resid 4.170945e-05
## ... Similar to previous best
## Run 5 stress 0.07584372
## ... New best solution
## ... Procrustes: rmse 1.826573e-05 max resid 4.824922e-05
## ... Similar to previous best
## Run 6 stress 0.07584373
## ... Procrustes: rmse 2.971308e-05 max resid 8.044785e-05
## ... Similar to previous best
## Run 7 stress 0.07584373
## ... Procrustes: rmse 1.989714e-05 max resid 5.451174e-05
## ... Similar to previous best
## Run 8 stress 0.07584372
## ... New best solution
## ... Procrustes: rmse 1.0176e-05 max resid 2.764527e-05
## ... Similar to previous best
## Run 9 stress 0.07584372
## ... Procrustes: rmse 3.286768e-05 max resid 9.139466e-05
## ... Similar to previous best
## Run 10 stress 0.07584372
## ... Procrustes: rmse 2.101646e-05 max resid 5.725167e-05
## ... Similar to previous best
## Run 11 stress 0.07584374
## ... Procrustes: rmse 7.456094e-05 max resid 0.0002070444
## ... Similar to previous best
## Run 12 stress 0.07584372
## ... Procrustes: rmse 2.228091e-05 max resid 6.193606e-05
## ... Similar to previous best
## Run 13 stress 0.07584372
## ... Procrustes: rmse 1.132283e-05 max resid 3.080314e-05
## ... Similar to previous best
## Run 14 stress 0.07584372
## ... Procrustes: rmse 8.995452e-06 max resid 2.050501e-05
## ... Similar to previous best
## Run 15 stress 0.07584376
## ... Procrustes: rmse 9.259999e-05 max resid 0.0002540213
## ... Similar to previous best
## Run 16 stress 0.07584372
## ... Procrustes: rmse 2.765577e-05 max resid 7.681595e-05
## ... Similar to previous best
## Run 17 stress 0.07584372
## ... Procrustes: rmse 9.148062e-06 max resid 2.268319e-05
## ... Similar to previous best
## Run 18 stress 0.1349903
## Run 19 stress 0.07584373
## ... Procrustes: rmse 1.750789e-05 max resid 3.643618e-05
## ... Similar to previous best
## Run 20 stress 0.07584372
## ... Procrustes: rmse 7.94749e-06 max resid 2.163178e-05
## ... Similar to previous best
## *** Best solution repeated 12 times
cat("OK: NMDS ejecutado. Stress = ", round(nmds$stress, 3), "\n", sep = "")
## OK: NMDS ejecutado. Stress = 0.076
# 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().