ANOVA di R: aov(), tabel ANOVA, dan Tukey HSD

Mengapa ANOVA, bukan setumpuk uji t

Uji t membandingkan dua rata-rata. Dengan tiga kelompok atau lebih, godaannya adalah menjalankan uji t untuk setiap pasangan, dan justru di situlah jebakannya. Setiap uji pada taraf 0,05 membawa risiko positif palsu 5%, dan risiko itu menumpuk: dengan 4 kelompok ada 6 uji berpasangan dan peluang sekitar 26% mendapat setidaknya satu hasil "signifikan" yang palsu; dengan 5 kelompok (10 uji), sekitar 40%. Kamu akan menciptakan temuan dari derau semata.

ANOVA (ANalysis Of VAriance, analisis varians) mengatasi ini dengan mengajukan satu pertanyaan dengan satu p-value: apakah semua rata-rata kelompok sama, atau setidaknya ada satu yang berbeda? Meski namanya begitu, ANOVA membandingkan rata-rata, hanya saja caranya lewat analisis varians: kalau rata-rata kelompok lebih menyebar daripada yang bisa dijelaskan oleh derau di dalam kelompok, ada sesuatu yang nyata sedang terjadi.

ANOVA satu arah dengan aov()

PlantGrowth memang dibuat untuk ini: berat tanaman pada satu kondisi kontrol dan dua kondisi perlakuan. Buat modelnya dengan aov(), lalu (ini penting) cetak tabelnya dengan summary():

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

Baca formulanya sebagai "apakah weight bergantung pada group?" Variabel pengelompokan harus berupa factor. PlantGrowth$group sudah berupa factor, tapi kalau kelompokmu dikodekan sebagai angka (dosis, ID batch), bungkus dulu: aov(y ~ factor(dose), ...). Kalau tidak, aov() diam-diam memasang garis regresi melalui kode kelompok alih-alih membandingkan rata-rata kelompok: modelnya salah, dan tidak ada pesan error.

Membaca tabel ANOVA

Tabelnya terdiri dari dua baris (faktor dan residual) dan lima kolom. Sel demi sel:

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df: derajat kebebasan. Baris faktor mendapat jumlah kelompok − 1 (3 kelompok → 2); baris residual mendapat jumlah observasi − jumlah kelompok (30 − 3 = 27). Cara cepat untuk memastikan R membaca desain yang kamu maksud.
• Sum Sq: variasi, dibagi menjadi dua tumpukan. Baris group adalah tumpukan antarkelompok: seberapa jauh rata-rata kelompok dari rata-rata keseluruhan. Residuals adalah tumpukan dalam kelompok: seberapa besar tanaman bervariasi di sekitar rata-rata kelompoknya sendiri. Jika dijumlahkan, keduanya sama dengan total variasi dalam data.
• Mean Sq: setiap Sum Sq dibagi Df-nya, sehingga tumpukan variasi berubah menjadi laju per derajat kebebasan yang bisa dibandingkan. Mean Sq residual (0,389) adalah tingkat derau.
• F value: rasionya, yaitu Mean Sq dari group dibagi Mean Sq dari Residuals (1,8832 / 0,3886 ≈ 4,85). Jika semua rata-rata kelompok benar-benar sama, rasio ini berkisar di sekitar 1. Makin besar F, makin jauh perbedaan antarkelompok melampaui apa yang bisa dijelaskan oleh derau.
• Pr(>F): p-value, yaitu peluang mendapatkan F sebesar ini jika ketiga rata-rata sebenarnya identik. Di sini 0,016, cukup kecil pada taraf konvensional 0,05 untuk menolak "semuanya sama". Seperti biasa, ini bukan peluang bahwa hipotesis nol benar, dan nilai ini tidak mengatakan apa pun tentang kelompok mana yang berbeda atau seberapa besar perbedaannya.

Poin terakhir itulah batasan pentingnya: F yang signifikan berarti "ada perbedaan, entah di mana." Tidak lebih.

Kelompok mana yang berbeda? TukeyHSD()

Pertanyaan lanjutannya membutuhkan uji post-hoc. Honest Significant Difference dari Tukey menguji setiap pasangan sambil menjaga tingkat error family-wise tetap 5% untuk semua perbandingan sekaligus:

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

Satu baris per pasangan, empat angka per baris:

• diff: selisih rata-rata yang diperkirakan (kelompok yang disebut kedua dikurangi yang pertama).
• lwr, upr: interval kepercayaan family-wise 95% untuk selisih itu.
• p adj: p-value yang sudah disesuaikan karena ada tiga perbandingan.

Untuk PlantGrowth: trt2-trt1 menunjukkan selisih sekitar 0,87 dengan p adj ≈ 0,012 dan interval yang tidak menyentuh nol: perlakuan 2 tumbuh lebih besar daripada perlakuan 1. Kedua baris perlakuan lawan kontrol (trt1-ctrl, trt2-ctrl) punya interval yang melintasi nol dan p adj jauh di atas 0,05: tidak ada perlakuan yang bisa dibedakan dari kontrol pada sampel ini. Aturan praktisnya sama seperti interval kepercayaan di mana pun: interval tidak mencakup nol ⇔ p adj di bawah 0,05.

Perhatikan bagaimana F keseluruhan berkata "ada perbedaan" sementara Tukey menemukannya tepat di satu pasangan. Struktur dua langkah itulah inti desainnya: satu uji keseluruhan yang jujur, lalu kerja detektif berpasangan yang dikoreksi dengan benar.

ANOVA dua arah: dua faktor dan interaksinya

Dengan dua variabel pengelompokan, satu pemanggilan menguji keduanya, plus apakah keduanya berinteraksi. ToothGrowth menyilangkan jenis suplemen (supp) dengan dosis (0,5; 1; dan 2 mg, berupa angka, jadi perlu factor()):

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) diuraikan menjadi tiga efek, masing-masing satu baris tabel:

• supp: dirata-ratakan atas semua dosis, apakah jenis suplemen berpengaruh? (Ya: p ≈ 0,0002.)
• factor(dose): dirata-ratakan atas semua suplemen, apakah dosis berpengaruh? (Sangat jelas: p-nya sangat kecil.)
• supp:factor(dose): interaksinya, yaitu apakah efek suplemen berubah tergantung dosisnya? Di sini p ≈ 0,022, jadi ya. Jus jeruk mengalahkan asam askorbat pada dosis rendah, tapi selisihnya hilang pada dosis 2 mg.

Interaksi yang signifikan adalah label peringatan bagi efek utama: "jenis suplemen berpengaruh" hanya benar secara rata-rata, dan rata-rata itu menyembunyikan cerita yang bergantung pada dosis. Saat interaksinya signifikan, jelaskan kombinasinya (dengan TukeyHSD(fit2, "supp:factor(dose)") atau tabel rata-rata kelompok lewat ringkasan per kelompok) alih-alih hanya melaporkan efek utama. Pakai + sebagai ganti * hanya kalau kamu sengaja ingin model tanpa interaksi.

Asumsi, dan alternatif berbasis peringkat

Matematika ANOVA bersandar pada tiga asumsi, diurutkan dari yang paling tidak bisa ditawar:

• Independensi: observasi tidak saling memengaruhi. Tidak ada yang bisa memperbaiki pelanggarannya setelah data terkumpul; ini sifat dari desain penelitian.
• Varians yang sama antarkelompok: Mean Sq residual adalah satu perkiraan derau gabungan, jadi tiap kelompok seharusnya kira-kira sama berderaunya. Cek dalam satu baris: bartlett.test(weight ~ group, data = PlantGrowth) (p-value yang besar berarti tidak ada bukti bahwa variansnya tidak sama).
• Residual yang kira-kira normal: paling tidak penting pada desain yang seimbang dengan ukuran kelompok yang memadai; lihat sekilas histogram atau boxplot residual lewat residuals(fit).

Saat datanya sangat miring, ordinal, atau penuh pencilan, alternatif berbasis peringkat untuk ANOVA satu arah adalah kruskal.test(weight ~ group, data = PlantGrowth), saudara versi banyak kelompok dari wilcox.test(). Cukup satu baris, dan uji ini menukar sedikit kekuatan uji dengan ketahanan.

Yang kamu bawa pulang

• ANOVA mengajukan satu pertanyaan tentang 3+ rata-rata kelompok dengan satu p-value: solusi untuk masalah banyaknya uji t.
• fit <- aov(y ~ group, data = df); summary(fit), dan variabel pengelompokan harus berupa factor.
• F value adalah sinyal antarkelompok dibagi derau dalam kelompok; Pr(>F) yang kecil berarti "ada perbedaan, entah di mana", tidak lebih.
• TukeyHSD(fit) menemukan pasangan mana yang berbeda, dengan koreksi perbandingan berganda yang sudah tertanam.
• a * b memodelkan dua faktor plus interaksinya; interaksi yang signifikan berarti efek utama saja tidak menceritakan seluruh ceritanya.
• Cek kesamaan varians (bartlett.test) dan normalitas residual; kruskal.test() adalah alternatif berbasis peringkat.

Berikutnya: dari membandingkan rata-rata kelompok ke memodelkan hubungan, yaitu regresi linear dengan lm().