Seandainya kita ingin melakukan percobaan BIB dengan sejumlah \(v\) perlakuan dan \(b\) kelompok. Namun setiap kelompok hanya memuat \(k\) perlakuan, sehingga setiap perlakuan diulang sebanyak \(r\) kali dalam percobaan, dengan banyaknya ulangan setiap pasangan perlakuan adalah \(\lambda\).
library(ibd)
set.seed(123)
des<-bibd(v=4, b=4, r=3, k=3, lambda=2)
## v: banyaknya perlakuan
## b: banyaknya kelompok/blok
## r: banyaknya ulangan untuk masing-masing perlakuan
## k: banyaknya perlakuan pada setiap blok
## lambda: banyaknya ulangan pada setiap pasangan perlakuan
des$design## [,1] [,2] [,3]
## Block-1 1 3 4
## Block-2 2 3 4
## Block-3 1 2 4
## Block-4 1 2 3Output di atas menunjukkan bahwa pada blok 1, digunakan perlakuan 1, 3, dan 4. Pada blok 2 digunakan perlakuan 2, 3, dan 4, dan seterusnya. Perlu diingat pula bahwa pengacakan perlu dilakukan dalam menempatkan urutan perlakuan di setiap blok.
Illustration from Montgomery (2013)
Silahkan perhatikan kembali ilustrasi yang dibahas pada modul praktikum STA 521 pertemuan ke-2, yang diambil dari Montgomery (2013). Perhatikan pula tabulasi data yang diinput pada program berikut ini.
gasoline<-read.csv("https://raw.githubusercontent.com/raoy/data/master/example.csv")
gasoline## Treatment Block Response
## 1 1 1 73
## 2 3 1 73
## 3 4 1 75
## 4 1 2 74
## 5 2 2 75
## 6 3 2 75
## 7 2 3 67
## 8 3 3 68
## 9 4 3 72
## 10 1 4 71
## 11 2 4 72
## 12 4 4 75xtabs(Response~Treatment+Block, data=gasoline)## Block
## Treatment 1 2 3 4
## 1 73 74 0 71
## 2 0 75 67 72
## 3 73 75 68 0
## 4 75 0 72 75Sebagai catatan, nilai 0 pada output tabel di atas menunjukkan bahwa tidak ada pengamatan pada pasangan perlakuan dan blok tersebut.
Uji hipotesis untuk pengaruh perlakuan dapat dilakukan dengan uji-F yang merupakan rasio antara kuadrat tengah galat dari perlakuan (adjusted) terhadap kuadrat tengah galat. Begitu pula halnya dengan uji hipotesis untuk pengaruh blok, dilakukan menggunakan kuadrat tengah galat dari blok (adjusted).
library(daewr)
mod1<-aov( Response ~ as.factor(Block)+as.factor(Treatment), data = gasoline)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Block) 3 55.00 18.333 28.20 0.00147 **
## as.factor(Treatment) 3 22.75 7.583 11.67 0.01074 *
## Residuals 5 3.25 0.650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Silahkan buka kembali catatan Anda pada praktikum minggu lalu, dan perhatikan apakah F-hitung yang diperoleh pada perhitungan manual sama dengan output di atas?
Perhatikan bahwa tabel anova di atas belum menampilkan nilai blok (adjusted). Oleh karenanya, kita perlu menambahkan perintah berikut ini.
drop1(mod1,test="F")## Single term deletions
##
## Model:
## Response ~ as.factor(Block) + as.factor(Treatment)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 3.250 -1.675
## as.factor(Block) 3 66.083 69.333 29.048 33.889 0.0009528 ***
## as.factor(Treatment) 3 22.750 26.000 17.278 11.667 0.0107387 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Silahkan periksa kembali perhitungan manual Anda minggu lalu, dan bandingkan dengan F-hitung pada tabel di atas.
Untuk memperoleh tabel anova dengan nilai yang sudah di-adjust pada blok dan perlakuan, kita dapat pula menggunakan fungsi aov.ibd seperti di bawah ini.
mod2<-ibd::aov.ibd(Response ~ as.factor(Block)+as.factor(Treatment), data = gasoline)
mod2$ANOVA.table## Anova Table (Type III tests)
##
## Response: Response
## Sum Sq Df F value Pr(>F)
## (Intercept) 8948.7 1 13767.198 8.528e-10 ***
## as.factor(Block) 66.1 3 33.889 0.0009528 ***
## as.factor(Treatment) 22.7 3 11.667 0.0107387 *
## Residuals 3.3 5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(lsmeans)
lsmeans(mod1,pairwise~Treatment,adjust=("tukey"))## $lsmeans
## Treatment lsmean SE df lower.CL upper.CL
## 1 71.4 0.487 5 70.1 72.6
## 2 71.6 0.487 5 70.4 72.9
## 3 72.0 0.487 5 70.7 73.3
## 4 75.0 0.487 5 73.7 76.3
##
## Results are averaged over the levels of: Block
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1 - 2 -0.250 0.698 5 -0.358 0.9825
## 1 - 3 -0.625 0.698 5 -0.895 0.8085
## 1 - 4 -3.625 0.698 5 -5.192 0.0130
## 2 - 3 -0.375 0.698 5 -0.537 0.9462
## 2 - 4 -3.375 0.698 5 -4.834 0.0175
## 3 - 4 -3.000 0.698 5 -4.297 0.0281
##
## Results are averaged over the levels of: Block
## P value adjustment: tukey method for comparing a family of 4 estimatesgasoline$Treatment<-as.factor(gasoline$Treatment)
c1<-c(1,1,-1,-1)
c2<-c(1,-1,0,0)
c3<-c(0,0,1,-1)
contr<-cbind(c1, c2, c3)
mod3<-aov( Response ~ as.factor(Block)+Treatment, contrasts=list(Treatment=contr), data = gasoline)
summary.aov(mod3, split=list(Treatment=list("1-2 Vs 3-4"=1,
"1 vs 2"=2,
"3 vs 4"=3)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Block) 3 55.00 18.333 28.205 0.00147 **
## Treatment 3 22.75 7.583 11.667 0.01074 *
## Treatment: 1-2 Vs 3-4 1 10.67 10.667 16.410 0.00982 **
## Treatment: 1 vs 2 1 0.08 0.083 0.128 0.73492
## Treatment: 3 vs 4 1 12.00 12.000 18.462 0.00774 **
## Residuals 5 3.25 0.650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Lawson, J. (2014). Design and Analysis of Experiments with R (Vol. 115). CRC press.
Meier, L. (2021, February 3). Chapter 9 incomplete block designs | ANOVA: A short intro using R. https://stat.ethz.ch/~meierluk/teaching/anova/incomplete-block-designs.html
Montgomery, D. C. (2013). Design and analysis of experiments. John wiley & sons.