Analisis Ragam 2 Arah
Irwan Lovadi
3 September 2020
Tulisan ini memuat panduan analisis ragam dua arah (ANOVA 2 jalur/two-way ANOVA) dengan menggunakan contoh fiktif. Sebelum menggunakan kode R di bawah ini, penting untuk diketahui bahwa sebelum melakukan analisis ragam dua arah, asumsi untuk analisis tersebut diperiksa terlebih dahulu. Asumsi tersebut yaitu kehomogenan ragam dan normalitas distribusi data. Jika data tidak terdistribusi normal, maka lakukan transformasi data terlebih dahulu. Panduan untuk analisis ragam dua arah menggunakan tiga paket R, yaitu car, agricolae dan ggpubr. Jika akan melakukan analisis ini di laptop atau komputer Anda, pastikan bahwa ketiga paket tersebut telah terpasang.
PENTING
Untuk menggunakan panduan ini, pertama jalankan program R Studio, lalu klik File - New File - R script. Setelah itu, ketik kode R di tiap blok abu-abu yang tertera di bawah ini. Untuk menjalankan kode di tiap baris, arahkan pointer di awal kode atau di akhir kode kemudian tekan tombol ctrl (control) bersamaan dengan menekan tombol enter. Alternatif lainnya klik run yang ada di bagian kanan atas panel source.
Hasil operasi dari kode R yang dijalankan akan muncul di bagian panel console (panel di bagian kiri bawah) atau panel fitur (panel di bagian kanan bawah) jika hasil berupa grafik. Jika pekerjaan Anda sudah selesai, klik File dan pilih Save. Simpan hasil pekerjaan Anda di folder yang Anda inginkan.
Memuat paket
library(car) ## Loading required package: carDatalibrary(agricolae)
library(ggpubr) ## Loading required package: ggplot2Menyusun dataset
Dataset dapat disusun langsung di R jika titik data tidak terlalu banyak dan variabel yang dianalisis sedikit. Walaupun demikian disarankan agar dataset dibuat terlebih dahulu di program Excel dan disimpan dalam format .csv.
Dalam panduan ini, dataset dibuat langsung di R dan kemudian ditransformasi ke dalam format csv dengan menggunakan R. Jika Anda sudah memiliki data yang telah disusun di Excel atau program sejenisnya dan disimpan dalam format .csv, maka lewati tahapan “Menyusun dataset” dan langsung mulai dari “Memuat data”.
Data yang akan dianalisis di panduan ini merupakan data buatan. Seorang peneliti akan membandingkan kerapatan kerang bivalvia di tiga lokasi yang berbeda (A, B, dan C) dan di tiga zona (pinggir, tengah, dalam)
# Menyusun komponen dataset
Lokasi <- c(rep("A", 15), rep("B", 15), rep("C", 15))
Zona <- c(rep("Pinggir", 5), rep("Tengah", 5), rep("Dalam", 5), rep("Pinggir", 5), rep("Tengah", 5), rep("Dalam", 5), rep("Pinggir", 5), rep("Tengah", 5), rep("Dalam", 5))
Kerapatan <- c(15,18,21,14,17,23,26,19,18,21,25,27,29,23,26,25,27,31,29,26,21,23,25,19,24,26,29,31,27,25,12,13,15,11,10,17,19,21,18,16,9,8,11,7,9)
# Atur Lokasi dan Zona sebagai faktor
Lokasi <- factor(Lokasi)
Zona <- factor(Zona)
# Menyusun dataset
data<- data.frame(Lokasi, Zona, Kerapatan)
# Mengubah dataset ke dalam format .csv dan menyimpannya di salah satu folder di komputer atau laptop
## Di panduan ini, folder penyimpanan dataset adalah Kode R
### folder path dari file .csv adalah C:/DATA/Untan/Lab Ekologi/Pelatihan R/Kode R/
write.csv(data, "C:/DATA/Untan/Lab Ekologi/Pelatihan R/Kode R/bivalvia.csv")Memuat dan mengatur data
# Pengaturan direktori/folder kerja
setwd("C:/DATA/Untan/Lab Ekologi/Pelatihan R/Kode R/")
df <- read.csv("C:/DATA/Untan/Lab Ekologi/Pelatihan R/Kode R/bivalvia.csv")
# Mengatur urutan Zona supaya berurutan sesuai dengan yang diinginkan
df$Zona <- ordered(df$Zona, levels=c("Pinggir", "Tengah", "Dalam"))Pemeriksaan varians (Uji Levene’s)
Bagian ini menunjukan kode untuk pemeriksaan kehomogenan ragam, salah satu asumsi yang harus diperiksa sebelum analisis ragam satu arah. Uji yang digunakan adalah Uji Levene’s. Data memenuhi asumsi jika nilai p > 0.05.
# leveneTest : kode untuk instruksi pemeriksaan kehomogenan ragam dengan Uji Levene.
leveneTest(Kerapatan ~ Lokasi*Zona, data = df)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 8 0.4073 0.9089
## 36Pemeriksaan normalitas (dengan Uji Shapiro dan pemeriksaan grafis)
Kode di bawah ini digunakan untuk pemeriksaan kenormalan distribusi data. Ada 2 metode yang dipakai dalam panduan ini, yaitu Uji Shapiro-Wilk dan Mencermati QQ-Plot. Jika akan menggunakan data lain dengan nama kelompok/perlakuan serta variabel yang berbeda, ganti kata Kerapatan dengan nama variabel dari data yang akan dibandingkan.
Data memenuhi asumsi: 1) jika nilai p > 0.05 pada Uji Shapiro atau 2) jika titik data berada di antara garis putus-putus pada QQ-plot.
# shapiro.test: kode untuk instruksi pemeriksaan normalitas data dengan Uji Shapiro.
model=lm(Kerapatan ~ Lokasi*Zona, df)
shapiro.test(model$residuals)##
## Shapiro-Wilk normality test
##
## data: model$residuals
## W = 0.96939, p-value = 0.2752# qqPlot: kode untuk menyusun Q-Q plot untuk mengetahui apakah data memenuhi asumsi atau tidak.
qqPlot(df$Kerapatan)
## [1] 44 42Factorial ANOVA Type I (untuk jumlah ulangan seragam)
Hasil pemeriksaan menunjukkan bahwa asumsi ANOVA terpenuhi. Untuk itu, kode R untuk analisis ragam 2 arah dapat dijalankan. Jika nilai p < 0,05 untuk tiap faktor utama, berarti bahwa ada perbedaan yang signifikan. Jika nilai p < 0,05 untuk Interaksi, uji lanjut dapat dilakukan.
#summary: kode instruksi untuk penyusunan rangkuman nilai hasil analisa data.
Kerapatan.mod = aov(Kerapatan ~ Lokasi*Zona, data = df)
summary(Kerapatan.mod)## Df Sum Sq Mean Sq F value Pr(>F)
## Lokasi 2 1268.8 634.4 114.65 2.43e-16 ***
## Zona 2 32.5 16.3 2.94 0.0657 .
## Lokasi:Zona 4 486.7 121.7 21.99 2.98e-09 ***
## Residuals 36 199.2 5.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Post Hoc: Tukey
Hasil analisis ragam satu arah di atas menunjukkan bahwa faktor Lokasi dan Interaksi signifikan. Untuk itu perlu dilakukan Uji Tukey, khususnya pada Interaksi. Dalam contoh data di atas, diketahui bahwa jumlah ulangan masing-masing taraf/perlakuan/kelompok sama, yaitu 5. Kode dibawah ini juga memungkinkan Uji Tukey dengan jumlah ulangan yang tidak sama. Caranya adalah dengan mengganti nilai “unbalanced” dari FALSE ke TRUE.
Catatan: Huruf yang berbeda pada hasil Uji Tukey menunjukkan adanya perbedaan yang signifikan.
# HSD.test: kode instruksi untuk uji lanjut - Tukey HSD)
(HSD.test(model, c("Lokasi", "Zona"), unbalanced = FALSE)) ## $statistics
## MSerror Df Mean CV MSD
## 5.533333 36 20.13333 11.68363 4.905172
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Lokasi:Zona 9 4.662789 0.05
##
## $means
## Kerapatan std r Min Max Q25 Q50 Q75
## A:Dalam 26.0 2.236068 5 23 29 25 26 27
## A:Pinggir 17.0 2.738613 5 14 21 15 17 18
## A:Tengah 21.4 3.209361 5 18 26 19 21 23
## B:Dalam 27.6 2.408319 5 25 31 26 27 29
## B:Pinggir 27.6 2.408319 5 25 31 26 27 29
## B:Tengah 22.4 2.408319 5 19 25 21 23 24
## C:Dalam 8.8 1.483240 5 7 11 8 9 9
## C:Pinggir 12.2 1.923538 5 10 15 11 12 13
## C:Tengah 18.2 1.923538 5 16 21 17 18 19
##
## $comparison
## NULL
##
## $groups
## Kerapatan groups
## B:Dalam 27.6 a
## B:Pinggir 27.6 a
## A:Dalam 26.0 ab
## B:Tengah 22.4 bc
## A:Tengah 21.4 bcd
## C:Tengah 18.2 cd
## A:Pinggir 17.0 de
## C:Pinggir 12.2 ef
## C:Dalam 8.8 f
##
## attr(,"class")
## [1] "group"Visualisasi hasil Uji Tukey
Selain penyajian dalam bentuk tabel, hasil ANOVA 2 Jalur juga dapat disajikan dalam bentuk grafik. Dalam tutorial ini, hasil uji lanjut disajikan dalam bentuk grafik plot galat (error plot).Modifikasi grafik untuk data yang lain dapat dilakukan dengan mengganti nama faktor yang signifikan dari ANOVA (“Lokasi), variabel yang diamati (”Kerapatan“), warna (”Zona"), nilai ylim, posisi label geom_text (nilai x dan y) dan label.
ggerrorplot(df,"Lokasi", "Kerapatan", color = "Zona", palette = "locuszoom", desc_stat = "mean_se", ylim=c(0, 30)) +
geom_text(x = 0.73, y = 19.5, label = "de", size = 3.5) +
geom_text(x = 1, y = 25, label = "bcd", size = 3.5) +
geom_text(x = 1.27, y = 29, label = "ab", size = 3.5) +
geom_text(x = 1.74, y = 31, label = "a", size = 3.5) +
geom_text(x = 2, y = 25.5, label = "bc", size = 3.5) +
geom_text(x = 2.27, y = 31, label = "a", size = 3.5) +
geom_text(x = 2.73, y = 15, label = "ef", size = 3.5) +
geom_text(x = 3, y = 21, label = "cd", size = 3.5) +
geom_text(x = 3.27, y = 11, label = "f", size = 3.5) +
geom_text(x = 3, y = 2, label= "ANOVA 2 Jalur, p < 0,05")