Alternativas a la ANOVA de una vía para heteroskedastic de datos

Question

Tengo los datos de los 3 grupos de la biomasa de algas ($A$, $B$, $C$) que contienen desigual tamaño de la muestra ($n_A=15$, $n_B=13$, $n_C=12$) y me gustaría comparar si estos grupos son de la misma población.

ANOVA de una vía sin duda sería el camino a seguir, sin embargo al llevar a cabo pruebas de normalidad en mis datos, heteroskedascity parece a la cuestión principal. Mis datos en bruto, sin ningún tipo de transformación, se elaboró una relación de varianzas ($F_{\max} = 19.1$) que es mucho mayor que el valor crítico ($F_{\rm crit} = 4.16$) y por lo tanto no puedo realizar el ANOVA de una vía.

También traté de transformación para normalizar mis datos. Incluso después de los ensayos de las distintas transformaciones (registro, raíz cuadrada, cuadrado), el menor $F_{\max}$ producido después de la transformación con un $\log_{10}$ transformación se $7.16$, lo que era aún más alta en comparación con $F_{\rm crit}$.

Alguien aquí puede aconsejarme sobre dónde ir desde aquí? No puedo pensar en otros métodos de transformación para normalizar los datos. ¿Existen alternativas a un ANOVA de una vía?

P. S.: mis datos en bruto son las siguientes:

A: 0.178 0.195 0.225 0.294 0.315 0.341 0.36  0.363 0.371 0.398 0.407 0.409 0.432 
   0.494 0.719
B: 0.11  0.111 0.204 0.416 0.417 0.441 0.492 0.965 1.113 1.19  1.233 1.505 1.897
C: 0.106 0.114 0.143 0.435 0.448 0.51  0.576 0.588 0.608 0.64  0.658 0.788 0.958

Answer 1

Hay un número de opciones disponibles cuando se trata con heteroscedastic de datos. Lamentablemente, ninguno de ellos está garantizado para siempre. Aquí son algunas de las opciones que estoy familiarizado con:

• transformaciones
• Welch ANOVA
• mínimos cuadrados ponderados
• regresión robusta
• heterocedasticidad errores estándar consistentes
• bootstrap
• Test de Kruskal-Wallis
• la regresión logística ordinal

---

Actualización: Aquí está una demostración en R de algunas de las formas de ajuste de un modelo lineal (es decir, un ANOVA o regresión) cuando haya heterocedasticidad / heterogeneidad de la varianza.

Vamos a empezar por echar un vistazo a sus datos. Para mayor comodidad, los tengo cargados en dos marcos de datos llamado my.data (que está estructurado como el anterior y con una columna de cada grupo) y stacked.data (que tiene dos columnas: values número y ind con el grupo indicador).

Podemos formalmente prueba de heterocedasticidad con el test de Levene prueba:

library(car)
leveneTest(values~ind, stacked.data)
# Levene's Test for Homogeneity of Variance (center = median)
#       Df F value   Pr(>F)   
# group  2  8.1269 0.001153 **
#       38                    

Por supuesto, usted tiene heterocedasticidad. Vamos a comprobar para ver lo que las varianzas de los grupos. Una regla de oro es que los modelos lineales son bastante robustos a la heterogeneidad de la varianza tan larga como la máxima varianza no es más que $4\!\times$ mayor que el mínimo de la varianza, por lo que nos vamos a encontrar que la relación así:

apply(my.data, 2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01734578 0.33182844 0.06673060 
var(my.data$B, na.rm=T) / var(my.data$A, na.rm=T)
# [1] 19.13021

Sus variaciones difieren sustancialmente, con la más grande, B, siendo $19\!\times$ el más pequeño, A. Esta es una problemática a nivel de heteroscedsaticity.

Había pensado usar transformaciones tales como el registro o la raíz cuadrada para estabilizar la varianza. Que va a trabajar, en algunos casos, pero de Box-Cox transformaciones de tipo estabilizar la varianza apretando los datos de forma asimétrica, ya sea apretando hacia abajo con los datos más alta de exprimido más, o apretando hacia arriba, con los datos más exprimido al máximo. Por lo tanto, usted necesita la varianza de los datos a cambiar con la media para que esto funcione de manera óptima. Los datos tienen una enorme diferencia en la varianza, pero una diferencia relativamente pequeña entre las medias y medianas, es decir, las distribuciones en su mayoría se superponen. Como un ejercicio de aprendizaje, podemos crear un cierto parallel.universe.data mediante la adición de $2.7$ todos B valores y $.7$ a C's para mostrar cómo funcionaría:

parallel.universe.data = with(my.data, data.frame(A=A, B=B+2.7, C=C+.7))
apply(parallel.universe.data,       2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01734578 0.33182844 0.06673060 
apply(log(parallel.universe.data),  2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.12750634 0.02631383 0.05240742 
apply(sqrt(parallel.universe.data), 2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01120956 0.02325107 0.01461479 
var(sqrt(parallel.universe.data$B), na.rm=T) / 
var(sqrt(parallel.universe.data$A), na.rm=T)
# [1] 2.074217

El uso de la transformación de raíz cuadrada estabiliza los datos muy bien. Usted puede ver la mejora en el universo paralelo de los datos aquí:

En lugar de sólo tratar las diferentes transformaciones, un enfoque más sistemático para optimizar el Box-Cox parámetro $\lambda$ (a pesar de que generalmente se recomienda a la vuelta que a la más cercana interpretables de transformación). En el caso de la raíz cuadrada, $\lambda = .5$, o el registro, $\lambda = 0$, son aceptables, a pesar de que ninguno de los dos funciona. Para el universo paralelo de los datos, la raíz cuadrada es mejor:

boxcox(values~ind, data=stacked.data,    na.action=na.omit)
boxcox(values~ind, data=stacked.pu.data, na.action=na.omit) 

Dado que este caso es un análisis de la VARIANZA (es decir, ninguna de las variables continuas), una manera de lidiar con la heterogeneidad es el uso de la corrección de Welch para el denominador grados de libertad en el $F$-prueba (n.b., df = 19.445, un valor fraccionario, en lugar de df = 38):

oneway.test(values~ind, data=stacked.data, na.action=na.omit, var.equal=FALSE)
#  One-way analysis of means (not assuming equal variances)
# 
# data:  values and ind
# F = 4.1769, num df = 2.000, denom df = 19.445, p-value = 0.03097

Un enfoque más general es el uso de mínimos cuadrados ponderados. Desde algunos grupos (B) se esparcen más, los datos de los grupos de proporcionar menos información acerca de la ubicación de la media de los datos en otros grupos. Podemos dejar el modelo de incorporar esta proporcionando un peso con cada punto de datos. Un sistema común es utilizar el recíproco del grupo de la varianza como el peso:

wl = 1 / apply(my.data, 2, function(x){ var(x, na.rm=T) })
stacked.data$w = with(stacked.data, ifelse(ind=="A", wl[1], 
                                           ifelse(ind=="B", wl[2], wl[3])))
w.mod = lm(values~ind, stacked.data, na.action=na.omit, weights=w)
anova(w.mod)
# Response: values
#           Df Sum Sq Mean Sq F value  Pr(>F)  
# ind        2   8.64  4.3201  4.3201 0.02039 *
# Residuals 38  38.00  1.0000                  

Este rendimientos ligeramente diferente $F$ $p$- valores de la ponderado análisis de la VARIANZA (4.5089, 0.01749), pero se ha abordado la heterogeneidad de bien:

Mínimos cuadrados ponderados no es una panacea, sin embargo. Una incómoda realidad es que solo es justo si los pesos están apenas a la derecha, lo cual significa, entre otras cosas, que son conocidos a priori. No abordar la no-normalidad (tales como la distorsión) o valores atípicos. El uso de los pesos estimados a partir de los datos se suelen funcionar bien, aunque, especialmente si usted tiene la cantidad suficiente de datos para la estimación de la varianza con razonable precisión (esto es análogo a la idea de utilizar un $z$-tabla en lugar de una $t$-tabla cuando ha $50$ o $100$ grados de libertad), los datos son lo suficientemente normal, y no parece tener valores atípicos. Desafortunadamente, usted tiene relativamente pocos datos (13 o 15 por grupo), algunos de sesgo y, posiblemente, algunos valores atípicos. No estoy seguro de que estos son lo suficientemente mal como para hacer una gran cosa, pero usted puede mezclar mínimos cuadrados ponderados con robusta métodos. En lugar de utilizar la varianza como su medida de propagación (que es sensible a los valores atípicos, especialmente con baja $N$), se podría utilizar el recíproco de la inter-cuartil rango (que no se ve afectado por hasta el 50% de los valores atípicos en cada grupo). Estos pesos pueden entonces ser combinado con regresión robusta, utilizando una función de pérdida como de Tukey bisquare:

1 / apply(my.data,       2, function(x){ var(x, na.rm=T) })
#         A         B         C 
# 57.650907  3.013606 14.985628 
1 / apply(my.data,       2, function(x){ IQR(x, na.rm=T) })
#        A        B        C 
# 9.661836 1.291990 4.878049 

rw = 1 / apply(my.data, 2, function(x){ IQR(x, na.rm=T) })
stacked.data$rw = with(stacked.data, ifelse(ind=="A", rw[1], 
                                            ifelse(ind=="B", rw[2], rw[3])))
library(robustbase)
w.r.mod = lmrob(values~ind, stacked.data, na.action=na.omit, weights=rw)
anova(w.r.mod, lmrob(values~1,stacked.data,na.action=na.omit,weights=rw), test="Wald")
# Robust Wald Test Table
# 
# Model 1: values ~ ind
# Model 2: values ~ 1
# Largest model fitted by lmrob(), i.e. SM
# 
#   pseudoDf Test.Stat Df Pr(>chisq)  
# 1       38                          
# 2       40    6.6016  2    0.03685 *

Los pesos aquí no son tan extremas. La predicción de grupo significa que difieren ligeramente (A: WLS 0.36673, robusto 0.35722; B: WLS 0.77646, robusto 0.70433; C: WLS 0.50554, robusto 0.51845), con los medios de B y C menos tirado por los valores extremos.

En econometría la de Huber-White ("sandwich") error estándar es muy popular. Como la corrección de Welch, esto no requiere conocer las desviaciones a-priori y no requiere la estimación de los pesos de sus datos y/o contingente en un modelo que puede no ser la correcta. Por otro lado, no sé cómo incorporar este con un análisis de la VARIANZA, lo que significa que usted consigue solamente para las pruebas de individuales ficticio códigos, que me parece menos útil en este caso, pero voy a demostrar que ellos de todos modos:

library(sandwich)
mod = lm(values~ind, stacked.data, na.action=na.omit)
sqrt(diag(vcovHC(mod)))
# (Intercept)        indB        indC 
#  0.03519921  0.16997457  0.08246131 
2*(1-pt(coef(mod) / sqrt(diag(vcovHC(mod))), df=38))
#  (Intercept)         indB         indC 
# 1.078249e-12 2.087484e-02 1.005212e-01 

La función vcovHC calcula un heteroscedasticicy constante de varianza-covarianza de la matriz para su betas (su maniquí de códigos), que es lo que las cartas en la llamada a la función. Para obtener los errores estándar, extrae la diagonal principal y tomar la raíz cuadrada. Para obtener $t$-pruebas para sus betas, dividir los coeficientes estimados por el SEs y comparar los resultados con las $t$-distribución (es decir, el $t$-distribución con tu residual grados de libertad).

Para R usuarios en concreto, @TomWenseleers notas en los comentarios de abajo que el ?Anova en función de la car paquete puede aceptar un white.adjust argumento a obtener un $p$-valor para el factor de uso de heterocedasticidad consistente errores.

Anova(mod, white.adjust=TRUE)
# Analysis of Deviance Table (Type II tests)
# 
# Response: values
#           Df      F  Pr(>F)  
# ind        2 3.9946 0.02663 *
# Residuals 38                 
# ---
# Signif. codes:  0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1

Usted puede tratar de conseguir una estimación empírica de lo que la actual distribución de muestreo de la estadística de prueba se ve como por bootstrapping. En primer lugar, crear un verdadero null haciendo que todo el grupo significa exactamente igual. Luego de remuestreo con reemplazo y calcular el estadístico de prueba ($F$) en cada bootsample para obtener una estimación empírica de la distribución de muestreo de $F$ bajo la nula con sus datos, independientemente de su estatus con respecto a la normalidad o la homogeneidad. La proporción de la distribución de muestreo que es tan extremo o más extremo que el de su observó estadístico de prueba es el $p$-valor:

mod    = lm(values~ind, stacked.data, na.action=na.omit)
F.stat = anova(mod)[1,4]
  # create null version of the data
nullA  = my.data$A - mean(my.data$A)
nullB  = my.data$B - mean(my.data$B, na.rm=T)
nullC  = my.data$C - mean(my.data$C, na.rm=T)
set.seed(1)
F.vect = vector(length=10000)
for(i in 1:10000){
  A = sample(na.omit(nullA), 15, replace=T)
  B = sample(na.omit(nullB), 13, replace=T)
  C = sample(na.omit(nullC), 13, replace=T)
  boot.dat  = stack(list(A=A, B=B, C=C))
  boot.mod  = lm(values~ind, boot.dat)
  F.vect[i] = anova(boot.mod)[1,4]
}
1-mean(F.stat>F.vect)
# [1] 0.0485

En algunas formas, el arranque es la máxima reducción de asunción enfoque para llevar a cabo un análisis de los parámetros (por ejemplo, los medios), pero asume que los datos son una buena representación de la población, lo que significa que tiene una razonable del tamaño de la muestra. Desde su $n$'s son pequeños, puede ser menos fiables. Probablemente la máxima protección contra los no-normalidad y la heterogeneidad es el uso de una prueba no paramétrica. El basic no paramétrico de la versión de un ANOVA es la prueba de Kruskal-Wallis:

kruskal.test(values~ind, stacked.data, na.action=na.omit)
#  Kruskal-Wallis rank sum test
# 
# data:  values by ind
# Kruskal-Wallis chi-squared = 5.7705, df = 2, p-value = 0.05584

Aunque el test de Kruskal-Wallis es sin duda la mejor protección contra los errores de tipo I, que sólo puede ser usada con una sola variable categórica (es decir, no continua predictores o diseños factoriales) y tiene menos poder de todas las estrategias analizadas. Otro enfoque no paramétrico es el uso de regresión logística ordinal. Esto parece extraño que un montón de gente, pero sólo es necesario asumir que su respuesta de datos contienen legítimo ordinal de la información, que por supuesto que lo hacen, o bien todos los otros de la estrategia anterior no es válido así:

library(rms)
olr.mod = orm(values~ind, stacked.data)
olr.mod
# Model Likelihood          Discrimination          Rank Discrim.    
# Ratio Test                 Indexes                Indexes       
# Obs            41    LR chi2      6.63    R2                  0.149    rho     0.365
# Unique Y       41    d.f.            2    g                   0.829
# Median Y    0.432    Pr(> chi2) 0.0363    gr                  2.292
# max |deriv| 2e-04    Score chi2   6.48    |Pr(Y>=median)-0.5| 0.179
#                      Pr(> chi2) 0.0391

Puede que no sea claro desde la salida, pero la prueba de que el modelo como un todo, que en este caso es la prueba de sus grupos, es la chi2 bajo Discrimination Indexes. Dos versiones están en la lista, una prueba de razón de verosimilitud y una puntuación de la prueba. La prueba de razón de verosimilitud, se considera generalmente la mejor. Se obtiene un $p$-valor de 0.0363.