# -----------------------------
# 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 33 36 18 3 7 3 2 1 2 4 0 0 12
## Intermedia_P1 7 17 3 7 10 18 13 11 8 8 14 8 0
## Tardia_P1 13 3 4 13 21 11 4 6 21 0 2 21 20
## Madura_P1 1 3 6 1 9 8 16 2 5 24 3 6 17
## Pionera_P2 8 18 14 35 0 1 2 10 8 0 5 1 2
## Intermedia_P2 5 7 15 20 9 19 2 7 2 20 8 0 1
## Tardia_P2 8 2 15 1 4 15 17 16 19 12 10 0 0
## Madura_P2 6 2 1 8 0 9 15 2 6 10 4 17 1
## Pionera_P3 25 23 12 1 2 8 3 0 1 1 9 1 5
## Intermedia_P3 37 4 12 16 15 5 11 5 9 2 1 11 2
## Tardia_P3 5 3 6 0 7 22 8 10 12 5 13 5 9
## Madura_P3 1 5 5 13 5 4 5 5 2 14 2 10 11
## Sp14 Sp15
## Pionera_P1 0 9
## Intermedia_P1 2 1
## Tardia_P1 3 4
## Madura_P1 9 4
## Pionera_P2 3 1
## Intermedia_P2 3 1
## Tardia_P2 8 4
## Madura_P2 6 5
## Pionera_P3 0 3
## Intermedia_P3 4 12
## Tardia_P3 4 4
## Madura_P3 6 21
# -----------------------------
# 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 1.99 0.821 0.799
## 2 Intermedia_P1 Intermedia 1 14 2.47 0.907 0.935
## 3 Tardia_P1 Tardia 1 14 2.38 0.893 0.902
## 4 Madura_P1 Madura 1 15 2.38 0.886 0.881
## 5 Pionera_P2 Pionera 2 13 2.05 0.827 0.801
## 6 Intermedia_P2 Intermedia 2 14 2.31 0.882 0.874
## 7 Tardia_P2 Tardia 2 13 2.37 0.897 0.924
## 8 Madura_P2 Madura 2 14 2.38 0.892 0.903
## 9 Pionera_P3 Pionera 3 13 2.04 0.831 0.797
## 10 Intermedia_P3 Intermedia 3 15 2.38 0.880 0.877
## 11 Tardia_P3 Tardia 3 14 2.48 0.903 0.941
## 12 Madura_P3 Madura 3 15 2.46 0.898 0.909
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.03 0.0367 0.826 0.00502
## 2 Inter… 14.3 0.577 2.38 0.0809 0.890 0.0149
## 3 Tardia 13.7 0.577 2.41 0.0635 0.897 0.00500
## 4 Madura 14.7 0.577 2.41 0.0457 0.892 0.00605
## # ℹ 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.1097516
## Run 1 stress 0.1097516
## ... Procrustes: rmse 1.313502e-06 max resid 2.805925e-06
## ... Similar to previous best
## Run 2 stress 0.1097516
## ... Procrustes: rmse 2.09734e-06 max resid 4.285974e-06
## ... Similar to previous best
## Run 3 stress 0.1097516
## ... Procrustes: rmse 1.191379e-06 max resid 2.674147e-06
## ... Similar to previous best
## Run 4 stress 0.1097516
## ... Procrustes: rmse 5.319909e-07 max resid 1.330504e-06
## ... Similar to previous best
## Run 5 stress 0.1097516
## ... New best solution
## ... Procrustes: rmse 5.826081e-07 max resid 1.17815e-06
## ... Similar to previous best
## Run 6 stress 0.1097516
## ... Procrustes: rmse 1.452468e-06 max resid 2.897869e-06
## ... Similar to previous best
## Run 7 stress 0.1097516
## ... Procrustes: rmse 2.411699e-06 max resid 4.486274e-06
## ... Similar to previous best
## Run 8 stress 0.1097516
## ... Procrustes: rmse 2.280474e-06 max resid 4.606083e-06
## ... Similar to previous best
## Run 9 stress 0.1573215
## Run 10 stress 0.1097516
## ... Procrustes: rmse 1.521999e-06 max resid 3.179239e-06
## ... Similar to previous best
## Run 11 stress 0.1097516
## ... Procrustes: rmse 2.355476e-06 max resid 4.711175e-06
## ... Similar to previous best
## Run 12 stress 0.1097516
## ... Procrustes: rmse 6.582424e-07 max resid 1.224711e-06
## ... Similar to previous best
## Run 13 stress 0.2910527
## Run 14 stress 0.1097516
## ... Procrustes: rmse 2.859965e-06 max resid 5.776692e-06
## ... Similar to previous best
## Run 15 stress 0.1097516
## ... Procrustes: rmse 1.068511e-06 max resid 2.103236e-06
## ... Similar to previous best
## Run 16 stress 0.1097516
## ... Procrustes: rmse 2.774962e-06 max resid 5.553992e-06
## ... Similar to previous best
## Run 17 stress 0.1097516
## ... Procrustes: rmse 2.417125e-06 max resid 4.84199e-06
## ... Similar to previous best
## Run 18 stress 0.1097516
## ... Procrustes: rmse 1.181561e-06 max resid 2.348921e-06
## ... Similar to previous best
## Run 19 stress 0.1097516
## ... Procrustes: rmse 1.545716e-06 max resid 3.112014e-06
## ... Similar to previous best
## Run 20 stress 0.1097516
## ... Procrustes: rmse 1.078569e-06 max resid 2.071395e-06
## ... Similar to previous best
## *** Best solution repeated 14 times
cat("OK: NMDS ejecutado. Stress = ", round(nmds$stress, 3), "\n", sep = "")
## OK: NMDS ejecutado. Stress = 0.11
# 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().
cat("\n==============================\n")
##
## ==============================
cat("FASE 4: Curvas rango-abundancia\n")
## FASE 4: Curvas rango-abundancia
cat("==============================\n")
## ==============================
rank_df <- abund %>%
group_by(Etapa) %>%
summarise(across(all_of(especies), mean), .groups = "drop") %>%
pivot_longer(-Etapa, names_to = "Especie", values_to = "Abundancia") %>%
group_by(Etapa) %>%
arrange(desc(Abundancia), .by_group = TRUE) %>%
mutate(Rango = row_number())
ggplot(rank_df, aes(x = Rango, y = Abundancia)) +
geom_line() +
facet_wrap(~Etapa, scales = "free_y") +
theme_minimal() +
labs(title = "Curvas rango-abundancia por etapa")
cat("\nResumen FASE 4:\n")
##
## Resumen FASE 4:
cat("- Se promediaron abundancias por especie dentro de cada etapa.\n")
## - Se promediaron abundancias por especie dentro de cada etapa.
cat("- Se ordenaron especies por dominancia y se graficaron curvas.\n")
## - Se ordenaron especies por dominancia y se graficaron curvas.
cat("Comandos clave: summarise(across()), pivot_longer(), arrange(), row_number(), facet_wrap().\n")
## Comandos clave: summarise(across()), pivot_longer(), arrange(), row_number(), facet_wrap().
cat("\n==============================\n")
##
## ==============================
cat("FASE 5: Interacciones (red didactica)\n")
## FASE 5: Interacciones (red didactica)
cat("==============================\n")
## ==============================
tipos <- tibble(
Tipo = c("Mutualismo", "Comensalismo", "Competencia", "Depredacion", "Amensalismo", "Neutralismo"),
efecto_A = c(+1, +1, -1, +1, 0, 0),
efecto_B = c(+1, 0, -1, -1, -1, 0)
)
interacciones <- tibble(
A = sample(especies, 30, replace = TRUE),
B = sample(especies, 30, replace = TRUE),
Tipo = sample(tipos$Tipo, 30, replace = TRUE,
prob = c(0.15, 0.15, 0.25, 0.20, 0.10, 0.15))
) %>%
filter(A != B) %>%
left_join(tipos, by = "Tipo")
balance <- bind_rows(
interacciones %>% transmute(Especie = A, Efecto = efecto_A),
interacciones %>% transmute(Especie = B, Efecto = efecto_B)
) %>%
group_by(Especie) %>%
summarise(Balance_neto = sum(Efecto), n_interacciones = n(), .groups = "drop") %>%
arrange(desc(Balance_neto))
balance
## # A tibble: 15 × 3
## Especie Balance_neto n_interacciones
## <chr> <dbl> <int>
## 1 Sp04 3 3
## 2 Sp06 2 2
## 3 Sp08 2 2
## 4 Sp10 1 4
## 5 Sp14 1 1
## 6 Sp05 0 4
## 7 Sp07 0 7
## 8 Sp09 0 5
## 9 Sp15 0 8
## 10 Sp01 -1 3
## 11 Sp12 -1 3
## 12 Sp02 -2 4
## 13 Sp03 -2 3
## 14 Sp13 -2 5
## 15 Sp11 -4 4
g_int <- igraph::graph_from_data_frame(interacciones %>% select(A, B, Tipo), directed = TRUE)
plot(g_int, vertex.size = 16, vertex.label.cex = 0.7,
main = "Red de interacciones (didactica)")
cat("\nResumen FASE 5:\n")
##
## Resumen FASE 5:
cat("- Se simularon interacciones entre especies y se asignaron efectos.\n")
## - Se simularon interacciones entre especies y se asignaron efectos.
cat("- Se estimo balance neto por especie.\n")
## - Se estimo balance neto por especie.
cat("- Se construyo y grafico una red dirigida.\n")
## - Se construyo y grafico una red dirigida.
cat("Comandos clave: sample(), filter(), left_join(), bind_rows(), igraph::graph_from_data_frame().\n")
## Comandos clave: sample(), filter(), left_join(), bind_rows(), igraph::graph_from_data_frame().
cat("\n==============================\n")
##
## ==============================
cat("FASE 6: Niveles troficos y piramide de energia\n")
## FASE 6: Niveles troficos y piramide de energia
cat("==============================\n")
## ==============================
niveles <- tibble(
Especie = especies,
Nivel = case_when(
Especie %in% paste0("Sp", sprintf("%02d", 1:5)) ~ "Productores",
Especie %in% paste0("Sp", sprintf("%02d", 6:10)) ~ "Herbivoros",
Especie %in% paste0("Sp", sprintf("%02d", 11:14)) ~ "Carnivoros",
TRUE ~ "Apex"
)
)
biomasa_nivel <- abund %>%
select(ID, Etapa, all_of(especies)) %>%
pivot_longer(all_of(especies), names_to = "Especie", values_to = "Abund") %>%
left_join(niveles, by = "Especie") %>%
group_by(Etapa, ID, Nivel) %>%
summarise(Biomasa = sum(Abund), .groups = "drop") %>%
group_by(Etapa, Nivel) %>%
summarise(Biomasa_prom = mean(Biomasa), .groups = "drop")
E0 <- 10000
ef <- 0.10
orden <- c("Productores", "Herbivoros", "Carnivoros", "Apex")
energia_base <- tibble(
Nivel = orden,
Energia_teorica = E0 * ef^(match(Nivel, orden)-1)
)
energia <- biomasa_nivel %>%
left_join(energia_base, by = "Nivel") %>%
group_by(Etapa) %>%
mutate(prop_biomasa = Biomasa_prom / sum(Biomasa_prom),
Energia_asignada = Energia_teorica * prop_biomasa) %>%
ungroup()
ggplot(energia, aes(x = Nivel, y = Energia_asignada)) +
geom_col() +
facet_wrap(~Etapa) +
theme_minimal() +
labs(title = "Piramide de energia (modelo simple) por etapa",
y = "Energia asignada (unidades relativas)")
cat("\nResumen FASE 6:\n")
##
## Resumen FASE 6:
cat("- Se asignaron niveles troficos a especies.\n")
## - Se asignaron niveles troficos a especies.
cat("- Se calculo biomasa (proxy) por nivel y etapa.\n")
## - Se calculo biomasa (proxy) por nivel y etapa.
cat("- Se asigno energia por nivel (transferencia ~10%) y se grafico.\n")
## - Se asigno energia por nivel (transferencia ~10%) y se grafico.
cat("Comandos clave: case_when(), pivot_longer(), summarise(), mutate(), geom_col().\n")
## Comandos clave: case_when(), pivot_longer(), summarise(), mutate(), geom_col().
cat("\n==============================\n")
##
## ==============================
cat("FASE 7 (Opcional): Nutrientes N y P + razon N:P\n")
## FASE 7 (Opcional): Nutrientes N y P + razon N:P
cat("==============================\n")
## ==============================
nut <- parcelas %>%
mutate(
N_mgL = case_when(
Etapa == "Pionera" ~ runif(n(), 0.8, 1.6),
Etapa == "Intermedia" ~ runif(n(), 0.6, 1.2),
Etapa == "Tardia" ~ runif(n(), 0.4, 1.0),
TRUE ~ runif(n(), 0.3, 0.8)
),
P_mgL = case_when(
Etapa == "Pionera" ~ runif(n(), 0.05, 0.12),
Etapa == "Intermedia" ~ runif(n(), 0.04, 0.10),
Etapa == "Tardia" ~ runif(n(), 0.03, 0.08),
TRUE ~ runif(n(), 0.02, 0.06)
)
) %>%
mutate(
N_molar = N_mgL / 14,
P_molar = P_mgL / 31,
NP_molar = N_molar / P_molar,
Limitacion = case_when(
NP_molar < 16 ~ "Limitado por N",
NP_molar > 16 ~ "Limitado por P",
TRUE ~ "Cercano a equilibrio"
)
)
nut %>% select(ID, Etapa, N_mgL, P_mgL, NP_molar, Limitacion)
## # A tibble: 12 × 6
## ID Etapa N_mgL P_mgL NP_molar Limitacion
## <chr> <fct> <dbl> <dbl> <dbl> <chr>
## 1 Pionera_P1 Pionera 1.47 0.0961 33.8 Limitado por P
## 2 Intermedia_P1 Intermedia 1.02 0.0638 35.3 Limitado por P
## 3 Tardia_P1 Tardia 0.872 0.0433 44.6 Limitado por P
## 4 Madura_P1 Madura 0.678 0.0485 30.9 Limitado por P
## 5 Pionera_P2 Pionera 1.48 0.0819 40.0 Limitado por P
## 6 Intermedia_P2 Intermedia 0.843 0.0634 29.5 Limitado por P
## 7 Tardia_P2 Tardia 0.949 0.0779 27.0 Limitado por P
## 8 Madura_P2 Madura 0.746 0.0578 28.6 Limitado por P
## 9 Pionera_P3 Pionera 1.55 0.0893 38.4 Limitado por P
## 10 Intermedia_P3 Intermedia 0.741 0.0500 32.8 Limitado por P
## 11 Tardia_P3 Tardia 0.688 0.0326 46.7 Limitado por P
## 12 Madura_P3 Madura 0.566 0.0269 46.7 Limitado por P
ggplot(nut, aes(x = Etapa, y = NP_molar)) +
geom_boxplot() +
geom_hline(yintercept = 16, linetype = "dashed") +
theme_minimal() +
labs(title = "Razon molar N:P por etapa (linea = 16)",
y = "N:P (molar)")
cat("\nResumen FASE 7:\n")
##
## Resumen FASE 7:
cat("- Se simularon concentraciones de N y P por etapa.\n")
## - Se simularon concentraciones de N y P por etapa.
cat("- Se calculo N:P y se interpreto limitacion relativa usando el umbral 16.\n")
## - Se calculo N:P y se interpreto limitacion relativa usando el umbral 16.
cat("Comandos clave: runif(), case_when(), mutate(), geom_hline().\n")
## Comandos clave: runif(), case_when(), mutate(), geom_hline().