Diseño, Análisis Estadístico y Selección Multivariada en Cannabis F2

Aplicación del Modelo Cotes & Ñústez para 30 Líneas Segregantes

Author

Maria Jose Gutierrez - Juan Jose Bejarano - Silvia Valeria Avila Alferez

Code

# ==========================================
# 0. CARGA DE LIBRERÍAS
# ==========================================
library(tidyverse)
library(agricolae)
library(ggrepel)
library(DT)
library(plotly)
library(knitr)

set.seed(123)

SECCIÓN 1: Contexto Genético y Diseño Experimental

Explicación Metodológica

En esta primera sección definimos la estructura experimental del ensayo. Debido a que las 30 líneas segregantes $F_2$ son materiales genéticos nuevos en fase de evaluación inicial, solo contamos con una repetición por genotipo ($N=1$). Para poder evaluar y comparar estos materiales de forma estadísticamente válida frente a las variaciones del entorno (microclima, gradientes de luz o fertirriego), implementamos un Diseño Aumentado de Bloques Completos al Azar (DBCA).

Se incluyen 3 testigos comerciales estables (BocadilloEstable, DurbanPoison y AfganoRef), los cuales sí se repiten a lo largo de los 6 bloques de la sala. A partir de esta disposición, simulamos 6 variables agronómicas y organolépticas divididas en dos categorías:

Code

diccionario_vars <- data.frame(
  Variable = c("rend_flor_calA", "rend_flor_calB", "potencia_thc", "tricomas_score", "dias_floracion", "incidencia_botrytis"),
  Nombre = c("Rendimiento Flor A", "Rendimiento Flor B", "Potencia THC", "Densidad Tricomas", "Días a Floración", "Incidencia Botrytis"),
  Unidad = c("Gramos (g)", "Gramos (g)", "Porcentaje (%)", "Escala 1-5", "Días", "Escala 1-5"),
  Objetivo = c("Maximizar", "Maximizar", "Maximizar", "Maximizar", "Minimizar", "Minimizar"),
  Peso_ISG = c("+15", "+5", "+12", "+8", "-5", "-5")
)

datatable(
  diccionario_vars,
  caption = "Resumen de Variables Evaluadas en la Matriz Multivariada",
  options = list(dom = 't', scrollX = TRUE),
  rownames = FALSE
)

Code

# ==========================================
# 1. DISEÑO AUMENTADO BAJO BLOQUES (MODELO ÑÚSTEZ)
# ==========================================
testigos <- c("BocadilloEstable", "DurbanPoison", "AfganoRef") # g = 3
tratamientos_nuevos <- paste0("F2_", sprintf("%02d", 1:30))   # N = 30

# Construcción formal del DBCA Aumentado (6 bloques)
disenio <- design.dau(trt1 = testigos, trt2 = tratamientos_nuevos, r = 6, serie = 2)
campo <- disenio$book

# Identificación de Testigos (check = 1) vs Nuevos (check = 0)
campo <- campo %>%
  mutate(
    block = as.factor(block),
    check = ifelse(trt %in% testigos, 1, 0)
  )

# ==========================================
# 2. GENERACIÓN DE LAS 6 VARIABLES EN EL DISEÑO
# ==========================================
campo <- campo %>%
  group_by(check) %>%
  mutate(
    # Directas (Buscamos maximizar)
    rend_flor_calA      = rnorm(n(), mean = ifelse(first(check) == 1, 400, 330), sd = ifelse(first(check) == 1, 12, 30)),
    rend_flor_calB      = rnorm(n(), mean = ifelse(first(check) == 1, 100, 85),  sd = ifelse(first(check) == 1, 5, 15)),
    potencia_thc        = rnorm(n(), mean = ifelse(first(check) == 1, 22, 19),   sd = ifelse(first(check) == 1, 1, 3)),
    tricomas_score      = sample(1:5, n(), replace = TRUE, prob = if (first(check) == 1) c(0.05, 0.1, 0.2, 0.35, 0.3) else c(0.2, 0.2, 0.2, 0.2, 0.2)),
    # Invertidas (Buscamos minimizar)
    dias_floracion      = rnorm(n(), mean = ifelse(first(check) == 1, 60, 65),   sd = ifelse(first(check) == 1, 2, 6)),
    incidencia_botrytis = sample(1:5, n(), replace = TRUE, prob = if (first(check) == 1) c(0.4, 0.3, 0.2, 0.1, 0) else c(0.2, 0.2, 0.2, 0.2, 0.2))
  ) %>%
  ungroup()

SECCIÓN 2: Factores de Corrección ($FC$) e Índices de Selección ($IS$)

Explicación Metodológica

Para evitar seleccionar “falsos positivos” (plantas que rinden mucho solo porque les tocó un bloque con mejores condiciones de luz o aireación), aplicamos el modelo de Cotes & Ñústez (2001):

Tabla 1 (DBCA Factor Corrección $FC$): Se calcula el promedio del rendimiento de los testigos en cada bloque respecto al promedio global de los testigos. La diferencia resultante es el Factor de Corrección ($FC$) de ese bloque. Si un bloque tuvo condiciones ambientales superiores, su $FC$ será positivo y se le restará al valor observado de las líneas $F_2$ situadas en dicho bloque para “limpiar” el efecto ambiental.

Code

# ==========================================
# TABLA 1: CROQUIS DE BLOQUES Y FACTOR DE CORRECCIÓN (FC)
# ==========================================
tabla_1_fc <- campo %>%
  filter(check == 1) %>%
  group_by(block, trt) %>%
  summarise(valor = mean(rend_flor_calA), .groups = "drop") %>%
  pivot_wider(names_from = trt, values_from = valor) %>%
  rowwise() %>%
  mutate(
    Total_Bloque = sum(c_across(-block)),
    Media_Bloque = mean(c_across(-c(block, Total_Bloque)))
  ) %>%
  ungroup() %>%
  mutate(
    Media_Global = mean(Media_Bloque),
    FC = round(Media_Bloque - Media_Global, 2)
  )

tabla_1_fc %>%
  mutate(across(where(is.numeric), ~ round(., 2))) %>%
  datatable(
    options = list(pageLength = 6, scrollX = TRUE),
    caption = "Tabla 1: DBCA Factor Corrección FC por Bloque"
  )

Tabla 2 (DBCA Índices de Selección $IS$): A través de un modelo lineal sobre los testigos, estimamos la Diferencia Mínima Significativa ($DMS$). Con este parámetro, se le asigna a cada variable un Índice de Selección individual ($IS = R + R'$), donde $R$ es el ranking numérico de la planta y $R'$ representa la distancia en múltiplos de $DMS$ frente al promedio de los testigos comerciales.

Code

# ==========================================
# FUNCIÓN ESTIMADORA DEL IS (COTES & ÑÚSTEZ)
# ==========================================
calcular_IS_variable <- function(df, var_name, alpha = 0.10, r = 6, g = 3) {
  fc_df <- df %>%
    filter(check == 1) %>%
    group_by(block) %>%
    summarise(mean_b = mean(.data[[var_name]])) %>%
    mutate(FC = mean_b - mean(mean_b))
  
  nuevos <- df %>%
    filter(check == 0) %>%
    left_join(fc_df %>% select(block, FC), by = "block") %>%
    mutate(valor_ajustado = .data[[var_name]] - FC)
  
  fit <- lm(as.formula(paste(var_name, "~ block + trt")), data = filter(df, check == 1))
  cme <- deviance(fit) / df.residual(fit)
  gle <- df.residual(fit)
  t_val <- qt(1 - (alpha/2), df = gle)
  dms <- t_val * sqrt(((r - 1)*(g - 1)*cme)/(r*g))
  prom_testigos <- mean(filter(df, check == 1)[[var_name]])
  
  res <- nuevos %>%
    arrange(valor_ajustado) %>%
    mutate(
      R = row_number(),
      dif = valor_ajustado - prom_testigos,
      R_prime = case_when(
        abs(dif) <= 1*dms ~ 1,
        abs(dif) <= 2*dms ~ 2,
        abs(dif) <= 3*dms ~ 3,
        TRUE ~ ceiling(abs(dif)/dms)
      ) * sign(dif),
      IS = R + R_prime
    )
  return(list(data = res %>% select(trt, block, valor_ajustado, IS), dms = dms, prom_t = prom_testigos))
}

# ==========================================
# TABLA 2: MATRIZ DE ÍNDICES DE SELECCIÓN (IS)
# ==========================================
res_v1 <- calcular_IS_variable(campo, "rend_flor_calA")
res_v2 <- calcular_IS_variable(campo, "rend_flor_calB")
res_v3 <- calcular_IS_variable(campo, "potencia_thc")
res_v4 <- calcular_IS_variable(campo, "tricomas_score")
res_v5 <- calcular_IS_variable(campo, "dias_floracion")
res_v6 <- calcular_IS_variable(campo, "incidencia_botrytis")

tabla_2_is <- res_v1$data %>% select(trt, block, rend_ajustado = valor_ajustado, IS_FlorA = IS) %>%
  left_join(res_v2$data %>% select(trt, block, IS_FlorB = IS), by = c("trt", "block")) %>%
  left_join(res_v3$data %>% select(trt, block, IS_THC = IS), by = c("trt", "block")) %>%
  left_join(res_v4$data %>% select(trt, block, IS_Tricomas = IS), by = c("trt", "block")) %>%
  left_join(res_v5$data %>% select(trt, block, IS_Dias = IS), by = c("trt", "block")) %>%
  left_join(res_v6$data %>% select(trt, block, IS_Botrytis = IS), by = c("trt", "block"))

tabla_2_is %>%
  mutate(across(where(is.numeric), ~ round(., 2))) %>%
  datatable(
    options = list(pageLength = 10, scrollX = TRUE),
    caption = "Tabla 2: DBCA Índices de Selección IS por Variable"
  )

SECCIÓN 3: Matriz Consolidada y Selección Final ($ISG$)

Explicación Metodológica

Con los índices $IS$ calculados para cada variable individual, consolidamos el Índice de Selección General ($ISG$). Este índice aplica una ponderación económica/agronómica según las prioridades del programa de fitomejoramiento:

\[\text{Puntos Positivos} = (15 \times IS_{\text{FlorA}}) + (5 \times IS_{\text{FlorB}}) + (12 \times IS_{\text{THC}}) + (8 \times IS_{\text{Tricomas}})\] \[\text{Puntos Negativos} = (5 \times IS_{\text{Días}}) + (5 \times IS_{\text{Botrytis}})\] \[\text{ISG} = \text{Puntos Positivos} - \text{Puntos Negativos}\]

Finalmente, se calcula el promedio global del $ISG$ en el ensayo. Aquellas líneas $F_2$ cuyo $ISG$ supere el promedio son marcadas como “SÍ (Clonar)” para avanzar a la sala de madres; las que quedan por debajo se marcan como “NO (Descartar)”.

Code

# ==========================================
#  TABLA 3: MATRIZ CONSOLIDADA ISG Y DECISIÓN DE SELECCIÓN
# ==========================================
tabla_3_isg <- tabla_2_is %>%
  mutate(
    Puntos_Positivos = (15 * IS_FlorA) + (5 * IS_FlorB) + (12 * IS_THC) + (8 * IS_Tricomas),
    Puntos_Negativos = (5 * IS_Dias) + (5 * IS_Botrytis),
    ISG = Puntos_Positivos - Puntos_Negativos,
    Seleccionado = ifelse(ISG > mean(ISG), "SÍ (Clonar)", "NO (Descartar)")
  ) %>%
  arrange(desc(ISG))

promedio_ISG <- mean(tabla_3_isg$ISG)

tabla_3_isg %>%
  mutate(across(where(is.numeric), ~ round(., 2))) %>%
  datatable(
    options = list(pageLength = 10, scrollX = TRUE),
    caption = "Tabla 3: DBCA Selección Final ISG"
  )

SECCIÓN 4: Análisis Gráfico y Evaluaciones

Explicación Metodológica

Esta sección traduce la matriz de datos numéricos a una serie de visualizaciones clave para la toma de decisiones:

Boxplot 1: Distribución del Rendimiento por Bloque

Evalúa la heterogeneidad ambiental observando cómo varió el rendimiento promedio entre los 6 bloques de la sala.

Code

p_box1 <- ggplot(campo, aes(x = block, y = rend_flor_calA, fill = block)) +
  geom_boxplot(alpha = 0.7, outlier.color = "red") +
  geom_jitter(aes(shape = factor(check)), width = 0.2, alpha = 0.6, size = 2) +
  scale_shape_manual(values = c("0" = 16, "1" = 17), labels = c("0" = "Línea F2", "1" = "Testigo")) +
  scale_fill_brewer(palette = "Set3") +
  labs(
    title = "Boxplot 1: Distribución del Rendimiento por Bloque",
    subtitle = "Muestra la variabilidad microclimática entre bloques en la sala (Puntos = Genotipos)",
    x = "Bloque Experimental",
    y = "Rendimiento Flor A (g)",
    shape = "Tipo de Genotipo",
    fill = "Bloque"
  ) +
  theme_minimal()

ggplotly(p_box1)

Boxplot 2: Comparativa Testigos vs. Población F2

Compara la distribución global de la población $F_2$ frente a los testigos comerciales en variables críticas, verificando si existe segregación transgresiva (plantas superiores a los mejores padres).

Code

df_box_comparacion <- campo %>%
  select(trt, check, Rendimiento = rend_flor_calA, THC = potencia_thc, Botrytis = incidencia_botrytis) %>%
  pivot_longer(cols = c("Rendimiento", "THC", "Botrytis"), names_to = "Variable", values_to = "Valor") %>%
  mutate(Tipo = ifelse(check == 1, "Testigos Comerciales", "Líneas F2 Candidateadas"))

p_box2 <- ggplot(df_box_comparacion, aes(x = Tipo, y = Valor, fill = Tipo)) +
  geom_boxplot(alpha = 0.6, outlier.shape = NA) +
  geom_jitter(width = 0.15, alpha = 0.5, size = 1.8) +
  facet_wrap(~ Variable, scales = "free_y") +
  scale_fill_manual(values = c("Testigos Comerciales" = "#1976d2", "Líneas F2 Candidateadas" = "#388e3c")) +
  labs(
    title = "Boxplot 2: Comparativa Testigos vs. Población F2",
    subtitle = "Evaluación del rango de segregación en variables clave de cannabis",
    x = "",
    y = "Valor de la Variable"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

p_box2

Boxplot 3: Separación de Grupos por Índice ISG

Muestra la separación clara y sin solapamiento de los puntajes $ISG$ entre el grupo seleccionado y el descartado.

Code

p_box3 <- ggplot(tabla_3_isg, aes(x = Seleccionado, y = ISG, fill = Seleccionado)) +
  geom_boxplot(alpha = 0.7, width = 0.4) +
  geom_jitter(width = 0.1, size = 3, alpha = 0.8) +
  geom_hline(yintercept = promedio_ISG, linetype = "dashed", color = "black") +
  scale_fill_manual(values = c("SÍ (Clonar)" = "#2e7d32", "NO (Descartar)" = "#c62828")) +
  labs(
    title = "Boxplot 3: Separación de Grupos por Índice ISG",
    subtitle = "Compara la dispersión del índice final entre plantas seleccionadas y descartadas",
    x = "Decisión de Selección",
    y = "Índice de Selección General (ISG)"
  ) +
  theme_light()

ggplotly(p_box3)

Gráfico 1: Índice de Selección General Multivariado (ISG - DBCA)

Muestra el ranking ordenado de genotipos según su índice $ISG$ respecto a la línea de corte del umbral promedio.

Code

p1 <- ggplot(tabla_3_isg, aes(x = reorder(trt, ISG), y = ISG, fill = Seleccionado)) +
  geom_col(show.legend = TRUE) +
  geom_hline(yintercept = promedio_ISG, linetype = "dashed", color = "black", size = 1) +
  scale_fill_manual(values = c("SÍ (Clonar)" = "#2e7d32", "NO (Descartar)" = "#c62828")) +
  labs(title = "Gráfico 1: Índice de Selección General Multivariado (ISG - DBCA)",
       subtitle = paste("Umbral Promedio del Ensayo =", round(promedio_ISG, 1)),
       x = "Genotipos F2 Cannabis", y = "Índice ISG") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 8))

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Code

ggplotly(p1)

Gráfico 2: Ajuste por Microclima de Bloque (Método Ñústez)

Pone en evidencia el efecto de la corrección por microclima (método Ñústez), comparando el rendimiento observado (crudo) contra el ajustado (genético real).

Code

df_g2 <- campo %>%
  filter(check == 0) %>%
  select(trt, block, Observado = rend_flor_calA) %>%
  left_join(tabla_2_is %>% select(trt, block, Ajustado = rend_ajustado), by = c("trt", "block")) %>%
  pivot_longer(cols = c("Observado", "Ajustado"), names_to = "Tipo", values_to = "Rendimiento")

p2 <- ggplot(df_g2, aes(x = reorder(trt, Rendimiento), y = Rendimiento, color = Tipo, group = Tipo)) +
  geom_point(size = 3) + geom_line(alpha = 0.5) +
  scale_color_manual(values = c("Observado" = "#fb8c00", "Ajustado" = "#1565c0")) +
  labs(title = "Gráfico 2: Ajuste por Microclima de Bloque (Método Ñústez)",
       x = "Genotipos F2 Cannabis", y = "Rendimiento Flor A (g)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 8))

ggplotly(p2)

Gráfico 3: Matriz de Selección Multivariada

Un plano cartesiano de dispersión multivariada que cruza el Rendimiento Ajustado vs. el $ISG$ Final para ubicar estratégicamente cada código de planta.

Code

p3 <- ggplot(tabla_3_isg, aes(x = rend_ajustado, y = ISG)) +
  geom_point(aes(color = Seleccionado), size = 4) +
  geom_text_repel(aes(label = trt), size = 3) +
  geom_vline(xintercept = res_v1$prom_t, linetype = "dotted", color = "gray30") +
  geom_hline(yintercept = promedio_ISG, linetype = "dashed", color = "black") +
  scale_color_manual(values = c("SÍ (Clonar)" = "#2e7d32", "NO (Descartar)" = "#757575")) +
  labs(title = "Gráfico 3: Matriz de Selección Multivariada",
       x = "Rendimiento Flor A Ajustado (g)", y = "Índice ISG") +
  theme_light()

p3

SECCIÓN 5: Resumen de Genotipos Seleccionados

Explicación Metodológica

Esta sección final consolida los resultados en una tabla ejecutiva que resume el Top 5 de líneas $F_2$ élite. Estas plantas han demostrado un rendimiento ajustado superior, alto contenido de cannabinoides y un puntaje de selección multivariado que justifica su paso a la Sala de Madres para clonación masiva y conservación genotipica.

Top 5 Genotipos F2 Seleccionados por ISG para Clonación
Genotipo F2	Rend. Ajustado (g)	IS Flor A	IS THC	Índice ISG Total
F2_05	396.2107	29	25	775
F2_19	342.4488	18	32	573
F2_16	338.7824	16	17	554
F2_03	369.0037	25	5	552
F2_28	328.8163	10	28	542

SECCIÓN 6: Discusión Crítica y Propuestas de Mejora

1. Discusión Crítica sobre el Método Estadístico Usado

La metodología de Cotes & Ñústez (2001) implementada a través del Diseño Aumentado (DBCA) ha demostrado ser una herramienta sumamente eficiente para programas de mejoramiento en fases tempranas ($F_2$). Permite evaluar una población segregante amplia sin requerir múltiples réplicas vegetativas por individuo, optimizando el área física y los costos operativos.

Sin embargo, el método presenta limitaciones metodológicas claras:

Falta de repeticiones en genotipos evaluados ($N=1$): Al no tener repeticiones de cada planta $F_2$, el cálculo del error experimental depende exclusivamente de los testigos comerciales. Si una planta $F_2$ sufre un microestrés localizado (daño mecánico, obstrucción de gotero), su dato no puede ser contrastado con otra réplica idéntica.
Ponderación arbitraria de la matriz de selección: Los pesos asignados ($15, 12, 8, 5$) dependen del criterio subjetivo del fitomejorador y no de un parámetro económico de mercado rigurosamente modelado.
Asunción de homogeneidad lineal del bloque: El Factor de Corrección ($FC$) asume que el gradiente ambiental afecta exactamente de la misma manera a todas las líneas, lo cual no siempre ocurre debido a la interacción genotipo-ambiente ($G \times A$).

2. Propuestas de Mejora para Próximas Generaciones ($F_3$ / Clones)

Con el fin de incrementar la ganancia genética por ciclo de selección, se recomiendan las siguientes mejoras metodológicas para las fases posteriores:

Modelos Mixtos y BLUPs (REML/BLUP): Evolucionar del modelo lineal clásico basado en $FC$ hacia un enfoque de Modelos Lineales Mixtos mediante el paquete lme4 en R. Tratar los efectos de bloque como aleatorios permite calcular los Best Linear Unbiased Predictions (BLUPs), ajustando los valores fenotípicos por el fenómeno de contracción (shrinkage), lo que reduce drásticamente la tasa de falsos positivos.
Índice Económico de Smith-Hazel: Sustituir los ponderadores fijos del $ISG$ por un Índice de Selección Económico de Smith-Hazel, derivado de la matriz de varianzas y covarianzas genéticas junto con el valor marginal real del mercado ($/g THC, $/g Flor Calidad A).
Incorporación de Sensores y Variabilidad Ambiental Continua (ANCOVA): Implementar redes de sensores de microclima (Déficit de Presión de Vapor - VPD, Parámetros fotosintéticos PAR por bloque) para utilizar el entorno como una covariable continua en un Análisis de Covarianza (ANCOVA), en lugar de bloques discretos estáticos.
Micro-replicación por Enraizamiento Temprano: En la fase $F_3$, tomar 3 esquejes por individuo antes de inducir la floración. Esto permitirá realizar un DBCA tradicional con réplicas clonales ($r=3$), eliminando la limitación de $N=1$ y permitiendo estimar la heredabilidad en sentido amplio ($H^2$) para cada variable agronómica.

Diseño, Análisis Estadístico y Selección Multivariada en Cannabis F2

SECCIÓN 1: Contexto Genético y Diseño Experimental

Explicación Metodológica

SECCIÓN 2: Factores de Corrección (\(FC\)) e Índices de Selección (\(IS\))

Explicación Metodológica

SECCIÓN 3: Matriz Consolidada y Selección Final (\(ISG\))

Explicación Metodológica

SECCIÓN 4: Análisis Gráfico y Evaluaciones

Explicación Metodológica

Boxplot 1: Distribución del Rendimiento por Bloque

Boxplot 2: Comparativa Testigos vs. Población F2

Boxplot 3: Separación de Grupos por Índice ISG

Gráfico 1: Índice de Selección General Multivariado (ISG - DBCA)

Gráfico 2: Ajuste por Microclima de Bloque (Método Ñústez)

Gráfico 3: Matriz de Selección Multivariada

SECCIÓN 5: Resumen de Genotipos Seleccionados

Explicación Metodológica

SECCIÓN 6: Discusión Crítica y Propuestas de Mejora

1. Discusión Crítica sobre el Método Estadístico Usado

2. Propuestas de Mejora para Próximas Generaciones (\(F_3\) / Clones)