##Analisis de varianza (ANOVA) # ANOVA

Ingreso de datos

UCS <- c(115, 125, 83, 105, 115, 96, 78, 160, 135, 100,
         91, 76, 57, 80, 110, 42, 76, 70, 95, 70,
         46, 91, 55, 46, 42, 87, 67, 67, 55, 80,
         80, 130, 92, 78, 53, 74, 68, 115, 78, 53
)
UCS
##  [1] 115 125  83 105 115  96  78 160 135 100  91  76  57  80 110  42  76  70  95
## [20]  70  46  91  55  46  42  87  67  67  55  80  80 130  92  78  53  74  68 115
## [39]  78  53
Bloques <- c(rep("Bq1",10),rep("Bq2",10),
                 rep("Bq3",10),rep("Bq4",10))
Bloques
##  [1] "Bq1" "Bq1" "Bq1" "Bq1" "Bq1" "Bq1" "Bq1" "Bq1" "Bq1" "Bq1" "Bq2" "Bq2"
## [13] "Bq2" "Bq2" "Bq2" "Bq2" "Bq2" "Bq2" "Bq2" "Bq2" "Bq3" "Bq3" "Bq3" "Bq3"
## [25] "Bq3" "Bq3" "Bq3" "Bq3" "Bq3" "Bq3" "Bq4" "Bq4" "Bq4" "Bq4" "Bq4" "Bq4"
## [37] "Bq4" "Bq4" "Bq4" "Bq4"

###Diagrama de caja o bigote

boxplot(UCS ~ Bloques, col=rainbow(4),
        main="Diagrama de Caja")

Prueba Anova usando R

\(H_0:\quad{\mu_1}={\mu_2}={\mu_3}={\mu_4}\)

\(H_1: \quad\text{Al menos uno es diferente}\)

Compilamos el modelo ANOVA mediante la función aov, así:

anova1 <- aov(UCS ~ Bloques)
summary(anova1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Bloques      3  12115    4038   8.503 0.000212 ***
## Residuals   36  17097     475                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova1 <- aov(UCS ~ Bloques) summary(anova1)

str(anova1)##str muestra la estructura completa

Si al detectar que existen diferencias significativas, se procede a realizar el análisis de Tukey ## Análisis de Tukey, este análisi de diferencias honestamente significativas.

Análsis de Tukey (HSD)

library(stats)
par(cex=0.8)##con cex se hace la letra más pequeña
intervalos <- TukeyHSD(anova1)
intervalos
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = UCS ~ Bloques)
## 
## $Bloques
##          diff        lwr        upr     p adj
## Bq2-Bq1 -34.5 -60.748029  -8.251971 0.0059294
## Bq3-Bq1 -47.6 -73.848029 -21.351971 0.0001216
## Bq4-Bq1 -29.1 -55.348029  -2.851971 0.0249272
## Bq3-Bq2 -13.1 -39.348029  13.148029 0.5416747
## Bq4-Bq2   5.4 -20.848029  31.648029 0.9448355
## Bq4-Bq3  18.5  -7.748029  44.748029 0.2467914
str(intervalos)
## List of 1
##  $ Bloques: num [1:6, 1:4] -34.5 -47.6 -29.1 -13.1 5.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:6] "Bq2-Bq1" "Bq3-Bq1" "Bq4-Bq1" "Bq3-Bq2" ...
##   .. ..$ : chr [1:4] "diff" "lwr" "upr" "p adj"
##  - attr(*, "class")= chr [1:2] "TukeyHSD" "multicomp"
##  - attr(*, "orig.call")= language aov(formula = UCS ~ Bloques)
##  - attr(*, "conf.level")= num 0.95
##  - attr(*, "ordered")= logi FALSE
plot(intervalos,col=rainbow(6))

Validación del modelo ANOVA

A partir de los residuos del modelo se comprobará, si el modelo ANOVA es adecuado. Los supuestos que se deben cumplir son tres: independencia, homocedasticidad y normalidad.

Análisis de residuos

residuos <- anova1$residuals
residuos
##     1     2     3     4     5     6     7     8     9    10    11    12    13 
##   3.8  13.8 -28.2  -6.2   3.8 -15.2 -33.2  48.8  23.8 -11.2  14.3  -0.7 -19.7 
##    14    15    16    17    18    19    20    21    22    23    24    25    26 
##   3.3  33.3 -34.7  -0.7  -6.7  18.3  -6.7 -17.6  27.4  -8.6 -17.6 -21.6  23.4 
##    27    28    29    30    31    32    33    34    35    36    37    38    39 
##   3.4   3.4  -8.6  16.4  -2.1  47.9   9.9  -4.1 -29.1  -8.1 -14.1  32.9  -4.1 
##    40 
## -29.1

Análisis de la normalidad de los residuos

par(mfrow=c(2,2), cex=0.7)

plot(residuos, pch=20, col=4, 
     main = "Diagrama de residuos",
     xlab="Índice")

abline(h=0, col=2)
abline(h=30, col=5)
abline(h=-30, col=5)
par(mfrow=c(2,2), cex=0.7)

plot(residuos, pch=20, col=4, 
     main = "Diagrama de residuos",
     xlab="Índice")

abline(h=0, col=2)
abline(h=30, col=5)
abline(h=-30, col=5)

hist(residuos,main="Histograma de los residuos")

qqnorm(residuos, pch=20, col=4)
qqline(residuos, col=2)

boxplot(residuos, horizontal = T, main="Diagrama de caja")

par(mfrow=c(1,1))

Test de Normalidad de los residuos

\(H_0:\) Los residuos siguen una distribución normal

\(H_1:\) Los residuos no siguen una distribución normal

shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.96659, p-value = 0.2791

El valor-p es 0,2791, mayor que el nivel de significancia por lo cual no se rechaza la hipótesis nula, eso indica que los residuos siguen una distribución normal.

Homocedasticidad de los residuos

\(H_O:\) La varianza de los residuos es constante

\(H_1:\) La varianza de los residuos no es constante

Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos, sin embargo se realiza la prueba Barlett.

Test de Bartlett

bartlett.test(residuos ~ Bloques)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuos by Bloques
## Bartlett's K-squared = 1.4513, df = 3, p-value = 0.6936

El test de Bartlett indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (la varianza es constante).

Otra prueba para verificar la homocedasticidad es la de Breush-Pagan.la cual se encuentra en la libreria lmtest.

Prueba de Breush-Pagan

\(H_0:\) La varianza de los residuos es constante (homocedasticidad)

\(H_1:\) La varianza de los residuos no es constante(heterocedasticidad)

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(formula = anova1)
## 
##  studentized Breusch-Pagan test
## 
## data:  anova1
## BP = 1.7998, df = 3, p-value = 0.615

El valor-p es mayor que el nivel de significancia por lo cual no se tiene evidencia estadistica suficiente para rechazar la hipótesis nula, eso indica que la varianza de los residuos se considera constante.

Independencia de los residuos

En la gráfica de los residuos de manera visual aparentemente no hay ninguna tendencia sin embargo se procede a realizar la prueba de indepedencia o de autocorrelación de Durbin-Watson con la función dwtest, así:

Prueba de Durbin-Watson

\(H_0:\) Los residuos son independientes (No están autocorrelacionados)

\(H_1:\) Los residuos no son independientes (Están autocorrelacionados)

dwtest(anova1)
## 
##  Durbin-Watson test
## 
## data:  anova1
## DW = 2.2201, p-value = 0.5867
## alternative hypothesis: true autocorrelation is greater than 0

dwtest(anova1)

El valor-p es mayor que el nivel de significancia por lo cual no se tiene evidencia estadística suficiente para rechazar la hipótesis nula,eso indica que los residuos son independientes.

Identificación de los valores Atípicos (outliers)

Primero se detecta si hay valores atipicos, para este caso se usará la función outlierTest de la librería car, así:

library(car)
## Loading required package: carData
outlierTest(model = anova1)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 8 2.531537           0.016002      0.64006

Si el valor de|r-Student| > 3 atípico, se trata de un valor atípico.

La observación 8 tiene un valor \(|r-Student|\) menor que 3, por lo cual no se considera atípico.

Si se detectaría que hay la presencia de valores atípico o outliers, se procede a verificar si éstos son influyentes, usando la función:influence.measures(model = anova1), por ejemplo:

summary(influence.measures(model = anova1))
## Potentially influential observations of
##   aov(formula = UCS ~ Bloques) :
## 
##    dfb.1_ dfb.BlB2 dfb.BlB3 dfb.BlB4 dffit cov.r   cook.d hat  
## 8   0.84  -0.60    -0.60    -0.60     0.84  0.63_*  0.15   0.10
## 32  0.00   0.00     0.00     0.58     0.83  0.65_*  0.15   0.10

Para validar si los valores atípicos son influyentes se emplean los siguientes párametros: > La distancia de Cook la cual no debe superar el valor de uno (cook.d>1), si supera este valor se trata de un “valor atípico influyente”

El segundo valor es el valor hat, que se calcula de la siguiente manera

\[hat > 2.5(p+1)/n \] donde : p número de variables independientes, n es número de observaciones

n=length(residuos)
n
## [1] 40
hat0 <- 2.5*(4+1)/n
hat0  
## [1] 0.3125
hatvalues(model = anova1)
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
## 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

Alli podemos observar que ninguna observación es outlier influyente

Visualizacion de los valores atípicos

influencePlot(model = anova1, main="Valores atípico influyentes")

##       StudRes Hat       CookD
## 2   0.6622717 0.1 0.012376440
## 4  -0.2960654 0.1 0.002498164
## 8   2.5315372 0.1 0.154766593
## 32  2.4765767 0.1 0.149110630

Al validar todas esas pruebas, podemos concluir que el modelo está muy bien definido.