İkili ve çoklu karşılaştırma yöntemleri için hangi paketleri ve nasıl kullanacağımız ile ilgili örnekler verilecektir. İlk olarak derste kullandığımız sürücü kursu verisini ele alalım.

Örnek 1:

A, B, C ve D biçiminde 4 farklı sürücü kursunun öğrencilerinin ehliyet sınavından aldığı notların ortalamalarına göre karşılaştırılmasını inceleyelim.

x<- c(70,65,75,72,74,68,50,52,60,62,54,55,80,82,77,88,83,85,90,92,82,86,88,84)
faktor<-c(rep("A",6),rep("B",6),rep("C",6),rep("D",6))
veri<-data.frame(x,faktor)
veri
##     x faktor
## 1  70      A
## 2  65      A
## 3  75      A
## 4  72      A
## 5  74      A
## 6  68      A
## 7  50      B
## 8  52      B
## 9  60      B
## 10 62      B
## 11 54      B
## 12 55      B
## 13 80      C
## 14 82      C
## 15 77      C
## 16 88      C
## 17 83      C
## 18 85      C
## 19 90      D
## 20 92      D
## 21 82      D
## 22 86      D
## 23 88      D
## 24 84      D
boxplot(x ~ faktor)

Bu 4 kurs arasında fark olup olmadığını ANOVA ile inceleyelim.

anova<-aov(x~faktor,data=veri)
summary(anova)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## faktor       3   3567  1189.2   73.78 5.43e-11 ***
## Residuals   20    322    16.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(F\) testinin sonucuna göre \(p-value<0.005\) olduğundan \(H_{0}:\mu _{1}=\mu _{2}=\mu _{3}=\mu _{4}\) hipotezi red edilir.

Varsayımların Kontrolü

Ayrıca, bu verinin ANOVA varsayımlarını sağladığını göstermiştik. Şimdi, tekrar gözden geçirelim. Artıklar üzerinden normallik testleri yapılır.

#normallik testleri
qqnorm(residuals(anova))

shapiro.test(residuals(anova))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(anova)
## W = 0.95437, p-value = 0.336
ks.test(residuals(anova),"pnorm",mean(residuals(anova)),sd(residuals(anova)))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  residuals(anova)
## D = 0.094792, p-value = 0.9823
## alternative hypothesis: two-sided
library(nortest)
lillie.test(residuals(anova))  #nortest package  "library(nortest)"```
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuals(anova)
## D = 0.094792, p-value = 0.8346

Yukarıda yapılan tüm test sonuçlarına göre \(p-value>0.005\) olduğundan artıkların Normal dağılıma sahiptirsıfır hipotezi kabul edilir.

Şimdi, faktör düzeylerinin homojen varyanslılığını kontrol edelim. Grafiksel yöntem ve test ile yapabiliriz.

plot(fitted.values(anova),residuals(anova))

plot(anova,1)

bartlett.test(x~faktor)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  x by faktor
## Bartlett's K-squared = 0.30707, df = 3, p-value = 0.9587
library(car)
leveneTest(x, factor(faktor)) 
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.1007 0.9587
##       20
leveneTest(x, factor(faktor),mean)
## Levene's Test for Homogeneity of Variance (center = mean)
##       Df F value Pr(>F)
## group  3  0.1917 0.9008
##       20

Tahmin edlen yanıt değişkeni ile onlara karşılık gelen artıklar kullanılarak çizilen bu grafiklerde noktalar iki doğru arasında kaldığından homojen varyanslılık varsayımı sağlanır.

Ayrıca, yapılan testler için \(p-value>0.005\) olduğundan homojen varyanslılık varsayımı sağlanır.

Bu durumda, faktör düzeylerinin ortalamaları arasındaki farklılıkların hangi düzey veya düzeylerden kaynaklandığını belirlemek için ikili ve çoklu karşılaştırma yöntemleri kullanılır.

İkili ve Çoklu Karşılaştırmalar

İkili ve çoklu karşılaştırmalarda “multcomp” paketindeki “glht” (General Linear Hypotheses) fonksiyonundan faydalanacağız. Ayrıca, “agricolae” paketinde faktör düzeylerini karşılaştırmak için bir çok test bulunmaktadır. Bu testler: LSD, Bonferroni and other p-adjust, HSD, Waller, Student Newman Keuls SNK, Duncan, REGW, Scheffe.

1. Lineer bağıntılar ve Dik lineer bağıntılar

\(H_{0}:\mu _{1}+\mu _{2}=\mu _{3}+\mu _{4}\) veya \(H_{0}:\mu _{1}+\mu _{2}-\mu _{3}-\mu _{4}=0\) hipotezini lineer bağıntılar yöntemi ile test edelim. “glht” fonksiyonu içinde bu ortalamalar ile ilgili faktörleri aşağıdaki gibi yazarak bu test yapılır.

library(multcomp)
result<-glht(anova,linfct = mcp(faktor = "A + B -C - D == 0"))
result
## 
##   General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Linear Hypotheses:
##                    Estimate
## A + B - C - D == 0   -43.33
summary(result) 
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = x ~ faktor, data = veri)
## 
## Linear Hypotheses:
##                    Estimate Std. Error t value Pr(>|t|)    
## A + B - C - D == 0  -43.333      3.278  -13.22 2.41e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Bu test hem \(t\) hem de \(F\) testi ile yapılabilir, buradaki sonuçta \(t\) testinin değeri verilir. Görüldüğü gibi \(p-value<0.005\) olduğundan \(H_{0}\) hipotezi reddedilir.

\(H_{0}:\mu _{1}-\mu _{2}=0\) gibi ikili karşılaştırmalar için de bu yöntem kullanılaabilir. Ayrıca, aynı anda birden fazla hipotezi de test edebiliriz.

Örneğin, \(H_{0}:\mu _{1}+\mu _{2}-\mu _{3}-\mu _{4}=0\) ve \(H_{0}:\mu _{1}-\mu _{2}=0\) hipotezlerini test edelim.

result1<-glht(anova,linfct = mcp(faktor = c("A + B -C - D == 0", "A - B == 0")))
result1
## 
##   General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Linear Hypotheses:
##                    Estimate
## A + B - C - D == 0   -43.33
## A - B == 0            15.17
summary(result1) 
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = x ~ faktor, data = veri)
## 
## Linear Hypotheses:
##                    Estimate Std. Error t value Pr(>|t|)    
## A + B - C - D == 0  -43.333      3.278 -13.220  < 1e-10 ***
## A - B == 0           15.167      2.318   6.544 4.47e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

\(H_{0}:\mu _{1}+\mu _{2}-\mu _{3}-\mu _{4}=0\) için \(p-value<1e-10\) ve \(H_{0}:\mu _{1}-\mu _{2}=0\) için \(p-value=4.47e-06\) olarak bulunmuştur. Her iki test içinde \(H_{0}\) hipotezleri reddedilir.

3 farklı hipotezi test edelim. \(H_{0}:\mu _{1}+\mu _{2}-\mu _{3}-\mu _{4}=0\) , \(H_{0}:\mu _{1}-\mu _{2}=0\) ve \(H_{0}:\mu _{3}-\mu _{4}=0\) hipotezlerini test edelim.

result2<-glht(anova,linfct = mcp(faktor = c("A + B -C - D == 0", "A - B == 0", "C - D ==0") ) )
result2
## 
##   General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Linear Hypotheses:
##                    Estimate
## A + B - C - D == 0   -43.33
## A - B == 0            15.17
## C - D == 0            -4.50
summary(result2) 
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = x ~ faktor, data = veri)
## 
## Linear Hypotheses:
##                    Estimate Std. Error t value Pr(>|t|)    
## A + B - C - D == 0  -43.333      3.278 -13.220   <0.001 ***
## A - B == 0           15.167      2.318   6.544   <0.001 ***
## C - D == 0           -4.500      2.318  -1.941    0.181    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Sonuca göre ilk iki test için \(H_{0}\) hipotezleri reddedilir. Ancak, üçünücü \(H_{0}\) hipotezi için \(p-value=0.181>0.005\) olduğundan \(H_{0}:\mu _{3}-\mu _{4}=0\) kabul edilir. Dolayısıyla, C ve D faktörlerinin ortlamaları arasında anlamlı bir fark yoktur.

2. Tukey Testi

Tukey testi R da iki farklı yol ile yapılabilir. Birincisi “stats” paketindeki “TukeyHSD” fonksiyonu ile, ikincisi ise “multcomp” paketindeki “glht” fonksiyonunu kullanmaktır. Tukey testi ile faktör düzeyleri arasında tüm mümkün ikili karşılaştırmalar yapılır. Bu örneğimiz için 4 faktör olduğundan \(\dbinom{4}{2}=6\) farklı ikili karşılaştırma yapılabilir.

TukeyHSD(anova)  
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = x ~ faktor, data = veri)
## 
## $faktor
##          diff        lwr       upr     p adj
## B-A -15.16667 -21.654056 -8.679277 0.0000125
## C-A  11.83333   5.345944 18.320723 0.0002920
## D-A  16.33333   9.845944 22.820723 0.0000044
## C-B  27.00000  20.512610 33.487390 0.0000000
## D-B  31.50000  25.012610 37.987390 0.0000000
## D-C   4.50000  -1.987390 10.987390 0.2430915
plot(TukeyHSD(anova))  #karşılaştırmalar için %95 güven aralıklarını çizdiririz

Sonuca göre son satırda D-C için \(p-value=0.2430915>0.005\) olduğundan bu faktörler ile ilgili hipotezimiz \(H_{0}:\mu _{3}-\mu _{4}=0\) hipotezi kabul edilir. Yukarıda bulduğumuz sonuç ile aynıdır. Bu hipotez haricinde diğer 5 karşılaştırma için \(H_{0}\) hipotezleri reddedilir.

Ayrıca, güven aralığı grafiğine bakarak hipotez testini yapabiliriz. Grafiğin en altında D-C için olan grafik \(0\) değerini içerdiğinden \(H_{0}\) hipotezi kabul edilir. Diğer grafiklerde böyle bir durum söz konusu değildir.

Ayrıca, Tukey testini aşağıdaki gibi de yapabiliriz. Aynı sonuçlar elde edilir.

comparisons1<-glht(anova, linfct = mcp(faktor = "Tukey"))   #"multcomp" paketi 
summary(comparisons1)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = x ~ faktor, data = veri)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)    
## B - A == 0  -15.167      2.318  -6.544   <0.001 ***
## C - A == 0   11.833      2.318   5.105   <0.001 ***
## D - A == 0   16.333      2.318   7.047   <0.001 ***
## C - B == 0   27.000      2.318  11.649   <0.001 ***
## D - B == 0   31.500      2.318  13.590   <0.001 ***
## D - C == 0    4.500      2.318   1.941    0.243    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(comparisons1)   # tüm karsılatırmalar icin güven araliklari
## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = x ~ faktor, data = veri)
## 
## Quantile = 2.7993
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##            Estimate lwr      upr     
## B - A == 0 -15.1667 -21.6549  -8.6785
## C - A == 0  11.8333   5.3451  18.3215
## D - A == 0  16.3333   9.8451  22.8215
## C - B == 0  27.0000  20.5118  33.4882
## D - B == 0  31.5000  25.0118  37.9882
## D - C == 0   4.5000  -1.9882  10.9882
plot(confint(comparisons1))  # guven araliklari icin grafik

3. Fisher LSD, Scheffe ve Duncan Testleri

Bu testler için “agricolae” paketini kullancağız.

Fisher LSD testi

Fisher LSD testi ile aşağıdaki gibi tüm mümkün ikili karşılaştırmalar yapılabilir.

library(agricolae) 
out<-LSD.test(anova,"faktor",group=FALSE)   ### Fisher LSD test
out
## $statistics
##    MSerror Df     Mean       CV  t.value      LSD
##   16.11667 20 73.91667 5.431193 2.085963 4.834857
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none faktor   4  0.05
## 
## $means
##          x      std r      LCL      UCL Min Max  Q25  Q50   Q75
## A 70.66667 3.777124 6 67.24791 74.08543  65  75 68.5 71.0 73.50
## B 55.50000 4.636809 6 52.08124 58.91876  50  62 52.5 54.5 58.75
## C 82.50000 3.834058 6 79.08124 85.91876  77  88 80.5 82.5 84.50
## D 87.00000 3.741657 6 83.58124 90.41876  82  92 84.5 87.0 89.50
## 
## $comparison
##       difference pvalue signif.        LCL         UCL
## A - B   15.16667 0.0000     ***  10.331809  20.0015242
## A - C  -11.83333 0.0001     *** -16.668191  -6.9984758
## A - D  -16.33333 0.0000     *** -21.168191 -11.4984758
## B - C  -27.00000 0.0000     *** -31.834857 -22.1651425
## B - D  -31.50000 0.0000     *** -36.334857 -26.6651425
## C - D   -4.50000 0.0664       .  -9.334857   0.3348575
## 
## $groups
## NULL
## 
## attr(,"class")
## [1] "group"
#out$comparison

Yukarıdaki sonuçlarda olduğu gibi C-D için \(p-value=0.0664>0.005\) olduğundan \(H_{0}:\mu _{3}-\mu _{4}=0\) hipotezi kabul edilir.

Not: Normal dağılıma sahip iki grubun karşılaştırması için \(t\) testini kullanırız. Örneğin, A ve B gruplarını karşılaştırmak için aşağıdaki gibi yaparız.

t.test(x[1:6],x[7:12],alternative = "two.sided",var.equal = TRUE)  # varyanslar eşit değil ise "var.equal = FALSE" kullanılır.
## 
##  Two Sample t-test
## 
## data:  x[1:6] and x[7:12]
## t = 6.2119, df = 10, p-value = 9.989e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   9.726579 20.606754
## sample estimates:
## mean of x mean of y 
##  70.66667  55.50000

Fisher LSD testinde yapılan \(t\) testidir. Ancak,bu iki yöntemdeki \(t\) testlerinin serbestlik dereceleri arasında fark vardır. Bu nedenle sonuçları arasında farklılıklar olabilir.

Soru Tüm ikili karşılaştırmaları \(t\) testi kullanarak (yani “t.test”) test ediniz. Sonuçlarınızı Fisher LSD testinin sonuçları ile karşılaştırınız.

Scheffe testi

out1<-scheffe.test(anova,"faktor",group=TRUE)  #Scheffe test 
out1
## $statistics
##    MSerror Df        F     Mean       CV  Scheffe CriticalDifference
##   16.11667 20 3.098391 73.91667 5.431193 3.048799           7.066522
## 
## $parameters
##      test name.t ntr alpha
##   Scheffe faktor   4  0.05
## 
## $means
##          x      std r Min Max  Q25  Q50   Q75
## A 70.66667 3.777124 6  65  75 68.5 71.0 73.50
## B 55.50000 4.636809 6  50  62 52.5 54.5 58.75
## C 82.50000 3.834058 6  77  88 80.5 82.5 84.50
## D 87.00000 3.741657 6  82  92 84.5 87.0 89.50
## 
## $comparison
## NULL
## 
## $groups
##          x groups
## D 87.00000      a
## C 82.50000      a
## A 70.66667      b
## B 55.50000      c
## 
## attr(,"class")
## [1] "group"
#out1$comparison

Bu sonuçta $gruoups altında aynı harfler ile gösterilenler arasında anlamlı bir fark bulunmamaktadır, yani C ve D arasında. Benzer oalrak Fisher için de bu yapılabilir.

out1$groups
##          x groups
## D 87.00000      a
## C 82.50000      a
## A 70.66667      b
## B 55.50000      c

Ayrıca, scheffe.test fonksiyonunda group=FALSE alarak test istatistiklerini ve ilgili \(p-value\) değerlerini buluruz.

out2<-scheffe.test(anova,"faktor",group=FALSE)  
out2
## $statistics
##    MSerror Df        F     Mean       CV  Scheffe CriticalDifference
##   16.11667 20 3.098391 73.91667 5.431193 3.048799           7.066522
## 
## $parameters
##      test name.t ntr alpha
##   Scheffe faktor   4  0.05
## 
## $means
##          x      std r Min Max  Q25  Q50   Q75
## A 70.66667 3.777124 6  65  75 68.5 71.0 73.50
## B 55.50000 4.636809 6  50  62 52.5 54.5 58.75
## C 82.50000 3.834058 6  77  88 80.5 82.5 84.50
## D 87.00000 3.741657 6  82  92 84.5 87.0 89.50
## 
## $comparison
##       Difference pvalue sig        LCL        UCL
## A - B   15.16667 0.0000 ***   8.100144  22.233189
## A - C  -11.83333 0.0007 *** -18.899856  -4.766811
## A - D  -16.33333 0.0000 *** -23.399856  -9.266811
## B - C  -27.00000 0.0000 *** -34.066522 -19.933478
## B - D  -31.50000 0.0000 *** -38.566522 -24.433478
## C - D   -4.50000 0.3160     -11.566522   2.566522
## 
## $groups
## NULL
## 
## attr(,"class")
## [1] "group"

Duncan testi

out3<-duncan.test(anova,"faktor",group=FALSE)  
out3
## $statistics
##    MSerror Df     Mean       CV
##   16.11667 20 73.91667 5.431193
## 
## $parameters
##     test name.t ntr alpha
##   Duncan faktor   4  0.05
## 
## $duncan
##      Table CriticalRange
## 2 2.949998      4.834857
## 3 3.096506      5.074976
## 4 3.189616      5.227576
## 
## $means
##          x      std r Min Max  Q25  Q50   Q75
## A 70.66667 3.777124 6  65  75 68.5 71.0 73.50
## B 55.50000 4.636809 6  50  62 52.5 54.5 58.75
## C 82.50000 3.834058 6  77  88 80.5 82.5 84.50
## D 87.00000 3.741657 6  82  92 84.5 87.0 89.50
## 
## $comparison
##       difference pvalue signif.        LCL         UCL
## A - B   15.16667 0.0000     ***  10.331809  20.0015239
## A - C  -11.83333 0.0001     *** -16.668191  -6.9984761
## A - D  -16.33333 0.0000     *** -21.408309 -11.2583576
## B - C  -27.00000 0.0000     *** -32.074976 -21.9250243
## B - D  -31.50000 0.0000     *** -36.727576 -26.2724237
## C - D   -4.50000 0.0664       .  -9.334857   0.3348572
## 
## $groups
## NULL
## 
## attr(,"class")
## [1] "group"

C-D için \(p-value=0.0664>0.005\) olduğundan \(H_{0}:\mu _{3}-\mu _{4}=0\) hipotezi kabul edilir.

Örnek 2:

(Montgomery 3.8.1 Bölüm Chocolate and Cardivascular Health verisini inceleyelim.)

100 gr Dark Chocolate (DC), 200 ml yağlı süt ile birlikte 100 gr Dark Chocolate (DC+MK) ve 200 gr milk Chocolate (MC) olmak üzere 3 farklı biçimde çikolata yemenin kalp sağlığı üzerindeki etkisi araştırılmak isteniyor.

Bu amaçla yaşa ortalamaları 31.2 ile 33.2, ağırlık ortalamaları 62.7kg ile 68.9 kg ve vücut kitle indeksi 21.5 ile 22.3 kg/m^2 olan 7 kadın ve 5 erkek ile bu deney gerçekleştiriliyor.

Her bir kişi farklı günlerde bu 3 farklı tipteki çikolatadan bir tanesini yedikten 1 saat sonra kan plazmasındaki antioksidan değerleri ölçülmektedir. Rasgele olarak bu 12 kişi her bir tipteki çikolataları yedikten sonra ölçülen antioksidan değerleri aşağıdaki gibi bulunmuştur.

DC<-c(118.8,    122.6,  115.6   ,113.6, 119.5   ,115.9, 115.8   ,115.1, 116.9,  115.4   ,115.6, 107.9)
DC_MK<-c(105.4, 101.1,102.7 ,97.1   ,101.9  ,98.9   ,100.0, 99.8    ,102.6, 100.9,  104.5   ,93.5)
MC<-c(102.1,    105.8   ,99.6   ,102.7, 98.8,   100.9   ,102.8, 98.7,   94.7,   97.8,   99.7,   98.6)

data<-data.frame(c(DC, DC_MK,MC),rep(c("DC","DC_MK","MC"),rep(12,3)))
colnames(data)=c("antioxidant","grup")
data
##    antioxidant  grup
## 1        118.8    DC
## 2        122.6    DC
## 3        115.6    DC
## 4        113.6    DC
## 5        119.5    DC
## 6        115.9    DC
## 7        115.8    DC
## 8        115.1    DC
## 9        116.9    DC
## 10       115.4    DC
## 11       115.6    DC
## 12       107.9    DC
## 13       105.4 DC_MK
## 14       101.1 DC_MK
## 15       102.7 DC_MK
## 16        97.1 DC_MK
## 17       101.9 DC_MK
## 18        98.9 DC_MK
## 19       100.0 DC_MK
## 20        99.8 DC_MK
## 21       102.6 DC_MK
## 22       100.9 DC_MK
## 23       104.5 DC_MK
## 24        93.5 DC_MK
## 25       102.1    MC
## 26       105.8    MC
## 27        99.6    MC
## 28       102.7    MC
## 29        98.8    MC
## 30       100.9    MC
## 31       102.8    MC
## 32        98.7    MC
## 33        94.7    MC
## 34        97.8    MC
## 35        99.7    MC
## 36        98.6    MC
str(data)
## 'data.frame':    36 obs. of  2 variables:
##  $ antioxidant: num  119 123 116 114 120 ...
##  $ grup       : Factor w/ 3 levels "DC","DC_MK","MC": 1 1 1 1 1 1 1 1 1 1 ...
tapply(data$antioxidant,data$grup,summary)  #gruplar icin tanıimlayici istatistikler
## $DC
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   107.9   115.3   115.7   116.1   117.4   122.6 
## 
## $DC_MK
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   93.50   99.58  101.00  100.70  102.62  105.40 
## 
## $MC
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   94.70   98.67   99.65  100.18  102.25  105.80
boxplot(data$antioxidant ~ data$grup)    #gruplar icin box plot

Bu 3 grubu karşılaştırmak için varyans analizini ANOVA uygulayalım.

anova<-aov(antioxidant~grup,data=data)
summary(anova)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## grup         2 1952.6   976.3   93.58 2.52e-14 ***
## Residuals   33  344.3    10.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA tablosuna göre \(F\) istatistiğinin değeri \(93.58\) olup \(p-value<0.05\) olduğundan

\(H_{0}:\mu _{DC}=\mu_{DC+MK}=\mu _{MC}\)

hipotezi red edilir. Böylece, en az iki grubun ortalamaları birbirinden farklıdır. Şimdi bu farkılığı oluşturan grupları belirleyelim. Bu işleme geçmeden önce varsayımların kontrolünü yapalım.

Varsayım Kontrolü

1. Normallik:

ANOVA modelinde oluşan artıklar üzerinden normallik testleri yapılır.

artıklar<-print(residuals(anova))
##          1          2          3          4          5          6          7 
##  2.7416667  6.5416667 -0.4583333 -2.4583333  3.4416667 -0.1583333 -0.2583333 
##          8          9         10         11         12         13         14 
## -0.9583333  0.8416667 -0.6583333 -0.4583333 -8.1583333  4.7000000  0.4000000 
##         15         16         17         18         19         20         21 
##  2.0000000 -3.6000000  1.2000000 -1.8000000 -0.7000000 -0.9000000  1.9000000 
##         22         23         24         25         26         27         28 
##  0.2000000  3.8000000 -7.2000000  1.9166667  5.6166667 -0.5833333  2.5166667 
##         29         30         31         32         33         34         35 
## -1.3833333  0.7166667  2.6166667 -1.4833333 -5.4833333 -2.3833333 -0.4833333 
##         36 
## -1.5833333

Soru: Artıklar nasıl hesaplanır, hatırlayınız. Sizde artıkları hesaplayarak yukarıdakiler ile karşılaştırınız.

qqnorm(residuals(anova))

shapiro.test(residuals(anova))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(anova)
## W = 0.96254, p-value = 0.2572
ks.test(residuals(anova),"pnorm",mean(residuals(anova)),sd(residuals(anova)))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  residuals(anova)
## D = 0.11635, p-value = 0.6714
## alternative hypothesis: two-sided

Ayrıca, Anderson-Darling ve Cramer-Von Mises testleri de kullanılabilir. Dağılım parametreleri veriden tahmin edildiği için “goftest” paketinde yer alan “ad.test” ve “cvm.test” fonksiyonları ile uyum iyiliği testini (Goodness-of-Fit Test) yapabiliriz.

library(goftest)
ad.test(artıklar,"pnorm",mean=mean(artıklar),sd=sd(artıklar),estimated=TRUE)
## 
##  Anderson-Darling test of goodness-of-fit
##  Braun's adjustment using 6 groups
##  Null hypothesis: Normal distribution
##  with parameters mean = -6.64013597193345e-17, sd = 3.13644855271205
##  Parameters assumed to have been estimated from data
## 
## data:  artıklar
## Anmax = 1.7644, p-value = 0.554

Benzer olarak, Cramer-Von Mises testi

cvm.test(artıklar,"pnorm",mean=mean(artıklar),sd=sd(artıklar),estimated=TRUE)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Braun's adjustment using 6 groups
##  Null hypothesis: Normal distribution
##  with parameters mean = -6.64013597193345e-17, sd = 3.13644855271205
##  Parameters assumed to have been estimated from data
## 
## data:  artıklar
## omega2max = 0.28485, p-value = 0.619

Yukarıda yapılan tüm testlerin sonucunda \(p-value>0.05\) olduğundan bu testlerdeki sıfır hipotezimiz

\(H_{0}:\) Artıklar Normal dağılıma sahiptir.

hipotezi kabul edilir.

2. Homoscedasticity

Bartlett ve Levene testleri ile kontrol edelim.

bartlett.test(data$antioxidant ~ data$grup)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  data$antioxidant by data$grup
## Bartlett's K-squared = 0.4247, df = 2, p-value = 0.8087
library(car)
leveneTest(data$antioxidant, factor(data$grup)) #medyana göre
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.0213  0.979
##       33
leveneTest(data$antioxidant, factor(data$grup),mean) #ortalamaya göre
## Levene's Test for Homogeneity of Variance (center = mean)
##       Df F value Pr(>F)
## group  2  0.0129 0.9871
##       33

Yukarıdaki test sonuçlarına göre \(p-value>0.05\) olduğundan gruplar arasında varyansların homojen olduğu hipotezi kabul edilir.

Not: Grupların örneklem varyansları:

 tapply(data$antioxidant,data$grup,sd)^2
##        DC     DC_MK        MC 
## 12.484470 10.465455  8.350606

İkili Karşılaştırmalar

Tüm mümkün ikli karrşılaştırmaları yapalım.

1. Fisher LSD test

library(agricolae) 
out<-LSD.test(anova,"grup",group=FALSE)   ### Fisher LSD test
out
## $statistics
##    MSerror Df     Mean       CV  t.value      LSD
##   10.43351 33 105.6472 3.057434 2.034515 2.682876
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none   grup   3  0.05
## 
## $means
##       antioxidant      std  r       LCL      UCL   Min   Max     Q25    Q50
## DC       116.0583 3.533337 12 114.16125 117.9554 107.9 122.6 115.325 115.70
## DC_MK    100.7000 3.235035 12  98.80292 102.5971  93.5 105.4  99.575 101.00
## MC       100.1833 2.889742 12  98.28625 102.0804  94.7 105.8  98.675  99.65
##           Q75
## DC    117.375
## DC_MK 102.625
## MC    102.250
## 
## $comparison
##            difference pvalue signif.       LCL       UCL
## DC - DC_MK 15.3583333 0.0000     *** 12.675458 18.041209
## DC - MC    15.8750000 0.0000     *** 13.192124 18.557876
## DC_MK - MC  0.5166667 0.6977         -2.166209  3.199542
## 
## $groups
## NULL
## 
## attr(,"class")
## [1] "group"

$comparison kısmında görüldüğü gibi DC ile DC_MK ve MC grupları (faktör düzeyleri) arasındaki karşılaştırmalar için \(p-value<0.05\) olduğundan bunlarla ilgili sıfır hipotezleri \(H_{0}:\mu _{DC}=\mu_{DC+MK}\) ve \(H_{0}:\mu _{DC}=\mu _{MC}\) hipotezleri red edilir.

Ancak, \(H_{0}:\mu _{DC+MK}=\mu_{MC}\) hipotezi kabul edilir. Çünkü, bu karşılaştırma için \(p-value=0.6977>0.05\) dır. Böylece, DC_MK ile MC grupları arasında anlamlı bir fark yoktur. Yani, Dark Chocolate ile süt içenler ile Milk Chocolate yiyen kişilerin kanlarındaki antioksidan değerleri bakımından anlamlı bir farklılık yoktur.

Aytıca, LSD.test fonksiyonunda group=TRUE seçersek $groups kısmında görüldüğü gibi DC_MK ve MC faktör düzeyleri aynı grupta olur, yani aralarında anlamlı bir fark yoktur.

out<-LSD.test(anova,"grup",group=TRUE)   ### Fisher LSD test
out
## $statistics
##    MSerror Df     Mean       CV  t.value      LSD
##   10.43351 33 105.6472 3.057434 2.034515 2.682876
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none   grup   3  0.05
## 
## $means
##       antioxidant      std  r       LCL      UCL   Min   Max     Q25    Q50
## DC       116.0583 3.533337 12 114.16125 117.9554 107.9 122.6 115.325 115.70
## DC_MK    100.7000 3.235035 12  98.80292 102.5971  93.5 105.4  99.575 101.00
## MC       100.1833 2.889742 12  98.28625 102.0804  94.7 105.8  98.675  99.65
##           Q75
## DC    117.375
## DC_MK 102.625
## MC    102.250
## 
## $comparison
## NULL
## 
## $groups
##       antioxidant groups
## DC       116.0583      a
## DC_MK    100.7000      b
## MC       100.1833      b
## 
## attr(,"class")
## [1] "group"

2. Tukey test

TukeyHSD(anova)  
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = antioxidant ~ grup, data = data)
## 
## $grup
##                 diff        lwr        upr     p adj
## DC_MK-DC -15.3583333 -18.594104 -12.122562 0.0000000
## MC-DC    -15.8750000 -19.110771 -12.639229 0.0000000
## MC-DC_MK  -0.5166667  -3.752438   2.719104 0.9190724
plot(TukeyHSD(anova))

Tukey testinin sonuç kısmındaki $grup dan görüldüğü gibi sadece MC ile DC_MK grupları arasında anlamlı bir fark yoktur. Bu sonuca iki farklı biçimde ulaşabilirsiniz.

  1. “p adj” değeri yani düzeltilmiş \(p-value\) kullanarak. Burada, MC ile DC_MK grupları için \(p-value=0.9190724>0.05\) olduğundan bu gruplar arasında anlamlı bir fark yoktur.

  2. İkili karşılaştırma için oluşturulan %95 lik güven aralığını kullanarak. Eğer bu güven aralığı \(0\) değerini içeriyorsa sıfır hipotezi \(H_{0}\) kabul edilir. Yukarıdaki sonuçta sadece MC ile DC_MK grupları için bu durum gerçekleşmiştir.

Sonuç olarak, bu örneğimiz için Fisher LSD ve Tukey testleri aynı sonuçları vermiştir.

Soru 1: Örnek 2 için diğer karşılaştırma testlerini kullanarak ikili karşılaştırmaları yapınız.

Soru 2: Montgomery, Design and Analysis of Experiments (9th Edition) kitabında, 3. bölüm sonunda bulunan 3.14, 3.15, 3.16 sorularını (sayfa 127) ilk olarak kağıda çözünüz ve sonrasında R programı ile de çözerek çözümlerinizi karşılaştırınız.