1. Introducción

1.1. ¿De qué trata este trabajo?

Cuando uno compra café casi siempre lo hace por costumbre, por el precio o simplemente porque es la marca que ya conoce, y pocas veces se pregunta si de verdad hay diferencias entre una marca y otra. Con esa duda decidimos comparar tres marcas de café: Café Córdoba, Sello Rojo y Medalla de Oro.

Para hacerlo realizamos una pequeña evaluación sensorial, que consiste en usar los sentidos (en este caso el gusto y el olfato) para medir las características de un alimento. Seis personas probaron cada café y le pusieron una nota de 1 a 10 al sabor y otra al aroma, donde 1 es muy malo y 10 es excelente. Lo que queremos saber es algo sencillo: ¿las personas notan diferencias reales entre las marcas, o las diferencias que aparecen en las notas son solo cuestión de azar?

1.2. ¿Cómo se organizaron los datos?

Tenemos tres marcas y cada una fue calificada por seis personas, así que en total son 18 notas de sabor y 18 notas de aroma. Como todas las marcas tienen la misma cantidad de notas, decimos que el diseño está balanceado, lo que facilita los análisis.

En todas las pruebas trabajamos con un nivel de significancia de 0.05 (5 %). Esto funciona como una regla de decisión: cada prueba nos entrega un número llamado valor p, y si ese valor es menor que 0.05 concluimos que sí hay una diferencia real; si es mayor, concluimos que no tenemos suficiente evidencia para afirmar que las marcas sean distintas.

1.3. ¿Cómo lo analizamos?

Dividimos el trabajo en dos partes que se complementan:

  1. En la Parte I usamos el ANOVA, que es la prueba más conocida para comparar promedios de varios grupos, y después comparamos las marcas de dos en dos.
  2. En la Parte II usamos la prueba de Kruskal-Wallis, que en lugar de trabajar con las notas directamente trabaja con el orden de las notas (de menor a mayor).

Decidimos usar las dos porque unas notas de 1 a 10, dadas por pocas personas, no siempre se comportan como el ANOVA necesita. Por eso primero revisamos si el ANOVA se puede usar y luego comparamos su resultado con el de la otra prueba, para ver si las dos nos cuentan lo mismo.

1.4. ¿Qué significa cada prueba?

Estadística descriptiva y diagrama de cajas

Antes de hacer cualquier prueba conviene mirar los datos. Para cada marca calculamos el promedio (la media), la mediana (la nota que queda en el medio cuando se ordenan todas) y la desviación estándar (que nos dice qué tanto se alejan las notas del promedio; si es pequeña, las personas estuvieron de acuerdo). El diagrama de cajas muestra todo eso en un dibujo: la línea del centro es la mediana, la caja agrupa la mitad central de las notas y los bigotes muestran hasta dónde llegan las demás.

ANOVA

El ANOVA compara los promedios de todas las marcas al mismo tiempo. Parte de una idea de base, llamada hipótesis nula, que dice que las marcas son iguales; la contraria, la hipótesis alternativa, dice que al menos una marca es distinta. Para decidir, compara cuánto varían las notas entre marcas con cuánto varían dentro de una misma marca (es decir, el desacuerdo normal entre las personas que probaron el mismo café). Si las diferencias entre marcas son mucho más grandes que ese desacuerdo normal, hay evidencia de que las marcas no son iguales. Eso sí, el ANOVA solo nos dice si hay diferencia, pero no entre cuáles marcas.

Supuestos del ANOVA

Para que los resultados del ANOVA sean confiables, los datos deben cumplir algunas condiciones. Las revisamos sobre los residuos, que son lo que le sobra o le falta a cada nota respecto al promedio de su marca:

  • Que los residuos sean normales, es decir, que se repartan de forma parecida a la famosa campana. Lo revisamos con las pruebas de Shapiro-Wilk y Kolmogorov-Smirnov y con el gráfico Q-Q, donde los puntos deben quedar cerca de la línea recta. Aquí nos conviene que el valor p sea mayor que 0.05, porque significa que no hay motivos para dudar de la normalidad.
  • Que la variación sea parecida en todas las marcas. Lo revisamos con las pruebas de Bartlett y Levene y con el gráfico de residuos contra valores ajustados, donde los puntos no deben formar ningún patrón. Igual que antes, queremos un valor p mayor que 0.05.
  • Que las notas sean independientes, o sea, que cada persona calificó por su cuenta sin influir en las demás. Esto se asume por la forma en que se recogieron los datos.

LSD y Tukey

Si el ANOVA nos dice que hay diferencias, la siguiente pregunta es entre cuáles marcas. Para eso sirven estas dos pruebas, que comparan las marcas de dos en dos. La LSD es más permisiva y detecta diferencias con más facilidad, pero se equivoca más cuando se hacen varias comparaciones. La de Tukey es más exigente y cuida mejor ese error, por eso suele considerarse más confiable.

Rangos

Consiste en ordenar todas las notas de menor a mayor y cambiar cada nota por su posición en la fila (su rango). Cuando dos o más notas son iguales, se reparte el promedio de las posiciones que ocupan. Después se suman y se promedian los rangos de cada marca: una marca con un rango promedio alto es una marca que tiende a recibir las notas más altas.

Prueba de Kruskal-Wallis

Es como el ANOVA, pero en lugar de usar las notas usa los rangos, y por eso no exige que los datos sean normales. Su idea de base (hipótesis nula) es que las tres marcas reciben notas parecidas, y la alternativa es que al menos una marca se diferencia de las otras. Si el valor p sale menor que 0.05, concluimos que sí hay diferencias.

Tamaño del efecto (épsilon cuadrado)

El valor p solo responde si hay diferencia, pero no qué tan grande es. Para eso calculamos el épsilon cuadrado, un número entre 0 y 1: mientras más cerca de 1, más pesa la marca en las notas que se dieron. Como guía, menos de 0.01 es insignificante, de 0.01 a 0.08 es pequeño, de 0.08 a 0.26 es moderado y desde 0.26 es grande.

Prueba de Dunn

Se usa solo cuando Kruskal-Wallis detecta diferencias, y su trabajo es decirnos entre cuáles marcas están, comparándolas de dos en dos. Como se hacen varias comparaciones, los valores p se corrigen con el método de Holm, para no equivocarnos por el hecho de comparar muchas veces.

2. Parte I: ANOVA y comparaciones múltiples

En esta parte se analiza cada atributo con el modelo lineal y el ANOVA, se hacen las comparaciones múltiples y se verifican los supuestos. Los bloques de código se presentan en el orden en que se ejecutan.

2.1. Preparación del entorno

Se cargan los paquetes necesarios: agricolae (para la prueba LSD) y car (para la prueba de Levene). También se definen los nombres de los tratamientos (las tres marcas), el número de tratamientos (t), el número de réplicas (r) y el total de observaciones (N).

library(agricolae)
## Warning: package 'agricolae' was built under R version 4.6.2
library(car)
## Cargando paquete requerido: carData
Tratamientos <- c("Café Córdoba", "Sello Rojo", "Medalla de Oro")
t <- 3  # número de tratamientos
r <- 6  # número de réplicas (datos balanceados)
N <- t * r

2.2. Análisis del sabor

2.2.1. Ingreso de los datos

Se guarda en el vector sabor las 18 calificaciones, ordenadas por marca (las primeras 6 son de Café Córdoba, las siguientes 6 de Sello Rojo y las últimas 6 de Medalla de Oro).

# 1. SABOR


# Datos de sabor (6 por café, en el orden de Tratamientos)
sabor <- c(8, 9, 7, 8, 9, 8,     # Café Córdoba
           9, 8, 9, 10, 8, 9,    # Sello Rojo
           7, 8, 6, 8, 7, 9)     # Medalla de Oro

2.2.2. Tabla de réplicas, totales y medias

Se organizan los datos en réplicas (r1 a r6): cada vector contiene una calificación por marca. Con cbind se forma una matriz donde cada fila es una marca y cada columna una réplica. Después se calcula el total de cada marca (suma de sus 6 calificaciones) y su media (total dividido entre r), y todo se presenta en una tabla resumen.

# Réplicas por tratamiento (fila = café, columna = réplica)
r1 <- c(8, 9, 7)
r2 <- c(9, 8, 8)
r3 <- c(7, 9, 6)
r4 <- c(8, 10, 8)
r5 <- c(9, 8, 7)
r6 <- c(8, 9, 9)
replicas <- cbind(r1, r2, r3, r4, r5, r6)
Total <- rowSums(replicas)
medias <- Total / r

Tabla_sabor <- data.frame(Tratamientos, r1, r2, r3, r4, r5, r6, Total, medias)
Tabla_sabor

2.2.3. Organización de los datos para el análisis

Se crea el factor cafes, que repite el nombre de cada marca 6 veces, y se une con las calificaciones en el data frame DATA_sabor. Este formato “largo” (una columna para el grupo y otra para la respuesta) es el que necesitan las funciones de modelo y de gráficos.

cafes <- factor(rep(Tratamientos, each = r), levels = Tratamientos)
DATA_sabor <- data.frame(cafes, y = sabor)
DATA_sabor

2.2.4. Diagrama de cajas

Muestra la distribución de las calificaciones de sabor de cada marca: la línea central es la mediana, la caja contiene el 50 % central de los datos y los bigotes muestran el rango.

boxplot(y ~ cafes, data = DATA_sabor, col = "green",
        main = "Calificación del sabor",
        xlab = "Tipo de café", ylab = "Sabor (escala de 1 a 10)")

¿Qué nos dice este gráfico? Cada caja representa una marca y muestra cómo se repartieron sus seis notas de sabor. En Café Córdoba las notas van de 7 a 9, la mitad central está entre 8 y 9 y la mediana es 8 (la raya queda en la parte baja de la caja). En Sello Rojo las notas van de 8 a 10, la caja va de 8 a 9 y la mediana es 9, así que es la marca que se ubica más arriba. En Medalla de Oro las notas van de 6 a 9, la caja va de 7 a 8 y la mediana es 7.5; es la que queda más abajo y la que tiene los bigotes más largos, es decir, la de notas más dispersas. No aparecen puntos sueltos fuera de las cajas, lo que significa que no hubo calificaciones atípicas.

En resumen, se ve una tendencia a favor de Sello Rojo y en contra de Medalla de Oro, pero las cajas se tocan o se solapan en varios valores (por ejemplo, todas incluyen notas de 8). Por eso la diferencia no es tan clara a simple vista y necesitamos las pruebas estadísticas para decidir si es real.

2.2.5. Modelo lineal y ANOVA

Con lm se ajusta el modelo en el que el sabor depende de la marca, y con aov se obtiene la tabla de análisis de varianza. summary(modelo_s) muestra los coeficientes (la media de la marca de referencia y la diferencia de las otras dos respecto a ella) y el coeficiente de determinación \(R^2\); summary(anova_s) muestra la tabla ANOVA con el estadístico \(F\) y su valor p.

modelo_s <- lm(y ~ cafes, data = DATA_sabor)
anova_s <- aov(modelo_s)
summary(modelo_s)
## 
## Call:
## lm(formula = y ~ cafes, data = DATA_sabor)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -1.5   -0.5    0.0    0.5    1.5 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           8.1667     0.3522  23.185 3.67e-13 ***
## cafesSello Rojo       0.6667     0.4981   1.338    0.201    
## cafesMedalla de Oro  -0.6667     0.4981  -1.338    0.201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8628 on 15 degrees of freedom
## Multiple R-squared:  0.3232, Adjusted R-squared:  0.233 
## F-statistic: 3.582 on 2 and 15 DF,  p-value: 0.05349
summary(anova_s)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## cafes        2  5.333  2.6667   3.582 0.0535 .
## Residuals   15 11.167  0.7444                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.2.6. Comparaciones múltiples (LSD y Tukey)

LSD.test agrupa las marcas con letras: las que comparten letra no difieren significativamente. TukeyHSD entrega, para cada par de marcas, la diferencia de medias, su intervalo de confianza del 95 % y el valor p ajustado; el gráfico plot(tukey_s) muestra esos intervalos, y los que no cruzan el cero indican diferencia significativa.

# Comparaciones múltiples
LSD.test(modelo_s, "cafes", group = TRUE, console = TRUE,
         main = "Sabor de los cafés")
## 
## Study: Sabor de los cafés
## 
## LSD t Test for y 
## 
## Mean Square Error:  0.7444444 
## 
## cafes,  means and individual ( 95 %) CI
## 
##                       y       std r        se      LCL      UCL Min Max  Q25
## Café Córdoba   8.166667 0.7527727 6 0.3522415 7.415882 8.917452   7   9 8.00
## Medalla de Oro 7.500000 1.0488088 6 0.3522415 6.749215 8.250785   6   9 7.00
## Sello Rojo     8.833333 0.7527727 6 0.3522415 8.082548 9.584118   8  10 8.25
##                Q50  Q75
## Café Córdoba   8.0 8.75
## Medalla de Oro 7.5 8.00
## Sello Rojo     9.0 9.00
## 
## Alpha: 0.05 ; DF Error: 15
## Critical Value of t: 2.13145 
## 
## least Significant Difference: 1.06177 
## 
## Treatments with the same letter are not significantly different.
## 
##                       y groups
## Sello Rojo     8.833333      a
## Café Córdoba   8.166667     ab
## Medalla de Oro 7.500000      b
tukey_s <- TukeyHSD(anova_s, "cafes", conf.level = 0.95)
plot(tukey_s)

tukey_s
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = modelo_s)
## 
## $cafes
##                                   diff        lwr         upr     p adj
## Sello Rojo-Café Córdoba      0.6666667 -0.6272494  1.96058271 0.3967618
## Medalla de Oro-Café Córdoba -0.6666667 -1.9605827  0.62724938 0.3967618
## Medalla de Oro-Sello Rojo   -1.3333333 -2.6272494 -0.03941729 0.0430246

¿Qué nos dice este gráfico? Cada línea horizontal es el intervalo de confianza del 95 % para la diferencia de promedios entre dos marcas (su centro es la diferencia que observamos), y la línea punteada vertical marca el 0, que significa “no hay diferencia”. La regla es sencilla: si el intervalo cruza el 0, esa diferencia no es significativa; si queda completamente a un lado, sí lo es.

Aquí las comparaciones de Sello Rojo con Café Córdoba (diferencia de +0.67 puntos, intervalo de -0.63 a 1.96) y de Medalla de Oro con Café Córdoba (-0.67, intervalo de -1.96 a 0.63) cruzan el 0, así que no hay diferencia significativa. La comparación de Medalla de Oro con Sello Rojo (-1.33, intervalo de -2.63 a -0.04) es la única que queda del lado negativo, pero su extremo derecho está casi pegado al 0. Es decir, es una diferencia que apenas alcanza a ser significativa, y por eso la tomamos con cautela.

2.2.7. Residuos y verificación de normalidad

Se extraen los residuos del modelo y se grafican en un gráfico Q-Q para verificar visualmente si siguen una distribución normal (los puntos deben seguir la línea).

# Verificación de supuestos
resid_s <- modelo_s$residuals
resid_s
##          1          2          3          4          5          6          7 
## -0.1666667  0.8333333 -1.1666667 -0.1666667  0.8333333 -0.1666667  0.1666667 
##          8          9         10         11         12         13         14 
## -0.8333333  0.1666667  1.1666667 -0.8333333  0.1666667 -0.5000000  0.5000000 
##         15         16         17         18 
## -1.5000000  0.5000000 -0.5000000  1.5000000
# Normalidad
qqnorm(resid_s)
qqline(resid_s)

¿Qué nos dice este gráfico? Este gráfico compara los residuos reales (eje vertical) con los valores que tendríamos si los residuos fueran perfectamente normales (eje horizontal). Si los puntos se pegan a la línea recta, los residuos se comportan como una normal y se cumple el supuesto. En nuestro caso los puntos siguen de cerca la línea, sin curvas marcadas ni puntos muy alejados en los extremos. Se ven pequeños “escalones” porque las notas son números enteros y varios residuos se repiten, y eso es normal. Esto concuerda con la prueba de Shapiro-Wilk del sabor (p = 0.986), que no encuentra motivos para dudar de la normalidad.

2.2.8. Homogeneidad de varianzas y gráficos de diagnóstico

El gráfico de residuos contra valores ajustados permite ver si la dispersión es parecida en todos los grupos. Luego par(mfrow = c(2, 2)) organiza en una sola ventana los cuatro gráficos de diagnóstico que genera plot(modelo_s), y par(mfrow = c(1, 1)) devuelve la ventana a su formato normal.

# Homogeneidad
plot(modelo_s$fitted.values, resid_s, main = "Residuals vs. Fitted (sabor)", pch = 20)

# Mezcla de gráficos
par(mfrow = c(2, 2))
plot(modelo_s)

par(mfrow = c(1, 1))

¿Qué nos dice este gráfico?

Gráfico de residuos contra valores ajustados (el primero). Como el modelo solo usa la marca, el valor ajustado de cada nota es el promedio de su marca, por eso los puntos aparecen en tres columnas verticales (en 7.5 para Medalla de Oro, 8.17 para Café Córdoba y 8.83 para Sello Rojo). Lo que miramos es que las tres columnas tengan una altura parecida y estén centradas en 0. La de Medalla de Oro es un poco más alta (de -1.5 a 1.5) que las otras dos (cerca de -1.2 a 1.2), pero no se forma un embudo ni ningún patrón, y las pruebas de Bartlett (p = 0.700) y Levene (p = 0.481) confirman que las varianzas pueden considerarse iguales.

Los cuatro gráficos de diagnóstico. Arriba a la izquierda está otra vez el de residuos contra ajustados, y la línea roja que lo acompaña debe ser casi plana y cercana a 0, como aquí. Arriba a la derecha está el gráfico Q-Q (con residuos estandarizados), que se lee igual que el anterior. Abajo a la izquierda está el Scale-Location, que sirve para revisar la homogeneidad de varianzas: se espera que la altura de las tres columnas sea parecida. Abajo a la derecha aparece el gráfico de residuos contra los niveles del factor; como todas las marcas tienen seis notas, todas las observaciones pesan lo mismo, y ninguna tiene una influencia especial sobre el modelo. Los números que acompañan a algunos puntos son simplemente el número de fila de las observaciones con mayor residuo; no son datos atípicos, solo los más alejados del promedio de su marca.

2.2.9. Pruebas formales de normalidad y homogeneidad

Se aplican Shapiro-Wilk y Kolmogorov-Smirnov (normalidad de los residuos) y Levene (igualdad de varianzas). Como las calificaciones son números enteros, hay valores repetidos, y R puede mostrar una advertencia en la prueba de Kolmogorov-Smirnov porque no admite empates de forma exacta; el resultado sigue siendo orientativo.

# Pruebas de normalidad y homogeneidad
shapiro.test(resid_s)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_s
## W = 0.98492, p-value = 0.9864
ks.test(resid_s, "pnorm", mean = mean(resid_s), sd = sd(resid_s))
## Warning in ks.test.default(resid_s, "pnorm", mean = mean(resid_s), sd =
## sd(resid_s)): ties should not be present for the one-sample Kolmogorov-Smirnov
## test
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resid_s
## D = 0.085202, p-value = 0.9994
## alternative hypothesis: two-sided
leveneTest(modelo_s$residuals ~ cafes, center = median, data = DATA_sabor)

2.3. Análisis del aroma

Se repite exactamente el mismo procedimiento de la sección 2.2, ahora con las calificaciones de aroma.

2.3.1. Ingreso de los datos

# 2. AROMA


aroma <- c(9, 8, 7, 8, 9, 8,     # Café Córdoba
           7, 8, 9, 6, 8, 7,     # Sello Rojo
           8, 7, 9, 8, 9, 6)     # Medalla de Oro

2.3.2. Tabla de réplicas, totales y medias

Se vuelven a definir r1 a r6 con los datos de aroma (una calificación por marca en cada réplica) y se calculan los totales y las medias de cada marca.

# Réplicas por tratamiento (fila = café, columna = réplica)
r1 <- c(9, 7, 8)
r2 <- c(8, 8, 7)
r3 <- c(7, 9, 9)
r4 <- c(8, 6, 8)
r5 <- c(9, 8, 9)
r6 <- c(8, 7, 6)
replicas <- cbind(r1, r2, r3, r4, r5, r6)
Total <- rowSums(replicas)
medias <- Total / r

Tabla_aroma <- data.frame(Tratamientos, r1, r2, r3, r4, r5, r6, Total, medias)
Tabla_aroma

2.3.3. Organización de los datos para el análisis

Se reutiliza el factor cafes creado antes y se une con las calificaciones de aroma en DATA_aroma.

DATA_aroma <- data.frame(cafes, y = aroma)
DATA_aroma

2.3.4. Diagrama de cajas

boxplot(y ~ cafes, data = DATA_aroma, col = "green",
        main = "Calificación del aroma",
        xlab = "Tipo de café", ylab = "Aroma (escala de 1 a 10)")

¿Qué nos dice este gráfico? En Café Córdoba las notas van de 7 a 9, la caja va de 8 a 9 y la mediana es 8. En Sello Rojo van de 6 a 9, la caja va de 7 a 8 y la mediana es 7.5. En Medalla de Oro también van de 6 a 9, pero la caja es más ancha (de 7 a 9), lo que indica que las personas se pusieron menos de acuerdo en el aroma de esta marca; su mediana es 8. No hay puntos atípicos.

Las tres cajas están muy encimadas y las medianas (8, 7.5 y 8) son casi iguales, así que a simple vista no se ve una marca claramente mejor que otra en aroma. Esto anticipa lo que confirmarán las pruebas.

2.3.5. Modelo lineal y ANOVA

modelo_a <- lm(y ~ cafes, data = DATA_aroma)
anova_a <- aov(modelo_a)
summary(modelo_a)
## 
## Call:
## lm(formula = y ~ cafes, data = DATA_aroma)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.833 -0.500  0.000  0.750  1.500 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           8.1667     0.4105  19.894 3.41e-12 ***
## cafesSello Rojo      -0.6667     0.5805  -1.148    0.269    
## cafesMedalla de Oro  -0.3333     0.5805  -0.574    0.574    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.006 on 15 degrees of freedom
## Multiple R-squared:  0.08081,    Adjusted R-squared:  -0.04175 
## F-statistic: 0.6593 on 2 and 15 DF,  p-value: 0.5316
summary(anova_a)
##             Df Sum Sq Mean Sq F value Pr(>F)
## cafes        2  1.333  0.6667   0.659  0.532
## Residuals   15 15.167  1.0111

2.3.6. Comparaciones múltiples (LSD y Tukey)

# Comparaciones múltiples
LSD.test(modelo_a, "cafes", group = TRUE, console = TRUE,
         main = "Aroma de los cafés")
## 
## Study: Aroma de los cafés
## 
## LSD t Test for y 
## 
## Mean Square Error:  1.011111 
## 
## cafes,  means and individual ( 95 %) CI
## 
##                       y       std r        se      LCL      UCL Min Max  Q25
## Café Córdoba   8.166667 0.7527727 6 0.4105101 7.291685 9.041648   7   9 8.00
## Medalla de Oro 7.833333 1.1690452 6 0.4105101 6.958352 8.708315   6   9 7.25
## Sello Rojo     7.500000 1.0488088 6 0.4105101 6.625018 8.374982   6   9 7.00
##                Q50  Q75
## Café Córdoba   8.0 8.75
## Medalla de Oro 8.0 8.75
## Sello Rojo     7.5 8.00
## 
## Alpha: 0.05 ; DF Error: 15
## Critical Value of t: 2.13145 
## 
## least Significant Difference: 1.237411 
## 
## Treatments with the same letter are not significantly different.
## 
##                       y groups
## Café Córdoba   8.166667      a
## Medalla de Oro 7.833333      a
## Sello Rojo     7.500000      a
tukey_a <- TukeyHSD(anova_a, "cafes", conf.level = 0.95)
plot(tukey_a)

tukey_a
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = modelo_a)
## 
## $cafes
##                                   diff       lwr       upr     p adj
## Sello Rojo-Café Córdoba     -0.6666667 -2.174625 0.8412918 0.5005539
## Medalla de Oro-Café Córdoba -0.3333333 -1.841292 1.1746252 0.8357270
## Medalla de Oro-Sello Rojo    0.3333333 -1.174625 1.8412918 0.8357270

¿Qué nos dice este gráfico? Igual que en el sabor, cada línea es el intervalo de confianza del 95 % para la diferencia entre dos marcas, y si cruza el 0 (línea punteada) la diferencia no es significativa. En el aroma las tres líneas cruzan el 0: Sello Rojo frente a Café Córdoba (-0.67, intervalo de -2.18 a 0.84), Medalla de Oro frente a Café Córdoba (-0.33, de -1.84 a 1.18) y Medalla de Oro frente a Sello Rojo (+0.33, de -1.18 a 1.84). Además, las líneas son largas y están centradas cerca del 0, lo que indica que las diferencias observadas son pequeñas frente a la variación normal entre catadores. Ninguna pareja de marcas se distingue en aroma.

2.3.7. Residuos y verificación de normalidad

# Verificación de supuestos
resid_a <- modelo_a$residuals
resid_a
##          1          2          3          4          5          6          7 
##  0.8333333 -0.1666667 -1.1666667 -0.1666667  0.8333333 -0.1666667 -0.5000000 
##          8          9         10         11         12         13         14 
##  0.5000000  1.5000000 -1.5000000  0.5000000 -0.5000000  0.1666667 -0.8333333 
##         15         16         17         18 
##  1.1666667  0.1666667  1.1666667 -1.8333333
# Normalidad
qqnorm(resid_a)
qqline(resid_a)

¿Qué nos dice este gráfico? Se interpreta igual que el del sabor: mientras más cerca de la línea estén los puntos, más normales son los residuos. Aquí los puntos siguen la línea de manera razonable, incluso en los extremos (el residuo más bajo, de -1.83, queda cerca de lo que se esperaría), y los pequeños escalones se deben a que las notas son enteros. Concuerda con la prueba de Shapiro-Wilk del aroma (p = 0.824): el supuesto de normalidad se cumple.

2.3.8. Homogeneidad de varianzas y gráficos de diagnóstico

# Homogeneidad
plot(modelo_a$fitted.values, resid_a, main = "Residuals vs. Fitted (aroma)", pch = 20)

# Mezcla de gráficos
par(mfrow = c(2, 2))
plot(modelo_a)

par(mfrow = c(1, 1))

¿Qué nos dice este gráfico?

Gráfico de residuos contra valores ajustados (el primero). Otra vez los puntos forman tres columnas verticales, una por marca (en 7.5 para Sello Rojo, 7.83 para Medalla de Oro y 8.17 para Café Córdoba). Las columnas de Sello Rojo (de -1.5 a 1.5) y de Medalla de Oro (de -1.83 a 1.17) son un poco más largas que la de Café Córdoba (de -1.17 a 0.83), es decir, hubo algo más de desacuerdo entre las personas en esas dos marcas. Aun así, no se forma un embudo ni ningún patrón, todo está centrado en 0 y las pruebas de Bartlett (p = 0.641) y Levene (p = 0.568) confirman que las varianzas pueden considerarse iguales.

Los cuatro gráficos de diagnóstico. Se leen igual que los del sabor: arriba a la izquierda, residuos contra ajustados con una línea roja casi plana; arriba a la derecha, el gráfico Q-Q de residuos estandarizados; abajo a la izquierda, el Scale-Location para revisar la homogeneidad de varianzas; y abajo a la derecha, los residuos contra los niveles del factor, donde ninguna observación tiene una influencia especial porque el diseño es balanceado. Ninguno de los cuatro muestra problemas serios, así que el modelo es adecuado para el aroma.

2.3.9. Pruebas formales de normalidad y homogeneidad

# Pruebas de normalidad y homogeneidad
shapiro.test(resid_a)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_a
## W = 0.97141, p-value = 0.8237
ks.test(resid_a, "pnorm", mean = mean(resid_a), sd = sd(resid_a))
## Warning in ks.test.default(resid_a, "pnorm", mean = mean(resid_a), sd =
## sd(resid_a)): ties should not be present for the one-sample Kolmogorov-Smirnov
## test
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resid_a
## D = 0.096636, p-value = 0.996
## alternative hypothesis: two-sided
leveneTest(modelo_a$residuals ~ cafes, center = median, data = DATA_aroma)

3. Parte II: Análisis no paramétrico del sabor

En esta parte se repite el análisis con un enfoque más completo, que incluye estadística descriptiva por grupo, la prueba de Bartlett, los rangos, la prueba de Kruskal-Wallis, el tamaño del efecto y la prueba de Dunn.

3.1. Ingreso de los datos

Se crea el vector sabor con las 18 calificaciones y el factor cafe con las tres marcas (cada una repetida 6 veces, con los niveles en un orden fijo). Con ambos se construye el data frame datos_sabor, que es la base de todo el análisis de esta parte.

sabor <- c(8, 9, 7, 8, 9, 8,
           9, 8, 9, 10, 8, 9,
           7, 8, 6, 8, 7, 9)

cafe <- factor(rep(c("Café Córdoba", "Sello Rojo",
                     "Medalla de Oro"), each = 6),
               levels = c("Café Córdoba", "Sello Rojo",
                          "Medalla de Oro"))

datos_sabor <- data.frame(cafe, sabor)

datos_sabor

3.2. Estadística descriptiva

La función aggregate calcula, para cada marca, el número de datos, la media, la mediana y la desviación estándar del sabor. Es la primera aproximación para ver qué marca tiene calificaciones más altas y cuál es más homogénea.

aggregate(sabor ~ cafe, data = datos_sabor,
          FUN = function(x) c(
            n = length(x),
            media = mean(x),
            mediana = median(x),
            desviacion = sd(x)
          ))

3.3. Gráfico de cajas

Representa la distribución de las calificaciones de cada marca, con un color distinto para cada una. Permite comparar medianas, dispersión y posibles valores atípicos.

boxplot(sabor ~ cafe, data = datos_sabor,
        col = c("#F4B183", "#9DC3E6", "#A9D18E"),
        ylab = "Sabor (escala de 1 a 10)",
        xlab = "Tipo de café",
        main = "Calificación del sabor")

¿Qué nos dice este gráfico? Es el mismo diagrama de la Parte I, ahora con un color para cada marca (naranja para Café Córdoba, azul para Sello Rojo y verde para Medalla de Oro), así que se lee igual. Sello Rojo es la marca con las notas más altas (mediana 9), Café Córdoba queda en el medio (mediana 8) y Medalla de Oro es la más baja (mediana 7.5) y la más dispersa. No hay datos atípicos. Hay una tendencia clara en el orden de las marcas, pero las cajas se solapan, por lo que falta ver si la diferencia es significativa.

3.4. ANOVA

Se ajusta el ANOVA de un factor con aov y se muestra la tabla con summary. Su resultado se contrasta más adelante con el de Kruskal-Wallis.

modelo_sabor <- aov(sabor ~ cafe, data = datos_sabor)

summary(modelo_sabor)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## cafe         2  5.333  2.6667   3.582 0.0535 .
## Residuals   15 11.167  0.7444                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.5. Verificación de supuestos

Se evalúan los supuestos del ANOVA sobre los residuos: la prueba de Shapiro-Wilk (normalidad), la prueba de Bartlett (igualdad de varianzas entre marcas) y el gráfico Q-Q con su línea de referencia.

shapiro.test(residuals(modelo_sabor))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo_sabor)
## W = 0.98492, p-value = 0.9864
bartlett.test(sabor ~ cafe, data = datos_sabor)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  sabor by cafe
## Bartlett's K-squared = 0.71314, df = 2, p-value = 0.7001
qqnorm(residuals(modelo_sabor))
qqline(residuals(modelo_sabor))

¿Qué nos dice este gráfico? Es el gráfico Q-Q de los residuos del ANOVA del sabor, y se lee igual que el de la Parte I: si los puntos se pegan a la línea, los residuos son normales. Los puntos siguen la línea sin desviaciones fuertes (con los pequeños escalones que producen las notas enteras), lo que concuerda con la prueba de Shapiro-Wilk (p = 0.986) que aparece justo arriba. La prueba de Bartlett (p = 0.700) indica que las varianzas son homogéneas, así que el ANOVA es una prueba válida para estos datos.

3.6. Rangos

Se ordenan las 18 calificaciones de menor a mayor y se asigna un rango a cada una (los empates reciben el promedio de sus posiciones). Luego se calcula, por marca, la suma y el promedio de los rangos, que son la base de la prueba de Kruskal-Wallis.

datos_sabor$rango <- rank(datos_sabor$sabor)

tapply(datos_sabor$rango, datos_sabor$cafe, sum)
##   Café Córdoba     Sello Rojo Medalla de Oro 
##           56.0           77.5           37.5
tapply(datos_sabor$rango, datos_sabor$cafe, mean)
##   Café Córdoba     Sello Rojo Medalla de Oro 
##       9.333333      12.916667       6.250000

3.7. Prueba de Kruskal-Wallis

Se contrasta si las tres marcas tienen la misma distribución de calificaciones de sabor. El resultado incluye el estadístico chi-cuadrado (\(H\)), los grados de libertad y el valor p.

kw_sabor <- kruskal.test(sabor ~ cafe, data = datos_sabor)

kw_sabor
## 
##  Kruskal-Wallis rank sum test
## 
## data:  sabor by cafe
## Kruskal-Wallis chi-squared = 5.1966, df = 2, p-value = 0.0744

3.8. Tamaño del efecto

Se calcula el épsilon cuadrado dividiendo el estadístico \(H\) entre \(N - 1\), donde \(N\) es el número total de observaciones. Indica qué tan grande es la diferencia entre las marcas, independientemente del valor p.

N_sabor <- nrow(datos_sabor)

epsilon2_sabor <- as.numeric(kw_sabor$statistic) /
                  (N_sabor - 1)

epsilon2_sabor
## [1] 0.3056827

3.9. Prueba de Dunn

El código decide automáticamente: si el valor p de Kruskal-Wallis es menor que 0.05, aplica la prueba de Dunn con ajuste de Holm para identificar entre cuáles marcas hay diferencias; si no, muestra un mensaje indicando que no hay diferencias significativas.

if (kw_sabor$p.value < 0.05) {
  FSA::dunnTest(sabor ~ cafe, data = datos_sabor,
                method = "holm")
} else {
  cat("No se encontraron diferencias significativas en sabor.")
}
## No se encontraron diferencias significativas en sabor.

4. Parte II: Análisis no paramétrico del aroma

Se repite el procedimiento de la sección 3 con las calificaciones de aroma.

4.1. Ingreso de los datos

Se crea el vector aroma y se une con el factor cafe ya definido en datos_aroma.

aroma <- c(9, 8, 7, 8, 9, 8,
           7, 8, 9, 6, 8, 7,
           8, 7, 9, 8, 9, 6)

datos_aroma <- data.frame(cafe, aroma)

datos_aroma

4.2. Estadística descriptiva

Número de datos, media, mediana y desviación estándar del aroma de cada marca.

aggregate(aroma ~ cafe, data = datos_aroma,
          FUN = function(x) c(
            n = length(x),
            media = mean(x),
            mediana = median(x),
            desviacion = sd(x)
          ))

4.3. Gráfico de cajas

boxplot(aroma ~ cafe, data = datos_aroma,
        col = c("#F4B183", "#9DC3E6", "#A9D18E"),
        ylab = "Aroma (escala de 1 a 10)",
        xlab = "Tipo de café",
        main = "Calificación del aroma")

¿Qué nos dice este gráfico? Es el mismo diagrama de la Parte I, con un color por marca (naranja para Café Córdoba, azul para Sello Rojo y verde para Medalla de Oro). Las tres cajas están muy encimadas y las medianas son casi iguales (8, 7.5 y 8), con una caja más ancha en Medalla de Oro. No hay datos atípicos. En el aroma no se ve una marca claramente por encima de las otras.

4.4. ANOVA

modelo_aroma <- aov(aroma ~ cafe, data = datos_aroma)

summary(modelo_aroma)
##             Df Sum Sq Mean Sq F value Pr(>F)
## cafe         2  1.333  0.6667   0.659  0.532
## Residuals   15 15.167  1.0111

4.5. Verificación de supuestos

shapiro.test(residuals(modelo_aroma))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo_aroma)
## W = 0.97141, p-value = 0.8237
bartlett.test(aroma ~ cafe, data = datos_aroma)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  aroma by cafe
## Bartlett's K-squared = 0.88828, df = 2, p-value = 0.6414
qqnorm(residuals(modelo_aroma))
qqline(residuals(modelo_aroma))

¿Qué nos dice este gráfico? Es el gráfico Q-Q de los residuos del ANOVA del aroma. Los puntos siguen razonablemente la línea, sin curvas ni valores muy alejados, lo que concuerda con la prueba de Shapiro-Wilk (p = 0.824). Junto con la prueba de Bartlett (p = 0.641), indica que se cumplen los supuestos del ANOVA para el aroma.

4.6. Rangos

datos_aroma$rango <- rank(datos_aroma$aroma)

tapply(datos_aroma$rango, datos_aroma$cafe, sum)
##   Café Córdoba     Sello Rojo Medalla de Oro 
##           66.5           46.5           58.0
tapply(datos_aroma$rango, datos_aroma$cafe, mean)
##   Café Córdoba     Sello Rojo Medalla de Oro 
##      11.083333       7.750000       9.666667

4.7. Prueba de Kruskal-Wallis

kw_aroma <- kruskal.test(aroma ~ cafe, data = datos_aroma)

kw_aroma
## 
##  Kruskal-Wallis rank sum test
## 
## data:  aroma by cafe
## Kruskal-Wallis chi-squared = 1.2946, df = 2, p-value = 0.5235

4.8. Tamaño del efecto

N_aroma <- nrow(datos_aroma)

epsilon2_aroma <- as.numeric(kw_aroma$statistic) /
                  (N_aroma - 1)

epsilon2_aroma
## [1] 0.07615268

4.9. Prueba de Dunn

if (kw_aroma$p.value < 0.05) {
  FSA::dunnTest(aroma ~ cafe, data = datos_aroma,
                method = "holm")
} else {
  cat("No se encontraron diferencias significativas en aroma.")
}
## No se encontraron diferencias significativas en aroma.

5. Conclusiones

Todas las pruebas se interpretaron con un nivel de significancia de \(\alpha = 0.05\).

5.1. Sabor

Descripción. Sello Rojo obtuvo la calificación media más alta (8.83; mediana 9), seguida de Café Córdoba (8.17; mediana 8) y Medalla de Oro (7.50; mediana 7.5), que además fue la marca con mayor dispersión (desviación estándar de 1.05, frente a 0.75 de las otras dos). Los rangos promedio siguen el mismo orden: Sello Rojo (12.92), Café Córdoba (9.33) y Medalla de Oro (6.25).

Supuestos. Los residuos no se apartan de la normalidad (Shapiro-Wilk: \(W = 0.985\), \(p = 0.986\)) y las varianzas son homogéneas (Bartlett: \(p = 0.700\); Levene: \(p = 0.481\)), por lo que el ANOVA es una prueba válida para estos datos.

Pruebas globales. El ANOVA arrojó \(F(2, 15) = 3.58\) con \(p = 0.053\), y la prueba de Kruskal-Wallis arrojó \(H = 5.20\) con 2 grados de libertad y \(p = 0.074\). En ambos casos el valor p es mayor que 0.05, de modo que no se rechaza la hipótesis nula: con este nivel de significancia no hay evidencia estadística suficiente de que las tres marcas difieran en sabor. Como Kruskal-Wallis no fue significativa, la prueba de Dunn no se ejecutó.

Tamaño del efecto. El épsilon cuadrado fue de 0.31 (y el \(R^2\) del modelo lineal, de 0.32), un efecto grande. Esto indica que la marca explica cerca de un tercio de la variación del sabor, y que la falta de significancia se debe en buena parte al tamaño reducido de la muestra (6 catadores por marca) más que a que las marcas sean idénticas.

Comparaciones múltiples. Las pruebas LSD y Tukey encontraron una diferencia entre Sello Rojo y Medalla de Oro (diferencia de 1.33 puntos; Tukey: \(p = 0.043\)), mientras que Café Córdoba no se distinguió de ninguna de las dos. Este resultado debe interpretarse con cautela: se obtuvo después de un ANOVA que no fue significativo (aunque estuvo muy cerca) y la prueba de Kruskal-Wallis no la confirma. Lo prudente es considerarlo una tendencia a favor de Sello Rojo, que valdría la pena confirmar con más catadores.

5.2. Aroma

Descripción. Las medias fueron muy parecidas: Café Córdoba 8.17 (mediana 8), Medalla de Oro 7.83 (mediana 8) y Sello Rojo 7.50 (mediana 7.5). Las diferencias entre marcas son de menos de un punto, y la dispersión de los datos (desviación estándar entre 0.75 y 1.17) es mayor que esas diferencias. Los rangos promedio fueron 11.08 (Café Córdoba), 9.67 (Medalla de Oro) y 7.75 (Sello Rojo).

Supuestos. Se cumplen la normalidad (Shapiro-Wilk: \(W = 0.971\), \(p = 0.824\)) y la homogeneidad de varianzas (Bartlett: \(p = 0.641\); Levene: \(p = 0.568\)).

Pruebas globales. El ANOVA dio \(F(2, 15) = 0.66\) con \(p = 0.532\) y Kruskal-Wallis dio \(H = 1.29\) con \(p = 0.523\). Ambas pruebas coinciden en que no se rechaza la hipótesis nula: no existen diferencias significativas de aroma entre las tres marcas, por lo que la prueba de Dunn tampoco se ejecutó.

Tamaño del efecto. El épsilon cuadrado fue de 0.08 (y el \(R^2\), de 0.08), un efecto pequeño. Las comparaciones múltiples (LSD y Tukey) tampoco encontraron ningún par de marcas distinto (todos los valores p de Tukey fueron mayores que 0.50).

5.3. Conclusión general

Con los datos obtenidos, los catadores no percibieron diferencias estadísticamente significativas entre Café Córdoba, Sello Rojo y Medalla de Oro, ni en el sabor ni en el aroma. En el aroma las tres marcas se comportan de manera prácticamente equivalente. En el sabor sí aparece una tendencia, con Sello Rojo mejor calificado y Medalla de Oro con las calificaciones más bajas, que está en el límite de la significancia y cuyo tamaño de efecto es grande.

Hay que tener en cuenta algunas limitaciones: la muestra es pequeña (6 calificaciones por marca), las calificaciones son una escala ordinal con muchos valores repetidos (empates), y la significancia es sensible a esos factores. Por eso el resultado del sabor debería verificarse con un panel más grande de catadores antes de afirmar que una marca es mejor que otra. Además, los dos enfoques usados (ANOVA y Kruskal-Wallis) coinciden en sus conclusiones, lo que da más confianza a los resultados.