DCA en R

Vamos a simular el experimento del fertilizante:

set.seed(123)

fertilizante_A <- c(118, 119, 120, 121, 122)
fertilizante_B <- c(121, 119, 120, 122, 118)
fertilizante_C <- c(135, 142, 138, 140, 145)

dca <- data.frame(
  rendimiento = c(fertilizante_A , fertilizante_B, fertilizante_C),
  grupo       = c(rep("A", 5), rep("B", 5), rep("C", 5))
)
dca
##    rendimiento grupo
## 1          118     A
## 2          119     A
## 3          120     A
## 4          121     A
## 5          122     A
## 6          121     B
## 7          119     B
## 8          120     B
## 9          122     B
## 10         118     B
## 11         135     C
## 12         142     C
## 13         138     C
## 14         140     C
## 15         145     C
str(dca)
## 'data.frame':    15 obs. of  2 variables:
##  $ rendimiento: num  118 119 120 121 122 121 119 120 122 118 ...
##  $ grupo      : chr  "A" "A" "A" "A" ...

H₀ (hipótesis nula): los datos SÍ son normales H₁ (hipótesis alterna): los datos NO son normales

Si el p-valor > 0.05 → no rechazamos H₀ → los datos son normales → podemos usar ANOVA.

shapiro.test(fertilizante_A)
## 
##  Shapiro-Wilk normality test
## 
## data:  fertilizante_A
## W = 0.98676, p-value = 0.9672
shapiro.test(fertilizante_B)
## 
##  Shapiro-Wilk normality test
## 
## data:  fertilizante_B
## W = 0.98676, p-value = 0.9672
shapiro.test(fertilizante_C)
## 
##  Shapiro-Wilk normality test
## 
## data:  fertilizante_C
## W = 0.99897, p-value = 0.9996

ANOVA

modelo_anova <- aov(rendimiento ~ grupo, data = dca)
summary(modelo_anova)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## grupo        2   1333   666.7   102.6 2.85e-08 ***
## Residuals   12     78     6.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cuando tenemos un p-valor por debajo de 0.05 rechazamos la hipotesis nula, es decir H₀ = los fertilizantes son IGUALES (no hay diferencia) H₁ = al menos un fertilizante es DIFERENTE como tenemos un p-valor brutalmente por debajo de 0.05 si tenemos el almenos un fertilizante diferente; pero cual es..

TukeyHSD(modelo_anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = rendimiento ~ grupo, data = dca)
## 
## $grupo
##             diff       lwr       upr p adj
## B-A 2.842171e-14 -4.301801  4.301801 1e+00
## C-A 2.000000e+01 15.698199 24.301801 1e-07
## C-B 2.000000e+01 15.698199 24.301801 1e-07

Kruskal-Wallis

Recodermos el flujo de decision:

¿Datos normales? SÍ → ANOVA NO → Kruskal-Wallis

Kruskal-Wallis hace lo mismo que ANOVA pero usando rangos, igual que Spearman usaba rangos en lugar de valores.

kruskal.test(rendimiento ~ grupo, data = dca)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  rendimiento by grupo
## Kruskal-Wallis chi-squared = 9.4595, df = 2, p-value = 0.008829

ANOVA detectó la diferencia con mucha más contundencia.

Potencia estadística = capacidad de detectar una diferencia real cuando realmente existe.

En este caso, con datos normales:

Cuando los datos NO son normales uso Kruskal-Wallis porque ANOVA asume normalidad y si ese supuesto se viola los resultados de ANOVA son poco confiables

En estadistica “poco confiable” significa que podria concluir que hay diferencias cuando no las hay, o que no las hay cuando si existen, eso es grave en un experimento real.

Flujo de decision completo

¿Los datos son normales? (Shapiro-Wilk) │ p > 0.05 p < 0.05 SÍ son normales NO son normales │ │ ANOVA Kruskal-Wallis │ │ Tukey HSD (post-hoc: (post-hoc) Dunn test)

Dos problemas de datos experimentales

Muy pocas muestras → vulnerable a outliers Demasiadas → costoso e innecesario

Concepto de potencia

En el concepto de potencia en estadistica radica en la capacidad para detectar diferecias reales. Con pocas muestras la potencia es baja, puedo tener diferencias reales entre fertilizantes y no detectarlas, con diferentes muestras la potencia sube.

library(pwr)

pwr.anova.test(
  k = 3,        # número de grupos
  f = 0.5,      # tamaño del efecto
  sig.level = 0.05,  # alpha
  power = 0.80       # potencia deseada
)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 3
##               n = 13.89521
##               f = 0.5
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

f = 0.5 es el tamaño del efecto, una medida estandarizada de que tan diferentes son las diferencias entre grupos.

f grande → diferencias GRANDES entre grupos → más fácil detectarlas → necesitas MENOS muestras f pequeño → diferencias PEQUEÑAS entre grupos → más difícil detectarlas → necesitas MÁS muestras

Lo vamos a pensar con los fertilizantes:

Si C produce 1000 kg más que A → diferencia obvia → pocas muestras bastan Si C produce 1 kg más que A → diferencia mínima → necesitas muchas muestras para estar seguro

En estadística los valores de f tienen una convención:

valores_f <- data.frame(  
  f = c(0.10 , 0.20, 0.40),
  tamaño = c("pequeño", "mediano", "grande")
)
valores_f
##     f  tamaño
## 1 0.1 pequeño
## 2 0.2 mediano
## 3 0.4  grande
library(pwr)
pwr.anova.test(
  k = 3,
  f = 0.5,
  sig.level = 0.05,
  power = 0.80
)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 3
##               n = 13.89521
##               f = 0.5
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

Ahora si la actualizo al 0.9 vamos a ver a si ya no son ~14 muestras

pwr.anova.test(
  k = 3,
  f = 0.5,
  sig.level = 0.05,
  power = 0.90
)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 3
##               n = 17.91287
##               f = 0.5
##       sig.level = 0.05
##           power = 0.9
## 
## NOTE: n is number in each group

entonces ya gracias a anova se cuantas muestras necesito para tener un intervalo de confianza que sea suficiente para llevar a cabo un experimento valido sin perdidas de dinero exageradas.

Flujo aprendido

  1. Calcular tamaño de muestra → pwr.anova.test() ↓
  2. Recolectar datos → data.frame() ↓
  3. Verificar normalidad → shapiro.test() ↓ Normal? SÍ → ANOVA → aov() → TukeyHSD() NO → Kruskal-Wallis → kruskal.test()