Alternativas al ANOVA unidireccional para datos heteroscedasticos

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

$F_{\max} = 19.1$$F_{\rm crit} = 4.16$

También intenté la transformación para normalizar mis datos. Incluso después de las pruebas de varias transformaciones (log, raíz cuadrada, cuadrada), la más baja producida después de la transformación con una transformación fue , que fue aún mayor en comparación con .$F_{\max}$$\log_{10}$$7.16$$F_{\rm crit}$

¿Alguien puede aconsejarme dónde ir desde aquí? No puedo pensar en otros métodos de transformación para normalizar por datos. ¿Hay alguna alternativa a un ANOVA unidireccional?

PD: mis datos en bruto están a continuación:

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

Hay una serie de opciones disponibles cuando se trata con datos heteroscedasticos. Desafortunadamente, ninguno de ellos está garantizado para funcionar siempre. Aquí hay algunas opciones con las que estoy familiarizado:

• transformaciones
• Welch ANOVA
• mínimos cuadrados ponderados
• regresión robusta
• heteroscedasticidad errores estándar consistentes
• oreja
• Prueba de Kruskal-Wallis
• regresión logística ordinal

Actualización: Aquí hay una demostración R de algunas formas de ajustar un modelo lineal (es decir, un ANOVA o una regresión) cuando tiene heterocedasticidad / heterogeneidad de varianza.

Comencemos por echar un vistazo a sus datos. Por conveniencia, los he cargado en dos marcos de datos llamados my.data(que está estructurado como arriba con una columna por grupo) y stacked.data(que tiene dos columnas: valuescon los números y indcon el indicador de grupo).

Podemos probar formalmente la heterocedasticidad con la prueba de Levene:

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                    

Efectivamente, tienes heterocedasticidad. Verificaremos cuáles son las variaciones de los grupos. Una regla general es que los modelos lineales son bastante robustos a la heterogeneidad de la varianza, siempre que la varianza máxima no sea más de mayor que la varianza mínima, por lo que también encontraremos esa relación: $4\!\times$

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 el mayor, Bsiendo los más pequeños,. Este es un nivel problemático de heteroscedsaticidad. $19\!\times$A

parallel.universe.data$2.7$B$.7$Cpara 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 esos datos bastante bien. Puede ver la mejora de los datos del universo paralelo aquí:

$\lambda$$\lambda = .5$$\lambda = 0$

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

$F$df = 19.445df = 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 utilizar mínimos cuadrados ponderados . Dado que algunos grupos ( B) se extienden más, los datos en esos grupos proporcionan menos información sobre la ubicación de la media que los datos en otros grupos. Podemos dejar que el modelo incorpore esto proporcionando un peso con cada punto de datos. Un sistema común es usar el recíproco de la varianza del grupo 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                  

$F$$p$4.50890.01749

$z$$t$$50$$100$$N$

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 extremos. Las medias de los grupos previstos difieren ligeramente ( A: WLS 0.36673, robusto 0.35722; B: WLS 0.77646, robusto 0.70433; C: WLS 0.50554, robusta 0.51845), con los medios de By Csiendo menos tirados por valores extremos.

En econometría, el error estándar de Huber-White ("sandwich") es muy popular. Al igual que la corrección de Welch, esto no requiere que conozca las variaciones a priori y no requiere que calcule los pesos de sus datos y / o contingente en un modelo que puede no ser correcto. Por otro lado, no sé cómo incorporar esto con un ANOVA, lo que significa que solo los obtienes para las pruebas de códigos ficticios individuales, lo que me parece menos útil en este caso, pero los demostraré 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 

vcovHC$t$$t$$t$

Rcarwhite.adjust$p$

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

$F$$F$$p$

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

$n$

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 la prueba de Kruskal-Wallis es definitivamente la mejor protección contra los errores de tipo I, solo se puede usar con una sola variable categórica (es decir, sin predictores continuos o diseños factoriales) y tiene el menor poder de todas las estrategias discutidas. Otro enfoque no paramétrico es utilizar la regresión logística ordinal . Esto parece extraño para muchas personas, pero solo debe suponer que sus datos de respuesta contienen información ordinal legítima, lo que seguramente hacen o de lo contrario cualquier otra estrategia anterior también es inválida:

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

chi2Discrimination Indexes$p$0.0363