Análisis de varianza de una vía

Dataset: \(\texttt{mollusc.txt}\)

Considere las mediciones de la tasa de impulso nervioso de 125 moluscos pertenecientes a 5 especies diferentes.

Columna 1: Nombre de la especie
Columna 2: Especie codificada (de 1 a 5)
Columna 3: Número de identificación dentro de cada especie
Columna 4: Tasa de impulso nervioso

La pregunta de interés es si la tasa de impulso nervioso depende de la especie del molusco.

1) Realice una exploración descriptiva y gráfica de los datos

mollusc = read.table("mollusc.txt", header = TRUE)
attach(mollusc)

head(mollusc)

##      Nombre     Especie Codificacion Codigo  Tasa
## 1 Ariolimax Columbianus            1      1 48.14
## 2 Ariolimax Columbianus            1      2 35.50
## 3 Ariolimax Columbianus            1      3 46.40
## 4 Ariolimax Columbianus            1      4 56.00
## 5 Ariolimax Columbianus            1      5 39.68
## 6 Ariolimax Columbianus            1      6 51.46

summary(mollusc)

##     Nombre            Especie           Codificacion       Codigo     
##  Length:125         Length:125         Min.   :1.000   Min.   : 1.00  
##  Class :character   Class :character   1st Qu.:1.000   1st Qu.: 7.00  
##  Mode  :character   Mode  :character   Median :2.000   Median :15.00  
##                                        Mean   :2.592   Mean   :18.89  
##                                        3rd Qu.:5.000   3rd Qu.:29.00  
##                                        Max.   :5.000   Max.   :54.00  
##       Tasa      
##  Min.   : 20.0  
##  1st Qu.: 46.4  
##  Median : 84.6  
##  Mean   :173.9  
##  3rd Qu.:325.0  
##  Max.   :800.0

pairs(mollusc[,c(3, 5)], panel = function(x, y) {
points(x, y)
lines(lowess(x, y), col = "red")
})

plot(Codificacion, Tasa)

1.1) ¿Cuántas observaciones existen para cada especie?

for (i in 1:5) {
  print(as.character(unique(mollusc$Codificacion[mollusc$Codificacion == i])))
  print(length(mollusc$Codificacion[mollusc$Codificacion == i]))
}

## [1] "1"
## [1] 54
## [1] "2"
## [1] 19
## [1] "3"
## [1] 10
## [1] "4"
## [1] 8
## [1] "5"
## [1] 34

1.2) ¿Cuánto es el promedio estimado de la tasa de impulso nervioso en cada especie?

for (i in 1:5) {
  print(as.character(unique(mollusc$Codificacion[mollusc$Codificacion == i])))
  print(mean(mollusc$Tasa[mollusc$Codificacion == i]))
}

## [1] "1"
## [1] 43.54963
## [1] "2"
## [1] 123.7421
## [1] "3"
## [1] 78.435
## [1] "4"
## [1] 203.5
## [1] "5"
## [1] 429.9676

1.3) Realice un diagrama de cajas por especie. ¿Sugiere el gráfico que la tasa de impulso nervioso depende de la especie del molusco?

boxplot(Tasa~Codificacion, data = mollusc)
for (i in 1:5) {
  abline(h = median(mollusc$Tasa[mollusc$Codificacion == i]), lty = 2, col = "red") 
}

Se puede intuir un efecto de dependencia entre la tasa de impulso nervioso y la especie del molusco.

2) Estime el modelo usando la función lm. Compare los coeficientes β estimados a los promedios obtenidos en el literal 1.2

mollusc$Codificacion = factor(mollusc$Codificacion)
modelo0 = lm(Tasa~Codificacion, data = mollusc)
summary(modelo0)

## 
## Call:
## lm(formula = Tasa ~ Codificacion, data = mollusc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -117.37  -22.34   -3.95   14.58  370.03 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     43.550      9.245   4.711 6.70e-06 ***
## Codificacion2   80.192     18.122   4.425 2.13e-05 ***
## Codificacion3   34.885     23.389   1.492    0.138    
## Codificacion4  159.950     25.738   6.215 7.71e-09 ***
## Codificacion5  386.418     14.874  25.980  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 67.94 on 120 degrees of freedom
## Multiple R-squared:  0.856,  Adjusted R-squared:  0.8512 
## F-statistic: 178.4 on 4 and 120 DF,  p-value: < 2.2e-16

data.frame(
  Medias_estimadas = aggregate(Tasa~Codificacion, data = mollusc, FUN = mean),
  B_estimados = modelo0$coefficients
)

##               Medias_estimadas.Codificacion Medias_estimadas.Tasa B_estimados
## (Intercept)                               1              43.54963    43.54963
## Codificacion2                             2             123.74211    80.19248
## Codificacion3                             3              78.43500    34.88537
## Codificacion4                             4             203.50000   159.95037
## Codificacion5                             5             429.96765   386.41802

Solo el estimador de \(\beta_0\) con el estimador de \(\mu_1\) coinciden, el resto de coeficientes de regresión difieren de las medias respectivas a cada clase.

3) Realice la prueba F del modelo ANOVA. ¿la tasa de impulso nervioso depende de la especie del molusco?

modelo_ANOVA = aov(Tasa~Codificacion, data = mollusc)
summary(modelo_ANOVA)

##               Df  Sum Sq Mean Sq F value Pr(>F)    
## Codificacion   4 3292873  823218   178.4 <2e-16 ***
## Residuals    120  553872    4616                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La clase de especie tiene un efecto en la tasa de impulso nervioso, debido a que al menos una media en las clases es diferente al resto.

4) Realice una comparación de pares en base al método de Tukey. ¿cuáles especies tiene una tasa de impulso nervioso distinta?

TukeyHSD(modelo_ANOVA)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Tasa ~ Codificacion, data = mollusc)
## 
## $Codificacion
##          diff          lwr       upr     p adj
## 2-1  80.19248   30.0005857 130.38437 0.0002047
## 3-1  34.88537  -29.8943258  99.66507 0.5699724
## 4-1 159.95037   88.6651049 231.23564 0.0000001
## 5-1 386.41802  345.2224208 427.61361 0.0000000
## 3-2 -45.30711 -118.8207573  28.20655 0.4335997
## 4-2  79.75789    0.4519236 159.06387 0.0479753
## 5-2 306.22554  252.3281840 360.12290 0.0000000
## 4-3 125.06500   35.8090743 214.32093 0.0015767
## 5-3 351.53265  283.8413744 419.22392 0.0000000
## 5-4 226.46765  152.5265325 300.40876 0.0000000

Con un nivel de significancia \(\alpha = 0.05\) las especie que tienen una tasa de impulso nervioso distinta son:

La especie 3 con la especie 1
La especie 3 con la especie 2

5) Compruebe los supuestos del modelo ANOVA. ¿Sus conclusiones de las preguntas 3 y 4 serían válidas?

# Supuesto de homocedasticidad.
plot(fitted.values(modelo_ANOVA),
rstandard(modelo_ANOVA),
xlab =" Fitted values",
ylab =" Studentized residuals",
main =" Residuals vs fitted values plot",
ylim = c(5, -5))
abline(h = 0, lty = "dashed", col = "blue")

Existen problemas de homocedasticidad debido a que hay distinta variabilidad en los niveles.

# Supuesto de independencia.
n_T = nrow(mollusc)
plot(rstandard(modelo_ANOVA)[-c(1)],
rstandard(modelo_ANOVA)[-c(n_T)],
xlab = "Studentized residuals at 1 lag",
ylab = "Studentized residuals",
main = "Sequence plot")
abline (a = 0, b = 1, lty = "dashed", col = "green")

# Supuesto de normalidad.
qqnorm(residuals(modelo_ANOVA))
qqline(residuals(modelo_ANOVA))

Se puede observar un fuerte problema de normalidad de los residuos, en el gráfico se evidencia un problema brusco de colas pesadas.

leveneTest(Tasa ~ Codificacion, data = mollusc)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   4  12.919 8.876e-09 ***
##       120                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Existen evidencia suficiente para decir que al menos existe un nivel con distinta variabilidad o varianza distinta.

6) Si los supuestos del modelo no se cumplen, realice la prueba no paramétrica de Kruskal-Wallis. ¿Llega a la misma conclusión de la pregunta 3?

kruskal.test(Tasa ~ Codificacion, data = mollusc)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Tasa by Codificacion
## Kruskal-Wallis chi-squared = 105.6, df = 4, p-value < 2.2e-16

Si se concluye igual que en la pregunta 3, de que al menos hay una media de un nivel que es diferente.}

7) Revise el siguiente tutorial []https://www.tutorialspoint.com/how-to-performpost- hoc-test-for-kruskal-wallis-in-r para realizar una prueba post-hoc no parámetrica basada en el test de Dunn. ¿Llega a la misma conclusión de la pregunta 4?

dunnTest(Tasa ~ Codificacion, data = mollusc)

## Dunn (1964) Kruskal-Wallis multiple comparison

##   p-values adjusted with the Holm method.

##    Comparison           Z      P.unadj        P.adj
## 1       1 - 2  -4.3451821 1.391602e-05 1.252442e-04
## 2       1 - 3  -2.4631774 1.377118e-02 6.885588e-02
## 3       2 - 3   0.7961653 4.259360e-01 4.259360e-01
## 4       1 - 4  -4.1255427 3.698615e-05 2.958892e-04
## 5       2 - 4  -0.9582823 3.379204e-01 6.758409e-01
## 6       3 - 4  -1.5071999 1.317594e-01 3.952783e-01
## 7       1 - 5 -10.0387941 1.029246e-23 1.029246e-22
## 8       2 - 5  -3.6265454 2.872383e-04 1.723430e-03
## 9       3 - 5  -3.7521858 1.752995e-04 1.227096e-03
## 10      4 - 5  -1.6156601 1.061678e-01 4.246713e-01

# media1 - media2
mean(mollusc$Tasa[mollusc$Codificacion == 1]) - mean(mollusc$Tasa[mollusc$Codificacion == 2])

## [1] -80.19248

8) ¿Qué método escoje para dar conclusiones, el paramétrico o el no paramétrico?

Debido a que no se cumplen algunos supuestos del modelo anova, lo mejor seria trabajar y dar conclusiones con métodos no paramétricos.