Bagaimana cara mendapatkan hasil tes post-hoc Tukey HSD dalam tabel yang menunjukkan pasangan yang dikelompokkan?

13

Saya akan senang untuk melakukan tes post-hoc TukeyHSD setelah Anova dua arah saya dengan R, memperoleh tabel berisi pasangan yang disortir yang dikelompokkan berdasarkan perbedaan yang signifikan. (Maaf tentang kata-katanya, saya masih baru dengan statistik.)

Saya ingin memiliki sesuatu seperti ini:

masukkan deskripsi gambar di sini

Jadi, dikelompokkan dengan bintang atau huruf.

Ada ide? Saya menguji fungsi HSD.test()dari agricolaepaket, tetapi tampaknya tidak menangani tabel dua arah.

terjepit
sumber

Jawaban:

22

The agricolae::HSD.testfungsi tidak tepat, tetapi Anda akan perlu untuk membiarkannya tahu bahwa Anda tertarik dalam istilah interaksi . Berikut ini adalah contoh dengan dataset Stata:

library(foreign)
yield <- read.dta("http://www.stata-press.com/data/r12/yield.dta")
tx <- with(yield, interaction(fertilizer, irrigation))
amod <- aov(yield ~ tx, data=yield)
library(agricolae)
HSD.test(amod, "tx", group=TRUE)

Ini memberikan hasil yang ditunjukkan di bawah ini:

Groups, Treatments and means
a        2.1     51.17547 
ab       4.1     50.7529 
abc      3.1     47.36229 
 bcd     1.1     45.81229 
  cd     5.1     44.55313 
   de    4.0     41.81757 
    ef   2.0     38.79482 
    ef   1.0     36.91257 
     f   3.0     36.34383 
     f   5.0     35.69507 

Mereka cocok dengan apa yang akan kita dapatkan dengan perintah berikut:

. webuse yield
. regress yield fertilizer##irrigation
. pwcompare fertilizer#irrigation, group mcompare(tukey)

-------------------------------------------------------
                      |                           Tukey
                      |     Margin   Std. Err.   Groups
----------------------+--------------------------------
fertilizer#irrigation |
                 1 0  |   36.91257   1.116571    AB    
                 1 1  |   45.81229   1.116571      CDE 
                 2 0  |   38.79482   1.116571    AB    
                 2 1  |   51.17547   1.116571         F
                 3 0  |   36.34383   1.116571    A     
                 3 1  |   47.36229   1.116571       DEF
                 4 0  |   41.81757   1.116571     BC   
                 4 1  |    50.7529   1.116571        EF
                 5 0  |   35.69507   1.116571    A     
                 5 1  |   44.55313   1.116571      CD  
-------------------------------------------------------
Note: Margins sharing a letter in the group label are
      not significantly different at the 5% level.

The multcomp paket juga menawarkan visualisasi simbolik ( 'kompak surat display', lihat Algoritma untuk Compact Surat Menampilkan: Perbandingan dan Evaluasi untuk rincian lebih lanjut) dari perbandingan berpasangan yang signifikan, meskipun tidak sekarang mereka dalam format tabel. Namun, ia memiliki metode merencanakan yang memungkinkan untuk dengan mudah menampilkan hasil menggunakan boxplots. Urutan presentasi dapat diubah juga (opsi decreasing=), dan memiliki lebih banyak opsi untuk beberapa perbandingan. Ada juga paket multcompView yang memperluas fungsionalitas tersebut.

Berikut adalah contoh yang sama dianalisis dengan glht:

library(multcomp)
tuk <- glht(amod, linfct = mcp(tx = "Tukey"))
summary(tuk)          # standard display
tuk.cld <- cld(tuk)   # letter-based display
opar <- par(mai=c(1,1,1.5,1))
plot(tuk.cld)
par(opar)

Perlakuan berbagi surat yang sama tidak berbeda secara signifikan, pada tingkat yang dipilih (default, 5%).

masukkan deskripsi gambar di sini

Kebetulan, ada proyek baru, saat ini di-host di R-Forge, yang terlihat menjanjikan: factorplot . Ini mencakup tampilan berbasis garis dan huruf, serta tinjauan umum matriks (melalui plot level) dari semua perbandingan berpasangan. Makalah kerja dapat ditemukan di sini: factorplot: Meningkatkan Presentasi Kontras Sederhana dalam GLMs

chl
sumber
Terima kasih banyak atas jawaban lengkap ini! Saya akan mencoba metode yang berbeda segera setelah saya mendapatkan beberapa menit. Bersulang!
stragu
Saya mencoba fungsi paket multcomp, taruh ketika saya menggunakan fungsi 'cld ()' saya mendapatkan kesalahan 'Kesalahan: sapply (split_names, length) == 2 tidak semua BENAR' Ada yang tahu kenapa?
Terjebak
1
@ chtfn Tampaknya ada masalah dengan label variabel. Sekilas melihat kode sumber menunjukkan bahwa pesan kesalahan ini berasal dari insert_absorb()yang mencoba untuk mengekstrak pasangan perawatan. Anda mungkin dapat mencoba mengubah pemisah yang Anda gunakan untuk mengkode level istilah interaksi Anda? Tanpa contoh kerja, sulit untuk mengatakan apa yang terjadi.
chl
Saya mengetahuinya: Saya memiliki poin dalam genotipe dan nama perawatan saya, dan karena qlht () menggunakan titik untuk membagi nama pasangan, itu aneh. Terima kasih banyak atas semua bantuan Anda, chl! :)
stragu
3
Saya melihat hari ini bahwa saya sekarang harus menambahkan console=TRUEdi HSD.test()dalam rangka untuk mendapatkan meja, dalam hal seseorang mencoba ini dan melihat tidak ada hasilnya. Mungkin pembaruan agricolae.
bertopang
2

Ada fungsi yang disebut TukeyHSDitu, menurut file bantuan, menghitung seperangkat interval kepercayaan pada perbedaan antara rata-rata tingkat faktor dengan probabilitas cakupan keluarga-bijaksana yang ditentukan. Interval didasarkan pada statistik rentang Studentized, metode "Perbedaan Jujur Signifikan" Tukey. Apakah ini melakukan apa yang Anda inginkan?

http://stat.ethz.ch/R-manual/R-patched/library/stats/html/TukeyHSD.html

smillig
sumber
Terima kasih atas tanggapan Anda. Ya, saya mencoba fungsi ini, tetapi memberi saya daftar perbandingan mentah. Yang saya inginkan adalah melihat mereka dikelompokkan seperti pada gambar dalam pertanyaan saya, untuk memiliki pandangan yang jelas tentang kelompok mana yang berbeda dengan kelompok mana, dan akhirnya menambahkan nama grup pada grafik saya (misalnya: a, ab, abc, bc , c)
bertengkar