datos$abund_total <- rowSums(especies)

ggplot(datos, aes(x=Parcela, y=abund_total)) +
  geom_boxplot(fill="skyblue") +
  theme_minimal() +
  labs(title="Abundancia total por parcela",
       x="Parcela",
       y="Total de individuos")

# Fase 2 — Diversidad alfa por parcela

library(ggplot2)

# asegurar que las abundancias sean numéricas
especies <- data.frame(lapply(especies, as.numeric))

# riqueza de especies
riqueza <- rowSums(especies > 0)

# total de individuos por parcela
N <- rowSums(especies)

# proporciones
p <- especies / N

# índice de Shannon
shannon <- -rowSums(p * log(p), na.rm = TRUE)

# índice de Simpson
simpson <- 1 - rowSums(p^2, na.rm = TRUE)

# tabla de resultados
div_alpha <- data.frame(riqueza, shannon, simpson)

# gráfico
ggplot(div_alpha, aes(x = riqueza, y = shannon)) +
  geom_point(color = "black", size = 3) +
  theme_minimal() +
  labs(title = "Diversidad alfa",
       x = "Riqueza de especies",
       y = "Índice de Shannon")

# Fase 3 — Diversidad beta (Bray-Curtis + NMDS con símbolos por etapa)

library(vegan)

# asegurar que especies sean numéricas
especies <- data.frame(lapply(especies, as.numeric))

# distancia Bray-Curtis
dist_bc <- vegdist(especies, method = "bray")

# NMDS
nmds <- metaMDS(dist_bc, k = 2, trymax = 100)
## Run 0 stress 0.2359893 
## Run 1 stress 0.2416836 
## Run 2 stress 0.24089 
## Run 3 stress 0.2401874 
## Run 4 stress 0.2456725 
## Run 5 stress 0.2384584 
## Run 6 stress 0.2420171 
## Run 7 stress 0.2497067 
## Run 8 stress 0.2423942 
## Run 9 stress 0.244837 
## Run 10 stress 0.2435351 
## Run 11 stress 0.2396573 
## Run 12 stress 0.2407004 
## Run 13 stress 0.2370278 
## Run 14 stress 0.2465409 
## Run 15 stress 0.2451368 
## Run 16 stress 0.2405292 
## Run 17 stress 0.2456837 
## Run 18 stress 0.2461727 
## Run 19 stress 0.2396027 
## Run 20 stress 0.2454913 
## Run 21 stress 0.2472796 
## Run 22 stress 0.2474893 
## Run 23 stress 0.2402888 
## Run 24 stress 0.2426131 
## Run 25 stress 0.2532829 
## Run 26 stress 0.2393044 
## Run 27 stress 0.2443609 
## Run 28 stress 0.2451614 
## Run 29 stress 0.2442635 
## Run 30 stress 0.2392652 
## Run 31 stress 0.2459351 
## Run 32 stress 0.240987 
## Run 33 stress 0.2412861 
## Run 34 stress 0.2442569 
## Run 35 stress 0.2476911 
## Run 36 stress 0.2449062 
## Run 37 stress 0.2387351 
## Run 38 stress 0.2419828 
## Run 39 stress 0.2473441 
## Run 40 stress 0.2417272 
## Run 41 stress 0.2454807 
## Run 42 stress 0.2407297 
## Run 43 stress 0.2440366 
## Run 44 stress 0.2418088 
## Run 45 stress 0.2450816 
## Run 46 stress 0.2411047 
## Run 47 stress 0.2404376 
## Run 48 stress 0.2417591 
## Run 49 stress 0.2441451 
## Run 50 stress 0.2439626 
## Run 51 stress 0.2425817 
## Run 52 stress 0.2410538 
## Run 53 stress 0.2420953 
## Run 54 stress 0.2431994 
## Run 55 stress 0.243366 
## Run 56 stress 0.2487255 
## Run 57 stress 0.2406456 
## Run 58 stress 0.236798 
## Run 59 stress 0.2436913 
## Run 60 stress 0.239722 
## Run 61 stress 0.237142 
## Run 62 stress 0.2505865 
## Run 63 stress 0.2389298 
## Run 64 stress 0.2485323 
## Run 65 stress 0.2443844 
## Run 66 stress 0.2446347 
## Run 67 stress 0.2401312 
## Run 68 stress 0.2442051 
## Run 69 stress 0.2437379 
## Run 70 stress 0.2463793 
## Run 71 stress 0.2367978 
## Run 72 stress 0.244891 
## Run 73 stress 0.2457738 
## Run 74 stress 0.2467106 
## Run 75 stress 0.2464753 
## Run 76 stress 0.2514881 
## Run 77 stress 0.2390793 
## Run 78 stress 0.2441824 
## Run 79 stress 0.2418518 
## Run 80 stress 0.2408587 
## Run 81 stress 0.2427618 
## Run 82 stress 0.243522 
## Run 83 stress 0.2427404 
## Run 84 stress 0.241556 
## Run 85 stress 0.2446409 
## Run 86 stress 0.2378338 
## Run 87 stress 0.2426781 
## Run 88 stress 0.2478905 
## Run 89 stress 0.2452677 
## Run 90 stress 0.2434193 
## Run 91 stress 0.2477175 
## Run 92 stress 0.2390524 
## Run 93 stress 0.244438 
## Run 94 stress 0.240365 
## Run 95 stress 0.2460099 
## Run 96 stress 0.2469245 
## Run 97 stress 0.2399522 
## Run 98 stress 0.2482688 
## Run 99 stress 0.237515 
## Run 100 stress 0.2405812 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     21: no. of iterations >= maxit
##     79: stress ratio > sratmax
# coordenadas
coord <- as.data.frame(scores(nmds))

# traer etapa y normalizar texto
coord$Etapa <- tolower(trimws(as.character(datos$Etapa)))

# crear vector de símbolos (por defecto círculo)
pch_vec <- rep(16, nrow(coord))        # ● pionera
pch_vec[coord$Etapa == "intermedia"] <- 17  # ▲
pch_vec[coord$Etapa == "tardia"] <- 15      # ■
pch_vec[coord$Etapa == "tardía"] <- 15      # ■ (con acento)
pch_vec[coord$Etapa == "madura"] <- 3       # +

# gráfico
plot(coord$NMDS1, coord$NMDS2,
     pch = pch_vec,
     col = "darkblue",
     xlab = "NMDS1",
     ylab = "NMDS2",
     main = "NMDS Bray-Curtis por etapa sucesional")

# leyenda
legend("topright",
       legend = c("Pionera","Intermedia","Tardía","Madura"),
       pch = c(16,17,15,3),
       col = "darkblue")

# Fase 4 — Curvas rango-abundancia por etapa sucesional

# asegurar que especies sean numéricas
especies <- data.frame(lapply(especies, as.numeric))

# obtener etapas únicas
etapas <- unique(datos$Etapa)

# hacer gráfico para cada etapa
for(e in etapas){
  
  # filtrar datos de la etapa
  sub_especies <- especies[datos$Etapa == e, ]
  
  # abundancia total por especie
  abund <- colSums(sub_especies)
  
  # ordenar abundancias
  rank_abund <- sort(abund, decreasing = TRUE)
  
  # gráfico
  plot(rank_abund,
       type = "b",
       pch = 19,
       col = "black",
       xlab = "Rango de especies",
       ylab = "Abundancia",
       main = paste("Curva rango-abundancia -", e))
}

# Fase 5 — Red de interacciones ecológicas con todas las especies

library(igraph)

# asegurar que especies sean numéricas
especies <- data.frame(lapply(especies, as.numeric))

# calcular abundancia total por especie
abund_especies <- colSums(especies)

# crear matriz de interacción entre todas las especies
matriz_red <- outer(abund_especies, abund_especies, FUN = function(x,y) (x+y)/2)

# convertir matriz en red
red <- graph_from_adjacency_matrix(matriz_red,
                                   mode = "undirected",
                                   weighted = TRUE,
                                   diag = FALSE)

# graficar red
plot(red,
     vertex.size = 15,
     vertex.color = "orange",
     vertex.label = names(abund_especies),
     vertex.label.color = "black",
     main = "Red ecológica de especies")

# Fase 6 — Pirámides de energía por etapa sucesional

par(mfrow=c(2,2))  # 4 gráficos en una sola figura

niveles <- c("Productores","Herbivoros","Carnivoros","Depredadores")

# Pionera
energia <- c(12000,1200,120,12)
bp <- barplot(energia,
              names.arg=niveles,
              col=c("darkgreen","yellowgreen","orange","red"),
              main="Pionera",
              ylab="Energia",
              ylim=c(0,13000))
text(bp, energia+200, energia)

# Intermedia
energia <- c(10000,1000,100,10)
bp <- barplot(energia,
              names.arg=niveles,
              col=c("darkgreen","yellowgreen","orange","red"),
              main="Intermedia",
              ylab="Energia",
              ylim=c(0,11000))
text(bp, energia+200, energia)

# Tardia
energia <- c(8000,800,80,8)
bp <- barplot(energia,
              names.arg=niveles,
              col=c("darkgreen","yellowgreen","orange","red"),
              main="Tardia",
              ylab="Energia",
              ylim=c(0,9000))
text(bp, energia+200, energia)

# Madura
energia <- c(6000,600,60,6)
bp <- barplot(energia,
              names.arg=niveles,
              col=c("darkgreen","yellowgreen","orange","red"),
              main="Madura",
              ylab="Energia",
              ylim=c(0,7000))
text(bp, energia+200, energia)

# Fase 7 — Relación N:P por etapa (igual estilo que Fase 1)

# etapas
etapas <- c("Pionera","Intermedia","Tardia","Madura")

# valores de ejemplo de relación N:P
relacion_np <- c(12, 10, 8, 6)

# gráfico tipo Fase 1
barplot(relacion_np,
        names.arg = etapas,
        col = "forestgreen",
        main = "Relación N:P por etapa sucesional",
        ylab = "Relación N:P")