Alternatywy dla jednokierunkowej ANOVA dla danych heteroskedastycznych

Mam dane z 3 grup biomasy alg ( , , ), które zawierają nierówne wielkości próbek ( n_A = 15 , n_B = 13 , n_C = 12 ) i chciałbym porównać, czy te grupy pochodzą z tej samej populacji.B C n A = 15 n B = 13 n C = $A$$B$$C$$n_A=15$$n_B=13$$n_C=12$

Jednokierunkowa ANOVA zdecydowanie byłaby właściwą drogą, jednak po przeprowadzeniu testów normalności na moich danych, heteroskedascity wydaje się głównym problemem. Moje nieprzetworzone dane, bez żadnej transformacji, wytworzyły współczynnik wariancji ( $F_{\max} = 19.1$ ), który jest znacznie wyższy niż wartość krytyczna ( $F_{\rm crit} = 4.16$ ) i dlatego nie mogę wykonać jednokierunkowej ANOVA .

Próbowałem także transformacji, aby znormalizować moje dane. Nawet po próbach różnych transformacji (log, pierwiastek kwadratowy, kwadrat), najniższy $F_{\max}$ uzyskany po transformacji z transformacją $\log_{10}$ wynosił $7.16$ , co było wciąż wyższe w porównaniu do $F_{\rm crit}$ .

Czy ktoś może mi doradzić, gdzie mam się udać? Nie mogę wymyślić innych metod transformacji normalizujących dane. Czy są jakieś alternatywy dla jednokierunkowej ANOVA?

PS: moje surowe dane są poniżej:

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

Istnieje wiele opcji dostępnych w przypadku danych heteroscedastycznych. Niestety, żadna z nich nie gwarantuje, że zawsze będzie działać. Oto kilka znanych mi opcji:

• transformacje
• Welch ANOVA
• ważone najmniejszych kwadratów
• solidna regresja
• heteroscedastyczność spójne błędy standardowe
• bootstrap
• Test Kruskala-Wallisa
• porządkowa regresja logistyczna

Aktualizacja: Oto demonstracja R niektórych sposobów dopasowania modelu liniowego (tj. ANOVA lub regresji), gdy masz heteroscedastyczność / heterogeniczność wariancji.

Zacznijmy od spojrzenia na twoje dane. Dla wygody mam je wczytane do dwóch ramek danych o nazwie my.data(które mają strukturę jak wyżej z jedną kolumną na grupę) i stacked.data(która ma dwie kolumny: valuesz liczbami i indze wskaźnikiem grupy).

Możemy formalnie przetestować heteroscedastyczność za pomocą testu Levene'a:

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                    

Rzeczywiście, masz heteroscedastyczność. Sprawdzimy, jakie są wariancje grup. Ogólna zasada jest taka, że ​​modele liniowe są dość odporne na niejednorodność wariancji, o ile maksymalna wariancja jest nie większa niż większa niż wariancja minimalna, więc znajdziemy również ten stosunek: $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

Twoje wariancje różnią się zasadniczo, z największym, Bbędąc najmniejsze, . Jest to problematyczny poziom heteroscedsatyczności. $19\!\times$A

Pomyślałeś, aby użyć transformacji, takich jak log lub pierwiastek kwadratowy, aby ustabilizować wariancję. W niektórych przypadkach będzie to działać, ale transformacje typu Box-Cox stabilizują wariancję, ściskając dane asymetrycznie, albo ściskając je w dół, gdy najwyższe dane są ściśnięte najbardziej, albo ściskając je w górę, gdy najniższe dane są ściśnięte najbardziej. Dlatego potrzebujesz wariancji danych, aby zmienić się ze średnią, aby działało to optymalnie. Twoje dane mają ogromną różnicę wariancji, ale stosunkowo niewielką różnicę między średnimi a medianami, tzn. Rozkłady w większości pokrywają się. Jako ćwiczenie możemy stworzyć niektóre parallel.universe.data, dodając do wszystkich wartości i do$2.7$B$.7$Caby pokazać, jak to będzie dział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

Zastosowanie transformacji pierwiastka kwadratowego dość dobrze stabilizuje te dane. Tutaj możesz zobaczyć poprawę danych równoległego wszechświata:

Zamiast tylko próbować różnych transformacji, bardziej systematycznym podejściem jest optymalizacja parametru Box-Cox (chociaż zwykle zaleca się zaokrąglenie tego do najbliższej możliwej do interpretacji transformacji). W twoim przypadku albo pierwiastek kwadratowy, , albo log, , są dopuszczalne, chociaż tak naprawdę nie działa. W przypadku danych z równoległego wszechświata pierwiastek kwadratowy jest najlepszy: λ = .5 λ = $\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) 

Ponieważ ten przypadek jest ANOVA (tj. Bez zmiennych ciągłych), jednym ze sposobów radzenia sobie z heterogenicznością jest zastosowanie korekcji Welcha do mianownika stopni swobody w teście (nb , zamiast wartości ułamkowej ): $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

Bardziej ogólnym podejściem jest stosowanie ważonych najmniejszych kwadratów . Ponieważ niektóre grupy ( B) rozprzestrzeniają się bardziej, dane w tych grupach dostarczają mniej informacji o lokalizacji średniej niż dane w innych grupach. Możemy pozwolić modelowi uwzględnić to, podając wagę dla każdego punktu danych. Powszechnym systemem jest stosowanie odwrotności wariancji grupowej jako wagi:

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                  

Daje to nieco inne wartości i niż nieważona ANOVA ( , ), ale dobrze rozwiązuje problem niejednorodności: $F$$p$4.50890.01749

Jednak najmniej ważone kwadraty nie są panaceum. Jednym niewygodnym faktem jest to, że jest to właściwe tylko wtedy, gdy wagi są odpowiednie, co oznacza między innymi, że są one znane z góry. Nie dotyczy to również nienormalności (takiej jak pochylenie) ani wartości odstających. Korzystanie z wagi szacunkowe dane często praca w porządku, chociaż, szczególnie jeśli masz wystarczająco dużo danych, aby oszacować wariancję z rozsądną dokładnością (ta jest analogiczna do idei używając -Tabela zamiast -Tabela gdy masz lubt 50 100 $z$$t$$50$$100$stopnie swobody), twoje dane są wystarczająco normalne i nie wydajesz się mieć żadnych wartości odstających. Niestety masz stosunkowo niewiele danych (13 lub 15 na grupę), niektóre wypaczenia i być może pewne wartości odstające. Nie jestem pewien, czy są one wystarczająco złe, aby zrobić z tego wielką sprawę, ale można mieszać ważone najmniejsze kwadraty z niezawodnymi metodami. Zamiast używać wariancji jako miary rozprzestrzeniania się (która jest wrażliwa na wartości odstające, szczególnie przy niskim ), możesz użyć odwrotności zakresu międzykwartylowego (na co nie ma wpływu nawet 50% wartości odstających w każdej grupie). Wagi te można następnie połączyć z solidną regresją przy użyciu innej funkcji utraty, takiej jak bisquare Tukeya: $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 *

Ciężary tutaj nie są tak ekstremalne. Przewidziane środki grupy różnią się nieznacznie ( A: WLS 0.36673, wytrzymałe 0.35722; B: WLS 0.77646, wytrzymałe 0.70433; C: WLS 0.50554, niezawodny 0.51845), przy czym środki Bi Csą mniej ciągnięte ekstremalnych wartości.

W ekonometrii błąd standardowy Hubera-White'a („sandwich”) jest bardzo popularny. Podobnie jak w przypadku korekcji Welcha, nie wymaga to uprzedniej znajomości odchyleń i nie wymaga oszacowania wag na podstawie danych i / lub uzależnienia od modelu, który może być nieprawidłowy. Z drugiej strony nie wiem, jak to włączyć z ANOVA, co oznacza, że ​​dostajesz je tylko do testów poszczególnych kodów fikcyjnych, co wydaje mi się mniej pomocne w tym przypadku, ale i tak je pokażę:

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 

Funkcja vcovHCoblicza heteroscedastyczną spójną macierz wariancji-kowariancji dla twoich bet (kodów fikcyjnych), co oznaczają litery w wywołaniu funkcji. Aby uzyskać standardowe błędy, wyodrębnij główną przekątną i weź pierwiastki kwadratowe. Aby uzyskać testy dla twoich bet, dzielisz swoje oszacowania współczynnika przez SE i porównujesz wyniki z odpowiednim rozkładem ( rozkład z resztkowym stopniem swobody). t $t$$t$$t$

W Rszczególności dla użytkowników @TomWenseleers zauważa w komentarzach poniżej, że funkcja ? Anova w carpakiecie może zaakceptować white.adjustargument, aby uzyskać wartość dla czynnika wykorzystującego błędy spójne z heteroscedastycznością. $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

Możesz spróbować uzyskać empiryczny szacunek tego, jak wygląda rzeczywisty rozkład próbkowania statystyki testowej przez ładowanie początkowe . Po pierwsze, tworzysz prawdziwy zerowy, czyniąc wszystkie środki grupy dokładnie równymi. Następnie ponownie próbkujesz z zastępowaniem i obliczasz swoją statystykę testową ( ) na każdej próbce, aby uzyskać empiryczną ocenę rozkładu próbkowania poniżej wartości zerowej z twoimi danymi, niezależnie od ich statusu w odniesieniu do normalności lub jednorodności. Proporcja tego rozkładu próbkowania, która jest tak ekstremalna lub bardziej ekstremalna niż obserwowana statystyka testu, to wartość : F $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

Pod pewnymi względami ładowanie początkowe jest ostatecznym podejściem z ograniczonym założeniem do przeprowadzania analizy parametrów (np. Średnich), ale zakłada, że ​​twoje dane stanowią dobrą reprezentację populacji, co oznacza, że ​​masz rozsądną wielkość próby. Ponieważ twoje są małe, może być mniej wiarygodne. Prawdopodobnie ostateczną ochroną przed nienormalnością i niejednorodnością jest zastosowanie testu nieparametrycznego. Podstawową nieparametryczną wersją ANOVA jest test Kruskala-Wallisa : $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

Chociaż test Kruskala-Wallisa jest zdecydowanie najlepszą ochroną przed błędami typu I, można go stosować tylko z jedną zmienną kategorialną (tj. Bez ciągłych predyktorów lub układów czynnikowych) i ma najmniejszą moc ze wszystkich omawianych strategii. Innym nieparametrycznym podejściem jest użycie porządkowej regresji logistycznej . Dla wielu osób wydaje się to dziwne, ale wystarczy założyć, że dane odpowiedzi zawierają prawidłowe informacje porządkowe, które z pewnością robią, w przeciwnym razie każda inna strategia powyżej również jest nieprawidłowa:

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

Wynik może nie być jasny, ale test modelu jako całości, który w tym przypadku jest testem twoich grup, jest chi2poniżej Discrimination Indexes. Wymieniono dwie wersje: test ilorazu wiarygodności i test punktowy. Test współczynnika prawdopodobieństwa jest zwykle uważany za najlepszy. Daje wartość wynoszącą . $p$0.0363