TALLER DE DISEÑOS EXPERIMENTALES - FACTORIAL SIMPLE EN ARREGLO EN BLOQUES

TEMA 3 SELECCIONADO: DISEÑO EN BLOQUES INCOMPLETOS (DBIB)

Caso: Rendimiento (t/ha) de 6 variedades de papa evaluadas en bloques

de 3 unidades experimentales cada uno (t = 6 tratamientos, k = 3)

Autor: NICOL ARIAS,ANDRES CUASQUER, MARIA PAULA PINILLA Y ANGELY PINZÓN

PUNTO 1. DEFINICION Y RAZON DE BLOQUEO (texto de referencia - ver informe) ## ============================================================================ # Se emplea un Diseno en Bloques Incompletos porque el numero de tratamientos # (6 variedades) es mayor que la capacidad de cada bloque (solo 3 unidades # experimentales por bloque). Los bloques se forman para controlar la # variabilidad natural del terreno: fertilidad, humedad y profundidad # efectiva del suelo. Al no caber las 6 variedades en un mismo bloque, no # existe ortogonalidad completa entre variedades y bloques, y las # comparaciones entre variedades usan informacion intra-bloque (mas precisa) # e inter-bloque (menos precisa). Este es el desarrollo estadistico completo # de ese diseno seleccionado.

============================================================================

PUNTO 2. MODELO ESTADISTICO (ya formalizado en el documento del grupo)

============================================================================

Y_ij = mu + T_i + beta_j + eps_ij

Y_ij : rendimiento (t/ha) de la variedad i en el bloque j

mu : media general del experimento

T_i : efecto fijo de la variedad i (i = 1,…,6)

beta_j: efecto de bloque j (j = 1,…,b); se asume FIJO porque los

bloques representan sectores especificos del lote experimental,

no una muestra aleatoria de una poblacion mas amplia de lotes.

eps_ij: error experimental, eps_ij ~ N(0, sigma^2) iid

Al ser bloque fijo, se usa aov() con F clasico y NO lmer(). Si en su caso

los bloques fueran una muestra aleatoria de fincas/lotes, este mismo

esqueleto de codigo cambiaria fitting con lmer(y ~ variedad + (1|block)).

## ---------------------------------------------------------------------------
## 0. PAQUETES
## ---------------------------------------------------------------------------
paquetes <- c("agricolae", "lme4", "car", "emmeans", "ggplot2", "dplyr", "multcomp")
instalar_faltantes <- paquetes[!(paquetes %in% installed.packages()[, "Package"])]
if (length(instalar_faltantes) > 0) install.packages(instalar_faltantes)

library(agricolae)  # design.bib(), BIB.test(), HSD.test()

## Warning: package 'agricolae' was built under R version 4.5.3

library(lme4)       # lmer() - por si se requiere version con bloque aleatorio

## Warning: package 'lme4' was built under R version 4.5.3

## Cargando paquete requerido: Matrix

library(car)         # leveneTest()

## Warning: package 'car' was built under R version 4.5.3

## Cargando paquete requerido: carData

## Warning: package 'carData' was built under R version 4.5.3

## Registered S3 method overwritten by 'car':
##   method           from
##   na.action.merMod lme4

library(emmeans)    # medias marginales ajustadas

## Warning: package 'emmeans' was built under R version 4.5.3

## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'

library(ggplot2)    # graficos
library(dplyr)       # manipulacion de datos

## 
## Adjuntando el paquete: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

—————————————————————————

PUNTO 3. GENERACION DE DATOS SIMULADOS EN R

—————————————————————————

set.seed(123)  # semilla fija para reproducibilidad
# 3.1 Definicion del factor de tratamiento A: variedad de papa (6 niveles)
variedades <- c("Pastusa", "Diacol_Capiro", "Criolla_Colombia",
                "ICA_Unica", "Superior", "Parda_Pastusa")
t <- length(variedades)   # numero de tratamientos = 6
k <- 3                    # tamano de bloque = unidades experimentales por bloque

# 3.2 Generacion del arreglo de campo (incidencia variedad-bloque).
# IMPORTANTE: el documento base de este grupo define los bloques como
# B1, B2, ..., B6 (j = 1,...,6), es decir b = 6 BLOQUES, no el diseno
# balanceado minimo de agricolae::design.bib() (que para t=6, k=3 arma
# automaticamente b=10, r=5). Con b=6 no es posible lograr un diseno
# perfectamente BALANCEADO (lambda constante) para 6 variedades en bloques
# de 3 -- por eso se construye a mano un diseno CICLICO: cada bloque
# contiene 3 variedades consecutivas (modulo 6). Este arreglo si cumple
# exactamente lo que pide el documento: b = 6 bloques, k = 3 unidades por
# bloque, r = 3 repeticiones por variedad, b*k = 6*3 = 18 unidades
# experimentales en total.
#   B1 = {v1,v2,v3}   B2 = {v2,v3,v4}   B3 = {v3,v4,v5}
#   B4 = {v4,v5,v6}   B5 = {v5,v6,v1}   B6 = {v6,v1,v2}
# (cada variedad aparece en exactamente 3 de los 6 bloques; el balance NO es
# perfecto -- algunos pares de variedades coinciden 2 veces, otros 1 vez y
# los pares "opuestos" (ej. v1-v4) no coinciden nunca -- pero es un diseno
# valido de bloques incompletos, coherente con el documento base).
b <- 6  # numero de bloques (B1...B6, segun el documento base)
bloques_incidencia <- list(
  B1 = variedades[c(1,2,3)],
  B2 = variedades[c(2,3,4)],
  B3 = variedades[c(3,4,5)],
  B4 = variedades[c(4,5,6)],
  B5 = variedades[c(5,6,1)],
  B6 = variedades[c(6,1,2)]
)

set.seed(123)  # para el orden aleatorio de las unidades dentro de cada bloque
book <- do.call(rbind, lapply(names(bloques_incidencia), function(nb) {
  vars_bloque <- sample(bloques_incidencia[[nb]])  # aleatoriza el orden dentro del bloque
  data.frame(block = nb, variedad = vars_bloque, stringsAsFactors = FALSE)
}))
book$plots <- 101:(100 + nrow(book))
book <- book[, c("plots", "block", "variedad")]

book$block    <- factor(book$block, levels = paste0("B", 1:6))
book$variedad <- factor(book$variedad, levels = variedades)

r <- as.numeric(table(book$variedad)[1])       # repeticiones por variedad (debe ser 3)
cat("Diseno en bloques incompletos generado -> t =", t, " k =", k, " b =", b, " r =", r, "\n")

## Diseno en bloques incompletos generado -> t = 6  k = 3  b = 6  r = 3

cat("Total de unidades experimentales = b * k =", b * k, "\n")

## Total de unidades experimentales = b * k = 18

# lambda NO es constante en este diseno (no es un BIBD balanceado perfecto);
# se reporta el rango observado para dejarlo documentado:
pares <- combn(variedades, 2)
lambda_pares <- apply(pares, 2, function(p) {
  sum(sapply(bloques_incidencia, function(bl) all(p %in% bl)))
})
cat("Rango de coincidencias por par de variedades (lambda no es constante): ",
    paste(range(lambda_pares), collapse = " a "), "\n")

## Rango de coincidencias por par de variedades (lambda no es constante):  0 a 2

# 3.3 Condiciones edaficas por bloque (fertilidad, humedad, profundidad
# efectiva del suelo -- las mismas que el documento base senala como razon
# de bloqueo). Se simulan primero estas variables y LUEGO se construye el
# efecto de bloque a partir de ellas, para que el bloqueo tenga sentido
# agronomico (bloques con mejores condiciones = mayor rendimiento base) en
# vez de ser un numero aleatorio sin relacion con la razon de bloqueo.
set.seed(456)
condiciones_edaficas <- data.frame(
  block       = factor(paste0("B", 1:6), levels = levels(book$block)),
  fertilidad  = round(rnorm(b, mean = 50, sd = 8), 1),   # indice de fertilidad (0-100)
  humedad     = round(rnorm(b, mean = 30, sd = 5), 1),    # % de humedad del suelo
  profundidad = round(rnorm(b, mean = 40, sd = 6), 1)     # profundidad efectiva (cm)
)

# Indice edafico compuesto: promedio de las 3 variables estandarizadas
# (cada una en su propia escala, por eso se estandarizan antes de combinar).
condiciones_edaficas$indice_edafico <- rowMeans(scale(condiciones_edaficas[, c("fertilidad","humedad","profundidad")]))

# 3.4 Grafico de condiciones edaficas por bloque (mapa de calor: color =
# nivel relativo de cada condicion entre bloques; numero = valor real
# simulado). Esta es la evidencia visual de POR QUE se bloqueo el terreno.
condiciones_long <- condiciones_edaficas %>%
  select(block, fertilidad, humedad, profundidad) %>%
  tidyr::pivot_longer(cols = c(fertilidad, humedad, profundidad),
                       names_to = "variable", values_to = "valor") %>%
  group_by(variable) %>%
  mutate(valor_z = as.numeric(scale(valor))) %>%
  ungroup() %>%
  mutate(variable = dplyr::recode(variable,
                            fertilidad  = "Fertilidad (indice)",
                            humedad     = "Humedad del suelo (%)",
                            profundidad = "Profundidad efectiva (cm)"))

p_edafico <- ggplot(condiciones_long, aes(x = block, y = variable, fill = valor_z)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(aes(label = valor), color = "black", size = 4) +
  scale_fill_gradient2(low = "firebrick", mid = "khaki", high = "forestgreen",
                        midpoint = 0, name = "Nivel\nrelativo") +
  labs(title = "Condiciones edaficas por bloque",
       subtitle = "Razon de bloqueo: fertilidad, humedad y profundidad efectiva del suelo",
       x = "Bloque", y = NULL) +
  theme_minimal()
print(p_edafico)

# 3.5 Simulacion del rendimiento con efectos conocidos + error aleatorio.
# Se fijan efectos "verdaderos" de variedad para poder verificar despues que
# el analisis los recupera correctamente.
mu <- 25  # rendimiento medio general del experimento (t/ha)

efecto_variedad <- c(Pastusa = 2.0, Diacol_Capiro = -1.0, Criolla_Colombia = 3.5,
                      ICA_Unica = 0.0, Superior = -2.5, Parda_Pastusa = 1.0)

# Efecto de bloque: ya NO es un numero aleatorio sin relacion con nada; se
# deriva directamente del indice edafico compuesto (fertilidad + humedad +
# profundidad) que se acaba de graficar arriba. Un bloque con mejores
# condiciones edaficas (indice mas alto) tiene un efecto positivo sobre el
# rendimiento; uno con peores condiciones, un efecto negativo.
efecto_bloque <- condiciones_edaficas$indice_edafico * 1.5
names(efecto_bloque) <- condiciones_edaficas$block

sigma_error <- 1.2  # desviacion estandar del error experimental

book$rendimiento <- mu +
  efecto_variedad[as.character(book$variedad)] +
  efecto_bloque[as.character(book$block)] +
  rnorm(nrow(book), mean = 0, sd = sigma_error)
book$rendimiento <- round(book$rendimiento, 2)

head(book)

##   plots block         variedad rendimiento
## 1   101    B1 Criolla_Colombia       30.91
## 2   102    B1          Pastusa       28.52
## 3   103    B1    Diacol_Capiro       23.10
## 4   104    B2 Criolla_Colombia       27.13
## 5   105    B2    Diacol_Capiro       22.98
## 6   106    B2        ICA_Unica       25.94

str(book)

## 'data.frame':    18 obs. of  4 variables:
##  $ plots      : int  101 102 103 104 105 106 107 108 109 110 ...
##  $ block      : Factor w/ 6 levels "B1","B2","B3",..: 1 1 1 2 2 2 3 3 3 4 ...
##  $ variedad   : Factor w/ 6 levels "Pastusa","Diacol_Capiro",..: 3 1 2 3 2 4 4 5 3 4 ...
##  $ rendimiento: num  30.9 28.5 23.1 27.1 23 ...

# 3.6 TABLA DE INCIDENCIA / LIBRO DE CAMPO: que variedad cae en cada bloque.
# Esta es la tabla que establece la estructura del diseno: filas = bloques,
# columnas = variedades; "X" indica que esa variedad SI esta presente en ese
# bloque (recordar: en un DBIB cada bloque solo contiene k = 3 de las 6
# variedades).
tabla_incidencia <- table(book$block, book$variedad)
tabla_incidencia_legible <- ifelse(tabla_incidencia > 0, "X", "")
cat("\n--- Tabla de incidencia: variedades presentes en cada bloque ---\n")

## 
## --- Tabla de incidencia: variedades presentes en cada bloque ---

print(tabla_incidencia_legible, quote = FALSE)

##     
##      Pastusa Diacol_Capiro Criolla_Colombia ICA_Unica Superior Parda_Pastusa
##   B1 X       X             X                                                
##   B2         X             X                X                               
##   B3                       X                X         X                     
##   B4                                        X         X        X            
##   B5 X                                                X        X            
##   B6 X       X                                                 X

# Version alternativa, mas parecida a un "libro de campo": una fila por
# bloque, con las 3 variedades que contiene listadas en columnas separadas.
libro_de_campo <- book %>%
  arrange(block) %>%
  group_by(block) %>%
  summarise(variedades_en_el_bloque = paste(variedad, collapse = " | "),
            .groups = "drop")
cat("\n--- Libro de campo: variedades asignadas a cada bloque ---\n")

## 
## --- Libro de campo: variedades asignadas a cada bloque ---

print(libro_de_campo, n = Inf)

## # A tibble: 6 × 2
##   block variedades_en_el_bloque                     
##   <fct> <chr>                                       
## 1 B1    Criolla_Colombia | Pastusa | Diacol_Capiro  
## 2 B2    Criolla_Colombia | Diacol_Capiro | ICA_Unica
## 3 B3    ICA_Unica | Superior | Criolla_Colombia     
## 4 B4    ICA_Unica | Superior | Parda_Pastusa        
## 5 B5    Parda_Pastusa | Superior | Pastusa          
## 6 B6    Diacol_Capiro | Pastusa | Parda_Pastusa

# Verificacion de balance: cada variedad debe repetirse r veces en total, y
# cada par de variedades debe coincidir en un bloque exactamente lambda veces.
cat("\nRepeticiones por variedad (debe ser r =", r, "para todas):\n")

## 
## Repeticiones por variedad (debe ser r = 3 para todas):

print(table(book$variedad))

## 
##          Pastusa    Diacol_Capiro Criolla_Colombia        ICA_Unica 
##                3                3                3                3 
##         Superior    Parda_Pastusa 
##                3                3

—————————————————————————

PUNTO 4. ANALISIS DESCRIPTIVO

—————————————————————————

# 4.1 Medias y desviaciones estandar por variedad (factor A)
resumen_variedad <- book %>%
  group_by(variedad) %>%
  summarise(media = mean(rendimiento),
            sd    = sd(rendimiento),
            n     = n(),
            .groups = "drop") %>%
  mutate(ee = sd / sqrt(n),                       # error estandar
         ic95 = qt(0.975, df = n - 1) * ee)        # semi-ancho IC 95%
print(resumen_variedad)

## # A tibble: 6 × 6
##   variedad         media    sd     n    ee  ic95
##   <fct>            <dbl> <dbl> <int> <dbl> <dbl>
## 1 Pastusa           26.7 2.46      3 1.42   6.12
## 2 Diacol_Capiro     23.6 0.978     3 0.564  2.43
## 3 Criolla_Colombia  28.7 1.98      3 1.14   4.92
## 4 ICA_Unica         25.2 0.704     3 0.406  1.75
## 5 Superior          23.1 1.75      3 1.01   4.35
## 6 Parda_Pastusa     25.9 2.26      3 1.31   5.61

# 4.2 Medias y desviaciones estandar por bloque
resumen_bloque <- book %>%
  group_by(block) %>%
  summarise(media = mean(rendimiento),
            sd    = sd(rendimiento),
            n     = n(),
            .groups = "drop")
print(resumen_bloque)

## # A tibble: 6 × 4
##   block media    sd     n
##   <fct> <dbl> <dbl> <int>
## 1 B1     27.5 4.00      3
## 2 B2     25.4 2.14      3
## 3 B3     25.6 2.11      3
## 4 B4     24.5 0.217     3
## 5 B5     23.2 1.88      3
## 6 B6     27.0 1.99      3

# 4.3 Grafico de medias por variedad con barras de error (IC 95%)
p_variedad <- ggplot(resumen_variedad, aes(x = reorder(variedad, -media), y = media)) +
  geom_col(fill = "#4C8C4A", width = 0.6) +
  geom_errorbar(aes(ymin = media - ic95, ymax = media + ic95), width = 0.2) +
  labs(title = "Rendimiento medio por variedad (IC 95%)",
       x = "Variedad", y = "Rendimiento (t/ha)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))
print(p_variedad)

# 4.4 Grafico de medias por bloque (evalua si el bloqueo fue efectivo:
# los bloques deben diferir claramente entre si)
p_bloque <- ggplot(resumen_bloque, aes(x = block, y = media)) +
  geom_col(fill = "#8C6D4C", width = 0.6) +
  labs(title = "Rendimiento medio por bloque",
       subtitle = "Diferencias marcadas entre bloques = bloqueo efectivo",
       x = "Bloque", y = "Rendimiento (t/ha)") +
  theme_minimal()
print(p_bloque)

# 4.5 Grafico de "perfiles" (= grafico de interaccion Variedad x Bloque).
# Este es EL MISMO grafico que se vuelve a usar mas abajo, en el Punto 5
# (Paso 1 del ANOVA), como la forma DESCRIPTIVA de observar la interaccion
# en un diseno donde el ANOVA no la puede testear formalmente (ver objeto
# p_interaccion). Se genera aqui tambien porque el taller lo pide como parte
# del analisis descriptivo.
p_perfil <- ggplot(book, aes(x = block, y = rendimiento, color = variedad, group = variedad)) +
  geom_point(size = 2) +
  geom_line(linewidth = 0.5, alpha = 0.7, na.rm = TRUE) +
  labs(title = "Perfiles de rendimiento por variedad a traves de los bloques",
       subtitle = "Lineas discontinuas = la variedad no esta presente en ese bloque",
       x = "Bloque", y = "Rendimiento (t/ha)", color = "Variedad") +
  theme_minimal()
print(p_perfil)

## ————————————————————————— ## PUNTO 5. ANALISIS DE VARIANZA (ANOVA) ## —————————————————————————

## --- Paso 1: modelo con interaccion Variedad x Bloque -----------------------
# En un DBIB cada combinacion variedad-bloque tiene, a lo sumo, UNA sola
# observacion (no hay repeticion dentro de celda) y ademas la mayoria de las
# combinaciones ni siquiera existen (variedad ausente de ese bloque). Esto
# significa que el termino de interaccion queda COMPLETAMENTE CONFUNDIDO con
# el error experimental: no quedan grados de libertad para estimarlo por
# separado. R lo evidencia generando NA en la fila de interaccion (terminos
# "aliased"/no estimables).
mod_interaccion <- aov(rendimiento ~ block + variedad + block:variedad, data = book)
cat("\n--- Intento de modelo con interaccion (para evidenciar la aliasing) ---\n")

## 
## --- Intento de modelo con interaccion (para evidenciar la aliasing) ---

print(summary(mod_interaccion))

##                Df Sum Sq Mean Sq
## block           5  37.86   7.573
## variedad        5  58.20  11.640
## block:variedad  7   6.95   0.994

# Interpretacion esperada: la fila block:variedad se queda con 0 grados de
# libertad de RESIDUALES (Residuals: 0 df), o R elimina coeficientes por
# singularidad ("... not defined because of singularities"). Aunque en la
# tabla aparezca una suma de cuadrados (SS) distinta de cero para
# block:variedad, ESO NO ES EVIDENCIA de interaccion real: sin grados de
# libertad de error no hay F ni p-valor calculable, por lo que esa fila NO
# se puede interpretar ni reportar como "significativa" o "no significativa".
# En los datos simulados de este script la interaccion es, ademas, realmente
# INEXISTENTE: el rendimiento se genero de forma puramente aditiva
# (mu + efecto_variedad + efecto_bloque + error, sin ningun termino cruzado
# variedad:bloque -- ver Punto 3.3). Cualquier numero que aparezca en esa
# fila del ANOVA es artefacto de la falta de grados de libertad, no una
# senal real de interaccion.
# CONCLUSION: en un diseno en bloques incompletos la aditividad
# bloque-tratamiento no es una hipotesis que se pueda contrastar con el ANOVA
# clasico; se ASUME como condicion de validez del diseno desde el momento en
# que se decide usarlo (ver "Punto de investigacion adicional" mas abajo).
# Por lo tanto se elimina el termino de interaccion y se pasa al modelo aditivo.

# Grafico DESCRIPTIVO de interaccion Variedad x Bloque (no es una prueba
# estadistica, es la unica forma valida de "ver" la interaccion en un DBIB
# ya que el ANOVA no la puede testear formalmente). Lineas paralelas =
# comportamiento aditivo esperado; lineas que se cruzan sugeririan
# interaccion, aunque sin soporte estadistico formal en este diseno.
p_interaccion <- ggplot(book, aes(x = block, y = rendimiento, color = variedad, group = variedad)) +
  geom_point(size = 2) +
  geom_line(linewidth = 0.5, alpha = 0.7, na.rm = TRUE) +
  labs(title = "Grafico de interaccion (descriptivo): Variedad x Bloque",
       subtitle = "No es una prueba formal: en un DBIB la interaccion no tiene grados de libertad de error para testearse",
       x = "Bloque", y = "Rendimiento (t/ha)", color = "Variedad") +
  theme_minimal()
print(p_interaccion)

## --- Paso 2: modelo final (aditivo, con analisis intra-bloque) --------------
# Para un DBIB la comparacion de variedades debe hacerse con las SUMAS DE
# CUADRADOS AJUSTADAS INTRA-BLOQUE (no las medias crudas), porque el diseno
# no es ortogonal. agricolae::BIB.test() implementa exactamente ese analisis:
# ajusta el efecto de variedad por bloque (y viceversa) usando la teoria de
# bloques incompletos balanceados.
modelo_bib <- agricolae::BIB.test(block = book$block,
                                   trt   = book$variedad,
                                   y     = book$rendimiento,
                                   test  = "tukey",
                                   group = TRUE)
# Nota: el argumento correcto en agricolae::BIB.test() se llama "test"
# (no "method"), y acepta "lsd","tukey","duncan","waller","snk". Tampoco
# existe un argumento "k": el tamano de bloque se calcula internamente a
# partir de la estructura de "block" y "trt".
cat("\n--- ANOVA ajustado del DBIB (intra-bloque, BIB.test) ---\n")

## 
## --- ANOVA ajustado del DBIB (intra-bloque, BIB.test) ---

print(modelo_bib$ANOVA)

## NULL

# Tabla ANOVA "clasica" equivalente para reportar en el informe (Tipo I:
# bloques ajustados primero, variedad ajustada por bloque = analisis
# intra-bloque aproximado). Se usa tambien para extraer los residuales para
# el chequeo de supuestos del Punto 6.
modelo_aditivo <- aov(rendimiento ~ block + variedad, data = book)
cat("\n--- Tabla ANOVA (Modelo aditivo: Bloque + Variedad) ---\n")

## 
## --- Tabla ANOVA (Modelo aditivo: Bloque + Variedad) ---

print(summary(modelo_aditivo))

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## block        5  37.86   7.573   7.622 0.00942 **
## variedad     5  58.20  11.640  11.716 0.00271 **
## Residuals    7   6.95   0.994                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

—————————————————————————

PUNTO 6. VERIFICACION DE SUPUESTOS (sobre el modelo final del Paso 2)

—————————————————————————

residuales <- residuals(modelo_aditivo)
ajustados  <- fitted(modelo_aditivo)

# 6.1 Normalidad de residuales
cat("\n--- Shapiro-Wilk (normalidad de residuales) ---\n")

## 
## --- Shapiro-Wilk (normalidad de residuales) ---

print(shapiro.test(residuales))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.93395, p-value = 0.2278

qqnorm(residuales, main = "Grafico Q-Q de residuales")
qqline(residuales, col = "red", lwd = 2)

# 6.2 Homogeneidad de varianzas entre variedades
cat("\n--- Bartlett (homogeneidad de varianzas) ---\n")

## 
## --- Bartlett (homogeneidad de varianzas) ---

print(bartlett.test(residuales ~ book$variedad))

## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuales by book$variedad
## Bartlett's K-squared = 2.9493, df = 5, p-value = 0.7078

cat("\n--- Levene (homogeneidad de varianzas, mas robusto que Bartlett) ---\n")

## 
## --- Levene (homogeneidad de varianzas, mas robusto que Bartlett) ---

print(car::leveneTest(residuales ~ book$variedad))

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  5  0.3203 0.8913
##       12

# 6.3 Independencia de residuales: residuales vs valores ajustados
df_diag <- data.frame(ajustados = ajustados, residuales = residuales)
p_diag <- ggplot(df_diag, aes(x = ajustados, y = residuales)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Residuales vs valores ajustados",
       x = "Valores ajustados", y = "Residuales") +
  theme_minimal()
print(p_diag)

# Interpretacion general (se completa segun los p-valores obtenidos al
# correr el script): p > 0.05 en Shapiro-Wilk indica que no hay evidencia
# para rechazar normalidad; p > 0.05 en Bartlett/Levene indica homogeneidad
# de varianzas; una nube de puntos sin patron en el grafico de residuales
# vs ajustados apoya independencia/homocedasticidad.


## ---------------------------------------------------------------------------
## PUNTO DE INVESTIGACION ADICIONAL: aditividad bloque x tratamiento en DBIB
## ---------------------------------------------------------------------------
cat("\n--- Punto de investigacion adicional (aditividad en DBIB) ---\n")

## 
## --- Punto de investigacion adicional (aditividad en DBIB) ---

cat("
En un DBIB, cada combinacion variedad-bloque aparece 0 o 1 vez, por lo que
la prueba clasica de un grado de libertad de Tukey para no-aditividad
(funcion agricolae::nonadditivity(), pensada para tablas de doble entrada
COMPLETAS con una observacion por celda) no es aplicable de forma directa,
porque requiere que todas las combinaciones existan.

Condiciones agronomicas donde SI podria esperarse una interaccion real
variedad x bloque en este caso: si los bloques representan gradientes
marcados de humedad o profundidad efectiva del suelo, y las variedades
difieren en su tolerancia a esas condiciones (p. ej. una variedad con
sistema radicular mas profundo que rinde mejor en suelos de baja humedad,
mientras otra depende de riego), el rendimiento relativo de las variedades
podria no ser consistente entre bloques (esto equivale a una interaccion
genotipo x ambiente).

Herramientas en R para explorar esto mas alla del ANOVA del Paso 1:
  - Modelos AMMI (Additive Main effects and Multiplicative Interaction),
    disponibles en agricolae::AMMI(), tratando cada bloque como un
    'ambiente': permiten visualizar y cuantificar interaccion
    genotipo-ambiente aun sin balance perfecto de celdas.
  - Diagnostico de Mandel (regresion de los efectos de tratamiento dentro
    de cada bloque contra el efecto de bloque) para detectar
    no-aditividad en tablas incompletas.
  - Transformaciones de Box-Cox (MASS::boxcox()) sobre la variable
    respuesta cuando se sospecha que la interaccion aparente se debe a
    una escala de medicion inadecuada, no a una interaccion biologica real.
Si se detecta interaccion real, el investigador debe: (1) reportar las
medias por variedad DENTRO de cada bloque en lugar de una unica media
general, (2) evaluar si conviene analizar los bloques como 'ambientes'
separados (enfoque multiambiente), y (3) revisar si la formacion de los
bloques fue realmente homogenea o si mezclo condiciones edaficas distintas.")

## 
## En un DBIB, cada combinacion variedad-bloque aparece 0 o 1 vez, por lo que
## la prueba clasica de un grado de libertad de Tukey para no-aditividad
## (funcion agricolae::nonadditivity(), pensada para tablas de doble entrada
## COMPLETAS con una observacion por celda) no es aplicable de forma directa,
## porque requiere que todas las combinaciones existan.
## 
## Condiciones agronomicas donde SI podria esperarse una interaccion real
## variedad x bloque en este caso: si los bloques representan gradientes
## marcados de humedad o profundidad efectiva del suelo, y las variedades
## difieren en su tolerancia a esas condiciones (p. ej. una variedad con
## sistema radicular mas profundo que rinde mejor en suelos de baja humedad,
## mientras otra depende de riego), el rendimiento relativo de las variedades
## podria no ser consistente entre bloques (esto equivale a una interaccion
## genotipo x ambiente).
## 
## Herramientas en R para explorar esto mas alla del ANOVA del Paso 1:
##   - Modelos AMMI (Additive Main effects and Multiplicative Interaction),
##     disponibles en agricolae::AMMI(), tratando cada bloque como un
##     'ambiente': permiten visualizar y cuantificar interaccion
##     genotipo-ambiente aun sin balance perfecto de celdas.
##   - Diagnostico de Mandel (regresion de los efectos de tratamiento dentro
##     de cada bloque contra el efecto de bloque) para detectar
##     no-aditividad en tablas incompletas.
##   - Transformaciones de Box-Cox (MASS::boxcox()) sobre la variable
##     respuesta cuando se sospecha que la interaccion aparente se debe a
##     una escala de medicion inadecuada, no a una interaccion biologica real.
## Si se detecta interaccion real, el investigador debe: (1) reportar las
## medias por variedad DENTRO de cada bloque en lugar de una unica media
## general, (2) evaluar si conviene analizar los bloques como 'ambientes'
## separados (enfoque multiambiente), y (3) revisar si la formacion de los
## bloques fue realmente homogenea o si mezclo condiciones edaficas distintas.

—————————————————————————

PUNTO 7. COMPARACION DE MEDIAS

—————————————————————————

# Justificacion: hay 6 variedades (multiples comparaciones posibles, C(6,2)=15
# pares) y NINGUNA de ellas es un testigo/control frente al cual comparar
# las demas (todas son variedades comerciales/lineas de interes en pie de
# igualdad) -> se descarta Dunnett. Con 15 comparaciones se prefiere
# controlar la tasa de error familiar -> se elige TUKEY HSD sobre Fisher-LSD
# (que no controla el error tipo I por comparacion multiple). Ademas
# BIB.test() con method = "tukey" ya calcula las medias AJUSTADAS
# intra-bloque (no las medias crudas), que es lo correcto en un diseno no
# ortogonal como este.
cat("\n--- Comparacion de medias ajustadas (Tukey, intra-bloque) ---\n")

## 
## --- Comparacion de medias ajustadas (Tukey, intra-bloque) ---

print(modelo_bib$means)

##                  book$rendimiento mean.adj        SE r       std   Min   Max
## Pastusa                  26.65333 26.48306 0.6325855 3 2.4643322 27.13 30.91
## Diacol_Capiro            23.60333 21.76639 0.6325855 3 0.9775650 22.98 24.73
## Criolla_Colombia         28.67667 28.65667 0.6325855 3 1.9813463 24.56 25.94
## ICA_Unica                25.17000 25.53722 0.6325855 3 0.7037755 24.63 28.55
## Superior                 23.06333 23.80806 0.6325855 3 1.7509236 23.86 28.52
## Parda_Pastusa            25.94000 26.85528 0.6325855 3 2.2603318 21.05 24.23
##                     Q25   Q50    Q75
## Pastusa          27.560 27.99 29.450
## Diacol_Capiro    23.040 23.10 23.915
## Criolla_Colombia 24.785 25.01 25.475
## ICA_Unica        24.635 24.64 26.595
## Superior         25.720 27.58 28.050
## Parda_Pastusa    22.480 23.91 24.070

print(modelo_bib$groups)

##                  book$rendimiento groups
## Criolla_Colombia         28.65667      a
## Parda_Pastusa            26.85528     ab
## Pastusa                  26.48306     ab
## ICA_Unica                25.53722     ab
## Superior                 23.80806     bc
## Diacol_Capiro            21.76639      c

# Grafico de barras con letras de agrupacion
grupos <- modelo_bib$groups
grupos$variedad <- rownames(grupos)
colnames(grupos)[1] <- "media_ajustada"

p_grupos <- ggplot(grupos, aes(x = reorder(variedad, -media_ajustada), y = media_ajustada)) +
  geom_col(fill = "#4C7A9C", width = 0.6) +
  geom_text(aes(label = groups, y = media_ajustada + 0.5), size = 5) +
  labs(title = "Medias ajustadas por variedad (Tukey HSD, intra-bloque)",
       x = "Variedad", y = "Rendimiento ajustado (t/ha)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))
print(p_grupos)

# --- Eficiencia relativa del bloqueo (ER) frente a un DCA ------------------
# Formula clasica (Cochran & Cox, 1957) para bloques incompletos balanceados:
#   ER = [ (r-1)*MSB + r*(k-1)*MSE ] / [ (r*k-1)*MSE ]
# donde MSB es el cuadrado medio de bloques SIN ajustar por tratamiento y
# MSE es el cuadrado medio del error intra-bloque del modelo aditivo.
mod_bloque_solo <- aov(rendimiento ~ block, data = book)
MSB <- summary(mod_bloque_solo)[[1]]["block", "Mean Sq"]
MSE <- summary(modelo_aditivo)[[1]]["Residuals", "Mean Sq"]
ER  <- ((r - 1) * MSB + r * (k - 1) * MSE) / ((r * k - 1) * MSE)
cat("\nEficiencia relativa del bloqueo (DBIB vs DCA hipotetico): ER =",
    round(ER, 2), "\n")

## 
## Eficiencia relativa del bloqueo (DBIB vs DCA hipotetico): ER = 2.66

cat(ifelse(ER > 1,
           "ER > 1: el bloqueo fue beneficial, redujo el error experimental.\n",
           "ER < 1: el bloqueo no aporto precision frente a un DCA (posible sobrebloqueo).\n"))

## ER > 1: el bloqueo fue beneficial, redujo el error experimental.

—————————————————————————

PUNTO 8. INTERPRETACION ESTADISTICA (parrafo generado con base en los

resultados obtenidos arriba; revisar y ajustar la redaccion final)

—————————————————————————

tabla_anova <- summary(modelo_aditivo)[[1]]
p_variedad_val <- tabla_anova["variedad", "Pr(>F)"]
p_bloque_val   <- tabla_anova["block", "Pr(>F)"]
alpha_sig <- ifelse(p_variedad_val < 0.01, "0.01",
                     ifelse(p_variedad_val < 0.05, "0.05", "no significativo (p > 0.05)"))

cat("\n--- PARRAFO DE INTERPRETACION ESTADISTICA (borrador automatico) ---\n")

## 
## --- PARRAFO DE INTERPRETACION ESTADISTICA (borrador automatico) ---

# Nota: el parrafo se construye como UNA sola linea logica (con "\n" como
# caracter de salto de linea dentro del texto, no como salto de linea real
# del archivo) para evitar el error "Incomplete expression" que aparece en
# la consola de R si se ejecuta o se selecciona el bloque por partes.
sig_variedad_txt <- ifelse(p_variedad_val < 0.05, sprintf("resultando significativo al nivel alpha = %s", alpha_sig), "resultando no significativo al nivel alpha = 0.05")
sig_bloque_txt   <- ifelse(p_bloque_val < 0.05, "el bloque SI explica", "el bloque NO explica claramente")
efecto_ER_txt    <- ifelse(ER > 1, "mejoro", "no mejoro")

parrafo_estadistico <- sprintf("El analisis de varianza intra-bloque del diseno en bloques incompletos balanceados mostro un efecto de variedad (factor A) con un valor p = %.4f, %s. El efecto de bloque presento un valor p = %.4f, lo que indica que %s una fuente de variacion relevante para el rendimiento, apoyando la decision de usar un diseno en bloques en lugar de un diseno completamente al azar. Dado que la interaccion variedad x bloque no pudo estimarse de forma independiente por la falta de grados de libertad propia de los disenos incompletos (ver Paso 1 del ANOVA), la validez de las comparaciones entre variedades se apoya en el supuesto de aditividad, discutido en el punto de investigacion adicional. La eficiencia relativa del bloqueo fue ER = %.2f, lo que indica que el bloqueo %s la precision del experimento frente a un diseno completamente al azar equivalente.\n", p_variedad_val, sig_variedad_txt, p_bloque_val, sig_bloque_txt, ER, efecto_ER_txt)
cat(parrafo_estadistico)

## El analisis de varianza intra-bloque del diseno en bloques incompletos balanceados mostro un efecto de variedad (factor A) con un valor p = 0.0027, resultando significativo al nivel alpha = 0.01. El efecto de bloque presento un valor p = 0.0094, lo que indica que el bloque SI explica una fuente de variacion relevante para el rendimiento, apoyando la decision de usar un diseno en bloques en lugar de un diseno completamente al azar. Dado que la interaccion variedad x bloque no pudo estimarse de forma independiente por la falta de grados de libertad propia de los disenos incompletos (ver Paso 1 del ANOVA), la validez de las comparaciones entre variedades se apoya en el supuesto de aditividad, discutido en el punto de investigacion adicional. La eficiencia relativa del bloqueo fue ER = 2.66, lo que indica que el bloqueo mejoro la precision del experimento frente a un diseno completamente al azar equivalente.

## ---------------------------------------------------------------------------
## PUNTO 9. INTERPRETACION AGRONOMICA  (parrafo generado con base en los
## resultados; revisar y ajustar la redaccion final)
## ---------------------------------------------------------------------------
mejor_variedad <- grupos$variedad[which.max(grupos$media_ajustada)]
rendimiento_max <- round(max(grupos$media_ajustada), 2)

cat("\n--- PARRAFO DE INTERPRETACION AGRONOMICA (borrador automatico) ---\n")

## 
## --- PARRAFO DE INTERPRETACION AGRONOMICA (borrador automatico) ---

parrafo_agronomico <- sprintf("Desde el punto de vista agronomico, la variedad %s presento el mayor rendimiento ajustado (%.2f t/ha), por lo que, de confirmarse con datos reales, seria la variedad recomendada para las condiciones edaficas del lote evaluado. El uso del diseno en bloques incompletos permitio controlar la heterogeneidad del terreno (fertilidad, humedad y profundidad efectiva del suelo) sin necesidad de que las 6 variedades coincidieran en cada bloque, lo cual habria sido logisticamente inviable con solo 3 unidades experimentales por bloque. La consistencia de la respuesta de las variedades a traves de los bloques no pudo verificarse formalmente por la falta de grados de libertad para la interaccion (Paso 1 del ANOVA); se recomienda, en un estudio real, aumentar el numero de bloques o repetir el ensayo en distintas localidades/campanas para evaluar la estabilidad de la variedad superior antes de una recomendacion de manejo definitiva.\n", mejor_variedad, rendimiento_max)
cat(parrafo_agronomico)

## Desde el punto de vista agronomico, la variedad Criolla_Colombia presento el mayor rendimiento ajustado (28.66 t/ha), por lo que, de confirmarse con datos reales, seria la variedad recomendada para las condiciones edaficas del lote evaluado. El uso del diseno en bloques incompletos permitio controlar la heterogeneidad del terreno (fertilidad, humedad y profundidad efectiva del suelo) sin necesidad de que las 6 variedades coincidieran en cada bloque, lo cual habria sido logisticamente inviable con solo 3 unidades experimentales por bloque. La consistencia de la respuesta de las variedades a traves de los bloques no pudo verificarse formalmente por la falta de grados de libertad para la interaccion (Paso 1 del ANOVA); se recomienda, en un estudio real, aumentar el numero de bloques o repetir el ensayo en distintas localidades/campanas para evaluar la estabilidad de la variedad superior antes de una recomendacion de manejo definitiva.