# ESTIMACIÓN DE TAC MEDIANTE MODELO LINEAL PONDERADO (WLS)
# Instalar paquetes requeridos si faltan
if(!require(tidyverse)) install.packages("tidyverse")
## Cargando paquete requerido: tidyverse
## Warning: package 'tidyverse' was built under R version 4.6.1
## Warning: package 'ggplot2' was built under R version 4.6.1
## Warning: package 'tibble' was built under R version 4.6.1
## Warning: package 'tidyr' was built under R version 4.6.1
## Warning: package 'readr' was built under R version 4.6.1
## Warning: package 'purrr' was built under R version 4.6.1
## Warning: package 'dplyr' was built under R version 4.6.1
## Warning: package 'stringr' was built under R version 4.6.1
## Warning: package 'forcats' was built under R version 4.6.1
## Warning: package 'lubridate' was built under R version 4.6.1
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
if(!require(emmeans)) install.packages("emmeans")
## Cargando paquete requerido: emmeans
## Warning: package 'emmeans' was built under R version 4.6.1
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
if(!require(multcompView)) install.packages("multcompView")
## Cargando paquete requerido: multcompView
## Warning: package 'multcompView' was built under R version 4.6.1
library(tidyverse)
library(emmeans)
library(multcompView)
# =========================================================================
# 1. IMPORTACIÓN Y CONFIGURACIÓN DE VARIABLES
# =========================================================================
datos_plantas <- read.csv("growth.csv", sep = ";", header = TRUE) %>%
rename(
Genotipo = any_of(c("Genotipo", "genotipo", "Variedad")),
Tratamiento = any_of(c("Tratamiento", "tratamiento", "Trat")),
Tiempo = any_of(c("Tiempo", "tiempo", "Dias")),
Biomasa = any_of(c("Biomasa", "biomasa", "Peso_Seco", "Peso"))
) %>%
drop_na(Genotipo, Tratamiento, Tiempo, Biomasa) %>%
mutate(
# Convertimos el Tiempo a factor (Categórica) como solicita el enfoque
Tiempo_Fact = as.factor(paste0("t", Tiempo)),
Genotipo = as.factor(Genotipo),
Tratamiento = as.factor(Tratamiento)
)
# Encontrar automáticamente delta t (t2 - t1 numérico)
tiempos_num <- sort(unique(datos_plantas$Tiempo))
dt <- tiempos_num[2] - tiempos_num[1]
# =========================================================================
# 2. CÁLCULO DE PESOS PONDERADOS PARA CORREGIR HETEROCEDASTICIDAD
# =========================================================================
# Calculamos la varianza de la biomasa para cada combinación en cada tiempo
varianzas_grupo <- datos_plantas %>%
group_by(Genotipo, Tratamiento, Tiempo_Fact) %>%
summarise(var_biomasa = var(Biomasa), .groups = 'drop')
# Acoplamos las varianzas a la base original y definimos el peso ponderado (1 / var)
datos_ponderados <- datos_plantas %>%
left_join(varianzas_grupo, by = c("Genotipo", "Tratamiento", "Tiempo_Fact")) %>%
mutate(
# Si la varianza es cero (raro), asignamos un peso neutro de 1
peso_wls = if_else(var_biomasa > 0, 1 / var_biomasa, 1)
)
# =========================================================================
# 3. AJUSTE DEL MODELO LINEAL PONDERADO (WLS) - CORREGIDO
# =========================================================================
# Usamos Tiempo_Fact (categórica) en el modelo lineal para cumplir tu enfoque.
# Mantenemos los pesos ponderados para controlar la varianza de t2.
modelo_wls <- lm(Biomasa ~ Genotipo * Tratamiento * Tiempo_Fact,
data = datos_ponderados,
weights = peso_wls)
print("=== RESUMEN DEL MODELO LINEAL PONDERADO (WLS) ===")
## [1] "=== RESUMEN DEL MODELO LINEAL PONDERADO (WLS) ==="
print(anova(modelo_wls))
## Analysis of Variance Table
##
## Response: Biomasa
## Df Sum Sq Mean Sq F value Pr(>F)
## Genotipo 3 120.475 40.158 40.1584 < 2.2e-16 ***
## Tratamiento 3 12.919 4.306 4.3064 0.0054448 **
## Tiempo_Fact 1 261.691 261.691 261.6912 < 2.2e-16 ***
## Genotipo:Tratamiento 9 29.512 3.279 3.2791 0.0008127 ***
## Genotipo:Tiempo_Fact 3 2.684 0.895 0.8948 0.4442676
## Tratamiento:Tiempo_Fact 3 6.274 2.091 2.0915 0.1015275
## Genotipo:Tratamiento:Tiempo_Fact 9 6.913 0.768 0.7681 0.6460655
## Residuals 283 283.000 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# =========================================================================
# 4. EXTRACCIÓN DE LA TAC Y COMPARACIONES POST-HOC (emmeans)
# =========================================================================
# 4A. Calculamos las medias marginales para cada combinación bajo el modelo WLS
medias_tiempo <- emmeans(modelo_wls, ~ Tiempo_Fact | Genotipo * Tratamiento)
# 4B. Calculamos la diferencia categórica (t2 - t1) para obtener el delta de biomasa
contrastes_tiempo <- pairs(medias_tiempo, reverse = TRUE)
# 4C. Convertimos el contraste a TAC dividiendo la tasa por el Delta t (dt)
grid_tac <- summary(contrastes_tiempo) %>%
mutate(
TAC_Promedio = estimate / dt,
Error_Estandar = SE / dt
) %>%
select(Genotipo, Tratamiento, TAC_Promedio, Error_Estandar)
print("=== ESTIMACIONES DE LA TAC (PENDIENTES DEL MODELO WLS) ===")
## [1] "=== ESTIMACIONES DE LA TAC (PENDIENTES DEL MODELO WLS) ==="
print(grid_tac, row.names = FALSE)
## Genotipo Tratamiento TAC_Promedio Error_Estandar
## CCN51 T1 0.16058333 0.07367467
## IMC67 T1 0.21462500 0.04583713
## TCS01 T1 0.29281667 0.04527419
## TCS19 T1 0.23100000 0.04915219
## CCN51 T2 0.14606667 0.03898202
## IMC67 T2 0.17170000 0.03622983
## TCS01 T2 0.20126667 0.03706014
## TCS19 T2 0.15094048 0.06521126
## CCN51 T3 0.14191667 0.05800835
## IMC67 T3 0.21361667 0.03618065
## TCS01 T3 0.16391667 0.04660867
## TCS19 T3 0.12333333 0.03589890
## CCN51 T4 0.17640000 0.04928972
## IMC67 T4 0.17905000 0.04137447
## TCS01 T4 0.09403333 0.05328677
## TCS19 T4 0.16414167 0.05029944
# =========================================================================
# 5. POST-HOC: TODOS LOS CONTRASTES ESPECÍFICOS POR GENOTIPO
# =========================================================================
# Calculamos las medias marginales de la TAC agrupadas por Genotipo y Tratamiento
medias_tac <- emmeans(modelo_wls, ~ Tratamiento | Genotipo)
## NOTE: Results may be misleading due to involvement in interactions
print("=== 5A. EJECUCIÓN DE TODOS LOS CONTRASTES ESPECÍFICOS DE HIPÓTESIS ===")
## [1] "=== 5A. EJECUCIÓN DE TODOS LOS CONTRASTES ESPECÍFICOS DE HIPÓTESIS ==="
# Definimos manualmente la lista expandida de contrastes requeridos
contrastes_especificos <- contrast(medias_tac, method = list(
"Efecto Sequía (T1 vs T3)" = c("T1" = 1, "T2" = 0, "T3" = -1, "T4" = 0),
"Efecto Sequía + Cd (T2 vs T4)" = c("T1" = 0, "T2" = 1, "T3" = 0, "T4" = -1),
"Efecto Cadmio (T1 vs T2)" = c("T1" = 1, "T2" = -1, "T3" = 0, "T4" = 0),
"Efecto Cadmio + Seq (T3 vs T4)"= c("T1" = 0, "T2" = 0, "T3" = 1, "T4" = -1)
), adjust = "bonferroni") # Ajuste estricto de Bonferroni para los 4 contrastes
# Convertimos a data.frame para extraer la significancia estadística de forma segura
df_contrastes <- as.data.frame(contrastes_especificos)
print(df_contrastes, row.names = FALSE)
## Genotipo = CCN51:
## contrast estimate SE df t.ratio p.value
## Efecto Sequía (T1 vs T3) -2.173000 2.813118 283 -0.772 1.0000
## Efecto Sequía + Cd (T2 vs T4) -1.323000 1.885250 283 -0.702 1.0000
## Efecto Cadmio (T1 vs T2) 7.404500 2.500560 283 2.961 0.0133
## Efecto Cadmio + Seq (T3 vs T4) 8.254500 2.283638 283 3.615 0.0014
##
## Genotipo = IMC67:
## contrast estimate SE df t.ratio p.value
## Efecto Sequía (T1 vs T3) -1.520250 1.751877 283 -0.868 1.0000
## Efecto Sequía + Cd (T2 vs T4) -1.310500 1.649849 283 -0.794 1.0000
## Efecto Cadmio (T1 vs T2) 2.680250 1.752792 283 1.529 0.5094
## Efecto Cadmio + Seq (T3 vs T4) 2.890000 1.648878 283 1.753 0.3229
##
## Genotipo = TCS01:
## contrast estimate SE df t.ratio p.value
## Efecto Sequía (T1 vs T3) 1.627000 1.949335 283 0.835 1.0000
## Efecto Sequía + Cd (T2 vs T4) 3.935000 1.947213 283 2.021 0.1770
## Efecto Cadmio (T1 vs T2) 0.068500 1.755245 283 0.039 1.0000
## Efecto Cadmio + Seq (T3 vs T4) 2.376500 2.123832 283 1.119 1.0000
##
## Genotipo = TCS19:
## contrast estimate SE df t.ratio p.value
## Efecto Sequía (T1 vs T3) -2.369000 1.825980 283 -1.297 0.7822
## Efecto Sequía + Cd (T2 vs T4) 2.484964 2.470686 283 1.006 1.0000
## Efecto Cadmio (T1 vs T2) -2.412214 2.449817 283 -0.985 1.0000
## Efecto Cadmio + Seq (T3 vs T4) 2.441750 1.853884 283 1.317 0.7555
##
## Results are averaged over the levels of: Tiempo_Fact
## P value adjustment: bonferroni method for 4 tests
# Convertimos los p-valores en asteriscos de significancia científica tradicional
df_asteriscos <- df_contrastes %>%
mutate(
Letra_Signif = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
TRUE ~ "ns"
),
# Homologamos exactamente el nombre del contraste para cruzarlo con el gráfico
Contraste_Nombre = case_when(
contrast == "Efecto Sequía (T1 vs T3)" ~ "1. Sequía sin Cadmio (T1 vs T3)",
contrast == "Efecto Sequía + Cd (T2 vs T4)" ~ "2. Sequía con Cadmio (T2 vs T4)",
contrast == "Efecto Cadmio (T1 vs T2)" ~ "3. Cadmio a 80% CC (T1 vs T2)",
contrast == "Efecto Cadmio + Seq (T3 vs T4)" ~ "4. Cadmio a 40% CC (T3 vs T4)"
)
)
# =========================================================================
# 6. CÁLCULO DE REDUCCIÓN PORCENTUAL PARA LOS 4 CONTRASTES
# =========================================================================
reduccion_porcentual <- grid_tac %>%
select(Genotipo, Tratamiento, TAC_Promedio) %>%
pivot_wider(names_from = Tratamiento, values_from = TAC_Promedio) %>%
mutate(
`1. Sequía sin Cadmio (T1 vs T3)` = ((T1 - T3) / T1) * 100,
`2. Sequía con Cadmio (T2 vs T4)` = ((T2 - T4) / T2) * 100,
`3. Cadmio a 80% CC (T1 vs T2)` = ((T1 - T2) / T1) * 100,
`4. Cadmio a 40% CC (T3 vs T4)` = ((T3 - T4) / T3) * 100
) %>%
select(Genotipo, starts_with("1"), starts_with("2"), starts_with("3"), starts_with("4")) %>%
pivot_longer(cols = -Genotipo, names_to = "Contraste_Nombre", values_to = "Reduccion_Porcentaje") %>%
# Acoplamos las etiquetas de significancia (asteriscos) calculadas en el Post-Hoc
left_join(df_asteriscos, by = c("Genotipo", "Contraste_Nombre"))
# =========================================================================
# 7. GRÁFICO FINAL EN CUADRÍCULA DE PANEL MATRICIAL (2x2)
# =========================================================================
ggplot(reduccion_porcentual, aes(x = Genotipo, y = Reduccion_Porcentaje, fill = Contraste_Nombre)) +
geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.6, color = "black", alpha = 0.85) +
# Añadimos los asteriscos controlando dinámicamente si la barra es positiva o negativa
geom_text(
aes(
label = Letra_Signif,
y = if_else(Reduccion_Porcentaje >= 0, Reduccion_Porcentaje + 3, Reduccion_Porcentaje - 5)
),
position = position_dodge(0.7),
vjust = if_else(reduccion_porcentual$Reduccion_Porcentaje >= 0, 0, 1),
fontface = "bold",
size = 3.8,
color = "black"
) +
geom_hline(yintercept = 0, linetype = "solid", color = "black", linewidth = 0.5) +
# Dividimos el lienzo en una matriz de 2x2 para evaluar los efectos de forma independiente
facet_wrap(~ Contraste_Nombre, nrow = 2, ncol = 2, scales = "free_y") +
labs(
title = "Análisis de Reducción de la TAC por Estresores Aislados e Interacciones",
subtitle = "Significancia de contrastes específicos calculada mediante WLS con ajuste Bonferroni\n(* p<0.05, ** p<0.01, *** p<0.001, ns = No Significativo)",
y = "Reducción de la Tasa de Crecimiento (%)",
x = "Genotipos Evaluados"
) +
scale_fill_manual(values = c(
"1. Sequía sin Cadmio (T1 vs T3)" = "#3a86ff", # Azul
"2. Sequía con Cadmio (T2 vs T4)" = "#8338ec", # Morado
"3. Cadmio a 80% CC (T1 vs T2)" = "#ff006e", # Rosa/Rojo
"4. Cadmio a 40% CC (T3 vs T4)" = "#fb5607" # Naranja
)) +
theme_bw() +
theme(
strip.text = element_text(face = "bold", size = 10),
axis.text.x = element_text(angle = 45, hjust = 1, face = "bold", size = 9),
legend.position = "none",
panel.grid.major.x = element_blank(),
plot.title = element_text(face = "bold", size = 13)
)
