# ---------------------------------------------
# FASE 2. Tema 2: Diversidad alfa (por parcela)
# ---------------------------------------------
riqueza <- vegan::specnumber(mat)
shannon <- vegan::diversity(mat, index = "shannon")
simpson <- vegan::diversity(mat, index = "simpson")
pielou <- shannon / log(pmax(riqueza, 1))
res_alfa <- abund %>%
select(ID, Etapa, Parcela) %>%
mutate(Riqueza = riqueza,
Shannon = shannon,
Simpson = simpson,
Pielou = pielou)
print(res_alfa)
## # A tibble: 12 × 7
## ID Etapa Parcela Riqueza Shannon Simpson Pielou
## <chr> <fct> <int> <int> <dbl> <dbl> <dbl>
## 1 Pionera_P1 Pionera 1 13 1.80 0.752 0.703
## 2 Intermedia_P1 Intermedia 1 12 2.23 0.878 0.899
## 3 Tardia_P1 Tardia 1 14 2.40 0.894 0.910
## 4 Madura_P1 Madura 1 15 2.42 0.897 0.894
## 5 Pionera_P2 Pionera 2 12 2.05 0.828 0.823
## 6 Intermedia_P2 Intermedia 2 14 2.37 0.891 0.900
## 7 Tardia_P2 Tardia 2 15 2.49 0.908 0.920
## 8 Madura_P2 Madura 2 15 2.34 0.879 0.862
## 9 Pionera_P3 Pionera 3 13 2.19 0.862 0.853
## 10 Intermedia_P3 Intermedia 3 13 2.14 0.850 0.834
## 11 Tardia_P3 Tardia 3 15 2.37 0.886 0.875
## 12 Madura_P3 Madura 3 13 2.24 0.872 0.875
res_etapa <- res_alfa %>%
group_by(Etapa) %>%
summarise(
Riqueza_prom = mean(Riqueza), Riqueza_sd = sd(Riqueza),
Shannon_prom = mean(Shannon), Shannon_sd = sd(Shannon),
Simpson_prom = mean(Simpson), Simpson_sd = sd(Simpson),
Pielou_prom = mean(Pielou), Pielou_sd = sd(Pielou),
.groups = "drop"
)
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… 12.7 0.577 2.01 0.194 0.814 0.0563
## 2 Inter… 13 1 2.25 0.119 0.873 0.0207
## 3 Tardia 14.7 0.577 2.42 0.0634 0.896 0.0113
## 4 Madura 14.3 1.15 2.33 0.0881 0.882 0.0128
## # ℹ 2 more variables: Pielou_prom <dbl>, Pielou_sd <dbl>
# ---------------------------------------------
# Grafica con colores
# ---------------------------------------------
ggplot(res_alfa, aes(x = Etapa, y = Shannon, fill = Etapa)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(aes(color = Etapa), width = 0.12, size = 2) +
theme_minimal() +
labs(title = "Diversidad (Shannon) por etapa de sucesion",
x = "Etapa",
y = "Indice de Shannon") +
scale_fill_manual(values = c("forestgreen","goldenrod","orange","darkgreen")) +
scale_color_manual(values = c("forestgreen","goldenrod","orange","darkgreen"))
# ---------------------------------------------
# FASE 3. Tema 3: Diversidad beta (comparación entre parcelas)
# ---------------------------------------------
# Distancia de Bray-Curtis entre parcelas
dist_bray <- vegan::vegdist(mat, method = "bray")
# NMDS (ordenamiento)
nmds <- vegan::metaMDS(mat, distance = "bray", k = 2, trymax = 100)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1246719
## Run 1 stress 0.1246719
## ... New best solution
## ... Procrustes: rmse 5.244923e-06 max resid 1.030692e-05
## ... Similar to previous best
## Run 2 stress 0.1246719
## ... Procrustes: rmse 7.71903e-06 max resid 1.441002e-05
## ... Similar to previous best
## Run 3 stress 0.1246722
## ... Procrustes: rmse 0.0001758076 max resid 0.0003351495
## ... Similar to previous best
## Run 4 stress 0.2109697
## Run 5 stress 0.2109699
## Run 6 stress 0.1365088
## Run 7 stress 0.1246719
## ... Procrustes: rmse 2.752051e-05 max resid 5.482544e-05
## ... Similar to previous best
## Run 8 stress 0.1246719
## ... Procrustes: rmse 1.820937e-06 max resid 4.579406e-06
## ... Similar to previous best
## Run 9 stress 0.2233244
## Run 10 stress 0.2153307
## Run 11 stress 0.1365094
## Run 12 stress 0.2424096
## Run 13 stress 0.1365091
## Run 14 stress 0.1246719
## ... Procrustes: rmse 4.883068e-06 max resid 9.47842e-06
## ... Similar to previous best
## Run 15 stress 0.1246719
## ... Procrustes: rmse 4.340509e-06 max resid 6.770627e-06
## ... Similar to previous best
## Run 16 stress 0.2295725
## Run 17 stress 0.1246719
## ... Procrustes: rmse 7.564766e-06 max resid 1.451213e-05
## ... Similar to previous best
## Run 18 stress 0.1246719
## ... Procrustes: rmse 1.316589e-05 max resid 2.721122e-05
## ... Similar to previous best
## Run 19 stress 0.1761933
## Run 20 stress 0.2095094
## *** Best solution repeated 9 times
# Coordenadas del NMDS
coords <- as.data.frame(nmds$points)
coords$ID <- rownames(coords)
# Unir con información de etapa
coords <- coords %>%
left_join(abund %>% select(ID, Etapa), by = "ID")
print(coords)
## MDS1 MDS2 ID Etapa
## 1 -0.35107390 -0.019733018 Pionera_P1 Pionera
## 2 -0.18462156 -0.280942268 Intermedia_P1 Intermedia
## 3 0.07188079 0.171968184 Tardia_P1 Tardia
## 4 0.12883693 0.063821723 Madura_P1 Madura
## 5 -0.17785474 0.279006477 Pionera_P2 Pionera
## 6 -0.00876143 0.002055373 Intermedia_P2 Intermedia
## 7 0.05498570 0.063650575 Tardia_P2 Tardia
## 8 0.24153095 0.095721912 Madura_P2 Madura
## 9 -0.03933602 -0.138381056 Pionera_P3 Pionera
## 10 -0.26668554 0.029250761 Intermedia_P3 Intermedia
## 11 0.33494075 -0.012544989 Tardia_P3 Tardia
## 12 0.19615807 -0.253873674 Madura_P3 Madura
# ---------------------------------------------
# Gráfica NMDS con colores por etapa
# ---------------------------------------------
ggplot(coords, aes(x = MDS1, y = MDS2, color = Etapa)) +
geom_point(size = 4) +
theme_minimal() +
labs(title = "Ordenamiento NMDS (diversidad beta)",
x = "NMDS1",
y = "NMDS2") +
scale_color_manual(values = c("forestgreen","goldenrod","orange","darkgreen"))
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, color = Etapa)) +
geom_line(size = 1.2) +
facet_wrap(~Etapa, scales = "free_y") +
theme_minimal() +
labs(title = "Curvas rango-abundancia por etapa",
x = "Rango de especie",
y = "Abundancia promedio") +
scale_color_manual(values = c("forestgreen","goldenrod","orange","darkgreen"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
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().
# 6. FASE 5 - Interacciones y red (didactica)
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 Sp12 2 5
## 2 Sp06 1 4
## 3 Sp07 1 4
## 4 Sp01 0 9
## 5 Sp03 0 3
## 6 Sp04 0 2
## 7 Sp09 0 1
## 8 Sp15 0 5
## 9 Sp02 -1 3
## 10 Sp05 -1 1
## 11 Sp11 -1 2
## 12 Sp14 -1 3
## 13 Sp10 -2 2
## 14 Sp13 -2 3
## 15 Sp08 -3 7
# Crear red
g_int <- igraph::graph_from_data_frame(interacciones %>% select(A, B, Tipo), directed = TRUE)
# Colores según balance
V(g_int)$color <- ifelse(balance$Balance_neto[match(V(g_int)$name, balance$Especie)] > 0,
"forestgreen",
ifelse(balance$Balance_neto[match(V(g_int)$name, balance$Especie)] < 0,
"orange",
"gold"))
# Graficar red con colores
plot(g_int,
vertex.size = 18,
vertex.label.cex = 0.8,
vertex.color = V(g_int)$color,
edge.arrow.size = 0.4,
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()
# -----------------------------
# Grafica con colores
# -----------------------------
ggplot(energia, aes(x = Nivel, y = Energia_asignada, fill = Nivel)) +
geom_col() +
facet_wrap(~Etapa) +
theme_minimal() +
labs(title = "Piramide de energia (modelo simple) por etapa",
x = "Nivel trofico",
y = "Energia asignada (unidades relativas)") +
scale_fill_manual(values = c("forestgreen","goldenrod","orange","red"))
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.21 0.0951 28.3 Limitado por P
## 2 Intermedia_P1 Intermedia 0.707 0.0500 31.3 Limitado por P
## 3 Tardia_P1 Tardia 0.912 0.0446 45.3 Limitado por P
## 4 Madura_P1 Madura 0.679 0.0235 64.0 Limitado por P
## 5 Pionera_P2 Pionera 1.19 0.0545 48.5 Limitado por P
## 6 Intermedia_P2 Intermedia 0.651 0.0974 14.8 Limitado por N
## 7 Tardia_P2 Tardia 0.628 0.0398 34.9 Limitado por P
## 8 Madura_P2 Madura 0.530 0.0300 39.2 Limitado por P
## 9 Pionera_P3 Pionera 0.811 0.0953 18.8 Limitado por P
## 10 Intermedia_P3 Intermedia 0.797 0.0543 32.5 Limitado por P
## 11 Tardia_P3 Tardia 0.736 0.0468 34.8 Limitado por P
## 12 Madura_P3 Madura 0.437 0.0551 17.5 Limitado por P
# -----------------------------
# Grafica con colores
# -----------------------------
ggplot(nut, aes(x = Etapa, y = NP_molar, fill = Etapa)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(aes(color = Etapa), width = 0.12, size = 2) +
geom_hline(yintercept = 16, linetype = "dashed", color = "red") +
theme_minimal() +
labs(title = "Razon molar N:P por etapa (linea = 16)",
x = "Etapa sucesional",
y = "N:P (molar)") +
scale_fill_manual(values = c("forestgreen","goldenrod","orange","darkgreen")) +
scale_color_manual(values = c("forestgreen","goldenrod","orange","darkgreen"))
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().