LibrerÃas requeridas
library(openxlsx)
library(readxl)
library(knitr)
library(kableExtra)
library(tidyverse)
library(car)
library(dplyr)
library(broom)
library(writexl)
library(lmtest)
library(patchwork)
library(camcorder)
library(showtext)
library(ggtext)
library(waffle)
Cargamos la base de datos original y además, de ser necesario, cargamos los agrupamientos por separado, considerando que cada agrupamiento son los distintos tratamientos que se aplican.
La base original consiste en tener la variable independiente, que corresponde a los tratamientos, mientras que la variable dependiente son los resultados obtenidos por cada longitud de onda o lambda.
Se propone realizar un análisis de varianzas (ANOVA) para revisar si se encuentran diferencias estadÃsticamente significativas dentro y entre los grupos.
Para llevar a cabo un ANOVA, se debe realizar una prueba de normalidad en los residuos y de homogeneidad de varianzas. Esto debe ser llevado a cabo por grupos debido a que puede que cada grupo tenga supuestos diferentes.
A continuación, se adjunta algoritmo para realizar las pruebas correspondientes para verificar supuestos. Este análisis, lo ideal serÃa llevarlo a cabo de manera general y además por agrupamiento, para esto es solo cambiar de base de datos, funcionan de la misma manera.
div21$Grupo = as.factor(div21$Grupo)
var_dependientes = setdiff(names(div21[,2:101]), c("Y.Values", "Grupo", "Method"))
# Función para evaluar normalidad y homogeneidad de varianzas
ajustar_modelo = function(var) {
formula = as.formula(paste(var, "~ Y.Values"))
modelo = lm(formula, data = div21)
residuos = residuals(modelo)
# Prueba de normalidad (Shapiro-Wilk)
p_normalidad = shapiro.test(residuos)$p.value
# Prueba de homogeneidad de varianza (Breusch-Pagan)
p_homogeneidad = bptest(modelo)$p.value
return(data.frame(variable = var, p_normalidad, p_homogeneidad))
}
En caso de cumplirse todos los supuestos, se lleva a cabo un ANOVA.
En caso de no cumplirse el supuesto de homogeneidad de varianzas, se utiliza un ANOVA Welch.
En caso de no cumplirse el supuesto de normalidad, se utiliza Kruskal-Wallis
# Función para realizar ANOVA o Kruskal-Wallis según normalidad
ajustar_anova_o_kruskal = function(var, p_normalidad, p_homogeneidad) {
formula = as.formula(paste(var, "~ SampleID"))
if (p_normalidad >= 0.05) { # Si los datos son normales
if (p_homogeneidad >= 0.05) {
resultado_anova = aov(formula, data = div21) # ANOVA clásico
resumen_anova = summary(resultado_anova)[[1]]
return(data.frame(
Variable = var,
Test_Usado = "ANOVA",
F_value = resumen_anova["F value"][1,1],
gl_Entre = resumen_anova["Df"][1,1],
gl_Dentro = resumen_anova["Df"][2,1],
p_value = resumen_anova["Pr(>F)"][1,1]
))
} else {
formula = as.formula(paste(var, "~ SampleID"))
resultado_welch = oneway.test(formula, data = div21, var.equal = FALSE) # ANOVA de Welch
return(data.frame(
Variable = var,
Test_Usado = "Welch-ANOVA",
F_value = resultado_welch$statistic,
gl_Entre = resultado_welch$parameter[1],
gl_Dentro = resultado_welch$parameter[2],
p_value = resultado_welch$p.value
))
}
} else { # Si los datos NO son normales → Kruskal-Wallis
formula = formula = as.formula(paste(var, "~ SampleID"))
resultado_kruskal = kruskal.test(formula, data = div21)
return(data.frame(
Variable = var,
Test_Usado = "Kruskal-Wallis",
F_value = NA, # Kruskal no tiene F-value
gl_Entre = resultado_kruskal$parameter,
gl_Dentro = NA,
p_value = resultado_kruskal$p.value
))
}
}
Al tener submuestras, se debe verificar que no hay variabilidad significativa dentro de estas, por lo que para lleva a cabo este análisis se realiza un análisis dentro de cada grupo, tomando en consideración los supuestos asociados dentro de estos. Considerando: Normalidad y homogeneidad de varianzas, se realiza anova. Normalidad y heterocedasticidad de varianzas, se realiza anova de Welch. No normalidad, se realiza Kruskal-Wallis.
# Función para realizar ANOVA o Kruskal-Wallis dentro de cada grupo
ajustar_anova_dentro_grupo = function(var) {
resultados = list()
for (grupo in unique(div21$SampleID)) {
subdata = subset(div21, SampleID == grupo)
# Verificar que haya más de un nivel en Y.Values dentro del grupo
if (length(unique(subdata$Grupo)) > 1) {
# Prueba de normalidad en cada grupo
p_normalidad = shapiro.test(subdata[[var]])$p.value
p_homogeneidad = bptest(lm(as.formula(paste(var, "~ Grupo")), data = subdata))$p.value
if (p_normalidad >= 0.05) {
if (p_homogeneidad >= 0.05) {
resultado_anova = aov(as.formula(paste(var, "~ Grupo")), data = subdata)
resumen_anova = summary(resultado_anova)[[1]]
resultados[[grupo]] = data.frame(
Variable = var,
Grupo = grupo,
Test_Usado = "ANOVA",
F_value = resumen_anova["F value"][1,1],
gl_Entre = resumen_anova["Df"][1,1],
gl_Dentro = resumen_anova["Df"][2,1],
p_value = resumen_anova["Pr(>F)"][1,1]
)
} else {
resultado_welch = oneway.test(as.formula(paste(var, "~ Grupo")), data = subdata, var.equal = FALSE)
resultados[[grupo]] = data.frame(
Variable = var,
Grupo = grupo,
Test_Usado = "Welch-ANOVA",
F_value = resultado_welch$statistic,
gl_Entre = resultado_welch$parameter[1],
gl_Dentro = resultado_welch$parameter[2],
p_value = resultado_welch$p.value
)
}
} else {
resultado_kruskal = kruskal.test(as.formula(paste(var, "~ Grupo")), data = subdata)
resultados[[grupo]] = data.frame(
Variable = var,
Grupo = grupo,
Test_Usado = "Kruskal-Wallis",
F_value = NA, # Kruskal no tiene F-value
gl_Entre = resultado_kruskal$parameter,
gl_Dentro = NA,
p_value = resultado_kruskal$p.value
)
}
} else {
resultados[[grupo]] = data.frame(
Variable = var,
Grupo = grupo,
Test_Usado = NA,
F_value = NA,
gl_Entre = NA,
gl_Dentro = NA,
p_value = NA
)
}
}
return(do.call(rbind, resultados))
}
# Evaluar normalidad y homogeneidad de varianza
resultados_modelo = do.call(rbind, lapply(var_dependientes, ajustar_modelo))
# Aplicar ANOVA o Kruskal-Wallis según normalidad
resultados_finales = list()
resultados_dentro_grupo = list()
for (var in var_dependientes) {
p_normalidad = resultados_modelo[resultados_modelo$variable == var, "p_normalidad"]
p_homogeneidad = resultados_modelo[resultados_modelo$variable == var, "p_homogeneidad"]
resultados_finales[[var]] = ajustar_anova_o_kruskal(var, p_normalidad, p_homogeneidad)
resultados_dentro_grupo[[var]] = ajustar_anova_dentro_grupo(var)
}
Para visualizar lo recolectado anteriormente, se rescata en objetos dichos resultados y se guardan en un archivo con formato .xlsx
# Combinar resultados y exportar a Excel
resultados_modelo = do.call(rbind, lapply(var_dependientes, ajustar_modelo))
resultados_df = do.call(rbind, resultados_finales)
resultados_dentro_df = do.call(rbind, resultados_dentro_grupo)
# write.xlsx(list(Entre_Grupos = resultados_df, Dentro_Grupos = resultados_dentro_df, Supuestos= resultados_modelo), "Resultados_Analisis_23.xlsx", rowNames = FALSE)
# Mostrar resultados en consola
# print(resultados_df)
# print(resultados_dentro_df)
Para analizar los anterior de manera más simple se realizan dos gráficos de waffle, que muestran si es que el análisis de varianzas es significativo o no lo es.
graf_waffle + graf_waffle1 + graf_waffle2 + graf_waffle3 +graf_waffle4+graf_waffle5+
graf_waffle6 + graf_waffle7 + graf_waffle8 + graf_waffle9 + graf_waffle10 +
graf_waffle11 + graf_waffle12 + graf_waffle13 + graf_waffle14 + graf_waffle15 +
graf_waffle16 + graf_waffle17 + graf_waffle18 + graf_waffle19 + graf_waffle20 +
plot_layout(nrow = 21)
graf_waffleD + graf_waffleD1 + graf_waffleD2 + graf_waffleD3 +graf_waffleD4 + graf_waffleD5+
graf_waffleD6 + graf_waffleD7 + graf_waffleD8 + graf_waffleD9 + graf_waffleD10 +
graf_waffleD11 +
plot_layout(nrow = 12)
graf_waffleD12 + graf_waffleD13 + graf_waffleD14 + graf_waffleD15 +
graf_waffleD16 + graf_waffleD17 + graf_waffleD18 + graf_waffleD19 + graf_waffleD20 +
plot_layout(nrow = 9)