Uji Dunnett dalam R mengembalikan nilai yang berbeda setiap kali

Saya menggunakan pustaka R 'multcomp' ( http://cran.r-project.org/web/packages/multcomp/ ) untuk menghitung tes Dunnett. Saya menggunakan skrip di bawah ini:

Group <- factor(c("A","A","B","B","B","C","C","C","D","D","D","E","E","F","F","F"))
Value <- c(5,5.09901951359278,4.69041575982343,4.58257569495584,4.79583152331272,5,5.09901951359278,4.24264068711928,5.09901951359278,5.19615242270663,4.58257569495584,6.16441400296898,6.85565460040104,7.68114574786861,7.07106781186548,6.48074069840786)
data <- data.frame(Group, Value)
aov <- aov(Value ~ Group, data)
summary(glht(aov, linfct=mcp(Group="Dunnett")))

Sekarang jika saya menjalankan skrip ini melalui Konsol R beberapa kali, saya mendapatkan hasil yang sedikit berbeda setiap kali. Ini salah satu contohnya:

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0 -0.35990    0.37009  -0.972  0.76545   
C - A == 0 -0.26896    0.37009  -0.727  0.90019   
D - A == 0 -0.09026    0.37009  -0.244  0.99894   
E - A == 0  1.46052    0.40541   3.603  0.01710 * 
F - A == 0  2.02814    0.37009   5.480  0.00104 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Dan ini satu lagi:

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)    
B - A == 0 -0.35990    0.37009  -0.972   0.7654    
C - A == 0 -0.26896    0.37009  -0.727   0.9001    
D - A == 0 -0.09026    0.37009  -0.244   0.9989    
E - A == 0  1.46052    0.40541   3.603   0.0173 *  
F - A == 0  2.02814    0.37009   5.480   <0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Seperti yang Anda lihat, dua hasil di atas berbeda sangat sedikit, tetapi cukup untuk memindahkan grup terakhir (F) dari dua bintang menjadi tiga bintang, yang menurut saya mengkhawatirkan.

Saya punya beberapa pertanyaan tentang ini:

1. Kenapa ini terjadi ?! Tentunya jika Anda memasukkan data yang sama di setiap kali Anda harus mengeluarkan data yang sama.
2. Apakah ada semacam nomor acak yang digunakan di suatu tempat dalam perhitungan Dunnett?
3. Apakah variasi kecil ini setiap kali sebenarnya merupakan masalah?

Saya menjawab dua pertanyaan pertama Anda bersama-sama melalui contoh.

library(multcomp)

Group <- factor(c("A","A","B","B","B","C","C","C","D","D","D","E","E","F","F","F"))
Value <- c(5,5.09901951359278,4.69041575982343,4.58257569495584,4.79583152331272,5,5.09901951359278,4.24264068711928,5.09901951359278,5.19615242270663,4.58257569495584,6.16441400296898,6.85565460040104,7.68114574786861,7.07106781186548,6.48074069840786)
data <- data.frame(Group, Value)

fit <- aov(Value ~ Group, data)

set.seed(20140123)
Dunnet <- glht(fit, linfct=mcp(Group="Dunnett"))
summary(Dunnet)

Hasil:

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0 -0.35990    0.37009  -0.972  0.76536   
C - A == 0 -0.26896    0.37009  -0.727  0.90012   
D - A == 0 -0.09026    0.37009  -0.244  0.99895   
E - A == 0  1.46052    0.40541   3.603  0.01794 * 
F - A == 0  2.02814    0.37009   5.480  0.00112 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Jalankan lagi (tanpa mengatur seed):

summary(Dunnet)

Hasil yang berbeda:

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0 -0.35990    0.37009  -0.972  0.76535   
C - A == 0 -0.26896    0.37009  -0.727  0.90020   
D - A == 0 -0.09026    0.37009  -0.244  0.99895   
E - A == 0  1.46052    0.40541   3.603  0.01767 * 
F - A == 0  2.02814    0.37009   5.480  0.00105 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Jalankan lagi (dengan set seed):

set.seed(20140123)
Dunnet <- glht(fit, linfct=mcp(Group="Dunnett"))
summary(Dunnet)

Hasil yang sama:

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0 -0.35990    0.37009  -0.972  0.76536   
C - A == 0 -0.26896    0.37009  -0.727  0.90012   
D - A == 0 -0.09026    0.37009  -0.244  0.99895   
E - A == 0  1.46052    0.40541   3.603  0.01794 * 
F - A == 0  2.02814    0.37009   5.480  0.00112 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Dengan mengatur seed sebelum setiap run, Anda mendapatkan hasil yang konsisten. Oleh karena itu nampak bahwa angka acak sedang digunakan dalam perhitungan nilai-p.

Apakah saya pikir variasi kecil ini merupakan masalah? Saya tidak begitu menyukainya, tetapi saya akan hidup dengannya. Menggunakan set seed akan membuat hasil Anda dapat direproduksi. Saya akan merekomendasikan untuk tidak memikirkan nilai-p dalam hal berapa banyak bintang di sebelahnya - lebih baik pilih a$alpha$itu bermakna dan bermanfaat bagi Anda. Saya mencoba untuk tidak terjebak dalam apa yang terjadi 5 atau 6 poin desimal kecuali jika proyek yang saya kerjakan benar-benar membutuhkan tingkat ketelitian. Dalam hal ini saya pikir sebagian besar orang akan setuju bahwa bahkan jika perhitungan p-value sedikit berubah interpretasi hasilnya akan sama.

Anda benar, ada generasi bilangan acak yang terlibat, dan itu membuat kalkulasi bervariasi dari run-to-run. Pelakunya sebenarnya bukan prosedur Dunnett, tetapi distribusi t multivariat yang diperlukan untuk penyesuaian satu langkah.

Kode berikut menunjukkan contoh penghitungan $P(X<0)$ dengan vektor 5 dimensi $X$ memiliki multivariat $T_5$ distribusi dengan korelasi yang dapat ditukar:

> library(mvtnorm)
> cr2 <- matrix(rep(0.3, 25), nr=5); diag(cr2) <- 1
> cr2
     [,1] [,2] [,3] [,4] [,5]
[1,]  1.0  0.3  0.3  0.3  0.3
[2,]  0.3  1.0  0.3  0.3  0.3
[3,]  0.3  0.3  1.0  0.3  0.3
[4,]  0.3  0.3  0.3  1.0  0.3
[5,]  0.3  0.3  0.3  0.3  1.0
> b <- pmvt(lower=rep(-Inf,5), upper=rep(0,5), delta=rep(0,5), df=5, corr=cr2)
> a <- pmvt(lower=rep(-Inf,5), upper=rep(0,5), delta=rep(0,5), df=5, corr=cr2)
> all.equal(a,b)
[1] "Attributes: < Component 1: Mean relative difference: 0.1527122 >"
[2] "Mean relative difference: 0.0003698006"     

Jika ini yang menjadi perhatian, panggil saja set.seeddengan argumen apa pun sebelum perhitungan untuk membuatnya benar-benar dapat diproduksi kembali.

Omong-omong, ada pengakuan dan kuantifikasi kesalahan dalam output dari glht:

> ss <- summary(glht(aov, linfct=mcp(Group="Dunnett")))
> attr(ss$test$pvalues, "error")
[1] 0.0006597562