# -----------------------------
# 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      23   20   10    1    1    9    0    0    1    1    3    2    9
## Intermedia_P1    5   19   10   11   28    1    7   14    0    0    0    4    2
## Tardia_P1        4    0    1   16   27   11    4    0   12    7   16    4    6
## Madura_P1        8    0    5   21   33    2    4    2    2    6    4    4   27
## Pionera_P2      13   20   28    6    9    1   10    0    5    1    1    0    1
## Intermedia_P2   14   11    9   19   18    9   19   11    3    1    1    2    0
## Tardia_P2        3    4    4   11    3    8   13   30   10    5   33    0    0
## Madura_P2        0    1    0    6   11   20    4   11   14    4   14    5    1
## Pionera_P3      29   22   35   14    5    1    0   19   11    0    0    2    5
## Intermedia_P3   15   14   21   17   12    5    3    6    5    6    8    5    3
## Tardia_P3        1   15    6    2    0    5    6   28   16    7    2    3    0
## Madura_P3        0   13    3    2    5    7   17    2   25    3    7    1    1
##               Sp14 Sp15
## Pionera_P1       0   10
## Intermedia_P1   14    4
## Tardia_P1       11   12
## Madura_P1        3   21
## Pionera_P2       0    0
## Intermedia_P2    5   13
## Tardia_P2       16    7
## Madura_P2        0    5
## Pionera_P3       2    0
## Intermedia_P3    3   15
## Tardia_P3        4    0
## Madura_P3        1    5
# -----------------------------
# 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 × 22
##    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…    23    20    10     1     1     9     0     0     1
##  2 Intermed…       1 Inte…     5    19    10    11    28     1     7    14     0
##  3 Tardia          1 Tard…     4     0     1    16    27    11     4     0    12
##  4 Madura          1 Madu…     8     0     5    21    33     2     4     2     2
##  5 Pionera         2 Pion…    13    20    28     6     9     1    10     0     5
##  6 Intermed…       2 Inte…    14    11     9    19    18     9    19    11     3
##  7 Tardia          2 Tard…     3     4     4    11     3     8    13    30    10
##  8 Madura          2 Madu…     0     1     0     6    11    20     4    11    14
##  9 Pionera         3 Pion…    29    22    35    14     5     1     0    19    11
## 10 Intermed…       3 Inte…    15    14    21    17    12     5     3     6     5
## 11 Tardia          3 Tard…     1    15     6     2     0     5     6    28    16
## 12 Madura          3 Madu…     0    13     3     2     5     7    17     2    25
## # ℹ 10 more variables: Sp10 <int>, Sp11 <int>, Sp12 <int>, Sp13 <int>,
## #   Sp14 <int>, Sp15 <int>, Riqueza <int>, Shannon <dbl>, Simpson <dbl>,
## #   Pielou <dbl>
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      0.577         2.00     0.0491        0.835     0.0109
## 2 Inter…         13.7      1.53          2.37     0.160         0.891     0.0207
## 3 Tardia         12.7      0.577         2.24     0.126         0.866     0.0247
## 4 Madura         13.3      1.15          2.20     0.0197        0.861     0.0124
## # ℹ 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.1382918 
## Run 1 stress 0.1601306 
## Run 2 stress 0.2614448 
## Run 3 stress 0.2359999 
## Run 4 stress 0.1540641 
## Run 5 stress 0.1545858 
## Run 6 stress 0.1382919 
## ... Procrustes: rmse 0.0002036543  max resid 0.000318869 
## ... Similar to previous best
## Run 7 stress 0.1879122 
## Run 8 stress 0.1382919 
## ... Procrustes: rmse 0.0002904816  max resid 0.0004725002 
## ... Similar to previous best
## Run 9 stress 0.2308853 
## Run 10 stress 0.1468114 
## Run 11 stress 0.3114879 
## Run 12 stress 0.1540641 
## Run 13 stress 0.182391 
## Run 14 stress 0.1545858 
## Run 15 stress 0.179486 
## Run 16 stress 0.1382918 
## ... Procrustes: rmse 0.0001804889  max resid 0.0002813324 
## ... Similar to previous best
## Run 17 stress 0.1539309 
## Run 18 stress 0.141615 
## Run 19 stress 0.1446524 
## Run 20 stress 0.1539309 
## *** Best solution repeated 3 times
cat("OK: NMDS ejecutado. Stress = ", round(nmds$stress, 3), "\n", sep = "")
## OK: NMDS ejecutado. Stress = 0.138
# 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 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.1382918 
## Run 1 stress 0.1912676 
## Run 2 stress 0.1382918 
## ... Procrustes: rmse 9.937033e-05  max resid 0.0001583856 
## ... Similar to previous best
## Run 3 stress 0.1382919 
## ... Procrustes: rmse 0.0002983314  max resid 0.0004950996 
## ... Similar to previous best
## Run 4 stress 0.1414167 
## Run 5 stress 0.1420803 
## Run 6 stress 0.1540641 
## Run 7 stress 0.1541543 
## Run 8 stress 0.2600689 
## Run 9 stress 0.1539309 
## Run 10 stress 0.1545858 
## Run 11 stress 0.1446524 
## Run 12 stress 0.1545858 
## Run 13 stress 0.1854356 
## Run 14 stress 0.1540641 
## Run 15 stress 0.1601306 
## Run 16 stress 0.1541543 
## Run 17 stress 0.190198 
## Run 18 stress 0.1382919 
## ... Procrustes: rmse 0.0002851825  max resid 0.0004676199 
## ... Similar to previous best
## Run 19 stress 0.1382919 
## ... Procrustes: rmse 0.0002890476  max resid 0.0004890907 
## ... Similar to previous best
## Run 20 stress 0.2822425 
## *** Best solution repeated 4 times
cat("OK: NMDS ejecutado. Stress = ", round(nmds$stress, 3), "\n", sep = "")
## OK: NMDS ejecutado. Stress = 0.138
# 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("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 Sp13               4               5
##  2 Sp02               2               4
##  3 Sp08               1               2
##  4 Sp12               1               5
##  5 Sp09               0               6
##  6 Sp14               0               3
##  7 Sp15               0               6
##  8 Sp03              -1               3
##  9 Sp04              -1               7
## 10 Sp05              -1               1
## 11 Sp07              -1               1
## 12 Sp01              -2               2
## 13 Sp10              -2               2
## 14 Sp06              -3               5
## 15 Sp11              -3               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.36  0.0553     54.5 Limitado por P
##  2 Intermedia_P1 Intermedia 1.16  0.0969     26.6 Limitado por P
##  3 Tardia_P1     Tardia     0.737 0.0352     46.4 Limitado por P
##  4 Madura_P1     Madura     0.697 0.0410     37.7 Limitado por P
##  5 Pionera_P2    Pionera    1.21  0.0693     38.7 Limitado por P
##  6 Intermedia_P2 Intermedia 1.05  0.0632     36.8 Limitado por P
##  7 Tardia_P2     Tardia     0.961 0.0466     45.6 Limitado por P
##  8 Madura_P2     Madura     0.691 0.0558     27.4 Limitado por P
##  9 Pionera_P3    Pionera    0.951 0.0946     22.3 Limitado por P
## 10 Intermedia_P3 Intermedia 0.820 0.0766     23.7 Limitado por P
## 11 Tardia_P3     Tardia     0.624 0.0554     24.9 Limitado por P
## 12 Madura_P3     Madura     0.371 0.0549     15.0 Limitado por N
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().