##1. Fundamentos del Modelo Binomial (Gianinetti, 2020)¿Por qué se utiliza este modelo?El análisis clásico de datos de germinación mediante un ANOVA tradicional presenta serios problemas metodológicos. Las semillas presentan una respuesta binaria dicotómica: germinan (éxito) o no germinan (fracaso). Al agruparlas en placas de Petri, la variable resultante es un conteo de éxitos obtenidos de un total fijo (\(n\)).Modelar esto bajo los supuestos del ANOVA (distribución Normal y varianza constante) es incorrecto porque:Heterocedasticidad intrínseca: En una distribución binomial, la varianza depende directamente de la media (\(\text{Varianza} = n \cdot p \cdot (1-p)\)). Placas con tasas de germinación cercanas al 50% exhibirán una varianza máxima, mientras que aquellas cercanas al 0% o 100% tendrán una varianza mínima.Inadecuación de transformaciones: Transformaciones clásicas como el arcoseno de la raíz cuadrada (\(\text{arcsin}(\sqrt{p})\)) fallan al estabilizar la varianza cuando existen valores extremos (0% o 100%) y pueden arrojar predicciones absurdas fuera del rango real (como porcentajes negativos o mayores al 100%).¿Cómo funciona el modelo?La solución óptima descrita por Gianinetti (2020) es el uso de un Modelo Lineal Generalizado Mixto (GzLMM) con familia Binomial y función de enlace Logit.La función Logit toma la probabilidad acotada de germinación (\(p\), que va de 0 a 1) y la transforma en un predictor lineal continuo que se extiende de \(-\infty\) a \(+\infty\):\[\text{logit}(p) = \ln\left(\frac{p}{1-p}\right)\]Adicionalmente, al tratarse de un modelo mixto, se añade un efecto aleatorio para identificar la variabilidad única de cada placa de Petri ((1 | ID_Placa)). Esto corrige la correlación ambiental que comparten las semillas dentro de un mismo recipiente (sobredispersión), evitando subestimar los errores estándar y reduciendo drásticamente la tasa de falsos positivos.

##2. Simulación de Datos (Diseño Experimental Controlado) Simulamos un escenario típico basado en el artículo: 3 tratamientos, 5 réplicas (placas) por tratamiento y 50 semillas por placa. Añadimos un componente de error aleatorio a nivel de placa para simular condiciones reales de laboratorio.

# Cargar librerías necesarias
library(lme4)
library(ggplot2)
library(car)

# Configuración de reproducibilidad
set.seed(42)
n_replicas        <- 5
semillas_por_placa <- 50

# Estructuración del diseño experimental
datos_germinacion <- expand.grid(
  Tratamiento = c("Control", "Tratamiento_A", "Tratamiento_B"),
  Placa       = paste0("Placa_", 1:n_replicas)
)
datos_germinacion$ID_Placa       <- factor(paste(datos_germinacion$Tratamiento,
                                                  datos_germinacion$Placa, sep = "_"))
datos_germinacion$Total_Semillas <- semillas_por_placa

# Probabilidades teóricas impuestas en la simulación
prob_base   <- c("Control" = 0.85, "Tratamiento_A" = 0.40, "Tratamiento_B" = 0.65)
error_placa <- rnorm(nrow(datos_germinacion), mean = 0, sd = 0.35)

# Generación aleatoria basada en la distribución binomial
datos_germinacion$Germinadas <- sapply(seq_len(nrow(datos_germinacion)), function(i) {
  p       <- prob_base[as.character(datos_germinacion$Tratamiento[i])]
  logit_p <- log(p / (1 - p)) + error_placa[i]
  rbinom(n = 1, size = semillas_por_placa, prob = 1 / (1 + exp(-logit_p)))
})

# Cálculo de proporciones observadas y fracasos
datos_germinacion$No_Germinadas <- datos_germinacion$Total_Semillas - datos_germinacion$Germinadas
datos_germinacion$Porcentaje_Obs <- (datos_germinacion$Germinadas / datos_germinacion$Total_Semillas) * 100

# Presentación de los datos simulados
knitr::kable(head(datos_germinacion, 6),
             caption = "Estructura de las primeras 6 filas del set de datos generado",
             digits  = 1)
Estructura de las primeras 6 filas del set de datos generado
Tratamiento Placa ID_Placa Total_Semillas Germinadas No_Germinadas Porcentaje_Obs
Control Placa_1 Control_Placa_1 50 44 6 88
Tratamiento_A Placa_1 Tratamiento_A_Placa_1 50 21 29 42
Tratamiento_B Placa_1 Tratamiento_B_Placa_1 50 35 15 70
Control Placa_2 Control_Placa_2 50 43 7 86
Tratamiento_A Placa_2 Tratamiento_A_Placa_2 50 13 37 26
Tratamiento_B Placa_2 Tratamiento_B_Placa_2 50 29 21 58 
# Ajuste del modelo lineal generalizado mixto
modelo_binom <- glmer(
  cbind(Germinadas, No_Germinadas) ~ Tratamiento + (1 | ID_Placa),
  data   = datos_germinacion,
  family = binomial(link = "logit")
)

summary(modelo_binom)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(Germinadas, No_Germinadas) ~ Tratamiento + (1 | ID_Placa)
##    Data: datos_germinacion
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##      91.0      93.9     -41.5      83.0        11 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2095 -0.5255 -0.1757  0.4761  1.6751 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID_Placa (Intercept) 0.07693  0.2774  
## Number of obs: 15, groups:  ID_Placa, 15
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                2.0991     0.2378   8.826  < 2e-16 ***
## TratamientoTratamiento_A  -2.4456     0.2984  -8.195 2.50e-16 ***
## TratamientoTratamiento_B  -1.2569     0.3013  -4.171 3.03e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtT_A
## TrtmntTrt_A -0.799       
## TrtmntTrt_B -0.786  0.627

##4.1 Porcentaje de Germinación Observado por Tratamiento Este primer gráfico permite visualizar la dispersión real de los datos crudos. Es posible notar cómo la variabilidad física se distribuye a lo largo de las réplicas debido al efecto de placa introducido en la simulación.

# Paleta de colores consistente para los reportes
paleta <- c("Control"       = "#2ecc71",
            "Tratamiento_A" = "#e74c3c",
            "Tratamiento_B" = "#3498db")

ggplot(datos_germinacion,
       aes(x = Tratamiento, y = Porcentaje_Obs, fill = Tratamiento)) +
  geom_boxplot(alpha = 0.65, color = "grey30", outlier.shape = NA) +
  geom_jitter(aes(color = Tratamiento), width = 0.12,
              size = 2.8, alpha = 0.9, show.legend = FALSE) +
  scale_fill_manual(values  = paleta) +
  scale_color_manual(values = paleta) +
  labs(
    title    = "Porcentaje de Germinación Observado por Tratamiento",
    subtitle = "Datos crudos (puntos) acompañados de diagramas de caja y bigotes",
    x        = "Tratamiento",
    y        = "% Semillas Germinadas",
    fill     = "Tratamiento"
  ) +
  theme_classic(base_size = 13) +
  theme(
    legend.position  = "top",
    plot.title       = element_text(face = "bold"),
    plot.subtitle    = element_text(color = "grey40")
  )

##4.2 Probabilidades Predichas por el Modelo (Efectos Fijos)Utilizando el Método Delta, extraemos los errores estándar combinados y aplicamos la función logística inversa (plogis) para retrotransformar los estimadores lineales a porcentajes puros de germinación. Esto asegura que las estimaciones y sus intervalos de confianza respeten estrictamente los límites biológicos (\(0\%\) a \(100\%\)).

# Extraer coeficientes de efectos fijos y matriz de varianza-covarianza
fe  <- fixef(modelo_binom)
vc  <- vcov(modelo_binom)

# Matriz de contraste para aislar cada tratamiento
L <- rbind(
  Control       = c(1, 0, 0),
  Tratamiento_A = c(1, 1, 0),
  Tratamiento_B = c(1, 0, 1)
)
rownames(L) <- c("Control", "Tratamiento_A", "Tratamiento_B")

pred_df <- data.frame(
  Tratamiento = rownames(L),
  logit_est   = as.vector(L %*% fe)
)

# Cálculo del error estándar (SE) aproximado mediante combinaciones lineales
pred_df$se <- sqrt(diag(L %*% vc %*% t(L)))

# Retrotransformación matemática a escala porcentual de 0 a 100
pred_df$Prob_Pred <- plogis(pred_df$logit_est) * 100
pred_df$IC_Inf    <- plogis(pred_df$logit_est - 1.96 * pred_df$se) * 100
pred_df$IC_Sup    <- plogis(pred_df$logit_est + 1.96 * pred_df$se) * 100

# Renderizado de las predicciones medias marginales
ggplot(pred_df, aes(x = Tratamiento, y = Prob_Pred, fill = Tratamiento)) +
  geom_col(alpha = 0.80, color = "grey30", width = 0.5) +
  geom_errorbar(aes(ymin = IC_Inf, ymax = IC_Sup),
                width = 0.18, linewidth = 1.1, color = "grey20") +
  geom_text(aes(label = paste0(round(Prob_Pred, 1), "%")),
            vjust = -0.6, size = 4.2, fontface = "bold") +
  scale_fill_manual(values = paleta) +
  labs(
    title    = "Probabilidades Predichas por el Modelo GzLMM",
    subtitle = "Medias estimadas e Intervalos de Confianza (95%) vía Método Delta",
    x        = "Tratamiento",
    y        = "% Germinación Predicha por el Modelo",
    fill     = "Tratamiento"
  ) +
  ylim(0, 105) +
  theme_classic(base_size = 13) +
  theme(
    legend.position = "none",
    plot.title      = element_text(face = "bold"),
    plot.subtitle   = element_text(color = "grey40")
  )

##4.3 Residuos de Pearson vs. Valores PredichosUn supuesto fundamental de los Modelos Lineales Generalizados es que los residuos ponderados (Pearson) deben comportarse de forma homogénea alrededor de cero. Si los puntos se localizan mayoritariamente dentro de la banda de \(\pm 2\), se confirma que el modelo no padece problemas severos de sobredispersión ni de datos atípicos fuera de control.

# Construcción del set de datos para diagnóstico
resid_df <- data.frame(
  Predichos   = fitted(modelo_binom),
  Residuos    = residuals(modelo_binom, type = "pearson"),
  Tratamiento = datos_germinacion$Tratamiento
)

# Renderizado de residuos de Pearson
ggplot(resid_df, aes(x = Predichos, y = Residuos, color = Tratamiento)) +
  geom_hline(yintercept = 0, linetype = "dashed",
             color = "red", linewidth = 0.8) +
  geom_hline(yintercept = c(-2, 2), linetype = "dotted",
             color = "orange", linewidth = 0.6) +
  geom_point(size = 3, alpha = 0.85) +
  scale_color_manual(values = paleta) +
  labs(
    title    = "Diagnóstico: Residuos de Pearson vs. Valores Predichos",
    subtitle = "Las líneas punteadas naranjas marcan el umbral óptimo de variabilidad (±2)",
    x        = "Valores Predichos (Probabilidad en escala 0-1)",
    y        = "Residuos estandarizados de Pearson",
    color    = "Tratamiento"
  ) +
  theme_classic(base_size = 13) +
  theme(
    legend.position = "top",
    plot.title      = element_text(face = "bold"),
    plot.subtitle   = element_text(color = "grey40")
  )

##5. Prueba de Hipótesis: Análisis de Devianza Tipo IIPara realizar pruebas de hipótesis globales sobre variables categóricas en modelos generalizados, Gianinetti (2020) descarta las sumas de cuadrados ordinarias y promueve las pruebas basadas en la comparación de ratios de verosimilitud o aproximaciones de Wald. Empleamos un Análisis de Devianza Tipo II fundamentado en la distribución Chi-cuadrado (\(\chi^2\)).

# Análisis de devianza utilizando la función Anova de la librería 'car'
resultado_anova <- Anova(modelo_binom, type = "II")

knitr::kable(resultado_anova,
             caption = "Tabla del Análisis de Devianza Tipo II (Prueba de Wald)",
             digits  = 4)
Tabla del Análisis de Devianza Tipo II (Prueba de Wald)
Chisq Df Pr(>Chisq)
Tratamiento 68.716 2 0

##Conclusión Analítica FinalAl examinar la tabla de Devianza Tipo II generada arriba, se evalúa el p-valor (Pr(>Chisq)). Si este indicador cae por debajo del nivel de significancia crítico (\(\alpha = 0.05\)), se rechaza formalmente la hipótesis nula, concluyendo con rigurosidad científica que los tratamientos aplicados ejercen un efecto real y diferenciable sobre la capacidad germinativa de las semillas evaluadas, validando el marco teórico propuesto por Gianinetti (2020).