Alternativas à ANOVA unidirecional para dados heterocedásticos

Eu tenho dados de 3 grupos de biomassa de algas ( , , ) que contêm tamanhos de amostra desiguais ( , , ) e gostaria de comparar se esses grupos são da mesma população.B C n A = 15 n B = 13 n C = $A$$B$$C$$n_A=15$$n_B=13$$n_C=12$

A ANOVA unidirecional seria definitivamente o caminho a seguir, no entanto, ao realizar testes de normalidade nos meus dados, a heteroskedascity parece ser o principal problema. Meus dados brutos, sem nenhuma transformação, produziram uma razão de variações ( ) que é muito maior que o valor crítico ( F _ {\ rm crit} = 4.16 ) e, portanto, não posso executar ANOVA unidirecional .$F_{\max} = 19.1$$F_{\rm crit} = 4.16$

Eu também tentei a transformação para normalizar meus dados. Mesmo após tentativas de várias transformações (log, raiz quadrada, quadrado), o menor $F_{\max}$ produzido após a transformação com uma transformação $\log_{10}$ foi de $7.16$ , o que foi ainda maior em comparação com $F_{\rm crit}$ .

Alguém aqui pode me aconselhar sobre onde ir daqui? Não consigo pensar em outros métodos de transformação para normalizar pelos dados. Existem alternativas para uma ANOVA unidirecional?

PS: meus dados brutos estão abaixo:

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

Há várias opções disponíveis ao lidar com dados heterocedásticos. Infelizmente, nenhum deles é garantido para sempre trabalhar. Aqui estão algumas opções que eu estou familiarizado:

• transformações
• Welch ANOVA
• mínimos quadrados ponderados
• regressão robusta
• erros padrão consistentes de heterocedasticidade
• bootstrap
• Teste de Kruskal-Wallis
• regressão logística ordinal

Atualização: Aqui está uma demonstração R de algumas maneiras de ajustar um modelo linear (isto é, uma ANOVA ou uma regressão) quando você tem heterocedasticidade / heterogeneidade de variação.

Vamos começar analisando seus dados. Por conveniência, eu os carrego em dois quadros de dados chamados my.data(que está estruturado como acima com uma coluna por grupo) e stacked.data(que tem duas colunas: valuescom os números e indcom o indicador de grupo).

Podemos testar formalmente a heterocedasticidade com o teste 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                    

Com certeza, você tem heterocedasticidade. Vamos verificar quais são as variações dos grupos. Uma regra prática é que os modelos lineares são bastante robustos à heterogeneidade da variação, desde que a variação máxima não seja superior a maior que a variação mínima, também encontraremos essa proporção: $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

Suas variações diferem substancialmente, com as maiores B, sendo o menor,. Este é um nível problemático de heteroscedsaticidade. $19\!\times$A

Você pensou em usar transformações como o log ou a raiz quadrada para estabilizar a variação. Isso funcionará em alguns casos, mas as transformações do tipo Box-Cox estabilizam a variação pressionando os dados assimetricamente, pressionando-os para baixo com os dados mais altos pressionados mais ou pressionando-os para cima com os dados mais baixos pressionados mais. Portanto, você precisa que a variação de seus dados mude com a média para que isso funcione de maneira ideal. Seus dados têm uma enorme diferença de variação, mas uma diferença relativamente pequena entre as médias e medianas, ou seja, as distribuições se sobrepõem principalmente. Como exercício de ensino, podemos criar alguns parallel.universe.dataadicionando a todos os valores e 0,7 a$2.7$B$.7$Cpara mostrar como funcionaria:

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

O uso da transformação de raiz quadrada estabiliza esses dados muito bem. Você pode ver a melhoria dos dados do universo paralelo aqui:

Em vez de apenas tentar transformações diferentes, uma abordagem mais sistemática é otimizar o parâmetro Box-Cox (embora seja geralmente recomendável arredondá-lo para a transformação interpretável mais próxima). No seu caso, a raiz quadrada, λ = 0,5 , ou o log, λ = 0 , são aceitáveis, embora nenhum deles realmente funcione. Para os dados do universo paralelo, a raiz quadrada é a melhor: $\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

Uma abordagem mais geral é usar mínimos quadrados ponderados . Como alguns grupos ( B) se espalham mais, os dados desses grupos fornecem menos informações sobre a localização da média do que os dados de outros grupos. Podemos deixar que o modelo incorpore isso, fornecendo um peso a cada ponto de dados. Um sistema comum é usar o recíproco da variação do grupo como o 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 *

Os pesos aqui não são tão extremos. Os meios previstos grupo diferem ligeiramente ( A: WLS 0.36673, robusto 0.35722; B: WLS 0.77646, robusto 0.70433; C: WLS 0.50554, robusto 0.51845), com os meios de Be Csendo menos puxado por valores extremos.

Na econometria, o erro padrão de Huber-White ("sanduíche") é muito popular. Como a correção de Welch, isso não exige que você conheça as variações a priori e não exige que você estime pesos de seus dados e / ou contingente em um modelo que pode não estar correto. Por outro lado, não sei como incorporar isso a uma ANOVA, o que significa que você os obtém apenas para os testes de códigos fictícios individuais, o que me parece menos útil nesse caso, mas os demonstrarei de qualquer maneira:

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

Embora o teste de Kruskal-Wallis seja definitivamente a melhor proteção contra erros do tipo I, ele só pode ser usado com uma única variável categórica (ou seja, sem preditores contínuos ou desenhos fatoriais) e possui o menor poder de todas as estratégias discutidas. Outra abordagem não paramétrica é usar a regressão logística ordinal . Isso parece estranho para muitas pessoas, mas você só precisa supor que seus dados de resposta contenham informações ordinais legítimas, o que elas certamente contêm ou todas as outras estratégias acima também são inválidas:

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