Anova una via

El nombre ANOVA significa Análisis de varianza y se usa porque la idea original es dividir la varianza general en la respuesta a la debida a cada uno de los factores y el error. Los predictores se denominan típicamente factores que tienen una cierta cantidad de niveles y los parámetros se denominan a menudo efectos y estos pueden ser fijos pero desconocidos, llamados modelos de efectos fijos,o aleatorios cuando los parámetros se toman como variables aleatorias.

Modelo de Anova una via

Suponga un factor \(\alpha\) con \(k\) niveles \(i=1,2,3,...,k\) que ocuure \(j=1,2,...,J_{i}\) observaciones por nivel. El modelo que representa el efecto del factor sobre la variable de interes \(Y\) esta dado por:
\[y_{ij}=\mu+\alpha_{i}+\varepsilon_{ij}\]

Ejemplo

El conjunto de datos “coagulation” que estan en R consta del registro 24 tiempos de coagulación sanguínea (n=24, Y=coagulación sanguínea). Se asignaron aleatoriamente 24 animales a cuatro dietas diferentes y se tomaron muestras en un orden aleatorio (diseño aleatorio). Estos datos provienen de Box, Hunter y cazador (1978)
Lo primero es observar los datos y como es su estructura que se encuentran disponibles en la libreria “faraway” y la librería “MASS” para funciones especiales

library(MASS)
library(faraway)
data(coagulation)
head(coagulation)
##   coag diet
## 1   62    A
## 2   60    A
## 3   63    A
## 4   59    A
## 5   63    B
## 6   67    B

Como observamos tenemos dos columnas una con la variable de interes y la otra (factor) indica los niveles del factor dieta que son A,B,C y D.
Como primer paso es visualizar la distribución de la coagulación en los diferentes grupos de trtamiento. Tenga en cuenta que si la columna de tratamientos no esta como factor es necesario transformarla

plot(coag ~ diet, data=coagulation,
     main="Distribución de tiempos de coagulación sanguínea")

Como podemos ver en el grafico boxplot:
1. Outliers se presentan en la dieta C
2. Sesgo: esto será evidente a partir de una forma asimétrica de las cajas. El cual esta evidente en la dieta C y leve en la dieta D
3. Variación desigual: esto será evidente a partir de tamaños de caja claramente desiguales, que para este caso es evidente entre las dietas B y C.
En general, se recomienda realizar las pruebas para detectar los outlier y para el sesgo, junto con una descripción estadistica de cada grupo (diet)

dietaA<- subset.data.frame(coagulation,diet=="A")
dietaB<- subset.data.frame(coagulation,diet=="B")
dietaC<- subset.data.frame(coagulation,diet=="C")
dietaD<- subset.data.frame(coagulation,diet=="D")
summary(dietaA$coag)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   59.00   59.75   61.00   61.00   62.25   63.00
summary(dietaB$coag)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   63.00   64.25   65.50   66.00   66.75   71.00
summary(dietaC$coag)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   66.00   67.25   68.00   68.00   68.00   71.00
summary(dietaD$coag)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   56.00   59.75   61.50   61.00   63.00   64.00

Como se observan en los datos estadísticos de la muestra hay una diferencia en el valor medio de coagulación en las dietas C(68) y D(61)lo cual debemos verificar con una prueba estadística de comparación

Pruebas para detectar Outlier

Los valores atípicos en un conjunto de datos pueden provenir de las siguientes fuentes posibles:
1. muestras de datos contaminados;
2. puntos de datos de diferentes poblaciones con métodos de muestreo incorrectos,
3. respuesta significativa subyacente al tratamiento (p. ej., variación biológica de las muestras),
4. error en la recopilación o el análisis de datos.
Exiten diferentes metodologías y pruebas para detectar Outliers y varios enfoques para detectar valores atípicos en R, desde técnicas simples como estadísticas descriptivas (rango intercuartil, histograma, diagrama de caja y percentiles) hasta técnicas más formales como el filtro Hampel, el Grubbs, el Dixon y el Pruebas de Rosner para valores atípicos.
Recuerde que la eliminación de los datos atipicos, generalmente requiere una reflexión y analisis cuidadoso por parte del investigador dado que algunos dominios, es común eliminar los valores atípicos, y esto a menudo ocurre debido a un mal funcionamiento del proceso. En otros campos, los valores atípicos se mantienen porque contienen información valiosa.

Criterio de rango intercuartílico (IQR)

El criterio IQR significa que todas las observaciones que esten fuera del IQR son considerados como posibles valores atípicos por R. En otras palabras, todas las observaciones fuera del intervalo:
\[q_{0.75}-1.5IQR;q_{0.75}+1.5IQR\] También es posible determinar cuales son los valores de los posibles valores atípicos en función del criterio IQR gracias a la función boxplot.stats()$out:

boxplot.stats(dietaA$coag)$out
## numeric(0)
boxplot.stats(dietaB$coag)$out
## numeric(0)
boxplot.stats(dietaC$coag)$out
## [1] 71
boxplot.stats(dietaD$coag)$out
## numeric(0)

Aqui, nos indica que solo hay un valor atipico y este esta en la dieta C, con un valor de coagulación de 71.Gracias a la función which() es posible extraer el número de fila correspondiente a estos valores atípicos y con esta información, ahora puede volver fácilmente a las filas específicas del conjunto de datos para verificarlas

Test de Dixon

La prueba Q de Dixon es una prueba basada en hipótesis que se utiliza para identificar un solo valor atípico (valor mínimo o máximo) en un conjunto de datos univariado.Por lo tanto, si se sospecha que hay más de un valor atípico, la prueba debe realizarse individualmente en estos valores atípicos sospechosos.Esta prueba es aplicable a un conjunto de datos de muestra pequeño (el tamaño de la muestra está entre 3 y 30) es decir para muestras pequeñas y cuando los datos se distribuyen normalmente. Aunque la prueba Q de Dixon asume la normalidad, es resistente a la desviación de la normalidad, en nuestro caso el conjunto de datos tiene 24 observaciones (utilice str()para verificar).
Entonces, \(H_{0}\) Maximo valor de coagulación en la dieta C no es un outlier

library(outliers)
prueba_outlier<- dixon.test(coagulation$coag,
                       opposite = F)
#use F para probar el mas alto
prueba_outlier
## 
##  Dixon test for outliers
## 
## data:  coagulation$coag
## Q = 0.25, p-value = 0.6961
## alternative hypothesis: lowest value 56 is an outlier
prueba_outlier1<- dixon.test(coagulation$coag,
                       opposite = TRUE)
#use F para probar el mas alto y True para el más bajo
prueba_outlier1
## 
##  Dixon test for outliers
## 
## data:  coagulation$coag
## Q = 0.25, p-value = 0.6961
## alternative hypothesis: highest value 71 is an outlier

Como vemos en ambas caso ni el dato más bajo (56) y el más alto (71) corresponden a outlier dado que \(p_{valor}>\alpha=0.05\)

Como podemos ver en el conjunto de datos el valor más bajo (56) es el dato en la posición siete (7) de la dieta D y el dato más alto en coagulación (7) esta en la fila (7) dietaB y (13) dieta C.

Prueba de Rosner

Una prueba más potente es el test de Rosner para valores atípicos y tiene las ventajas de que:

  1. Se utiliza para detectar varios valores atípicos a la vez (a diferencia de la prueba de Grubbs y Dixon, que debe realizarse de forma iterativa para detectar múltiples valores atípicos), y
  2. Está diseñado para evitar el problema del enmascaramiento, donde un valor atípico que tiene un valor cercano a otro puede pasar desapercibido.

A diferencia de la prueba de Dixon y de Grubbs,la prueba de Rosner es más apropiada cuando el tamaño de la muestra es grande.

Asimetria

En las distribuciones simétricas los parámetros media, mediana y moda coinciden. Para este ejemplo, usaremos la libreria moments en R. Existen diferentes coefcientes de asimetria como:
-Fisher
-Pearson

Coeficiente de asimetría de Pearson: Mide la desviación de la simetría, expresada la diferencia entre la media y la mediana con respecto a la desviación estándar del grupo de mediciones la fórmula es: Sólo se puede utilizar en distribuciones uniformes, unimodales y moderadamente asimétricas.

library(moments) 
 library(modeest)
## 
## Attaching package: 'modeest'
## The following object is masked from 'package:moments':
## 
##     skewness
print(skewness(coagulation$coag)) 
## Warning: encountered a tie, and the difference between minimal and 
##                    maximal value is > length('x') * 'tie.limit'
## the distribution could be multimodal
## [1] 0.008797178
mlv(coagulation$coag, method = "mfv")
## [1] 63
hist(coagulation$coag)

Como vemos el coeficiente de asimetria es positivo y cercano a cero (0.008797178).La asimetría positiva, significa que la cola de la distribución de coagulación se alarga (a la derecha) para valores superiores a la media de coagulación y la moda es 63 de coagulación pero segun el resultado, esta indicando que es multimodal, lo que debe revisarse.

Estimación de Parámetros del modelo

Los efectos se pueden estimar utilizando el método de mínimos cuadrados y la primera prueba de interés es si hay una diferencia en los niveles del factor, para esto realizamos la siguiente prueba:
\[H_{0}:\alpha_{i}=0\] Y se utiliza la misma prueba F que en regresión lineal. El resultado de esta prueba será el mismo sin importar qué codificación/restricción usemos. Si No se rechaza la hipotesis nula, entonces hemos terminado el proceso de ANOVA no obstante, esta sujeto a una investigación de transformación y valores atípicos. Si rechazamos la hipotesis nula, debemos investigar qué niveles difieren (es decir, analisis POST- HOC)

Finalmente, debemos verificar que se dan las condiciones para realizar una prueba parametrica para comparar los niveles de coagulación a los diferentes niveles del factor dieta:
-Normalidad: cada muestra se toma de una población normalmente distribuida.
-Independencia de la muestra: cada muestra se ha extraído independientemente de las otras muestras. -Igualdad de varianza: que la varianza de los datos en los diferentes grupos debe ser la misma.

Teniendo en cuenta lo anterior (muestra de n=24 y de \(n_{i}=6\) para cada grupo de dieta) se aplicará una prueba No parametrica para las comparaciones.

Comparaciones con Pruebas No Parámetricas

Los métodos no paramétricos son útiles cuando no se cumple el supuesto de normalidad y el tamaño de la muestra es pequeño. Sin embargo, las pruebas no paramétricas no están completamente libres de supuestos acerca de los datos. Por ejemplo, es fundamental presuponer que las observaciones de las muestras son independientes y provienen de la misma distribución.
La prueba que se aplicará es el test “Prueba de Kruskal-Wallis” que es el equivalente de la ANOVA de un solo factor para datos no pareados.
Sin embargo, las pruebas no paramétricas no están completamente libres de supuestos acerca de los datos. Por ejemplo, es fundamental presuponer que las observaciones de las muestras son independientes y provienen de la misma distribución. Además, en los diseños de dos muestras, se requiere el supuesto de igualdad de forma y dispersión (Homocedasticidad)

Prueba Kruskal Wallis

El test de Kruskal-Wallis, se trata de una extensión del test de Mann-Whitney para más de dos grupos. Se trata por lo tanto de un test que emplea rangos para contrastar la hipótesis de que \(k\) muestras han sido obtenidas de una misma población.Bajo ciertas simplificaciones puede considerarse que el test de Kruskal-Wallis compara las medianas.

library(HH)
## Loading required package: lattice
## 
## Attaching package: 'lattice'
## The following object is masked from 'package:faraway':
## 
##     melanoma
## Loading required package: grid
## Loading required package: latticeExtra
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## 
## Attaching package: 'survival'
## The following objects are masked from 'package:faraway':
## 
##     rats, solder
## Loading required package: TH.data
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## Loading required package: gridExtra
## 
## Attaching package: 'HH'
## The following objects are masked from 'package:faraway':
## 
##     logit, vif
hov(coagulation$coag ~ coagulation$diet)
## 
##  hov: Brown-Forsyth
## 
## data:  coagulation$coag
## F = 0.64919, df:coagulation$diet = 3, df:Residuals = 20, p-value =
## 0.5926
## alternative hypothesis: variances are not identical

Como se observa el pvalor>nivel de significancia y entonces no se rechaza la hipotesis nula:

$H_{0}: {1}^{2}={2}{2}={3}^{2}={4}{2} $, luego procedemos a realizar la prueba de comparación

kruskal.test(coag ~ diet, data = coagulation)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  coag by diet
## Kruskal-Wallis chi-squared = 17.015, df = 3, p-value = 0.0007016

El test encuentra significancia en la diferencia de al menos dos grupos, porque p-value = 0.0007016. Para identificar donde se encuentran estas dierencias tenemos que proceder con el equivalente no parametrico de las t-Student pareadas:

pairwise.wilcox.test(coagulation$coag, 
                coagulation$diet, 
                p.adjust.method = "bonf",
                paired = F)
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  coagulation$coag and coagulation$diet 
## 
##   A     B     C    
## B 0.112 -     -    
## C 0.078 0.630 -    
## D 1.000 0.038 0.014
## 
## P value adjustment method: bonferroni

Como se puede observar la diferencia significativa en los niveles de coagulación presentan diferencias significativas entre la dieta C con la dieta D y entre la dieta D y B (\(p<0.05\)).

library(ggpubr)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:latticeExtra':
## 
##     layer
ggplot(coagulation, aes(x = diet, y = coag, fill = diet))+
  geom_boxplot()+
  geom_point()+
  ylim(0,71)+
  stat_compare_means(method = "kruskal",label.y =  71)+
  stat_compare_means(method = "wilcox.test", ref.group = "Control", label = "p.signif", label.y = 71)+
  theme_classic()
## Warning: Computation failed in `stat_compare_means()`:
## valor ausente donde TRUE/FALSE es necesario

Referencias

Herrera,Eddy.(2021). Notas de clase: Bioestadistica Aplicada con R.