library(ggplot2)
library(dplyr)
library(gridExtra)

Diseño Experimental para el Estudio de Maíz con Fertilización Nitrogenada y Densidad de Siembra

El maíz es uno de los cultivos más importantes a nivel mundial, y su rendimiento está influenciado por factores como la fertilización nitrogenada y la densidad de siembra. Este estudio evaluó el efecto de tres niveles de nitrógeno (0, 100 y 200 kg N/ha) y dos densidades de siembra (60,000 y 80,000 plantas/ha) en el rendimiento de grano de maíz bajo condiciones de riego. Se utilizó un diseño experimental de bloques completos al azar con arreglo factorial (3 × 2) y cuatro repeticiones. Los resultados mostraron que la interacción entre nitrógeno y densidad de siembra fue significativa (p < 0.05). La mayor producción (12.5 t/ha) se obtuvo con 200 kg N/ha y 80,000 plantas/ha. Se concluye que la fertilización nitrogenada y la densidad de siembra deben optimizarse conjuntamente para maximizar el rendimiento del maíz. ahora hazlo con este ejercicio

set.seed(31415)
n_repeticiones <- 4  # Cuatro repeticiones por tratamiento
niveles_N <- c(0, 100, 200)  # kg N/ha
densidades <- c(60000, 80000)  # plantas/ha
datos_maiz <- expand.grid(
  Repeticion = 1:n_repeticiones,
  Nitrogeno = niveles_N,
  Densidad = densidades
)
rendimiento_base <- c(
  6.2, 5.8,  # 0 kg N/ha
  9.1, 10.3,  # 100 kg N/ha
  11.4, 12.5  # 200 kg N/ha
)
datos_maiz$Rendimiento <- rep(rendimiento_base, each = n_repeticiones) + 
  rnorm(nrow(datos_maiz), mean = 0, sd = 0.5)
datos_maiz$Nitrogeno <- factor(datos_maiz$Nitrogeno)
datos_maiz$Densidad <- factor(datos_maiz$Densidad)
modelo_maiz <- aov(Rendimiento ~ Nitrogeno * Densidad, data = datos_maiz)
summary(modelo_maiz)
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## Nitrogeno           2  30.55   15.28   36.19 4.93e-07 ***
## Densidad            1 118.54  118.54  280.81 1.99e-12 ***
## Nitrogeno:Densidad  2   2.57    1.28    3.04   0.0729 .  
## Residuals          18   7.60    0.42                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
if (summary(modelo_maiz)[[1]]$"Pr(>F)"[3] < 0.05) {
  tukey_result <- TukeyHSD(modelo_maiz)
  print(tukey_result)
}
interaction_plot <- ggplot(datos_maiz, aes(x = Nitrogeno, y = Rendimiento, 
                                       color = Densidad, group = Densidad)) +
  stat_summary(fun = mean, geom = "line", linewidth = 1) +  # Corrección aquí
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  labs(title = "Interacción Nitrógeno-Densidad en rendimiento de maíz",
       x = "Nitrógeno (kg/ha)",
       y = "Rendimiento (t/ha)",
       color = "Densidad\n(plantas/ha)") +
  scale_color_manual(values = c("#1b9e77", "#d95f02")) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")
treatment_boxplot <- ggplot(datos_maiz, 
                           aes(x = interaction(Nitrogeno, Densidad), 
                           y = Rendimiento, 
                           fill = interaction(Nitrogeno, Densidad))) +
  geom_boxplot() +
  labs(title = "Rendimiento por combinación de tratamientos",
       x = "Combinación Nitrógeno-Densidad",
       y = "Rendimiento (t/ha)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none") +
  scale_fill_brewer(palette = "Set2")
grid.arrange(interaction_plot, treatment_boxplot, ncol = 1)

# Prueba Z para densidad 80000
grupo_200_80k <- subset(datos_maiz, Nitrogeno == "200" & Densidad == "80000")$Rendimiento
grupo_100_80k <- subset(datos_maiz, Nitrogeno == "100" & Densidad == "80000")$Rendimiento

z_80k <- (mean(grupo_200_80k) - mean(grupo_100_80k)) / 
         sqrt((sd(grupo_200_80k)^2 / length(grupo_200_80k)) + (sd(grupo_100_80k)^2 / length(grupo_100_80k)))

p_80k <- 2 * (1 - pnorm(abs(z_80k)))

cat("Densidad 80,000 plantas/ha:\n")
## Densidad 80,000 plantas/ha:
cat("Z =", round(z_80k, 2), " - p =", round(p_80k, 4), "\n\n")
## Z = 3.09  - p = 0.002
# Prueba Z para densidad 60000
grupo_200_60k <- subset(datos_maiz, Nitrogeno == "200" & Densidad == "60000")$Rendimiento
grupo_100_60k <- subset(datos_maiz, Nitrogeno == "100" & Densidad == "60000")$Rendimiento

z_60k <- (mean(grupo_200_60k) - mean(grupo_100_60k)) / 
         sqrt((sd(grupo_200_60k)^2 / length(grupo_200_60k)) + (sd(grupo_100_60k)^2 / length(grupo_100_60k)))

p_60k <- 2 * (1 - pnorm(abs(z_60k)))

cat("Densidad 60,000 plantas/ha:\n")
## Densidad 60,000 plantas/ha:
cat("Z =", round(z_60k, 2), " - p =", round(p_60k, 4), "\n")
## Z = 7.68  - p = 0
# Preparar datos solo para las comparaciones deseadas
datos_z <- datos_maiz %>%
  filter((Nitrogeno %in% c("100", "200")) & (Densidad %in% c("60000", "80000"))) %>%
  mutate(Tratamiento = paste0(Nitrogeno, "N-", Densidad))

# Convertir a factor ordenado
datos_z$Tratamiento <- factor(datos_z$Tratamiento, 
                              levels = c("100N-60000", "200N-60000", 
                                         "100N-80000", "200N-80000"))

# Calcular valores p por densidad
calcular_z_p <- function(df) {
  grupo_200 <- df %>% filter(Nitrogeno == "200") %>% pull(Rendimiento)
  grupo_100 <- df %>% filter(Nitrogeno == "100") %>% pull(Rendimiento)
  
  z <- (mean(grupo_200) - mean(grupo_100)) / 
       sqrt((sd(grupo_200)^2 / length(grupo_200)) + (sd(grupo_100)^2 / length(grupo_100)))
  
  p <- 2 * (1 - pnorm(abs(z)))
  return(round(p, 4))
}

# Calcular p para cada densidad
p_60000 <- calcular_z_p(filter(datos_z, Densidad == "60000"))
p_80000 <- calcular_z_p(filter(datos_z, Densidad == "80000"))

# Crear etiquetas para cada faceta
etiquetas <- c(
  `60000` = paste0("Densidad: 60,000 plantas/ha\np = ", p_60000),
  `80000` = paste0("Densidad: 80,000 plantas/ha\np = ", p_80000)
)

# Crear gráfico facetado
ggplot(datos_z, aes(x = Tratamiento, y = Rendimiento, fill = Tratamiento)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7) +
  geom_jitter(width = 0.1, color = "black", size = 2) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "red") +
  facet_wrap(~Densidad, labeller = as_labeller(etiquetas)) +
  labs(title = "Comparación de rendimiento entre 100 y 200 kg N/ha",
       x = "Tratamiento",
       y = "Rendimiento (t/ha)") +
  scale_fill_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c")) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        strip.text = element_text(face = "bold", size = 11))