Análisis de la Varianza No Paramétrico

Facultad de Ciencias Agrarias - Universidad Nacional de Catamarca

---

1 Prueba de Kruskal - Wallis

• La prueba de Kruskal–Wallis es una alternativa no paramétrica al análisis de varianza (ANOVA) de una vía (DCA).
• Se emplea cuando los datos no cumplen los supuestos clásicos del ANOVA, es decir, normalidad de las distribuciones y homogeneidad de varianzas.
• Este test compara k grupos independientes (con \(k > 2\)) y permite evaluar si provienen de la misma población o si existen diferencias significativas entre ellos.
• A diferencia del ANOVA, se basa en rangos en lugar de los valores originales de las observaciones.

1.1 Fundamento metodológico

El procedimiento transforma los datos en rangos:

1. Se ordenan todas las observaciones en forma conjunta.

2. Se asignan rangos consecutivos (1 al valor más bajo, \(N\) al valor más alto).

3. En caso de empates, se asigna el rango promedio.

4. Finalmente, se comparan las sumas de rangos promedio de cada grupo.

5. Si las distribuciones de los grupos son similares, los rangos estarán distribuidos de manera equilibrada; si no, se observarán diferencias sistemáticas en los rangos.

1.2 Fórmula del estadístico

El estadístico de Kruskal–Wallis se denota como “\(𝐻\)”.

\[ H = \frac{12}{N(N+1)} \sum_{i=1}^{k} \frac{R_i^2}{n_i} - 3(N+1) \]

donde:

\[\begin{align*} k & = \text{número de grupos o tratamientos} \\ n_i & = \text{número de observaciones en el grupo } i \\ N & = \sum_{i=1}^{k} n_i \quad \text{(total de observaciones)} \\ R_i & = \text{suma de los rangos asignados a las observaciones del grupo } i \end{align*}\]

1.3 Hipótesis estadística

\(𝐻_0\): todas las poblaciones tienen la misma distribución (no hay diferencias entre tratamientos).

\(𝐻_1\): al menos una población difiere.

1.4 Análisis estadístico

1.4.1 Instalamos y cargamos paquetes

# install.packages("readxl")
# install.packages("tidyverse")
# install.packages("car")
# install.packages("agricolae")

library(tidyverse)
library(readxl)
library(car)
library(agricolae)

1.4.2 La base de datos

Esta base forma parte del paquete básico de R y contiene datos de recuentos de insectos vivos en unidades experimentales agrícolas tratadas con diferentes insecticidas.

InsectSprays

##    count spray
## 1     10     A
## 2      7     A
## 3     20     A
## 4     14     A
## 5     14     A
## 6     12     A
## 7     10     A
## 8     23     A
## 9     17     A
## 10    20     A
## 11    14     A
## 12    13     A
## 13    11     B
## 14    17     B
## 15    21     B
## 16    11     B
## 17    16     B
## 18    14     B
## 19    17     B
## 20    17     B
## 21    19     B
## 22    21     B
## 23     7     B
## 24    13     B
## 25     0     C
## 26     1     C
## 27     7     C
## 28     2     C
## 29     3     C
## 30     1     C
## 31     2     C
## 32     1     C
## 33     3     C
## 34     0     C
## 35     1     C
## 36     4     C
## 37     3     D
## 38     5     D
## 39    12     D
## 40     6     D
## 41     4     D
## 42     3     D
## 43     5     D
## 44     5     D
## 45     5     D
## 46     5     D
## 47     2     D
## 48     4     D
## 49     3     E
## 50     5     E
## 51     3     E
## 52     5     E
## 53     3     E
## 54     6     E
## 55     1     E
## 56     1     E
## 57     3     E
## 58     2     E
## 59     6     E
## 60     4     E
## 61    11     F
## 62     9     F
## 63    15     F
## 64    22     F
## 65    15     F
## 66    16     F
## 67    13     F
## 68    10     F
## 69    26     F
## 70    26     F
## 71    24     F
## 72    13     F

str(InsectSprays)

## 'data.frame':    72 obs. of  2 variables:
##  $ count: num  10 7 20 14 14 12 10 23 17 20 ...
##  $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...

unique(InsectSprays$spray)

## [1] A B C D E F
## Levels: A B C D E F

1.4.3 Visualización

# Gráfico 
ggplot(InsectSprays, aes(x = spray, y = count, fill = spray)) +
  geom_boxplot(alpha = 0.6) +
  geom_jitter(width = 0.1, size = 2, alpha = 0.4) +
  labs(
    title = "Número de insectos vivos según tipo de spray",
    x = "Spray",
    y = "Cantidad de insectos") +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none")

1.4.4 ANOVA

# ANOVA
anova_mod <- aov(count ~ spray, data = InsectSprays)
summary(anova_mod)

##             Df Sum Sq Mean Sq F value Pr(>F)    
## spray        5   2669   533.8    34.7 <2e-16 ***
## Residuals   66   1015    15.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.4.5 Supuestos del modelo

1.4.5.1 Diagnóstico gráfico

# Gráfico de diagnóstico
plot(anova_mod, which = 1) 

plot(anova_mod, which = 2) 

1.4.5.2 Pruebas formales para supuestos

# Shapiro-Wilk para normalidad
shapiro.test(residuals(anova_mod))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(anova_mod)
## W = 0.96006, p-value = 0.02226

# Test de Bartlett para homocedasticidad
bartlett.test(count ~ spray, data = InsectSprays)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  count by spray
## Bartlett's K-squared = 25.96, df = 5, p-value = 9.085e-05

# test de Levene para homocedasticidad (requiere paquete car)
leveneTest(count ~ spray, data = InsectSprays)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  5  3.8214 0.004223 **
##       66                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.4.6 Prueba de Kruskal-Wallis

Es importante destacar que la prueba no identifica qué grupos difieren entre sí. Para ello se requieren análisis post-hoc.

# Opción 1
kruskal.test(count ~ spray, data = InsectSprays)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  count by spray
## Kruskal-Wallis chi-squared = 54.691, df = 5, p-value = 1.511e-10

1.4.7 La función kruskal() del paquete agricolae

La función kruskal() constituye una extensión de la prueba de Kruskal–Wallis, con la ventaja de que incorpora de manera automática comparaciones múltiples entre tratamientos. Esto la convierte en una herramienta valiosa, ya que permite no solo detectar diferencias globales entre tratamientos, sino también identificar cuáles tratamientos son estadísticamente equivalentes o distintos, sin necesidad de realizar pasos adicionales, para ello utiliza un procedimiento de comparaciones múltiples adaptado a los rangos.

# Opción 2
# La función kruskal() incorpora la prueba de comparación múltiple. 
resultados_KW <- kruskal(InsectSprays$count, InsectSprays$spray)
resultados_KW

## $statistics
##      Chisq Df      p.chisq  t.value      MSD
##   54.69134  5 1.510845e-10 1.996564 8.462804
## 
## $parameters
##             test p.ajusted             name.t ntr alpha
##   Kruskal-Wallis      none InsectSprays$spray   6  0.05
## 
## $means
##   InsectSprays.count     rank      std  r Min Max   Q25  Q50   Q75
## A          14.500000 52.16667 4.719399 12   7  23 11.50 14.0 17.75
## B          15.333333 54.83333 4.271115 12   7  21 12.50 16.5 17.50
## C           2.083333 11.45833 1.975225 12   0   7  1.00  1.5  3.00
## D           4.916667 25.58333 2.503028 12   2  12  3.75  5.0  5.00
## E           3.500000 19.33333 1.732051 12   1   6  2.75  3.0  5.00
## F          16.666667 55.62500 6.213378 12   9  26 12.50 15.0 22.50
## 
## $comparison
## NULL
## 
## $groups
##   InsectSprays$count groups
## F           55.62500      a
## B           54.83333      a
## A           52.16667      a
## D           25.58333      b
## E           19.33333     bc
## C           11.45833      c
## 
## attr(,"class")
## [1] "group"

1.4.8 Interpretación de resultados

$statistics:

• Chi-sq = 54.69 con 5 grados de libertad (Df = 5).
• p.chisq = 1.51 × 10⁻¹⁰, menor que 𝛼 = 0.05. Se rechaza \(H_0\).

Existen diferencias significativas entre los tratamientos (sprays) en cuanto al número de insectos. No todos los insecticidas tienen el mismo efecto.

$means: Medias y rangos promedio

Por cada tratamiento (A–F) se listan:

• Promedio de insectos (ej. C = 2.08, F = 16.66).
• Rango promedio (el valor clave en Kruskal–Wallis, ej. C = 11.45, F = 55.62).
• Otros descriptivos: desviación estándar, n (r), mínimo, máximo, cuartiles.

Los rangos promedio son los que se comparan en el test. Un rango mayor → más insectos presentes.

$groups:Grupos de significancia

1.4.9 Conclusión

La prueba de Kruskal–Wallis evidenció diferencias altamente significativas entre los tratamientos insecticidas (χ² = 54.69; p < 0.001).

• Los sprays C, D y E redujeron significativamente la cantidad de insectos en comparación con A, B y F.
• El spray C fue el más eficaz, presentando los menores conteos de insectos.
• A, B y F resultaron ineficaces, sin diferencias significativas entre ellos.

2 Prueba de Friedman

A continuación realizaremos el análisis no paramétrico para el caso en que el diseño tenga un factor a poner a prueba con tres o más tratamientos, estructura de parcela (DBCA) y no se cumplen los supuestos requeridos para el ANOVA.

2.1 Hipótesis estadística

H₀: Las medianas de todos los tratamientos son iguales (No hay diferencias entre los tratamientos).

H₁: Al menos una de las medianas de los tratamientos es diferente (Al menos un tratamiento difiere de los otros).

Si rechazamos la hipótesis nula la pregunta es: ¿cual Mediana diferente?

2.2 Estadístico

\[ H = \frac{12}{b \, k \, (k+1)} \sum_{i=1}^{k} R_i^2 - 3 \, b \, (k+1) \]

donde:

• \(b\): número de bloques
  
• \(k\): número de tratamientos
  
• \(R_i\): suma de rangos del tratamiento \(i\)

2.3 Caso práctico

En un estudio sobre la producción de madera en Pinus spp. se plantaron árboles en un campo ubicado en una zona montañosa con 3 sectores con relieve claramente diferenciados: ladera oriental, occidental y un valle. Allí se plantaron los árboles de la especie y se le aplicaron cuatro tratamientos: control(sin micorrizas), micorriza A, micorriza B y A+B. El resultado, expresado en peso total fué:

Tabla 1. Datos de Pinus spp.

Tratamientos	Oriental	Occidental	Valle
Control	450	320	500
Micorriza A	560	480	660
Micorriza B	500	390	490
Micorriza A+B	660	560	690

Los ingenieros forestales desean conocer si el peso totalsegún los tratamientos es igual, o si en algún grupo difiere. Por lo tanto se plantea:

\(H_0\): No hay diferencias entre tratamientos (todas las medianas son iguales).
\(H_1\): Al menos un tratamiento difiere.

2.4 Cargamos la base de datos

PINOS <- read_excel("pinos.xlsx")

2.4.1 Transformación de variables

str(PINOS)

## tibble [12 × 3] (S3: tbl_df/tbl/data.frame)
##  $ UBICACION   : chr [1:12] "oriental" "occidental" "valle" "oriental" ...
##  $ TRATAMIENTOS: chr [1:12] "control" "control" "control" "micorriza_A" ...
##  $ PESO        : num [1:12] 450 320 500 560 480 660 500 390 490 660 ...

PINOS <- PINOS %>% 
  mutate( 
    UBICACION = factor(UBICACION),
    TRATAMIENTOS = factor(TRATAMIENTOS))

str(PINOS)

## tibble [12 × 3] (S3: tbl_df/tbl/data.frame)
##  $ UBICACION   : Factor w/ 3 levels "occidental","oriental",..: 2 1 3 2 1 3 2 1 3 2 ...
##  $ TRATAMIENTOS: Factor w/ 4 levels "control","micorriza_A",..: 1 1 1 2 2 2 4 4 4 3 ...
##  $ PESO        : num [1:12] 450 320 500 560 480 660 500 390 490 660 ...

2.5 Visualización

ggplot(PINOS, aes(TRATAMIENTOS, PESO)) +
  geom_boxplot() + 
  coord_flip() + # cambiamos a posición horizontal las cajas
  scale_y_continuous(limits = c(300, 700)) +
  theme_light()

2.6 La función friedman del paquete agricolae

library(agricolae)
# Usamos la función del paquete agricolae::friedman
FRIEDMAN <- friedman(PINOS$UBICACION, PINOS$TRATAMIENTOS,PINOS$PESO, alpha = 0.05, group = TRUE)
FRIEDMAN

## $statistics
##   Chisq Df    p.chisq    F DFerror         p.F  t.value     LSD
##     8.2  3 0.04205418 20.5       6 0.001484194 2.446912 2.82545
## 
## $parameters
##       test             name.t ntr alpha
##   Friedman PINOS$TRATAMIENTOS   4  0.05
## 
## $means
##               PINOS.PESO rankSum      std r Min Max Q25 Q50 Q75
## control         423.3333       4 92.91573 3 320 500 385 450 475
## micorriza_A     566.6667       9 90.18500 3 480 660 520 560 610
## micorriza_AyB   636.6667      12 68.06859 3 560 690 610 660 675
## micorriza_B     460.0000       5 60.82763 3 390 500 440 490 495
## 
## $comparison
## NULL
## 
## $groups
##               Sum of ranks groups
## micorriza_AyB           12      a
## micorriza_A              9      b
## micorriza_B              5      c
## control                  4      c
## 
## attr(,"class")
## [1] "group"

• El valor del estadístico \(H\) se compara con una \(\chi^2\) con \(k-1\) gl.
• Si el valor-p < 0.05, se rechaza \(H_0\).
• En este caso, el paquete agricolae devuelve directamente los grupos de tratamientos formados por comparación múltiple.

A MODO DE CIERRE…

La culminación de esta Diplomatura en Bioestadística con Software R, representa no solo la adquisición de herramientas estadísticas aplicadas al análisis de datos en las ciencias biológicas, agronómicas y de la salud, sino también la consolidación de una comunidad académica orientada al intercambio, la colaboración y el aprendizaje permanente.

Confiamos en que esta experiencia constituya el punto de partida de nuevos proyectos y aplicaciones profesionales que contribuyan al desarrollo científico, social y personal de las y los participantes.