El Análisis de varianza es una de las herramientas más usadas y en general bien usada. A continuación desarrollo un ejemplo incorporando la libreria ANALITICA y solo se hará una tarea procedimental.

Modelo Diseño Completamente al Azar (DCA)

Diseño Completamente al Azar (DCA)

el modelo que se plantea es lineal y aditivo, del tipo

\(Y_{ij} = \mu + \tau_i + \varepsilon_{ij} \\\) donde \(\\\) \(\mu =\) media general \(\\\) \(\tau_i =\) effecto del tratamiento \(\\\) \(\varepsilon_{ij} =\) residuales \(\\\)

Librerias

¿Que librerias se deben ocupar?

library(tidyverse) 
library(Analitica)
library(nortest)
library(knitr) # para generar tablas
library(kableExtra)

Hipótesis Investigación

se debe establecer cual es la meta por alcanzar

\(H_i\) existe al menos un grupo que presenta mejores características de VD

Hipótesis estadística

\(H_0\) todos los grupos presentan iguales características \(\\\) \(H_1\) existe al menos un grupo que presenta mejores características de VD \(\\\)

esto se traduce en

\(H_0: \mu_1 =\mu_2=\mu_3 \\\) \(H_i: \mu_i \ne \mu_j\)

Analisis de Supuestos

S0 Calidad de Datos

Analizamos comportamientos anómalos y/o presencia de Valores Atípicos

VD = corresponde a la Variable Dependiente o Variable de Respuesta o de Medida

grupo = tratamientos.

descripYG(datos,VD,grupo)

##   Group  n    Mean Median        SD Kurtosis    Skewness          CV    Min
## 1     1 10 185.259 185.49 0.8474072 2.857842 -1.05079154 0.004574175 183.52
## 2     2 10 260.937 260.84 0.9599196 2.308984  0.53183290 0.003678741 259.78
## 3     3 10 150.129 150.16 1.9977401 2.404458 -0.06583392 0.013306824 146.83
##      Max      P25      P75    IQR
## 1 186.16 185.0125 185.8475 0.8350
## 2 262.68 260.2450 261.2875 1.0425
## 3 153.58 149.2300 151.4100 2.1800

a través del diagrama de caja, se observa la inexistencia de valores atípicos, en tanto los indices tienen comportamiento esperables. por tanto se puede continuar con el análisis con todos los datos.

S1 Normalidad

La hipótesis en cuestion es \(H_{0S1}: X \sim N(\mu,\sigma^2)\)

para el análisis de este supuesto es necesario presentar el modelo de análisis, el cual corresponde a un proceso paramétrico ANOVA

NOTA: recuerde que para que los resultados sean válidos, el df (degree free) siempre debe equivaler al numero de grupos menos 1, para este caso la variable grupo es númerica, por tanto debe transformarla a factor o aplicar la funcion as.factor

Forma incorrecta

mod1<-aov(VD~grupo,data = datos)
summary(mod1)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## grupo        1   6171    6171   2.978 0.0954 .
## Residuals   28  58012    2072                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Forma correcta

mod1<-aov(VD~as.factor(grupo),data = datos)
summary(mod1)
##                  Df Sum Sq Mean Sq F value Pr(>F)    
## as.factor(grupo)  2  64132   32066   17085 <2e-16 ***
## Residuals        27     51       2                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

luego se procede a la realización de algún test de normalidad (solo aplique uno), se aconseja revisar las características de los test (Shapiro-Wilk, Shapiro-Francia, Kolmogorov-Smirnov, Lillie, Pearson, Andersen, Cramer, Jarque_Bera con corrección)

lillie.test(mod1$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  mod1$residuals
## D = 0.098494, p-value = 0.6455
sf.test(mod1$residuals)
## 
##  Shapiro-Francia normality test
## 
## data:  mod1$residuals
## W = 0.96013, p-value = 0.2662
JBGTest(mod1$residuals)
## 
## Jarque-Bera (Glinskiy)
## 
## Variant: JB(Classic)
## 
## Statistic: 1.1664 , df = 2 , p-value = 0.5581
shapiro.test(mod1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.97286, p-value = 0.6202

De acuerdo al test, el p-valor 0.0708, es una valor mayor que 0.05, luego los residuos presentan comportamiento Normal.

S2 Homocedasticidad

La hipótesis en cuestión es \(H_{0S2}: Homocedasticidad\)

respecto de este supuesto, se debe tener claro número de replicas por grupo o tamaño de cada grupo (aparece en S0). todos presentan mismo tamaño y como el test de normalidad (aunque la significancia pueda ser similar debe aplicarse el adecuado, entre Levene, Fligner-Killen, Brown-Forsythe).

Escoger Levene o Barttlet depende de si los grupos tienen o no la misma cantidad de datos, cuando tienen el mismo n por grupo es mejor Bartllet (además de cumplir con el supuesto de normalidad), sino Levene

summary(Levene.Test(VD~grupo,data=datos),center="mean") # usar center="mean", cuando se cumple normalidad, sino usar "median"
## 
## --- Homoscedasticity Test Summary ---
## 
## Method applied         : Levene (median) - global by cells 
## F Statistic            : 3.510619 
## Degrees of freedom     : 2 (between), 27 (within)
## p-value                : 0.04413513 * 
## Decision (alpha = 0.05): Heteroscedastic (cells) 
## ----------------------------------------
summary(BartlettTest(VD~grupo,data=datos)) # usar cuando existe normalidad
## 
## --- Homoscedasticity Test Summary ---
## 
## Method applied         : Bartlett 
## Chi-squared Statistic  : 7.8705 
## Degrees of freedom     : 2 
## p-value                : 0.0195 * 
## Decision (alpha = 0.05): Heterocedastic 
## ----------------------------------------

esta versión de Test, nos presenta que el estadístico es significativo, por tanto se rechaza \(H_{0S2}: Homocedasticidad\)

Analisis sobre Hi

de acuerdo al Supuesto 1, se debe realizar un procedimiento Paramétrico, en este caso ANOVA, el cual fue provisto anteriormente.

mod1<-aov(VD~as.factor(grupo),data = datos) # es funcional con Homocedasticidad o heterocedasticidad, pero con Normalidad
summary(mod1)
##                  Df Sum Sq Mean Sq F value Pr(>F)    
## as.factor(grupo)  2  64132   32066   17085 <2e-16 ***
## Residuals        27     51       2                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El p-valor resulta menor que \(\alpha\), por tanto se rechaza la hipótesis que plantea que los grupos poseen iguales características. Por ende, se debe buscar entre que grupos se provoca la diferencia

en caso de no normalidad de los residuos aplicar Kruskal wallis

kruskal.test(VD~as.factor(grupo),data = datos)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  VD by as.factor(grupo)
## Kruskal-Wallis chi-squared = 25.806, df = 2, p-value = 2.49e-06

PostHoc

Tests PostHoc según comportamiento de los Supuestos (paquete Analitica)
Normales
No normales
Homocedásticos Heterocedástico Homocedásticos Heterocedástico
Tukey - TukeyTest() T2 - T2Test() MannWhitney - MWTest() Brunner Munzel - BMTest()
Scheffe - ScheffeTest() T3 - T3Test() Dunn - DunnTest() Bonferroni NP - BonferroniNPTest()
Holm - HolmTest() Games Howell - GHTest() Conover/Iman - ConoverTest() Jonckheere-Terpstra - JT_Test
Bonferroni - BonferroniTest() DSCF - DSCFTest()
Duncan - DuncanTest() Nemenyi - NemenyiTest()
Gabriel - GabrielTest()
SNK - SNKTest()
LSD - LSDTest()
Sidak - SidakTest()

El test adecuado debe analizarlo, hay algunos más conservadores, y otros que dependen de la cantidad de comparaciones, además deben seguir el cumplimiento de los supuestos. Dado que hay normalidad y heterocedásticidad. Aplicare el test T2 (alternativas T3, Games-Howell)

CASO S1 Aprobado, S2 Rechazado

res<-T2Test(mod1)
summary(res)
## =====================================
##   Multiple Comparison Method Summary
## =====================================
## Method used: Tamhane T2 
## 
## >> Group means:
##       1       2       3 
## 185.259 260.937 150.129 
## 
## >> Order of means (from highest to lowest):
## [1] "2" "1" "3"
## 
## >> Pairwise comparisons:
##    Comparacion Diferencia  t_value    gl p_value Significancia
## 1        1 - 2     75.678 186.8995 17.73       0           ***
## 11       1 - 3     35.130  51.1930 12.14       0           ***
## 2        2 - 3    110.808 158.0970 12.95       0           ***

Como se puede apreciar, la comparación entre cada grupo presenta diferencias significativas (p-valor menor que \(\alpha\)). Así mismo se destaca en la gráfica que muestra los promedios, que cada grupo posee diferentes letras indicando que no hay grupos que se asemejen.

plot(res)

CASO S1 Aprobado, S2 Aprobado

en el caso de aplicar Paramétrico Homocedástico Supuesto de Normalidad aprobado (ejemplo)

res<-ScheffeTest(mod1)
summary(res)
## =====================================
##   Multiple Comparison Method Summary
## =====================================
## Method used: Scheffe 
## 
## >> Group means:
##       1       2       3 
## 185.259 260.937 150.129 
## 
## >> Order of means (from highest to lowest):
## [1] "2" "1" "3"
## 
## >> Pairwise comparisons:
##   Comparacion Diferencia      SE2     F_obs Valor_Critico p_value Significancia
## 1       1 - 2     75.678 0.375367  7628.740        1.5868       0           ***
## 2       1 - 3     35.130 0.375367  1643.879        1.5868       0           ***
## 3       2 - 3    110.808 0.375367 16355.196        1.5868       0           ***
plot(res)

CASO S1 RECHAZADO

en este caso debe usar Kruskal Wallis,

modk<- kruskal.test(VD~grupo,data=datos)
modk
## 
##  Kruskal-Wallis rank sum test
## 
## data:  VD by grupo
## Kruskal-Wallis chi-squared = 25.806, df = 2, p-value = 2.49e-06

Como p valor menor que alfa, esto implica que existen grupos con diferencias, busquemos donde se encuentra

Caso Homocedastico

res<-DunnTest(VD~as.factor(grupo),data=datos)
summary(res)
## =====================================
##   Multiple Comparison Method Summary
## =====================================
## Method used: Dunn (no parametrico) 
## 
## >> Group means:
##    1    2    3 
## 15.5 25.5  5.5 
## 
## >> Order of means (from highest to lowest):
## [1] "2" "1" "3"
## 
## >> Pairwise comparisons:
##   Comparacion    z p_value p_ajustada Significancia
## 1       1 - 2 2.54  0.0111     0.0222             *
## 2       1 - 3 2.54  0.0111     0.0222             *
## 3       2 - 3 5.08  0.0000     0.0000           ***
plot(res)

Caso Heterocedastico

ud debe escoger uno de los que aparece en la tabla

res<-BonferroniNPTest(VD~as.factor(grupo),data=datos)
summary(res)
## =====================================
##   Multiple Comparison Method Summary
## =====================================
## Method used: Bonferroni (non-parametric) 
## 
## >> Group means:
##    1    2    3 
## 15.5 25.5  5.5 
## 
## >> Order of means (from highest to lowest):
## [1] "2" "1" "3"
## 
## >> Pairwise comparisons:
##   Comparacion   W p_value p_ajustada Significancia
## 1       1 - 2   0   2e-04      5e-04           ***
## 2       1 - 3 100   2e-04      5e-04           ***
## 3       2 - 3 100   2e-04      5e-04           ***
plot(res)