Alternatif untuk ANOVA satu arah untuk data heteroskedastik

36

Saya memiliki data dari 3 kelompok biomassa alga ( A , ,B ) yang berisi ukuran sampel yang tidak sama ( n A = 15 , n B = 13 , n C = 12 ) dan saya ingin membandingkan jika kelompok-kelompok ini berasal dari populasi yang sama .CnA=15nB=13nC=12

ANOVA satu arah pasti akan menjadi cara untuk pergi, namun setelah melakukan tes normalitas pada data saya, heteroskedascity tampaknya menjadi masalah utama. Data mentah saya, tanpa transformasi apapun, menghasilkan rasio varians ( ) yang sangat jauh lebih tinggi dari nilai kritis ( F c r i t = 4.16 ) dan karena itu saya tidak dapat melakukan ANOVA satu arah.Fmax=19.1Fcrit=4.16

Saya juga mencoba transformasi untuk menormalkan data saya. Bahkan setelah percobaan berbagai transformasi (log, akar kuadrat, kuadrat), terendah yang dihasilkan setelah transformasi dengan transformasi log 10 adalah 7,16 , yang masih lebih tinggi dibandingkan dengan F c r i t .Fmaxlog107.16Fcrit

Adakah yang bisa menyarankan saya ke mana harus pergi dari sini? Saya tidak bisa memikirkan metode transformasi lain untuk dinormalkan dengan data. Apakah ada alternatif selain ANOVA satu arah?

PS: data mentah saya di bawah:

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
Rick L.
sumber
2
Uji F jelas menunjukkan bahwa kelompok-kelompok tersebut bukan dari populasi yang sama. Langkah selanjutnya adalah menafsirkan hasil ini, dimulai dengan visualisasi yang jelas dan deskripsi kuantitatif dari data yang diuraikan oleh kelompok.
whuber
Ada juga tes permutasi dalam paket koin. Untuk ANOVA itu akan menjadi fungsi oneway_test. Contoh Quick-R: statmethods.net/stats/resampling.html
user16263

Jawaban:

73

Ada sejumlah opsi yang tersedia saat berurusan dengan data heteroskedastik. Sayangnya, tidak satu pun dari mereka dijamin untuk selalu bekerja. Berikut adalah beberapa opsi yang saya kenal:

  • transformasi
  • Welch ANOVA
  • kuadrat terkecil tertimbang
  • regresi yang kuat
  • heteroscedasticity konsisten kesalahan standar
  • bootstrap
  • Tes Kruskal-Wallis
  • regresi logistik ordinal

Pembaruan: Ini adalah demonstrasi di R beberapa cara pemasangan model linier (yaitu, ANOVA atau regresi) ketika Anda memiliki heteroskedastisitas / heterogenitas varians.

Mari kita mulai dengan melihat data Anda. Untuk kenyamanan, saya minta mereka memuatnya ke dalam dua frame data yang disebut my.data(yang terstruktur seperti di atas dengan satu kolom per grup) dan stacked.data(yang memiliki dua kolom: valuesdengan angka dan inddengan indikator grup).

Kami dapat secara resmi menguji heteroskedastisitas dengan uji 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                    

Benar saja, Anda memiliki heteroskedastisitas. Kami akan memeriksa untuk melihat apa varian dari grup tersebut. Aturan praktisnya adalah bahwa model linear cukup kuat untuk heterogenitas varians selama varians maksimum tidak lebih dari lebih besar dari varians minimum, jadi kami juga akan menemukan rasio itu: 4×

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

Varians Anda berbeda secara substansial, dengan yang terbesar B, yaitu yang terkecil,. Ini adalah tingkat heteroscedsaticity yang bermasalah. 19×A

parallel.universe.data2.7B.7CUntuk menunjukkan cara kerjanya:

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

Menggunakan transformasi akar kuadrat menstabilkan data tersebut dengan cukup baik. Anda dapat melihat peningkatan untuk data paralel semesta di sini:

masukkan deskripsi gambar di sini

λλ=.5λ=0

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

masukkan deskripsi gambar di sini

Fdf = 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

Pendekatan yang lebih umum adalah dengan menggunakan kuadrat terkecil tertimbang . Karena beberapa kelompok ( B) menyebar lebih banyak, data dalam kelompok tersebut memberikan lebih sedikit informasi tentang lokasi rata-rata daripada data dalam kelompok lain. Kita dapat membiarkan model memasukkan ini dengan memberikan bobot pada setiap titik data. Sistem yang umum adalah menggunakan kebalikan dari varians grup sebagai bobot:

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                  

Fhal4.50890.01749

masukkan deskripsi gambar di sini

zt50100N

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 *

Bobot di sini tidak ekstrem. Berarti kelompok diprediksi sedikit berbeda ( A: WLS 0.36673, kuat 0.35722; B: WLS 0.77646, kuat 0.70433; C: WLS 0.50554, kuat 0.51845), dengan cara Bdan Cyang kurang ditarik oleh nilai-nilai ekstrim.

Dalam ekonometrik kesalahan standar Huber-White ("sandwich") sangat populer. Seperti koreksi Welch, ini tidak mengharuskan Anda untuk mengetahui varian a-priori dan tidak mengharuskan Anda memperkirakan bobot dari data Anda dan / atau bergantung pada model yang mungkin tidak benar. Di sisi lain, saya tidak tahu bagaimana menggabungkan ini dengan ANOVA, yang berarti bahwa Anda hanya mendapatkannya untuk tes kode boneka individu, yang menurut saya kurang membantu dalam kasus ini, tetapi saya akan tetap menunjukkannya:

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 

vcovHCttt

Rcarwhite.adjusthal

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

FFhal

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

Meskipun uji Kruskal-Wallis jelas merupakan perlindungan terbaik terhadap kesalahan tipe I, tes ini hanya dapat digunakan dengan variabel kategori tunggal (yaitu, tidak ada prediktor berkelanjutan atau desain faktorial) dan memiliki kekuatan paling sedikit dari semua strategi yang dibahas. Pendekatan non-parametrik lain adalah dengan menggunakan regresi logistik ordinal . Ini terlihat aneh bagi banyak orang, tetapi Anda hanya perlu berasumsi bahwa data respons Anda mengandung informasi ordinal yang sah, yang pasti mereka lakukan atau strategi lain di atas tidak valid:

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 Indexeshal0.0363

gung - Reinstate Monica
sumber
1
Cintai jawaban Anda! Sekarang ada sedikit rumah untuk menangani data saya! :) Dari semua pendekatan ini, yang akan sangat Anda rekomendasikan untuk jenis data saya, terutama dalam hal menyediakan kekuatan statistik yang cukup? Saya tidak terlalu suka menggunakan pendekatan non-parametrik yaitu tes Kruskal-Wallis karena itu berarti saya cenderung melakukan kesalahan Tipe I? Koreksi saya jika saya salah.
Rick L.
2
Saya akan mengerjakan beberapa materi lagi di sini ketika saya memiliki kesempatan, tetapi tidak, tes Kruskal-Wallis mungkin memberikan perlindungan terbaik Anda terhadap kesalahan tipe I. OTOH, itu mungkin bukan opsi paling kuat yang Anda miliki.
gung - Reinstate Monica
4
Perhatikan bahwa di perpustakaan carada juga opsi untuk mengatur white.adjust=Tuntuk berurusan dengan heteroskedacity menggunakan White-adjusted heteroscedasticity standard error yang diperbaiki
Tom Wenseleers
3
Ya itu hanya menyebutkan untuk lm's, tetapi juga tampaknya bekerja untuk aov' s (pilihan untuk white.adjustyang white.adjust=c(FALSE, TRUE, "hc3", "hc0", "hc1", "hc2", "hc4")- untuk info lebih lanjut lihat ?hccm)
Tom Wenseleers
1
@Nakx, tidak, bukan yang saya tahu. Ini adalah aturan praktis; itu tidak mungkin di makalah dalam literatur statistik. Mungkin di beberapa buku pengantar.
gung - Reinstate Monica